51Nod 1238 - 最小公倍数之和 V3（毒瘤数学+杜教筛）

题目

戳这里

推导

${\sum }_{i=1}^n{\sum }_{j=1}^{n}lcm(i,j)$

$={\sum }_{i=1}^n{\sum }_{j=1}^n\frac{ij}{gcd(i,j)}$

$={\sum }_{i=1}^n{d}^{-1}{\sum }_{i=1}^n{\sum }_{j=1}^{n}ij[gcd(i,j)==d]$

$={\sum }_{i=1}^{n}d{\sum }_{i=1}^{\lfloor \frac{n}{d}\rfloor }{\sum }_{j=1}^{\lfloor \frac{n}{d}\rfloor }ij[gcd(i,j)==1]$

$={\sum }_{i=1}^{n}d{\sum }_{i=1}^{\lfloor \frac{n}{d}\rfloor }i{\sum }_{j=1}^{\lfloor \frac{n}{d}\rfloor }j[gcd(i,j)==1]$

$={\sum }_{i=1}^{n}d(2{\sum }_{i=1}^{\lfloor \frac{n}{d}\rfloor }i{\sum }_{j=1}^{i}j[gcd(i,j)==1]-1)$

$={\sum }_{i=1}^{n}d(2{\sum }_{i=1}^{\lfloor \frac{n}{d}\rfloor }i\frac{i\phi (i)+[i==1]}{2}-1)$

$={\sum }_{i=1}^{n}d{\sum }_{i=1}^{\lfloor \frac{n}{d}\rfloor }{i}^2\phi (i)$

子问题：

求${\sum }_{i=1}^n{i}^2\phi (i)$
令$f(i)=i^2\phi (i)$
使用狄利克雷卷积，卷一个$g(i)=i^2$
那么：

${\sum }_{i=1}^n(f\ast g)(i)$

$={\sum }_{i=1}^n{\sum }_{d|i}^{}f(d)g(\frac{i}{d})$

$={\sum }_{i=1}^n{\sum }_{d|i}^{}{d}^2\phi (d)(\frac{i}{d}{)}^2$

$={\sum }_{i=1}^n{i}^2{\sum }_{d|i}^{}\phi (d)$

$={\sum }_{i=1}^n{i}^3$

$=\frac{n^2(n+1{)}^2}{4}$

又因为：

${\sum }_{i=1}^n{\sum }_{d|i}^{}{d}^2\phi (d)(\frac{i}{d}{)}^2$

$={\sum }_{i=1}^n{i}^2{\sum }_{d=1}^{\lfloor \frac{n}{i}\rfloor }{d}^2\phi (d)$

$={\sum }_{i=2}^n{i}^2{\sum }_{d=1}^{\lfloor \frac{n}{i}\rfloor }{d}^2\phi (d)+{\sum }_{i=1}^n{i}^2\phi (i)$

$=\frac{n^2(n+1{)}^2}{4}$

所以：

${\sum }_{i=1}^n{i}^2\phi (i)=\frac{n^2(n+1{)}^2}{4}-{\sum }_{i=2}^n{i}^2{\sum }_{d=1}^{\lfloor \frac{n}{i}\rfloor }{d}^2\phi (d)$

使用杜教筛将时间复杂度降到$O(n^{\frac{2}{3}})$

数学太难了QAQ

代码：

#include<cstdio>
#include<cmath>
#include<cstring>
#include<map>
#include<algorithm>

#define maxn 5000000
#define INF 0x3f3f3f3f
#define MOD 1000000007
#define two 500000004
#define six 166666668

using namespace std;

inline long long getint()
{
	long long num=0,flag=1;char c;
	while((c=getchar())<'0'||c>'9')if(c=='-')flag=-1;
	while(c>='0'&&c<='9')num=num*10+c-48,c=getchar();
	return num*flag;
}

long long n;
bool not_prime[maxn+5];
int prime[maxn+5],cnt;
long long phi[maxn+5];
map<long long,long long>M;

inline void init()
{
	phi[1]=1;
	for(int i=2;i<=maxn;i++)
	{
		if(!not_prime[i])prime[++cnt]=i,phi[i]=i-1;
		for(int j=1;j<=cnt&&i*prime[j]<=maxn;j++)
		{
			not_prime[i*prime[j]]=1;
			if(i%prime[j])phi[i*prime[j]]=phi[i]*phi[prime[j]];
			else{phi[i*prime[j]]=phi[i]*prime[j];break;}
		}
	}
	for(int i=1;i<=maxn;i++)(phi[i]*=1ll*i*i%MOD)%=MOD;
	for(int i=1;i<=maxn;i++)(phi[i]+=phi[i-1])%=MOD;
}

inline long long getsqr(long long x)
{return x%MOD*((x+1)%MOD)%MOD*((2*x+1)%MOD)%MOD*six%MOD;}

inline long long solve(long long x)
{
	if(x<=maxn)return phi[x];
	if(M.count(x))return M[x];
	long long sum=x%MOD*((x+1)%MOD)%MOD*two%MOD;
	(sum*=sum)%=MOD;
	for(long long i=2,j;i<=x;i=j+1)
	{
		j=x/(x/i);
		(sum-=(getsqr(j)-getsqr(i-1))%MOD*solve(x/i)%MOD)%=MOD;
		(sum+=MOD)%=MOD;
	}
	return M[x]=sum;
}

int main()
{
	init();
	n=getint();
	long long sum=0;
	for(long long i=1,j;i<=n;i=j+1)
	{
		j=n/(n/i);
		(sum+=1ll*(j+i)%MOD*(j-i+1)%MOD*two%MOD*solve(n/i)%MOD)%=MOD;
	}
	printf("%lld\n",sum);
}