51Nod 1238 最小公倍数之和V3

题目传送门

分析：

现在我们需要求：

\(~~~~\sum_{i=1}^{n}\sum_{j=1}^{n}lcm(i,j)\)

\(=\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{i ~\cdot ~j}{gcd(i,j)}\)

\(=\sum_{d=1}^{n}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}i\cdot j \cdot [gcd(i,j)=1]\)

这一步推导可以看做是在枚举最小公因数，然后找除去最小公因数后的互质的两个数

又因为一个有趣的结论：

\(~~~~\sum_{i=1}^{n}i\cdot [gcd(i,n)=1]=\frac{n~\cdot ~\varphi(n)}{2}\)

对于一个数x，如果它和n互质，那么n-x也会和n互质，满足对称性

可以用类似等差数列求和的方法求和

我们继续推式子：

\(~~~~\sum_{d=1}^{n}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}i\cdot j \cdot [gcd(i,j)=1]\)

\(=\sum_{d=1}^{n}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}i^2\varphi(i)\)

前半截可以O(sqrt(n))算

我们现在要考虑如何快速求：

\(~~~~\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}i^2\varphi(i)\)

首先知道一个公式：

\(~~~~\sum_{i|n}\varphi(i)=n\)

然后：

\(~~~~\sum_{i=1}^{n}\sum_{d|i}\varphi(d)=\frac{n\cdot (n+1)}{2}\)

所以：

\(~~~~\sum_{i=1}^{n}\sum_{d|i}\varphi(d)\cdot i^2=\frac{n^2\cdot (n+1)^2}{4}\)

\(~~~~\sum_{i=1}^{n}\sum_{d|i}\varphi(d)\cdot d^2\cdot (\frac{i}{d})^2=\frac{n^2\cdot (n+1)^2}{4}\)

变换一下：

\(~~~~\sum_{i=1}^{n}\sum_{d=1}^{\lfloor\frac{n}{i}\rfloor}\varphi(d)\cdot d^2\cdot i^2=\frac{n^2\cdot (n+1)^2}{4}\)

\(~~~~\sum_{i=1}^{n}sum(\lfloor\frac{n}{i}\rfloor)\cdot i^2=\frac{n^2\cdot (n+1)^2}{4}\)

将i=1分离出来，得到：

\(~~~~sum(n)=\frac{n^2\cdot (n+1)^2}{4}-\sum_{i=2}^{n}sum(\lfloor\frac{n}{i}\rfloor)\cdot i^2\)

然后也可以O(sqrt(n))算了

整体复杂度就是O(n^(2/3))

辣鸡式子推半天，写的时候有几个点没开long long调半天。。

#include<cstdio>

#include<cstring>

#include<cmath>

#include<algorithm>

#include<vector>

#include<iostream>

#include<map>

#include<string>

#define maxn 5000005

#define INF 0x3f3f3f3f

#define MOD 1000000007

#define inv2 500000004

#define inv4 250000002

#define inv6 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;

long long pri[maxn],cnt,np[maxn],phi[maxn];

long long sum[maxn];

map<long long,long long>M;

long long ans;

inline void init()

{

	phi[1]=1;

	for(int i=2;i<maxn;i++)

	{

		if(!np[i])pri[++cnt]=i,phi[i]=i-1;

		for(int j=1;j<=cnt&&i*pri[j]<maxn;j++)

		{

			np[i*pri[j]]=1;

			if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];break;}

			phi[i*pri[j]]=phi[i]*(pri[j]-1);

		}

	}

	for(int i=1;i<maxn;i++)sum[i]=(phi[i]*i%MOD*i%MOD+sum[i-1])%MOD;

}

inline long long cal(long long x)

{return x*(x+1)%MOD*(2*x+1)%MOD*inv6%MOD;}

inline long long getans(long long x)

{

	if(x<maxn)return sum[x];

	if(M.count(x))return M[x];

	long long tmp=x%MOD;

	long long num=tmp*tmp%MOD*(tmp+1)%MOD*(tmp+1)%MOD*inv4%MOD;

	for(long long i=2,j;i<=x;i=j+1)

	{

		j=x/(x/i);

		num+=MOD-(cal(j%MOD)-cal((i-1)%MOD)+MOD)%MOD*getans(x/i)%MOD;

		num%=MOD;

	}

	return M[x]=num;

}

int main()

{

	N=getint();![](https://img2018.cnblogs.com/blog/1891964/202001/1891964-20200113163955776-377036624.jpg)

	init();

	for(long long i=1,j;i<=N;i=j+1)

	{

		j=N/(N/i);

		ans+=(i+j)%MOD*((j-i+1)%MOD)%MOD*inv2%MOD*getans(N/i)%MOD;

		ans%=MOD;

	}

	printf("%lld\n",ans);

}

巴特西