HDU 6053 莫比乌斯反演

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

题解

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

$$sum = \prod_{i = 1}^n a[i]$$

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

$$F(n) = \sum_{n|d} f(d)$$
反演得:
$$f(n) = \sum_{n|d}u(d / n) * F(d)$$
所以
$$f(1) = \sum_{i = 1}^{min(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;
}