版权声明:本文为博主原创文章,未经博主允许必须转载。 https://blog.csdn.net/C20181220_xiang_m_y/article/details/85038010
题目描述:
求
i=1∑nj=1∑nσ1(i∗j)
n≤109,mod 109+7
题目分析:
有个结论:
σk(i∗j)=x∣i∑y∣j∑[(x,y)==1] xk∗jk/yk
记
σk(n)为
n的约数的
k次幂之和,质因数分解
n,
n=p1a1p2a2...ptat
σk(n)=(1+p11k+p12k+⋯+p1a1k) ∗(1+p21k+p22k+⋯+p2a2k) … ∗(1+pt1k+pt2k+⋯+ptatk)
考虑一个质因子
p,它对
σk(n)的贡献为
1+pk+p2k+...+pak
那么看
σk(i∗j),考虑它们的公共质因子
p,在i中有a次,在j中有b次,造成的贡献就可以表示为
x=0∑ay=0∑b[(px,py)==1] pxk∗pbk/pyk
根据乘法原理,有
σk(i∗j)=x1=0∑a1y1=0∑b1[(px1,py1)==1] px1k∗pb1k/py1k ∗x2=0∑a2y2=0∑b2[(px2,py2)==1] px2k∗pb2k/py2k… ∗xt=0∑atyt=0∑bt[(pxt,pyt)==1] pxtk∗pbtk/pytk
只有当所有
(px,py)==1都满足时才会对答案造成贡献,所以可变形为:
σk(i∗j)=x∣i∑y∣j∑[(x,y)==1] xk∗jk/yk
ok,正文开始
i=1∑nj=1∑nσ1(i∗j)=i=1∑nj=1∑nx∣i∑y∣j∑[(x,y)==1] x∗j/y=k=1∑nμ(k)x=1∑kn⌊kxn⌋y=1∑knj=1∑kynxk∗kyj/ky=k=1∑nμ(k)kx=1∑kn⌊kxn⌋xj=1∑knj⌊kjn⌋=k=1∑nμ(k)k⎝⎜⎛x=1∑kn⌊kxn⌋x⎠⎟⎞2=k=1∑nμ(k)k⎝⎜⎛i=1∑knσ1(i)⎠⎟⎞2
线性筛
106,
μ(k)k标准的杜教筛,
σ1(n)的前缀和可以n
分块暴算
时间复杂度
O(n32)
- 其实有种更简单清晰的推法,第一步先不要把x,y提到最前面
=i=1∑nj=1∑nx∣i∑y∣j∑[(x,y)==1] x∗j/y=k=1∑nμ(k)i=1∑knj=1∑knx∣i∑y∣i∑xk∗jk/yk=k=1∑nμ(k)ki=1∑knj=1∑knx∣i∑xy∣j∑yj=k=1∑nμ(k)ki=1∑knj=1∑knσ1(i)∗σ1(j)=k=1∑nμ(k)k⎝⎛i=1∑knσ1(i)⎠⎞2
#include<cstdio>
#include<map>
using namespace std;
const int N = 1000000, mod = 1e9+7, inv2 = 5e8+4;
int p[N/10],sm[N+5],sd[N+5],mp[N+5];
bool v[N+5];
map<int,int>F;
void Prime()
{
sm[1]=sd[1]=1;int cnt=0;
for(int i=2;i<=N;i++)
{
if(!v[i]) p[++cnt]=i,sm[i]=-i,sd[i]=mp[i]=i+1;
for(int j=1,k;j<=cnt&&(k=p[j]*i)<=N;j++)
{
v[k]=1;
if(i%p[j]==0) {sm[k]=0,mp[k]=mp[i]*p[j]+1,sd[k]=sd[i]/mp[i]*mp[k];break;}
sm[k]=sm[i]*sm[p[j]];
sd[k]=sd[i]*sd[p[j]];
mp[k]=p[j]+1;
}
}
for(int i=1;i<=N;i++) sm[i]=(sm[i]+sm[i-1])%mod,sd[i]=(sd[i]+sd[i-1])%mod;
}
int Summu(int n)
{
if(n<=N) return sm[n];
if(F.count(n)) return F[n];
int ret=1;
for(int i=2,j;i<=n;i=j+1)
{
j=n/(n/i);
ret=(ret-1ll*(i+j)*(j-i+1)%mod*inv2%mod*Summu(n/i))%mod;
}
return F[n]=ret;
}
int calc(int n)
{
if(n<=N) return sd[n];
int ret=0;
for(int i=1,j;i<=n;i=j+1)
{
j=n/(n/i);
ret=(ret-1ll*(i+j)*(j-i+1)%mod*inv2%mod*(n/i))%mod;
}
return ret;
}
int main()
{
Prime();
int n,ans=0;
scanf("%d",&n);
for(int i=1,j,tmp,pre=0,tt;i<=n;i=j+1)
{
j=n/(n/i);
ans=(ans+1ll*((tmp=Summu(j))-pre)*(tt=calc(n/i))%mod*tt)%mod;
pre=tmp;
}
printf("%d",(ans+mod)%mod);
}