题意
d(k)表示k的所有约数的和。d(6) = 1 + 2 + 3 + 6 = 12。
定义S(N) = ∑1<=i<=N ∑1<=j<=N d(i*j)。
例如:S(3) = d(1) + d(2) + d(3) + d(2) + d(4) + d(6) + d(3) + d(6) + d(9) = 59,S(1000) = 563576517282。
给出正整数N,求S(N),由于结果可能会很大,输出Mod 1000000007(10^9 + 7)的结果。
2 <= N <= 10^9
分析
设
σ(i)
表示i的约数和。
这题要用到一个结论,就是
σ(ij)=∑p|i∑q|j[(p,q)=1]pjq
。
证明的话,可以考虑每个素因子分开考虑。设 r 为素数,
那么有
ans=∑i=1n∑j=1n∑p|i∑q|j[(p,q)=1]pjq
反演一下可以得到
ans=∑d=1nμ(d)∑i=1n∑j=1n∑p|i∑q|j[d|(p,q)]pjq
=∑d=1nμ(d)∑d|p∑d|qpq∑p|i∑q|jj
=∑d=1nμ(d)∑d|p∑d|qpq⌊np⌋q⌊nq⌋(⌊nq⌋+1)2
=∑d=1nμ(d)∑d|pp⌊np⌋∑d|q⌊nq⌋(⌊nq⌋+1)2
=∑d=1nμ(d)d(∑p=1⌊nd⌋p⌊ndp⌋)(∑q=1⌊nq⌋⌊ndq⌋(⌊ndq⌋+1)2)
这已经是一个分块的式子了,考虑把下取整里面的东西换一下。
∑i=1ni⌊ni⌋
∑i=1n⌊ni⌋(⌊ni⌋+1)2
也就是说我们要考虑如何快速求这两个东西。
其中第一个式子的下取整可以看做n以内有多少个数是i的倍数,那么可以得到
∑i=1ni⌊ni⌋=∑i=1nσ(i)
对第二个式子进行变形:
∑i=1n⌊ni⌋(⌊ni⌋+1)2=∑i=1n∑j=1⌊ni⌋j=∑i=1ni⌊ni⌋
也就是说上面两个式子的值居然是相等的!
那么就有
ans=∑d=1nμ(d)d(∑i=1⌊nd⌋σ(i))2
前面那部分的前缀和可以用杜教筛来求,后面那部分就小的预处理大的分块求即可。
代码
#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<map>
using namespace std;
typedef long long LL;
const int N=2000005;
const int MOD=1000000007;
int n,prime[N],yue[N],low[N],s[N],tot;
bool not_prime[N];
map<int,int> w;
void get_prime(int n)
{
yue[1]=s[1]=1;
for (int i=2;i<=n;i++)
{
if (!not_prime[i]) prime[++tot]=i,yue[i]=i+1,s[i]=-i,low[i]=i;
for (int j=1;j<=tot&&i*prime[j]<=n;j++)
{
not_prime[i*prime[j]]=1;
if (i%prime[j]==0)
{
low[i*prime[j]]=low[i]*prime[j];
yue[i*prime[j]]=(LL)yue[i]*prime[j]%MOD+yue[i/low[i]];
yue[i*prime[j]]-=yue[i*prime[j]]>=MOD?MOD:0;
break;
}
low[i*prime[j]]=prime[j];
yue[i*prime[j]]=(LL)yue[i]*yue[prime[j]]%MOD;
s[i*prime[j]]=s[i]*s[prime[j]];
}
}
for (int i=1;i<=n;i++) s[i]+=s[i]<0?MOD:0,s[i]+=s[i-1],s[i]-=s[i]>=MOD?MOD:0;
for (int i=1;i<=n;i++) yue[i]+=yue[i-1],yue[i]-=yue[i]>=MOD?MOD:0;
}
int get_s(int n)
{
if (n<=2000000) return s[n];
if (w[n]) return w[n];
int ans=1;
for (int i=2,last;i<=n;i=last+1)
{
last=n/(n/i);
ans+=MOD-(LL)((LL)last*(last+1)/2-(LL)i*(i-1)/2)%MOD*get_s(n/i)%MOD;
ans-=ans>=MOD?MOD:0;
}
return w[n]=ans;
}
int get_yue(int n)
{
if (n<=2000000) return yue[n];
int ans=0;
for (int i=1,last;i<=n;i=last+1)
{
last=n/(n/i);
ans+=(LL)((LL)last*(last+1)/2-(LL)i*(i-1)/2)%MOD*(n/i)%MOD;
ans-=ans>=MOD?MOD:0;
}
return ans;
}
int solve(int n)
{
int ans=0;
for (int i=1,last;i<=n;i=last+1)
{
last=n/(n/i);int w=get_yue(n/i);
ans+=(LL)(get_s(last)+MOD-get_s(i-1))*w%MOD*w%MOD;
ans-=ans>=MOD?MOD:0;
}
return ans;
}
int main()
{
get_prime(2000000);
scanf("%d",&n);
printf("%d",solve(n));
return 0;
}