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 ≤ $10^{6}$). The second line contains n integers a1, a2,…,an (1 ≤ ai ≤ $10^6$).

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 double
const 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;
}