ZOJ 3868 GCD Expectation (容斥＋莫比烏斯反演)

GCD Expectation

---

Time Limit: 4 Seconds     Memory Limit: 262144 KB

---

Edward has a set of n integers {a₁,a₂,...,a_n}. He randomly picks a nonempty subset {x₁,x₂,…,x_m} (each nonempty subset has equal probability to be picked), and would like to know the expectation of [gcd(x₁,x₂,…,x_m)]^k.

Note that gcd(x₁,x₂,…,x_m) is the greatest common divisor of {x₁,x₂,…,x_m}.

Input

There are multiple test cases. The first line of input contains an integerT indicating the number of test cases. For each test case:

The first line contains two integers n,k (1 ≤ n, k ≤ 10⁶). The second line containsn integers a₁, a₂,…,a_n (1 ≤a_i ≤ 10⁶).

The sum of values max{a_i} for all the test cases does not exceed 2000000.

Output

For each case, if the expectation is E, output a single integer denotesE · (2ⁿ - 1) modulo 998244353.

Sample Input

1
5 1
1 2 3 4 5

Sample Output

42

---

Author: LIN, Xi
Source: The 15th Zhejiang University Programming Contest

題目連結：http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemId=5480

題目大意：給一個集合，{xi}為它的一個非空子集，設E為[gcd(x₁,x₂,…,x_m)]^k   的期望，求E*(2^n - 1) mod 998244353

題目分析：首先一個有n個元素的集合的非空子集個數為2^n - 1，所以E的分母就是2^n - 1了，因此我們要求的只是E的分子,

設F(x)為gcd(xi) = x的個數，那麼ans = (1^k) * F(1) + (2^k) * F(2) + ... + (ma^k) * F(ma)

下面的問題就是如何快速的計算F(x)了，對於一個集合，先計算出x的倍數的個數，nlogn即可，然後就是基礎的容斥，假設現在要求gcd為1的，那就減去gcd為2的，gcd為3的，注意到6同時是2和3的倍數，也就是6的倍數被減了兩次，所以要加上gcd為6的，前面的係數剛好是數字對應的莫比烏斯函式，看到這題很多用dp來容斥的，其實本質和莫比烏斯函式一樣，但是莫比烏斯函式寫起來真的很簡單，2333333

#include <cstdio>
#include <cstring>
#include <algorithm>
#define ll long long
using namespace std;
int const MOD = 998244353;
int const MAX = 1e6 + 5;
ll two[MAX];
int p[MAX], mob[MAX], num[MAX], cnt[MAX];
bool noprime[MAX];
int n, k, ma, pnum;

void Mobius()
{
    pnum = 0;
    mob[1] = 1;
    for(int i = 2; i < MAX; i++)
    {
        if(!noprime[i])
        {
            p[pnum ++] = i;
            mob[i] = -1;
        }
        for(int j = 0; j < pnum && i * p[j] < MAX; j++)
        {
            noprime[i * p[j]] = true;
            if(i % p[j] == 0)
            {
                mob[i * p[j]] = 0;
                break;
            }
            mob[i * p[j]] = -mob[i];
        }
    }
}

ll qpow(ll x, ll n)
{
    ll res = 1;
    while(n != 0)
    {
        if(n & 1)
            res = (res * x) % MOD;
        x = (x * x) % MOD;
        n >>= 1;
    }
    return res;
}

void pre()
{
    Mobius();   
    two[0] = 1;
    for(int i = 1; i < MAX; i++)
        two[i] = two[i - 1] * 2ll % MOD;
}

int main()
{
    pre();
    int T;
    scanf("%d", &T);
    while(T --)
    {
        memset(num, 0, sizeof(num));
        memset(cnt, 0, sizeof(cnt));
        ma = 0;
        int tmp;
        scanf("%d %d", &n, &k);
        for(int i = 0; i < n; i++)
        {
            scanf("%d", &tmp);
            cnt[tmp] ++;
            ma = max(ma, tmp);
        }
        for(int i = 1; i <= ma; i++)
            for(int j = i; j <= ma; j += i)
                num[i] += cnt[j];       //求i的倍數的個數
        ll ans = 0;
        for(int i = 1; i <= ma; i++)    //列舉gcd
        {
            ll sum = 0;
            for(int j = i; j <= ma; j += i)  //容斥
                sum = (MOD + sum % MOD + mob[j / i] * (two[num[j]] - 1) % MOD) % MOD;
            ans = (MOD + ans % MOD + (sum * qpow(i, k)) % MOD) % MOD;
        }
        printf("%lld\n", ans);
    }
}