ZOJ3868:GCD Expectation

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 integer T 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 contains n 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 denotes E · (2ⁿ - 1) modulo 998244353.

Sample Input

15 11 2 3 4 5

Sample Output

42

对于N个数的序列，所有非空子集中，其期望是GCD的k次方

输出期望乘以（2^N-1）的值

题目中1的概率是26/31，2的概率是2/32,3,4,5的概率是1/32

期望则是42/32,所以答案为42，也就是说我们的目标是求出期望的分子部分即可

对于N的序列，肯定有2^N-1个非空子集，其中其最大的GCD不会大于原序列的max，那么我们用数组fun来记录其期望

例如题目中的，期望为1的有26个，期望为2的有2个，期望为3,4,5的都只有1个

我们可以拆分来算，首先对于1，期望为1,1的倍数有5个，那么这五个的全部非空子集为2^5-1种，得到S=(2^5-1)*1;

对于2,2的期望应该是2，但是在期望为1的时候所有的子集中，我们重复计算了2的期望，多以我们应该减去重复计算的期望数，现在2的期望应该作1算，那么对于2的倍数，有两个，2,4，其组成的非空子集有2^2-1个，所以得到S+=(2^2-1)*1

对于3,4,5同理；

#include <iostream>#include <stdio.h>#include <string.h>#include <stack>#include <queue>#include <map>#include <set>#include <vector>#include <math.h>#include <algorithm>using namespace std;#define ls 2*i#define rs 2*i+1#define up(i,x,y) for(i=x;i<=y;i++)#define down(i,x,y) for(i=x;i>=y;i--)#define mem(a,x) memset(a,x,sizeof(a))#define w(a) while(a)#define LL long longconst double pi = acos(-1.0);#define Len 1000005#define mod 998244353const LL inf = 1<<30;LL t,n,k;LL a[Len];LL two[Len],fun[Len],cnt[Len],vis[Len],maxn;LL power(LL x, LL y){    LL ans = 1;    w(y)    {        if(y&1)            ans=(ans*x)%mod;        x=(x*x)%mod;        y/=2;    }    return ans;}int main(){    LL i,j;    scanf("%lld",&t);    two[0] = 1;    up(i,1,Len-1)    two[i] = (two[i-1]*2)%mod;    w(t--)    {        mem(cnt,0);        mem(vis,0);        scanf("%lld%lld",&n,&k);        maxn = 0;        up(i,0,n-1)        {            scanf("%lld",&a[i]);            if(!vis[a[i]])            {                vis[a[i]] = 1;                cnt[a[i]] = 1;            }            else                cnt[a[i]]++;            maxn = max(maxn,a[i]);        }        fun[1] = 1;        up(i,2,maxn)        fun[i] = power(i,k);        up(i,1,maxn)        {            for(j = i+i; j<=maxn; j+=i)                fun[j]=(fun[j]-fun[i])%mod;        }        LL ans = (two[n]-1)*fun[1]%mod;        up(i,2,maxn)        {            LL cc = 0;            for(j = i; j<=maxn; j+=i)            {                if(vis[j]) cc+=cnt[j];            }            LL tem = (two[cc]-1)*fun[i]%mod;            ans = (ans+tem)%mod;        }        printf("%lld\n",(ans+mod)%mod);    }    return 0;}