逆天题无链接

题目描述

题解

我们将集合内元素取模 $n$ 过后，可以发现这是 $m$ 组 $0\sim n-1$ 的可重集。我们设 $F={\prod }_{i=0}^{n-1}(1+x^i)$（这里以及后面的多项式乘法都默认是长度为 $n$ 的循环卷积），那么有 $f_i=[x^i]F^m$。

再看这个平方，相当于是找一个子集加到 $i$，再找另一个子集减到 $0$。所以我们再设 $G={\prod }_{i=0}^{n-1}(1+x^{-i\mathrm{mod}n})$，那么答案即为 $[x^0](FG{)}^m$。

由于把 $0\sim n-1$ 的集合在 $\mathrm{mod}n$ 意义下取相反数仍然是 $0\sim n-1$，所以可以很快得到 F=G，最终得到答案的表达式：
$$Ans=[x^0]F^{2m}=[x^0]{\left(\prod _{i=0}^{n-1}(1+x^i)\right)}^{2m}$$
因为涉及到循环卷积，为了正常地进行数学推导，我们仍然要硬着头皮使用复数单位根（在模意义下单位根可能不存在，这意味着最终式子里如果带单位根则无法计算）：
对 $F$ 进行DFT，得到每处的点值
$$\begin{gathered} \widehat{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} \end{gathered}$$其中 $d=\gcd (i,n)$。
对于里面的式子，我们会发现它很像因式分解 $z^n-1={\prod }_{k=0}^{n-1}(z-{\omega }_n^k)$，我们带入 $z=-1$ 可得 $\widehat{f_i}={\left((-1{)}^{n/d}\left((-1{)}^{n/d}-1\right)\right)}^{2dm}={\left(2(n/d\mathrm{mod}2)\right)}^{2dm}$。

然后做IDFT，只需要求第0项：
$$Ans=\frac{1}{n}\sum _{i=0}^{n-1}\widehat{f_i}{\omega }_n^0=\frac{1}{n}\sum _{i=0}^{n-1}\widehat{f_i}$$从上面的推导我们可以得到一个重要结论，即 $d$ 相同的 $i$ 对应的 $\widehat{f_i}$ 是一样的，所以我们可以直接枚举 $d$ 计算答案：
$$Ans=\frac{1}{n}\sum _{d|n}\widehat{f_d}\phi (n/d)$$用 Pollard-Rho 求出 $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;
}