题面

对于 $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
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#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;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161