【NOI模拟赛】防AK题（生成函数，单位根，Pollard-Rho）

题面

对于 $r=0,1,\cdots ,n-1$，设 { 1 , 2 , ⋯ , n m } \{1,2,⋯,nm\} {1,2,⋯,nm} 中有 $f_r$ 个子集满足子集内元素之和 $\mathrm{mod}n=r$, 求出
$$\sum _{r=0}^{n-1}{f}_r^2\mathrm{mod}998244353.$$

一共 $T$ 组数据。

$T\le 20$，$n,m\le 10^{18}.$

题解

先写出 $f$ 的生成函数：
$$F={\left(\prod _{i=0}^n(1+x^i)\right)}^m\mathrm{mod}(x^n-1)$$

然后我们会发现一个性质：$f_i=f_{n-i}$ ，这个生成函数是对称的！

所以 ${\sum }_{r=0}^{n-1}{f}_r^2=[[x^0]]F^2$ ，把 $F$ 和自己循环卷积，刚好可以把每一项的平方加到常数项上。所以，我们只需要求如下函数的常数项：
$$F^{\prime }={\left(\prod _{i=0}^n(1+x^i)\right)}^{2m}\mathrm{mod}(x^n-1)$$

求循环卷积还有个方法，是用单位根，快速傅里叶变换，我们把正变换写出来：
$${\hat{f}}_k={\left(\prod _{i=0}^n(1+{\omega }_n^{ik})\right)}^{2m}$$

令 $d=\gcd (k,n)$ ，则
$${\hat{f}}_k={\left(\prod _{i=0}^n(1+{\omega }_{n/d}^{i(k/d)})\right)}^{2m}={\left(\prod _{i=0}^{n/d}(1+{\omega }_{n/d}^i)\right)}^{2dm}$$

接下来利用分圆多项式：$x^n-1=\prod (x-{\omega }_n^i)$ ，代入 $x=-1$，可以得出结论：
$${\hat{f}}_k=2^{2dm}(\frac{n}{d}\mathrm{mod}2)$$

故 ${\hat{f}}_k$ 只与 $\gcd (k,n)$ 有关，

然后再做逆变换，逆变换就是正变换的单位根取逆再除以 $n$ ，但我们求的是常数项，直接把每一项加起来再除以 $n$ 就行
$$ans=\frac{1}{n}\sum _{i=0}^{n-1}{\hat{f}}_i$$

我们用 Pollard-Rho 分解，然后枚举 $d=\gcd (i,n)$ ，计算 ${\hat{f}}_d$ 和 $\varphi (n/d)$ 乘起来。

计算 2 的幂可以预处理光速幂。

时间复杂度 $O(Td(n)\omega (n))$ ，众所周知，大整数的因数个数远小于根号。

CODE

#include<map>
#include<set>
#include<cmath>
#include<ctime>
#include<queue>
#include<stack>
#include<random>
#include<bitset>
#include<vector>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<unordered_map>
#pragma GCC optimize(2)
using namespace std;
#define MAXN 1000005
#define LL long long
#define ULL unsigned long long
#define ENDL putchar('\n')
#define DB double
#define lowbit(x) (-(x) & (x))
#define FI first
#define SE second
#define PR pair<int,int>
#define UIN unsigned int
#define SQ 317
int xchar() {
	static const int maxn = 1000000;
	static char b[maxn];
	static int pos = 0,len = 0;
	if(pos == len) pos = 0,len = fread(b,1,maxn,stdin);
	if(pos == len) return -1;
	return b[pos ++];
}
// #define getchar() xchar()
inline LL read() {
	LL f = 1,x = 0;int s = getchar();
	while(s < '0' || s > '9') {if(s<0)return -1;if(s=='-')f=-f;s = getchar();}
	while(s >= '0' && s <= '9') {x = (x<<1) + (x<<3) + (s^48);s = getchar();}
	return f*x;
}
void putpos(LL x) {if(!x)return ;putpos(x/10);putchar((x%10)^48);}
inline void putnum(LL x) {
	if(!x) {putchar('0');return ;}
	if(x<0) putchar('-'),x = -x;
	return putpos(x);
}
inline void AIput(LL x,int c) {putnum(x);putchar(c);}

const int MOD = 998244353;
int n,m,s,o,k;
LL N,M;
int pw[1<<20|5],PW[1<<20|5];
inline int pow2(LL n) {n%=(MOD-1);return PW[n>>20]*1ll*pw[n&1048575]%MOD;}
inline void MD(int &x) {if(x>=MOD)x-=MOD;}
inline int qkpow(int a,int b) {
	int res = 1;
	while(b > 0) {
		if(b & 1) res = res *1ll* a % MOD;
		a = a *1ll* a % MOD; b >>= 1;
	} return res;
}
LL safemul(LL a,LL b,LL MOD) {
	__int128 res = a; res *= b; res %= MOD;
	return res;
}
LL qkpow(LL a,LL b,LL MOD) {
	LL res = 1;
	while(b>0) {
		if(b & 1ll) res = safemul(res,a,MOD);
		a = safemul(a,a,MOD);
		b >>= 1;
	}
	return res % MOD;
}
mt19937 g(114514);
LL Rand() {return g();}
bool Miller_rabbin(LL p) {
	if(p == 2 || p == 3 || p == 5) return 1;
	if(p <= 1) return 0;
	LL nm = p-1;
	nm /= lowbit(nm);
	for(int id = 1;id <= 10;id ++) {
		LL a = Rand() % (p-1) + 1;
		if(qkpow(a,p-1,p) != 1ll) return 0;
		a = qkpow(a,nm,p);
		LL pre = p-1;
		while(a != 1) {
			pre = a;
			a = safemul(a,a,p);
		}
		if(pre != p-1) return 0;
	}
	return 1;
}
LL gcd(LL a,LL b) {return b==0 ? a:gcd(b,a%b);}
LL pd_rho(LL N) {
	if(N == 1 || Miller_rabbin(N)) return N;
	if(!(N&1)) return 2;
	while(1) {
		LL c = Rand()%(N-1)+1;
		LL a = c,b = (safemul(a,a,N) + c) % N;
		while(a != b) {
			LL dt = a-b; if(dt<0) dt=-dt;
			LL dis = gcd(dt,N);
			if(dis > 1) return dis;
			a = (safemul(a,a,N) + c) % N;
			b = (safemul(b,b,N) + c) % N;
			b = (safemul(b,b,N) + c) % N;
		}
	} return -1;
}
int F(LL x) {
	if((N/x)&1) return pow2((x%(MOD-1))*(M%(MOD-1)));
	return 0;
}
LL p[MAXN];
int cnp,w[MAXN],ans;
map<LL,int> mp;
void div(LL n) {
	if(n < 2) return ;
	LL p = pd_rho(n);
	if(p == n) {mp[n] ++; return ;}
	div(p); div(n/p);
	return ;
}
void dfs(int x,LL s,LL phi) {
	if(x > cnp) {
		MD(ans += phi%MOD*1ll*F(N/s)%MOD);
		return ;
	}
	LL pp = 1,ph = 1;
	for(int i = 0;i <= w[x];i ++) {
		dfs(x+1,s*pp,phi*ph);
		pp *= p[x];
		if(!i) ph *= (p[x]-1);
		else ph *= p[x];
	} return ;
}
int main() {
	pw[0] = 1; PW[0] = 1;
	for(int i = 1;i <= (1<<20);i ++) MD(pw[i] = pw[i-1]<<1);
	for(int i = 1;i <= (1<<20);i ++) PW[i] = PW[i-1] *1ll* pw[1<<20] % MOD;
	int T_ = read();
	while(T_ --) {
		N = read();
		M = read()<<1;
		ans = 0; cnp = 0;
		mp.clear();
		div(N);
		for(auto i = mp.begin();i != mp.end();i ++) {
			p[++ cnp] = i->FI;
			w[cnp] = i->SE;
		}
		dfs(1,1,1);
		ans = ans *1ll* qkpow(N%MOD,MOD-2) % MOD;
		AIput(ans,'\n');
	}
	return 0;
}