P3768 简单的数学题
题目描述
输入一个整数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$。
输入输出样例
说明
对于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_{i=1}^nsum_{j=1}^nijgcd(i,j)\
=sum_{i=1}^nsum_{j=1}^nij sum_{d|iwedge d|j}varphi(d)\
=sum_{d=1}^nvarphi(d)sum_{d|i}sum_{d|j}ij\
=sum_{d=1}^nvarphi(d)d^2(sum_{i=1}^{lfloorfrac nd
floor}i)^2\
=sum_{d=1}^nvarphi(d)d^2S(lfloorfrac nd
floor,3)
]
其中(S)表示自然数幂和。所以只要能求出(f(n)=varphi(n)n^2)的前缀和(F(n)=sum_{i=1}^nf(i))，这个式子就能整数分块做。
考虑(f(n)=varphi(n)n^2)它的形式，(f=varphicdot id^2)，符合杜教筛处理的形式。于是构造(g(n)=n^2)，则
[(f*g)(n)=sum_{d|n}varphi(d)d^2(frac nd)^2=n^3
]
而(g,f*g)的前缀和就是(S(n,2),S(n,3))，所以这个问题就解决了。最后列一下求和式：
[F(n)=sum_{i=1}^n(f*g)(i)-sum_{i=2}^ng(i)F(lfloor frac ni
floor)\
=sum_{i=1}^ni^3-sum_{i=2}^ni^2F(lfloor frac ni
floor)
]
时间复杂度(O(n^frac 23+n^frac 12))。话说这题数据真的水，我减号写成加号都有60分。
#include<bits/stdc++.h>
#define il inline
#define co const
template<class T>T read(){
T data=0,w=1;char ch=getchar();
for(;!isdigit(ch);ch=getchar())if(ch=='-') w=-w;
for(;isdigit(ch);ch=getchar()) data=data*10+ch-'0';
return data*w;
}
template<class T>il T read(T&x) {return x=read<T>();}
typedef long long LL;
using namespace std;
co int N=4641589;
int mod,i6;
int pri[N],tot,phi[N];
int add(int a,int b){
return (a+=b)>=mod?a-mod:a;
}
int mul(int a,int b){
return (LL)a*b%mod;
}
int fpow(int a,int b){
int ans=1;
for(;b;b>>=1,a=mul(a,a))
if(b&1) ans=mul(ans,a);
return ans;
}
void init(){
i6=fpow(6,mod-2);
pri[1]=phi[1]=1;
for(int i=2;i<N;++i){
if(!pri[i]) pri[++tot]=i,phi[i]=i-1;
for(int j=1;j<=tot&&i*pri[j]<N;++j){
pri[i*pri[j]]=1,phi[i*pri[j]]=phi[i]*phi[pri[j]];
if(i%pri[j]==0){
phi[i*pri[j]]=phi[i]*pri[j];
break;
}
}
}
for(int i=2;i<N;++i) phi[i]=add(phi[i-1],mul(phi[i],mul(i,i)));
}
int sum_s2(int n){
return mul(i6,mul(n,mul(n+1,2*n+1)));
}
int sum_s3(int n){
return n&1?mul(mul(n,(n+1)/2),mul(n,(n+1)/2)):mul(mul(n/2,n+1),mul(n/2,n+1));
}
map<LL,int> sf;
int sum_f(LL n){
if(n<N) return phi[n];
if(sf.count(n)) return sf[n];
int ans=sum_s3(n%mod);
for(LL l=2,r;l<=n;l=r+1){
r=n/(n/l);
ans=add(ans,mod-mul(add(sum_s2(r%mod),mod-sum_s2((l-1)%mod)),sum_f(n/l)));
}
return sf[n]=ans;
}
int solve(LL n){
int ans=0;
for(LL l=1,r;l<=n;l=r+1){
r=n/(n/l);
ans=add(ans,mul(add(sum_f(r),mod-sum_f(l-1)),sum_s3(n/l%mod)));
}
return ans;
}
int main(){
cerr<<(sizeof(pri)+sizeof(phi))/1024.0/1024<<endl;
read(mod);
init();
printf("%d
",solve(read<LL>()));
return 0;
}