ZOJ 3868 GCD Expectation DP

dp[i] 表示公约数为i时有多少种组合

先预处理一遍dp[i]这是的dp[i]表示含有公约数i或者i的倍数的组合有多少个

再倒着dp dp[i] - = Sigma(dp[j]) (j是i的倍数 2i,3i,4i.....)

结果既为 Sigma[ dp[i]*pow(i,k) ]

GCD Expectation

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

1
5 1
1 2 3 4 5

Sample Output

42

Source: The 15th Zhejiang University Programming Contest

Submit    Status

/* ***********************************************
Author        :CKboss
Created Time  :2015年04月13日 星期一 09时29分19秒
File Name     :I.cpp
************************************************ */

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <string>
#include <cmath>
#include <cstdlib>
#include <vector>
#include <queue>
#include <set>
#include <map>

using namespace std;

typedef long long int LL;
const LL mod=998244353LL;
const int maxn=2001000;

LL power(LL x,int n)
{
	LL e=1LL;
	while(n) { if(n&1) e=(e*x)%mod; x=(x*x)%mod; n/=2; }
	return e%mod;
}

int n,k,mx;
int vis[maxn];
LL two[maxn],dp[maxn];

int main()
{
    //freopen("in.txt","r",stdin);
    //freopen("out.txt","w",stdout);

	two[0]=1LL;
	for(int i=1;i<maxn;i++) { two[i]=(two[i-1]*2LL)%mod; }

	int T_T;
	scanf("%d",&T_T);
	while(T_T--)
	{
		scanf("%d%d",&n,&k);
		memset(vis,0,sizeof(vis)); memset(dp,0,sizeof(dp)); mx=0;
		for(int i=0;i<n;i++) { int x; scanf("%d",&x); mx=max(mx,x); vis[x]++; }
		dp[1]=(two[n]-1+mod)%mod;
		for(int i=2;i<=mx;i++)
		{
			int temp=0;
			for(int j=i;j<=mx;j+=i) if(vis[j]) temp+=vis[j];
			dp[i]=(two[temp]-1+mod)%mod;
		}
		for(int i=mx;i>=1;i--)
		{
			for(int j=i+i;j<=mx;j+=i)
				dp[i]=(dp[i]-dp[j]+mod)%mod;
		}
		LL sum=0;
		for(int i=1;i<=mx;i++)
			sum=(sum+(dp[i]*power(i,k))%mod)%mod;
		cout<<sum<<endl;
	}
    return 0;
}