我们把答案写成
2∑ni=1∑ij=1iσ(ij)−∑ni=1σ(i2)
,对两部分分别求和。
后一部分比较简单,
σ(i2)
是积性函数,且
σ((pxq⋅p)2)=(σ((px)2)+p2x+1+p2x+2)⋅σ(q2)
,可以线性筛。
考虑前一部分,首先我们知道,
σ(mn)=∑i|n∑j|mije(gcd(ni,j))
。那么
====∑i=1ni∑j=1iσ(ij)∑i=1ni∑j=1i∑a|ia∑b|jjb∑k|gcd(a,b)μ(k)∑k=1nμ(k)∑i=1⌊nk⌋ik∑a|iak∑j=1i∑b|jjb∑k=1nμ(k)k2∑i=1⌊nk⌋g(i)∑K=1n∑k|Kμ(k)k2g(Kk)
在求出 g 函数之后,我们就可以
最后的复杂度是 O(nlogn+q) 。
#include<cstdio>
#include<algorithm>
using namespace std;
#define LL long long
const int p=1000000007,maxn=1000000;
int inv[maxn+10],prm[maxn+10],mu[maxn+10],mnp[maxn+10],mnpw[maxn+10],
f[maxn+10],g[maxn+10],h[maxn+10],
ssigma2[maxn+10],sigma[maxn+10],
n,tot;
int inc(int x,int y)
{
x+=y;
return x>=p?x-p:x;
}
int dec(int x,int y)
{
x-=y;
return x<0?x+p:x;
}
int pow(int b,int k)
{
int ret=1;
for (;k;k>>=1,b=(LL)b*b%p)
if (k&1) ret=(LL)ret*b%p;
return ret;
}
void init()
{
int x;
inv[1]=mu[1]=sigma[1]=ssigma2[1]=1;
for (int i=2;i<=maxn;i++)
{
if (!mnp[i])
{
prm[++tot]=i;
mnp[i]=mnpw[i]=i;
mu[i]=p-1;
inv[i]=pow(i,p-2);
sigma[i]=i+1;
ssigma2[i]=inc(inc((LL)i*i%p,i),1);
}
for (int j=1;j<=tot&&(LL)i*prm[j]<=maxn;j++)
{
x=i*prm[j];
inv[x]=(LL)inv[i]*inv[prm[j]]%p;
if (i%prm[j])
{
mnp[x]=mnpw[x]=prm[j];
mu[x]=dec(0,mu[i]);
sigma[x]=(LL)sigma[i]*sigma[prm[j]]%p;
ssigma2[x]=(LL)ssigma2[i]*ssigma2[prm[j]]%p;
}
else
{
mnp[x]=prm[j];
mnpw[x]=mnpw[i]*prm[j];
mu[x]=0;
sigma[x]=inc(sigma[i],(LL)prm[j]*mnpw[i]%p*sigma[i/mnpw[i]]%p);
ssigma2[x]=inc(ssigma2[i],(LL)inc((LL)mnpw[i]*mnpw[i]%p*mnp[i]%p,(LL)mnpw[i]*mnpw[i]%p*mnp[i]%p*mnp[i]%p)*ssigma2[i/mnpw[i]]%p);
break;
}
}
}
for (int i=1;i<=maxn;i++)
{
ssigma2[i]=inc(ssigma2[i-1],(LL)i*ssigma2[i]%p);
h[i]=inc(h[i-1],sigma[i]);
g[i]=(LL)i*h[i]%p*sigma[i]%p;
for (int j=i;j<=maxn;j+=i) f[j]=inc(f[j],(LL)g[i]*mu[j/i]%p*(j/i)%p*(j/i)%p);
f[i]=inc(f[i-1],f[i]);
}
}
int main()
{
init();
int T;
scanf("%d",&T);
for (int K=1;K<=T;K++)
{
scanf("%d",&n);
printf("Case #%d: %d\n",K,dec(inc(f[n],f[n]),ssigma2[n]));
}
}