HDU 6053 莫比乌斯反演

那道题想到枚举1e5以内的质因子, 但是会重复
然后队友提醒了mobius反演, 果然还是太菜了啊

题解

首先, 题目提到任意区间满足条件, 也就是gcd(b[1], b[2], b[3]…b[n]) >= 2 就行了, 很容易得出总的b数量为

sum=∏i=1na[i]

而其中不满足条件的就是gcd(b[1], b[2], b[3]…b[n]) = 1的数量
我们定义
F(n)为gcd为n的倍数的b数量, f(n)为gcd为n的b数量
则有
F(n)=∏i=1a[i]/n

F(n)=∑n|df(d)

反演得:

f(n)=∑n|du(d/n)∗F(d)

所以

f(1)=∑i=1min(a[1],a[2]...a[n])u(i)∗F(i)

直接求会超时, 注意到是整数除法, 读入的时候预处理分块用快速幂可以降低时间复杂度

然后sum-f(1)即可

code:

#include <iostream>#include <cstdio>#include <cstring>using namespace std;typedef long long ll;const int N = 100010;const ll mod = 1e9 + 7;int mobius[N];int prime[N];int a[N];bool vis[N];int cnt;int pre[N << 1];int n;inline ll min(ll a, ll b){    return a <= b ? a : b;}inline ll max(ll a, ll b){    return a >= b ? a : b;}ll fastPowMod(ll a, ll b, ll mod){    ll res = 1 % mod;    while(b)    {        if(b & 1) res = res * a % mod;        a = a * a % mod;        b >>= 1;    }    return res;}void getMobius(){    cnt = 0;    memset(vis, 0, sizeof vis);    memset(mobius, 0, sizeof mobius);    mobius[1] = 1;    for(int i = 2; i <= 100000; ++i)    {        if(!vis[i])        {           prime[++cnt] = i;           mobius[i] = -1;        }        for(int j = 1; j <= cnt && i * prime[j] <= 100000; ++j)        {            vis[i * prime[j]] = 1;            if(i % prime[j] == 0)            {                mobius[i * prime[j]] = 0;                break;            }            mobius[i * prime[j]] = -mobius[i];        }    }}int main(){    getMobius();    int t, kas = 0;    scanf("%d", &t);    while(t--)    {        scanf("%d", &n);        memset(pre, 0, sizeof pre);        ll mn = N, mx = 0;        ll a_multi = 1;        for(int i = 1; i <= n; ++i)        {            scanf("%d", a + i);            ++pre[a[i]];            a_multi *= a[i];            a_multi %= mod;            mn = min(mn, a[i]);            mx = max(mx, a[i]);        }        for(int i = 1; i < (N << 1); ++i) pre[i] += pre[i - 1];//前缀和用于分块        ll rs = 0;        for(int i = 1; i <= mn; ++i)        {            ll r = 1;            for(ll j = 1; j <= mx / i; ++j)            {                r = r * fastPowMod(j, pre[j * i + i - 1] - pre[j * i - 1], mod) % mod;//            }            r = (r * mobius[i] % mod + mod) % mod;            rs = (rs + r) % mod;        }        printf("Case #%d: %I64d\n", ++kas, ((a_multi - rs) % mod + mod) % mod);    }    return 0;}