CF839 D 莫比乌斯反演

题目大意

定义一个牛逼序列，内含若干个下标，记作(i_1,i_2,i_3...i_k),要求(i_1 < i_2 < i_3...<i_k)。同时(gcd(a_{i_1},a_{i_2},...a_{i_k}) > 1)。这样的一个牛逼序列，会产生其(gcd*k)的牛逼值。给定一个序列({a})。求他所有可能的牛逼序列的牛逼值的总和。

思路

首先不难想到一个递推：(f[i])表示选([1,i])的元素并且强制选择(a_i)。接着会发现事实上转移是不能进行的，因为这个划分不能控制最大公约数的能不能取到某个特定的值。把(i)维度换成约数事实上也是一样的，因为根本不能从现有信息推知选择某些数得到的最大公约数能否是一个特定的值而非别的。

考虑划分答案(f[n] = sumlimits_{S} |S| * n [gcd(S) == n])其中(S)表示选出的集合。直接求比较困难，考虑套路反演：(g[n] = sumlimits_{S}|S|[n | gcd(S)]),满足(g[n] = sumlimits_{n | d}f(d) / d)。反演拿到(f[n] = n * sumlimits_{n | d}u(d/n)*g(d))。考虑求(g[n])，根据定义，需要求所有(gcd)是(n)的倍数的长度，这样的集合可以通过排列组合求：因为(gcd)是(n)的倍数，所以不管怎么选，只能从是(n)的倍数的集合中选数，设是(n)的倍数的数有(cnt[n])个，那么(g[n] = sumlimits_{L = 1}^{cnt[d]}L*C_{cnt[d]}^L = cnt[d] * (2^{cnt[d] - 1}))。故(f[n] = n * sumlimits_{n | d}u(d / n) * cnt[d] * 2 ^{cnt[d] - 1})。换枚举对象，就有(f[n] = n * sumlimits_{i = 1}^{N / n}u(i)*cnt[i*d]*2^{cnt[i*d] - 1})。

又已知答案是(sumlimits_{i = 2}^N f[i] = sumlimits_{i = 2}^Nisumlimits_{j = 1}^{N / n}u(j)*cnt[i * j] * 2^{cnt[i * d] - 1} = sumlimits_{k = 2}^N cnt[k] * 2^{cnt[k] - 1}sumlimits_{d geq 2,d|k}d*u(k/d))。对于最后一部分，若加上(d=1)的部分，则求和答案为(phi(k))，但是由于(dgeq 2)，所以答案要去掉(d=1)的部分，由于(d=1)总是有(d|k)，所以这部分求和可以直接表示成(phi(k) - u(k))。所以答案(sumlimits_{k=2}^N cnt[k] * 2^{cnt[k] - 1]}*(phi(k) - u(k)))。

代码

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define forn(i,x,n) for(int i = x;i <= n;++i)
#define forr(i,x,n) for(int i = n;i >= x;--i)
#define Angel_Dust ios::sync_with_stdio(0);cin.tie(0)

const int N = 1e6+7,MOD = 1e9+7;
int primes[N],cnt,a[N];
ll mobius[N],phi[N],dcnt[N],pw[N];
bool st[N];

void init()
{
    mobius[1] = 1;phi[1] = 1;
    for(int i = 2;i < N;++i)
    {
        if(!st[i])  primes[cnt++] = i,mobius[i] = -1,phi[i] = i - 1;
        for(int j = 0;j < cnt && primes[j] * i < N;++j)
        {
            int t = primes[j] * i;
            st[t] = 1;
            if(i % primes[j] == 0)
            {
                phi[t] = phi[i] * primes[j];
                mobius[t] = 0;
                break;
            }
            mobius[t] = mobius[i] * -1;
            phi[t] = phi[i] * (primes[j] - 1);
        }
    }
    pw[0] = 1;
    forn(i,1,N - 1) pw[i] = pw[i - 1] * 2 % MOD;
}

int main()
{
//    freopen("D:\OJCodes\TTT\in.txt","r",stdin);
//    freopen("D:\OJCodes\TTT\out.txt","w",stdout);

    init();
    int n;scanf("%d",&n);
    forn(i,1,n)
    {
        scanf("%d",&a[i]);
        for(int j = 1;j <= a[i] / j;++j)
        {
            if(a[i] % j)    continue;
            ++dcnt[j];
            if(j != a[i] / j)   ++dcnt[a[i] / j];
        }
    }
    ll res = 0;
    forn(k,2,1000000)
    {
        ll nxt = dcnt[k];
        nxt = nxt * pw[dcnt[k] - 1] % MOD;
        nxt = nxt * (phi[k] - mobius[k]) % MOD;
        nxt = ((nxt) % MOD + MOD) % MOD;
        res = (res + nxt) % MOD;
    }
    printf("%lld
",res);
//    fclose(stdin);
//    fclose(stdout);
    return 0;
}