题目大意
求
∑i=1n∑j=1nijgcd(i,j)
对质数p取膜的结果
n<=10^10 5*10^8<=p<=1.1*10^9
传送门https://www.luogu.org/problem/show?pid=3768
分析
看到这种两层求和带gcd的式子 可以考虑直接用莫比乌斯反演求解。
但注意到i j的上界都是n,所以用
ϕ
推导公式也可以解决这题。(感觉可能会简便一些)
∑i=1n∑j=1nijgcd(i,j)
=∑d=1nd3∑i=1⌊nd⌋∑j=1⌊nd⌋ij[i⊥j]
=∑d=1nd3G(⌊nd⌋)
右边的式子是关于 ⌊nd⌋ 的 可以看成一个函数
如果能够快速计算出G函数 那么这个式子可以分块快速求出
G(n)=∑i=1n∑j=1nij[i⊥j]
=∑i=1ni∑j=1nj[i⊥j]
=2∑i=1ni∑j=1ij[i⊥j]−1
=2∑i=1ni(iϕ(i)+[i==1]2)−1
=∑i=1n(i2ϕ(i)+[i==1])−1
=∑i=1ni2ϕ(i)
(嗯 看起来很简单吧..)
G(n)显然是一个积性函数的前缀和
如果预处理出G函数 那么可以设计一个O(n)的算法 可以通过60%的数据。
记
G(n)=∑ni=1f(i)
考虑将
ϕ
通过狄利克雷卷积消去
已知
ϕ∗Id0=Id1
假设
f∗h=T
T(n)=∑d|nf(d)h(nd)=∑d|nd2ϕ(d)h(nd)=∑d|n(d2h(nd))ϕ(d)
当
h=Id2
时
d2可以与h(nd)配对成n2后面的∑phi(d)也可以合并成一个n
卷积即可得到
T=Id3
那么就可以用杜教筛快速计算出G(n)了
G(n)=n3−∑i=2ni2G(⌊ni⌋)
总时间复杂度 O(n23)
#include<stdio.h>
#include<math.h>
#include<string.h>
#include<stdlib.h>
#include<map>
#include<set>
#include<algorithm>
#include<queue>
#include<time.h>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fd(i,a,b) for(int i=a;i>=b;i--)
#define fe(x) for(int ii=be[x];ii;ii=e[ii].ne)
using namespace std;
const int mn=20000000;
int p,pri[mn/10],pt,mm;
long long n,p2,p6,ans,g[mn+20],mg[100000];
bool vis[100000],isp[mn+20];
inline long long qpow(int a,int n){
long long x=a,ans=1;
while (n){
if (n&1) ans=ans*x%p;
x=x*x%p;
n>>=1;
}
return ans;
}
inline long long sum2(long long n){
n%=p;
return (n+n+1)%p*(n+1)%p*n%p*p6%p;
}
inline long long sum3(long long n){
n%=p;
n=n*(n+1)%p*p2%p;
return n*n%p;
}
inline long long G(long long x){
if (x<=mm) return g[x];
int ms=n/x;
if (vis[ms]) return mg[ms];
vis[ms]=1;
long long i=2,j,te;
long long &ans=mg[ms];
ans=sum3(x);
while (i<=x){
te=x/i;
j=x/te;
ans-=(sum2(j)-sum2(i-1))*G(te)%p;
i=j+1;
}
return ans%=p;
}
int main(){
scanf("%d%lld",&p,&n);
mm=mn;if (n<mm) mm=n;
g[1]=1;
fo(i,2,mm){
if (!isp[i]){
g[i]=i-1;
pri[++pt]=i;
}
fo(j,1,pt){
long long te=i*(long long)pri[j];
if (te>mn) break;
isp[te]=1;
if (i%pri[j]==0){
g[te]=g[i]*pri[j];
break;
}
g[te]=g[i]*(pri[j]-1);
}
}
fo(i,2,mm) g[i]=(g[i-1]+g[i]*i%p*i)%p;
p2=qpow(2,p-2);p6=qpow(6,p-2);
long long i=1,j,te;
while (i<=n){
te=n/i;
j=n/te;
ans+=(sum3(j)-sum3(i-1))*G(te)%p;
i=j+1;
}
ans=(ans%p+p)%p;
printf("%lld\n",ans);
return 0;
}