ACM学习历程—ZOJ 3868 GCD Expectation（莫比乌斯 || 容斥原理）

Description

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

15 11 2 3 4 5

Sample Output

这个题目和之前做的莫比乌斯很像。

首先d | gcd的情况的方案数比较好做。

这个很容易想到是2^num[d]-1。其中num[d]表示数列里面d的倍数的个数。

然后就是无脑莫比乌斯了，

不过第一次我把快速幂放在了第二层for循环里面T了一发，导致重复计算了quickPow(d, k)。

代码：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <set>
#include <queue>
#include <vector>
#define LL long long
#define MOD 998244353

using namespace std;

const int maxN = 1e6+10;
int n, k, len;
int pow2[maxN], num[maxN];
int prime[maxN], u[maxN];
bool vis[maxN];

//莫比乌斯反演
//F(n)和f(n)是定义在非负整数集合上的两个函数
//并且满足条件F(n) = sum(f(d)) (d|n)
//那么f(n) = sum(u(d)*F(n/d)) (d|n)
//case 1: if d = 1, u(d) = 1
//case 2: if d = p1*p2...pk, (pi均为互异素数), u(d) = (-1)^k
//case 3: else, u(d) = 0;
//性质1: sum(u(d)) (d|n) = 1 (n == 1) or 0 (n > 1)
//性质2: sum(u(d)/d) (d|n) = phi(n)/n
//另一种形式:F(d) = sum(f(n)) (d|n) => f(d) = sum(u(n/d)*F(n)) (d|n)
//线性筛选求莫比乌斯反演函数代码
void mobius()
{
    memset(vis, false,sizeof(vis));
    u[1] = 1;
    int cnt = 0;
    for(int i = 2; i < maxN; i++)
    {
        if(!vis[i])
        {
            prime[cnt++] = i;
            u[i] = -1;
        }
        for(int j = 0; j < cnt && i*prime[j] < maxN; j++)
        {
            vis[i*prime[j]] = true;
            if(i%prime[j])
                u[i*prime[j]] = -u[i];
            else
            {
                u[i*prime[j]] = 0;
                break;
            }
        }
    }
}

void init()
{
    mobius();
    pow2[0] = 1;
    for (int i = 1; i < maxN; ++i)
        pow2[i] = 2*pow2[i-1]%MOD;
}

void input()
{
    int tmp;
    len = 0;
    scanf("%d%d", &n, &k);
    memset(num, 0, sizeof(num));
    for (int i = 0; i < n; ++i)
    {
        scanf("%d", &tmp);
        len = max(len, tmp);
        num[tmp]++;
    }
    for (int i = 1; i <= len; ++i)
    {
        for (int j = i*2; j <= len; j += i)
            num[i] += num[j];
    }
}

//快速幂m^n
LL quickPow(LL x, int n)
{
    if (n == 0)
        return 1;
    if (x == 0)
        return 0;
    LL a = 1;
    while (n)
    {
        a *= n&1 ? x : 1;
        a %= MOD;
        n >>= 1 ;
        x *= x;
        x %= MOD;
    }
    return a;
}

void work()
{
    int ans = 0, val;
    for (int i = 1; i <= len; ++i)
    {
        val = quickPow(i, k);
        for (int j = i; j <= len; j += i)
        {
            ans += ((LL)u[j/i]*val%MOD)*(pow2[num[j]]-1)%MOD;
            ans = (ans%MOD+MOD)%MOD;
        }
    }
    printf("%d
", ans);
}

int main()
{
    //freopen("test.in", "r", stdin);
    init();
    int T;
    scanf("%d", &T);
    for (int times = 1; times <= T; ++times)
    {
        input();
        work();
    }
    return 0;
}