ZOJ 3868 GCD Expectation 莫比乌斯反演

GCD Expectation

Time Limit: 4 Seconds Memory Limit: 262144 KB

Edward has a set of n integers {a1, a2,…,an}. He randomly picks a nonempty subset {x1, x2,…,xm} (each nonempty subset has equal probability to be picked), and would like to know the expectation of [gcd(x1, x2,…,xm)]k.

Note that gcd(x1, x2,…,xm) is the greatest common divisor of {x1, x2,…,xm}.

Input

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

The first line contains two integers n, k (1 ≤ n, k ≤ 106). The second line contains n integers a1, a2,…,an (1 ≤ ai ≤ 106).

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

Output

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

Sample Input

1
5 1
1 2 3 4 5

Sample Output

42

题意：

从N个数中取任意多个数进行一下gcd，然后求情况的gcd的k此方的和是多少。

思路：

我们能通过枚举gcd的值进行计算，那么我们只要快速的知道有多少种取法满足这个gcd就行了。我们一开始可以先预处理出F（x），表示能被x整除的数的个数，那么我们只要在这些数里面任意取一下数，那么他们的gcd一定能被x整除， 于是我们就有了F(x) gcd能被x整除的情况数。 而我们需要求的是gcd=d的情况数，那么我们能用莫比乌斯反演来做，本质就是容斥。复杂度是 max{a[i]} / 1 + max{a[i]} / 2 + max{a[i]} / 3 + ….. 就是max{a[i]} * log(max{a[i]}) 啦～

代码：

#include <cstdio>#include <iostream>#include <cmath>#include <algorithm>#include <cstring>#include <queue>#include <vector>#define rep(i,a,b) for(int i=(a);i<(b);++i)#define rrep(i,b,a) for(int i=(b);i>=(a);--i)#define clr(a,x) memset(a,(x),sizeof(a))#define ll long long#define lson l, m, rt<<1#define rson m+1,r,rt<<1|1#define mp make_pair#define ld long doubleconst int maxn = 1000000 + 5;const int mod = 998244353;ll a[maxn];int N, K;ll powerOfTwo[maxn];std::vector<int> mobius;int u[maxn];bool isprime[maxn];ll qpow(ll base,ll p){    ll result = 1;    while (p != 0) {        if (p & 1) result = result * base % mod;        base = base * base % mod;        p >>= 1;    }    return result;}void read_int(ll & x){    char ch = getchar();    while (ch < '0' || ch > '9') ch = getchar();    x = ch - '0'; ch = getchar();    while ('0' <= ch && ch <= '9') {        x = 10 * x + ch - '0';        ch = getchar();    }}void pre_init(){    powerOfTwo[0] = 1;    rep(i,1,maxn)        powerOfTwo[i] = powerOfTwo[i-1] * 2 % mod;    rep(i,1,maxn) isprime[i] = true, u[i] = 1;    rep(i,2,maxn) if (isprime[i]) {        u[i] = -1;        for(int j = i + i; j < maxn; j += i) {            isprime[j] = false;            u[j] = -u[j];        }        for(ll j = (ll)i * i; j < maxn; j += (ll)i * i) {            u[j] = 0;        }    }    rep(i,1,maxn) if (u[i] != 0) {        mobius.emplace_back(i);       // printf("%d\n",i);    }   // std::cout << (double)mobius.size() / maxn << std::endl;}ll F[maxn];void solve(){    clr(F,0);    int maxval = 0;    rep(i,0,N) {        ++F[a[i]];        if (a[i] > maxval) maxval = a[i];    }    rep(i,1,maxval+1) {        for(int j = i + i; j <= maxval; j += i)            F[i] += F[j];    }    rep(i,1,maxval+1) {        F[i] = (powerOfTwo[F[i]] - 1 + mod) % mod;    }    ll ans = 0;    rep(d,1,maxval+1) {        ll sum = 0;        for(auto x : mobius) {            if ((ll)x * d > maxval) break;            sum += (ll)u[x] * F[x * d];            if (sum >= mod) sum -= mod;            if (sum < 0) sum += mod;        }        ans = (ans + sum * qpow(d, K) % mod) % mod;    }    if (ans < 0) ans += mod;    printf("%lld\n",ans);}int main(){    //Getinput(); return 0;    #ifdef ACM        freopen("in.txt","r",stdin);        //freopen("data.in","r",stdin);        //freopen("data.out","w",stdout);    #endif // ACM    pre_init();    int T; std::cin >> T;    rep(cas,1,T+1) {        scanf("%d%d",&N,&K);        std::for_each(a,a+N,read_int);        solve();    }    return 0;}