传送门
惯用套路
[f(d)=dsum_{i=1}^{n}sum_{j=1}^{n}ij[gcd(i,j)==d]\
]
然后简单的莫比乌斯反演一下，得到
[f(d)=d^3sum_{T=1}^{n/d}mu(T)T^2sum_{i=1}^{n/Td}sum_{j=1}^{n/Td}ij\
]
然后知道答案式
[ans=sum_{d=1}^{n}f(d)=sum_{d=1}^{n}d^3sum_{T=1}^{n/d}mu(T)T^2sum_{i=1}^{n/Td}sum_{j=1}^{n/Td}ij\
]
这个时候去算，复杂度是(O(n))的，无法通过
然后设(k=Td)
[ans=sum_{k=1}^{n}sum_{i=1}^{n/k}sum_{j=1}^{n/k}ijsum_{d|k}mu(frac{k}{d})(frac{k}{d})^2d^3\
ans=sum_{k=1}^{n}sum_{i=1}^{n/k}sum_{j=1}^{n/k}ijk^2sum_{d|k}mu(frac{k}{d})d\
]
然后现在只要能够筛出(k^2sum_{d|k}mu(frac{k}{d})d)，就可以(O(sqrt{n}))解决了
我们设(g(k)=k^2sum_{d|k}mu(frac{k}{d})d)
现在考虑杜教筛，我们有一个很经典的式子
[mu*id=varphi\
g(k)=k^2varphi(k)
]
则我们设(S(n)=sum_{i=1}^{n}g(i))，套上杜教筛的式子
[S(n)h(1)=sum_{i=1}^{n}(h*g)-sum_{i=2}^{n}S(lfloorfrac{n}{i}
floor)h(i)
]
发现
[h*g=sum_{d|n}h(lfloorfrac{n}{d}
floor)g(d)=sum_{d|n}h(lfloorfrac{n}{d}
floor)d^2varphi(d)\
]
考虑如何把(varphi(d))消去
因为
[sum_{d|n}varphi(d)=n
]
则
[h*g=sum_{d|n}h(lfloorfrac{n}{d}
floor)d^2n\
]
然后我们可以考虑将(h(i))设为(i^2)
则
[h*g=n^3\
S(n)=sum_{i=1}^{n}i^3-sum_{i=2}^{n}i^2S(lfloorfrac{n}{i}
floor)
]
这个就已经可以处理了
然后对于这个答案式太复杂了
[ans=sum_{k=1}^{n}sum_{i=1}^{n/k}sum_{j=1}^{n/k}ijk^2sum_{d|k}mu(frac{k}{d})d\
]
我们设(s(x)=sum_{i=1}^{x}sum_{j=1}^{x}ij)，这个是可以用等差数列解决的
然后答案式就简洁多了
[ans=sum_{k=1}^{n}s(lfloorfrac{n}{k}
floor)g(k)
]
代码：
#include<cstdio>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<tr1/unordered_map>
using namespace std;
void read(int &x)
{
char ch;bool ok;
for(ok=0,ch=getchar();!isdigit(ch);ch=getchar())if(ch=='-')ok=1;
for(x=0;isdigit(ch);x=x*10+ch-'0',ch=getchar());if(ok)x=-x;
}
#define rg register
#define ll long long
const int maxn=3e6+10;
int mod,ans,phi[maxn],pri[maxn],tot,inv2,inv6,f[maxn];
ll n;bool vis[maxn];
tr1::unordered_map<ll,int>mp;
int mul(int x,int y){return 1ll*x*y-1ll*x*y/mod*mod;}
int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
int del(int x,int y){return x-y<0?x-y+mod:x-y;}
int mi(int a,int b){
int ans=1;
while(b){
if(b&1)ans=mul(ans,a);
b>>=1,a=mul(a,a);
}
return ans;
}
void prepare()
{
phi[1]=1;
for(rg int i=2;i<=3e6;i++){
if(!vis[i])pri[++tot]=i,phi[i]=i-1;
for(rg int j=1;j<=tot&&pri[j]*i<=3e6;j++){
vis[pri[j]*i]=1;
if(!(i%pri[j])){phi[pri[j]*i]=mul(pri[j],phi[i]);break;}
else phi[pri[j]*i]=mul(phi[pri[j]],phi[i]);
}
}
for(rg int i=1;i<=3e6;i++)f[i]=mul(mul(i,i),phi[i]);
for(rg int i=1;i<=3e6;i++)f[i]=add(f[i],f[i-1]);
}
int s(int n) {return mul(mul(mul(n*2+1,n+1),n),inv6);}
int square(ll l,ll r){return del(s(r%mod),s((l-1)%mod));}
int sum(ll n){return mul(mul((n+1)%mod,n%mod),inv2);}
int sum1(ll n){return mul(mul(n%mod,(n+1)%mod),inv2);}
int get_phi(ll n){
if(n<=3e6)return f[n];
if(mp[n])return mp[n];
int ans=mul(sum1(n),sum1(n));
for(rg ll i=2,j;i<=n;i=j+1)j=n/(n/i),ans=del(ans,mul(square(i,j),get_phi(n/i)));
return mp[n]=ans;
}
int main()
{
read(mod),scanf("%lld",&n);
prepare(),inv2=mi(2,mod-2),inv6=mi(6,mod-2);
for(rg ll i=1,j,t;i<=n;i=j+1){
j=n/(n/i),t=mul(sum(n/i),sum(n/i));
ans=add(ans,mul(t,del(get_phi(j),get_phi(i-1))));
}
printf("%d
",ans);
}