n<=1e9 给定一个模数
现场有13人AC..我..好菜啊
对原式化简
∑d=1nd∑i=1n∑j=1i∑k=1i[gcd(i,j,k)==d]∑d=1nd∑i=1n∑j=1i∑k=1i[gcd(i,j,k)==d]
∑d=1nd∑i=1⌊nd⌋∑j=1i∑k=1i[gcd(i,j,k)==1]∑d=1nd∑i=1⌊nd⌋∑j=1i∑k=1i[gcd(i,j,k)==1]
∑d=1nd∑i=1⌊nd⌋∑j=1i∑k=1iε(gcd(i,j,k))∑d=1nd∑i=1⌊nd⌋∑j=1i∑k=1iε(gcd(i,j,k))
∑d=1nd∑i=1⌊nd⌋∑j=1i∑k=1i∑d′|gcd(i,j,k)μ(d′)∑d=1nd∑i=1⌊nd⌋∑j=1i∑k=1i∑d′|gcd(i,j,k)μ(d′)
∑d=1nd∑d′=1⌊nd⌋μ(d′)∑i=1⌊ndd′⌋∑j=1i∑k=1i1∑d=1nd∑d′=1⌊nd⌋μ(d′)∑i=1⌊ndd′⌋∑j=1i∑k=1i1
设D为d∗d′d∗d′观察到后面为∑i=1⌊ndd′⌋i2∑i=1⌊ndd′⌋i2
将前半部分替换卷积 后半部分替换公式
∑D=1n∑d′|DDd′μ(d′)⌊ndd′⌋∗(⌊ndd′⌋+1)∗(2∗⌊ndd′⌋)6∑D=1n∑d′|DDd′μ(d′)⌊ndd′⌋∗(⌊ndd′⌋+1)∗(2∗⌊ndd′⌋)6
∑D=1nφ(D)⌊ndd′⌋∗(⌊ndd′⌋+1)∗(2∗⌊ndd′⌋)6∑D=1nφ(D)⌊ndd′⌋∗(⌊ndd′⌋+1)∗(2∗⌊ndd′⌋)6
因为后半部分可以考虑到分块计算所以只需要前半部分前缀和即可 那么杜教筛 预处理欧拉函数的前 2/3使得复杂度降为n23n23
顺带学习了发杜教筛 下面乘号部分代表卷积符号
考虑两个积性函数的卷积(f∗g)(n)=h(f∗g)(n)=h如何求积性函数前缀和
考虑φ∗1=idφ∗1=id
顾前缀和可以表示为∑i=1ni=∑i=1n∑d|i1(d)∗φ(nd)∑i=1ni=∑i=1n∑d|i1(d)∗φ(nd)常见套路
∑i=1ni=∑d=1n∑i=1⌊nd⌋φ(i)∑i=1ni=∑d=1n∑i=1⌊nd⌋φ(i)
设H[i]H[i]表示φφ的前缀和
那么显然n∗(n+1)2=∑d=1nH(nd)n∗(n+1)2=∑d=1nH(nd)
n∗(n+1)2−∑d=2nH(nd)=H(n)n∗(n+1)2−∑d=2nH(nd)=H(n)好了剩下部分可以递归搞了
#include<cstdio>
#include<algorithm>
#define ll long long
#define N 1100000
using namespace std;
int prime[N/10],tot,mod,inv6,n;
bool not_prime[N];ll s[N],ans,phi[N];
inline int ksm(int x,int t){
int tmp=1;for (;t;x=(ll)x*x%mod,t>>=1) if (t&1) tmp=(ll)tmp*x%mod;return tmp;
}
inline ll calc_phi(ll x){
if (x<=1e6) return phi[x];if (s[n/x]) return s[n/x];
ll tmp=x*x+x>>1;int last=1;
for (int i=2;i<=x;i=last+1){
last=x/(x/i);tmp-=calc_phi(x/i)*(last-i+1);
}return s[n/x]=tmp;
}
inline int calc(int x){
return (ll)x*(x+1)%mod*(x<<1|1)%mod*inv6%mod;
}
int main(){
freopen("sum.in","r",stdin);
for (int i=2;i<=1e6;++i){
if (!not_prime[i]) prime[++tot]=i,phi[i]=i-1;
for (int j=1;prime[j]*i<=1e6;++j){
not_prime[i*prime[j]]=1;
if (i%prime[j]==0){
phi[i*prime[j]]=prime[j]*phi[i];break;
}else phi[i*prime[j]]=phi[i]*phi[prime[j]];
}
}scanf("%d%d",&n,&mod);int last=1;inv6=ksm(6,mod-2);
//for (int i=1;i<=10;++i) printf("%d ",phi[i]);puts("");
phi[1]=1;for (int i=2;i<=1e6;++i) phi[i]+=phi[i-1];
for (int i=1;i<=n;i=last+1){
last=n/(n/i);ll ans1=(calc_phi(last)-calc_phi(i-1))%mod,ans2=calc(n/i);
ans+=ans1*ans2%mod;ans%=mod;
// printf("%lld %lld %lld %d\n",ans,ans1,ans2,i);
}printf("%lld\n",ans);
return 0;
}