逆天题无链接
题目描述
题解
我们将集合内元素取模 n n n 过后,可以发现这是 m m m 组 0 ∼ n − 1 0\sim n-1 0∼n−1 的可重集。我们设 F = ∏ i = 0 n − 1 ( 1 + x i ) F=\prod_{i=0}^{n-1}(1+x^i) F=∏i=0n−1(1+xi)(这里以及后面的多项式乘法都默认是长度为 n n n 的循环卷积),那么有 f i = [ x i ] F m f_i=[x^i]F^m fi=[xi]Fm。
再看这个平方,相当于是找一个子集加到 i i i,再找另一个子集减到 0 0 0。所以我们再设 G = ∏ i = 0 n − 1 ( 1 + x − i m o d n ) G=\prod_{i=0}^{n-1}(1+x^{-i\bmod n}) G=∏i=0n−1(1+x−imodn),那么答案即为 [ x 0 ] ( F G ) m [x^0](FG)^m [x0](FG)m。
由于把
0
∼
n
−
1
0\sim n-1
0∼n−1 的集合在
m
o
d
n
\bmod \,n
modn 意义下取相反数仍然是
0
∼
n
−
1
0\sim n-1
0∼n−1,所以可以很快得到 F=G,最终得到答案的表达式:
A
n
s
=
[
x
0
]
F
2
m
=
[
x
0
]
(
∏
i
=
0
n
−
1
(
1
+
x
i
)
)
2
m
Ans=[x^0]F^{2m}=[x^0]\left(\prod_{i=0}^{n-1}(1+x^i)\right)^{2m}
Ans=[x0]F2m=[x0](i=0∏n−1(1+xi))2m
因为涉及到循环卷积,为了正常地进行数学推导,我们仍然要硬着头皮使用复数单位根(在模意义下单位根可能不存在,这意味着最终式子里如果带单位根则无法计算):
对
F
F
F 进行DFT,得到每处的点值
f
i
^
=
(
∏
k
=
0
n
−
1
(
1
+
ω
n
i
k
)
)
2
m
=
(
∏
k
=
0
n
/
d
−
1
(
1
+
ω
n
/
d
k
d
)
)
2
d
m
\hat{f_i}=\left(\prod_{k=0}^{n-1}(1+\omega_n^{ik})\right)^{2m}\\ =\left(\prod_{k=0}^{n/d-1}(1+\omega_{n/d}^{kd})\right)^{2dm}\\
fi^=(k=0∏n−1(1+ωnik))2m=⎝
⎛k=0∏n/d−1(1+ωn/dkd)⎠
⎞2dm其中
d
=
gcd
(
i
,
n
)
d=\gcd(i,n)
d=gcd(i,n)。
对于里面的式子,我们会发现它很像因式分解
z
n
−
1
=
∏
k
=
0
n
−
1
(
z
−
ω
n
k
)
z^n-1=\prod_{k=0}^{n-1}(z-\omega_{n}^{k})
zn−1=∏k=0n−1(z−ωnk),我们带入
z
=
−
1
z=-1
z=−1 可得
f
i
^
=
(
(
−
1
)
n
/
d
(
(
−
1
)
n
/
d
−
1
)
)
2
d
m
=
(
2
(
n
/
d
m
o
d
2
)
)
2
d
m
\hat{f_i}=\left((-1)^{n/d}\left((-1)^{n/d}-1\right)\right)^{2dm}=\left(2(n/d\bmod 2)\right)^{2dm}
fi^=((−1)n/d((−1)n/d−1))2dm=(2(n/dmod2))2dm。
然后做IDFT,只需要求第0项:
A
n
s
=
1
n
∑
i
=
0
n
−
1
f
i
^
ω
n
0
=
1
n
∑
i
=
0
n
−
1
f
i
^
Ans=\frac{1}{n}\sum_{i=0}^{n-1}\hat{f_i}\omega_n^0=\frac{1}{n}\sum_{i=0}^{n-1}\hat{f_i}
Ans=n1i=0∑n−1fi^ωn0=n1i=0∑n−1fi^从上面的推导我们可以得到一个重要结论,即
d
d
d 相同的
i
i
i 对应的
f
i
^
\hat{f_i}
fi^ 是一样的,所以我们可以直接枚举
d
d
d 计算答案:
A
n
s
=
1
n
∑
d
∣
n
f
d
^
φ
(
n
/
d
)
Ans=\frac{1}{n}\sum_{d|n}\hat{f_d}\varphi(n/d)
Ans=n1d∣n∑fd^φ(n/d)用 Pollard-Rho 求出
n
n
n 的唯一分解后,再枚举因子计算即可。欧拉函数值可以在枚举的同时顺便求得。
复杂度不会算,反正跑得出奇地快就对了。
代码
#include<bits/stdc++.h>//JZM yyds!!
#define ll long long
#define lll __int128
#define uns unsigned
#define fi first
#define se second
#define IF (it->fi)
#define IS (it->se)
#define END putchar('\n')
#define lowbit(x) ((x)&-(x))
#define inline jzmyyds
using namespace std;
const int MAXN=1e6+5;
const ll INF=1e17;
ll read(){
ll x=0;bool f=1;char s=getchar();
while((s<'0'||s>'9')&&s>0){if(s=='-')f^=1;s=getchar();}
while(s>='0'&&s<='9')x=(x<<1)+(x<<3)+(s^48),s=getchar();
return f?x:-x;
}
int ptf[50],lpt;
void print(ll x,char c='\n'){
if(x<0)putchar('-'),x=-x;
ptf[lpt=1]=x%10;
while(x>9)x/=10,ptf[++lpt]=x%10;
while(lpt>0)putchar(ptf[lpt--]^48);
if(c>0)putchar(c);
}
const ll MOD=998244353;
ll ksm(ll a,ll b,const ll&mo=MOD){
ll res=1;
for(;b;b>>=1,a=a*a%mo)if(b&1)res=res*a%mo;
return res;
}
const lll E=1;
ll ksml(ll a,ll b,const ll&mo=MOD){
ll res=1;
for(;b;b>>=1,a=E*a*a%mo)if(b&1)res=E*res*a%mo;
return res;
}
mt19937_64 rng(1145141);
bool MillerRabin(const ll&p){
if(p==2||p==3||p==5||p==7||p==11)return 1;
if((~p&1)||!(p%3)||!(p%5)||!(p%7)||!(p%11))return 0;
ll k=p-1,m=0;
while(~k&1)k>>=1,m++;
for(int cnt=35;cnt--;){
ll a=ksml(rng()%(p-3)+2,k,p);
if(a==1||a==p-1)continue;
for(int t=0,tg=0;t<m;t++){
a=E*a*a%p;
if(!tg&&a==1)return 0;
if(a==p-1)tg=1;
}if(a!=1)return 0;
}return 1;
}
ll gcd(ll a,ll b){for(;b;a%=b,swap(a,b));return a;}
ll PollardRho(const ll&n){
while(1){
ll d=1,x=rng()%(n-1)+1,y=x,c=rng()%(n-1)+1;
for(int len=1,lim=127;;len<<=1,x=y,lim=127){
bool hc=0;
for(int cnt=len;cnt--;){
y=(E*y*y+c)%n,d=E*d*(x-y+n)%n,lim--;
if(x==y||!d){hc=1;break;}
if(!lim){
d=gcd(d,n),lim=127;
if(d>1)return d;
}
}
if(hc)break;
d=gcd(d,n);
if(d>1)return d;
}
}return 1919810;
}
unordered_map<ll,int>mp;
void findz(ll n,int cnt){
if(n==1)return;
if(MillerRabin(n)){mp[n]+=cnt;return;}
ll d=PollardRho(n);int cg=0;
while(n%d==0)n/=d,cg++;
findz(d,cg*cnt),findz(n,cnt);
}
ll n,m,a[105],ans;
int tot,b[105];
ll F(ll d){
return ksm(2,d%(MOD-1)*((m<<1)%(MOD-1))%(MOD-1));
}
void dfs(int x,ll d,ll phi){
if(x>tot){(ans+=phi%MOD*F(d)%MOD)%=MOD;return;}
if(a[x]==2){dfs(x+1,d<<b[x],phi);return;}
phi*=a[x]-1;
for(int i=1;i<b[x];i++)phi*=a[x];
for(int i=0;i<=b[x];i++)
dfs(x+1,d,phi),d*=a[x],phi/=(i==b[x]-1?a[x]-1:a[x]);
}
int main()
{
freopen("ntt.in","r",stdin);
freopen("ntt.out","w",stdout);
for(int NND=read();NND--;){
n=read(),m=read(),mp.clear();
findz(n,1),tot=0;
for(auto&x:mp)a[++tot]=x.fi,b[tot]=x.se;
ans=0,dfs(1,1,1);
print(ans*ksm(n%MOD,MOD-2)%MOD);
}
return 0;
}