Newcoder 148 E.Rikka with Equation（数论+莫比乌斯反演）

Description

对于一个长度为$n$的正整数序列$A$和一个正整数$m$，定义$f(A,m)$为满足同余方程 $$\sum _{i=1}^n{A}_i{x}_i\equiv 0(modm)$$ 的向量$x_i\in [0,m]$的个数。现在给出一个长度为$n$的序列$B$和一个正整数$M$，记$B$的$N=2^n-1$个非空子序列为$A_1,...,A_N$，对于每个$m\in [1,M]$，求$\sum _{i=1}^{N}f(A_i,m)$

Input

第一行一整数$T$表示用例组数，每组用例首先输入两个整数$n,M$，之后输入$n$个正整数$B_i$

$(1\le T\le 3,1\le n,M\le 10^5,1\le B_i\le 10^5)$

Output

令$w_m=\sum _{i=1}^{N}f(A_i,m)mod998244353$，输出${\oplus }_{i=1}^M{w}_i$，其中$\oplus$为异或操作

Sample Input

2 5 5 1 2 3 4 5 10 10 1 2 3 4 5 6 7 8 9 10

Sample Output

1079 933958261

Solution

首先考虑求$f(A,m)$，任意固定$x_1,...,x_{n-1}$，假设$A_i{x}_i\equiv y_i(modm)$，显然存在$k$使得 $$y_1+...+y_{n-1}=k\cdot gcd(A_1,A_2,...,A_{n-1},m)$$ 原方程成立等价于求满足 $$k\cdot gcd(A_1,...,A_{n-1},m)+A_n\cdot x_n\equiv 0(modm)$$ 的$x_n$个数，显然有$gcd(gcd(A_1,m),...,gcd(A_n,m),m)$个解，也即$gcd(A,m)$个解，故 $$f(A,m)=m^{n-1}\cdot gcd(A,m)$$ 假设$B$序列中为$d$倍数的数字有$n_d$个，记$f(d)$为$B$的子序列中满足$gcd(A,m)=d$的$m^{|A|-1}$求和，$F(d)$为$B$的子序列中满足$d|gcd(A,m)$的$m^{|A|-1}$求和，则有 $$F(d)=\sum _{d|i}f(i),F(d)=\sum _{i=1}^{n_d}{C}_{n_d}^i\cdot m^{i-1}=\frac{(m+1{)}^{n_d}-1}{m}$$ 由莫比乌斯反演 $$f(d)=\sum _{d|i}\mu (\frac{i}{d})F(i)$$ 故有 $$w_m=\sum _{i=1}^{N}f(A_i,m)=\sum _{d|m}f(d)\cdot d=\sum _{i|m}F(i)\sum _{d|i}d\cdot \mu (\frac{i}{d})=\sum _{d|m}\phi (d)\cdot F(d)$$ 预处理$\phi ,F$即可$O(mlogm)$得到$w_1,...,w_m$

Code

#include<cstdio> #include<cstring> using namespace std; typedef long long ll; #define maxn 100005 int euler[maxn],prime[maxn],res; void get_euler(int n=1e5) { euler[1]=1; res=0; for(int i=2;i<=n;i++) { if(!euler[i])euler[i]=i-1,prime[res++]=i; for(int j=0;j<res&&prime[j]*i<=n;j++) { if(i%prime[j]) euler[prime[j]*i]=euler[i]*(prime[j]-1); else { euler[prime[j]*i]=euler[i]*prime[j]; break; } } } } #define mod 998244353 int mul(int x,int y) { ll z=1ll*x*y; return z-z/mod*mod; } int add(int x,int y) { x+=y; if(x>=mod)x-=mod; return x; } int Pow(int x,ll y) { int ans=1; while(y) { if(y&1)ans=mul(ans,x); x=mul(x,x); y>>=1; } return ans; } int T,n,m,b[maxn],a[maxn]; void deal(int x) { for(int i=1;i*i<=x;i++) if(x%i==0) { a[i]++; if(i*i!=x)a[x/i]++; } } int inv[maxn],fact[maxn],sum[maxn],f[maxn]; void init(int n=1e5) { inv[1]=1; for(int i=2;i<=n;i++)inv[i]=mul(mod-mod/i,inv[mod%i]); sum[0]=1; for(int i=1;i<=n;i++)sum[i]=mul(sum[i-1],inv[i]); fact[0]=1; for(int i=1;i<=n;i++)fact[i]=mul(i,fact[i-1]); } int C(int n,int m) { return mul(fact[n],mul(sum[m],sum[n-m])); } int main() { init(); get_euler(); scanf("%d",&T); while(T--) { scanf("%d%d",&n,&m); memset(a,0,sizeof(a)); for(int i=1;i<=n;i++) { scanf("%d",&b[i]); deal(b[i]); } memset(f,0,sizeof(f)); for(int i=1;i<=m;i++) for(int j=i;j<=m;j+=i) f[j]=add(f[j],mul(euler[i],mul(add(Pow(j+1,a[i]),mod-1),inv[j]))); int ans=0; for(int i=1;i<=m;i++)ans^=f[i]; printf("%d\n",ans); } return 0; }