题目描述
由于出题人懒得写背景了，题目还是简单一点好。
输入一个整数n和一个整数p，你需要求出$(sum_{i=1}^nsum_{j=1}^n ijgcd(i,j))~mod~p$，其中gcd(a,b)表示a与b的最大公约数。
刚才题面打错了，已修改
输入输出格式
输入格式：一行两个整数p、n。
输出格式：一行一个整数$(sum_{i=1}^nsum_{j=1}^n ijgcd(i,j))~mod~p$。
输入输出样例
998244353 2000
883968974
说明
对于20%的数据，$n leq 1000$。
对于30%的数据，$n leq 5000$。
对于60%的数据，$n leq 10^6$，时限1s。
对于另外20%的数据，$n leq 10^9$，时限3s。
对于最后20%的数据，$n leq 10^{10}$，时限6s。
对于100%的数据，$5 imes 10^8 leq p leq 1.1 imes 10^9$且p为质数。
$$sum_{d=1}^ndsum_{i=1}^nsum_{j=1}^nij[gcd(i,j)==d]$$
$$sum_{d=1}^nd^3sum_{i=1}^{n/d}sum_{j=1}^{n/d}ij[gcd(i,j)==1]$$
莫比乌斯反演得
$$sum_{d=1}^nd^3sum_{i=1}^{n/d}mu(i)i^2(1+2+...[frac{n}{id}])^2$$
设$sum(x)=1+2+3+...x$
令$T=id$
$$sum_{T=1}^nsum(frac{n}{T})^2sum_{d|T}d^3frac{T}{d}^2mu(frac{T}{d}$$
$$sum_{T=1}^nsum(frac{n}{T})^2T^2sum_{d|T}dmu(frac{T}{d})$$
但是后面的$$T^2sum_{d|T}dmu(frac{T}{d})$$没法O(n)求前缀
这时候就要用到O(n$frac{2}{3}$)的杜教筛
首先$$(id*mu)(i)=varphi(i)$$
$mu$和$id(x)=x$的狄利克雷卷积等于$varphi(i)$
即$$sum_{d|T}dmu(frac{T}{d})=varphi(i)$$
$$T^2sum_{d|T}dmu(frac{T}{d})=T^2varphi(T)$$
令$$f(i)=i^2varphi(i)$$
$$S(n)=sum_{i=1}^nf(i)$$
根据狄利克雷卷积：
egin{aligned} sum_{i = 1}^N (f * g)(i) & = sum_{i = 1}^N sum_{d mid i} g(d) fleft(frac i d ight) \ & = sum_{d = 1}^N g(d) sum_{1 leq i leq N, d | i} fleft(frac i d ight) \ & = sum_{d = 1}^N g(d) sum_{1 leq i leq lfloor frac n d floor} f(i) \ & = sum_{d = 1}^N g(d) Sleft(leftlfloor frac n d ight floor ight) end{aligned}
于是有$$g(1)S(n)=sum_{i=1}^n(g*f)(i)-sum_{i=2}^ng(i)S(frac{n}{i})$$
构造经可能能约掉d2的g(x),发现g(x)=x2时最好
因为$$sum_{d|i}varphi(d)=i$$
$$(g*f)(i)=sum_{d|i}f(d)g(frac{i}{d})=sum_{d|i}d^2varphi(d)frac{i}{d}^2$$ $$=i^2sum_{d|i}varphi(d)=i^3$$
所以 $$S(n)=sum_{i=1}^ni^3-sum_{i=2}^ni^2S(frac{n}{i})$$
先预处理出S(x)的前N位(N值自定，一般为maxn$frac{2}{3}$)
否则就内部分块递归
外部直接数论分块
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<cmath> 5 #include<algorithm> 6 #include<map> 7 using namespace std; 8 typedef long long lol; 9 const lol N=8000000; 10 lol Mod,tot,n; 11 lol ans,phi[N+5]; 12 lol prime[N+5],inv2,inv6,MAX; 13 bool vis[N+5]; 14 map<lol,lol> M; 15 lol qpow(lol x,lol y) 16 { 17 lol res=1; 18 while (y) 19 { 20 if (y&1) res=res*x%Mod; 21 x=x*x%Mod; 22 y/=2; 23 } 24 return res; 25 } 26 void pre() 27 {lol i,j; 28 phi[1]=1; 29 for (i=2;i<=MAX;i++) 30 { 31 if (vis[i]==0) 32 { 33 tot++; 34 prime[tot]=i; 35 phi[i]=(i-1)%Mod; 36 } 37 for (j=1;j<=tot;j++) 38 { 39 if (i*prime[j]>MAX) break; 40 vis[i*prime[j]]=1; 41 if (i%prime[j]==0) 42 { 43 phi[i*prime[j]]=phi[i]*prime[j]%Mod; 44 break; 45 } 46 else phi[i*prime[j]]=phi[i]*(prime[j]-1)%Mod; 47 } 48 } 49 for (i=1;i<=MAX;i++) 50 phi[i]=(i*i%Mod*phi[i]%Mod+phi[i-1])%Mod; 51 } 52 lol query(lol x) 53 { 54 if (x<=N) return phi[x]; 55 if (M[x]) return M[x]; 56 lol s=0,i,pos,sum,p=x%Mod; 57 s=p*(p+1)%Mod*inv2%Mod; 58 s=s*s%Mod; 59 lol as=0; 60 for (i=2;i<=x;i=pos+1) 61 { 62 pos=x/(x/i); 63 p=pos%Mod; 64 sum=(p+1)*p%Mod*(2*p+1)%Mod*inv6%Mod; 65 p=(i-1)%Mod; 66 sum-=(p+1)*p%Mod*(2*p+1)%Mod*inv6%Mod; 67 sum=(sum+Mod)%Mod; 68 sum=sum*query(x/i)%Mod; 69 as+=sum;as%=Mod; 70 } 71 return M[x]=(s-as+Mod)%Mod; 72 } 73 int main() 74 {lol i,pos,x,sum; 75 cin>>Mod>>n; 76 MAX=min(N,n); 77 pre(); 78 inv2=qpow(2,Mod-2); 79 inv6=qpow(6,Mod-2); 80 ans=0; 81 for (i=1;i<=n;i=pos+1) 82 { 83 pos=(n/(n/i)); 84 x=n/i;x%=Mod; 85 sum=(x+1)*x%Mod*inv2%Mod; 86 sum=sum*sum%Mod; 87 ans=ans+(sum*((query(pos)-query(i-1)+Mod)%Mod))%Mod; 88 ans%=Mod; 89 } 90 printf("%lld",ans); 91 }