【2023ICPC网络赛I 】E. Magical Pair

当时在做洛谷U389682 最大公约数合并的时候我就想到把每个质因子分解出来然后跑高维前缀和，但是那一道题不是用这个方法，所有我也一直在思考这种做法是不是真的有用。因为昨天通过2024上海大学生程序设计竞赛I-六元组计数这道题我了解到了不少关于原根的性质，所以想着回来做去年网络赛的题目。因为我当时完全不了解原根，因此做不了这个题目，更看不懂题解，但是我现在已经大概掌握原根的知识，所以感觉做这道题还算比较轻松，而且这道题里面刚好就用到了曾经想到的质因数分解+高维前缀和，感觉十分有趣，于是写博客记录。

因为原根的题目一般都喜欢把0先处理掉，然后在处理不为0的情况。
很显然，如果想让左右两边为0，只需要满足$n|x,n|y$，方案数为$(n-1{)}^2$。

然后我的想法是先列出一个$(n-1)\ast (n-1)$的表，$1<=x<=n-1,1<=y<=n-1,a_{i,j}=i^j$。实际上，表上每一个位置都表示$n\ast (n-1)$个数，因为实际上它的行坐标可以加上若干个$n$，列坐标可以加上若干个$n-1$。然后就可以发现任意两个值相同的位置（可以一样）都恰好对应了一个解（我也不清楚如果以前没有看过题解能不能想到这一步）。比如两个坐标分别为$(x_1,y_1),(x_2,y_2)$，如果我们想形成一个解，那么就必须要让$x_1+k_{1}n=y_2+k_2(n-1),y_1+k_3(n-1)=x_2+k_{4}n$，这个根据扩展欧几里得可以知道如果想找到另外一个解，那么$k_1,k_4$都必须改变$n-1$，$k_2,k_3$都必须改变$n$，因此两个相同值的位置恰好对应一个解。

因此，我们就必须算出每个数出现的次数，答案就是这个出现次数的平方和。

根据原根的性质，和$n-1$最大公因数相同的数出现次数一样， 然后我们就把和$n-1$最大公因数相同的数全部在一个组，然后考虑组合组之间的影响（详见2024上海大学生程序设计竞赛I-六元组计数&原根知识详解。

我们就考虑每个集合的最小的数，那么只有它的因子所在的集合会出现这个数，实际上就是一个周期出现一次,$gcd(x,n-1)=y$的数有$y$个周期，所以贡献就是。

$(n-1{)}^2+{\sum }_{i|n-1}\varphi (\frac{n-1}{i})({\sum }_{j|i}\varphi (\frac{n-1}{j})j{)}^2$

我们可以预处理出$2\ast 10^7$以内的质数（实测表明n-1的所有素因子都不会超过4*10^14，所有我只用处理出这么多），然后分解完质因数就用高位前缀和把${\sum }_{j|i}\varphi (\frac{n-1}{j})j$计算出来即可。

#include<bits/stdc++.h>
#define rep(i,x,y) for(int i=x;i<=y;i++)
#define dwn(i,x,y) for(int i=x;i>=y;i--)
#define ll long long
using namespace std;
template<typename T>inline void qr(T &x){x=0;int f=0;char s=getchar();while(!isdigit(s))f|=s=='-',s=getchar();while(isdigit(s))x=x*10+s-48,s=getchar();x=f?-x:x;
}
int cc=0,buf[31];
template<typename T>inline void qw(T x){if(x<0)putchar('-'),x=-x;do{buf[++cc]=int(x%10);x/=10;}while(x);while(cc)putchar(buf[cc--]+'0');
}
const int N=2e7+10,M=2e5+10,mod=998244353;
int cnt;ll p[N];bool v[N];
ll n,ans;
struct node{ll x,y;
}a[210];int m;ll val[210];
ll phi[M],lim,dp[M],num[M];
void solve(){qr(n);if(n==2){puts("2");return;}n--;ll nn=n;ans=(n%mod)*(n%mod)%mod;m=0;rep(i,1,cnt){if(p[i]*p[i]>n)break;if(n%p[i]==0){a[++m].x=p[i];a[m].y=0;while(n%p[i]==0)a[m].y++,n/=p[i];}}if(n>1)a[++m]=(node){n,1};val[1]=1;rep(i,2,m)val[i]=val[i-1]*(a[i-1].y+1);lim=val[m]*(a[m].y+1);rep(i,0,lim-1){if(!i){phi[i]=1;num[i]=1;continue;}bool bk=0;rep(j,1,m){ll t=i/val[j]%(a[j].y+1);if(t){if(t==1)phi[i]=phi[i-val[j]]*(a[j].x-1);else phi[i]=phi[i-val[j]]*a[j].x;num[i]=num[i-val[j]]*a[j].x;break;}}}rep(i,0,lim-1)dp[i]=(phi[lim-1-i]%mod)*(num[i]%mod)%mod;rep(i,1,m)rep(j,val[i],lim-1)if(j/val[i]%(a[i].y+1))(dp[j]+=dp[j-val[i]])%=mod;rep(i,0,lim-1){dp[i]=dp[i]*dp[i]%mod;(ans+=(phi[lim-1-i]%mod)*dp[i]%mod)%=mod;}cout<<ans<<endl;
}
int main(){rep(i,2,20000000){if(!v[i])v[i]=1,p[++cnt]=i;for(int j=1;j<=cnt&&i*p[j]<=20000000ll;j++){v[i*p[j]]=1;if(i%p[j]==0)break;}}int tt;qr(tt);while(tt--)solve();return 0;
}