51nod 1238 最小公倍数之和 V3 题解

博客观赏效果更佳

题意简述

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

$n\le 10^{10}$ （简短题意，恶臭规模）

思路

总的来说，就是先把式子变一变，然后巧妙的用欧拉函数 $\phi$ 加上杜教筛和整除分块来算。

变形

首先我们知道 $\mathrm{lcm}(a,b)=\mathrm{lcm}(b,a)$

所以，当 $i\ne j$ 时候，我们计算 $j<i$ 的情况，然后 $\times 2$ 即可

对于 $i=j$ 的情况特判，这时候产生的贡献和为 $\frac{n(n+1)}{2}$

所以，答案化为

$\sum _{i=1}^n\sum _{j=1}^{i-1}\mathrm{lcm}(i,j)+\frac{n(n+1)}{2}$

我们都知道 $\mathrm{lcm}(i,j)=\frac{ij}{\gcd (i,j)}$。这个再换一下即可。

开始硬♂上

开局一公式，结论全靠推。让我们开心的推倒式子吧~

首先枚举 $gcd$，然后枚举互质的倍数 $i,j$。

原式

$=\sum _{d=1}^n\sum _{i=1}^{n/d}\sum _{j=1}^i[\gcd (i,j)=1]ijd$
$=\sum _{i=1}^n\sum _{j=1}^i[\gcd (i,j)=1]ij\sum _{d=1}^{n/i}d$ （枚举一对互质的 $i,j$ ，计算有多少个 $d$ 算到了这对 $i,j$，加起来）
设 $S_k(n)=\sum _{i=1}^n{i}^k$ 。显然， $k=1,2,3$ 时这个式子都能 $O(1)$ 算。
接下来变成
$=\sum _{i=1}^n\sum _{j=1}^i[\gcd (i,j)=1]ij{S}_1(n/i)$ （恭喜你少了一个 $d$）
考虑如何求 $\sum _{j=1}^i[\gcd (i,j)=1]j$，也就是求 $i$ 以内和 $i$ 互质的数的和。

{% fold %}

如何求 i 以内和 i 互质的数的和（会的跳过）

首先数量一共有 $\phi (i)$ 个。
然后我们打表发现，这东西和等差数列有一个很类似的性质：我们把 $[1,i]$ 中和 $i$ 互质的 $j$ 从小到大排一排，换一行，再从大到小排一排，对应位置和相等（都是 $i$）！
根据高斯当年想出来的办法，和就是（每个对应位置的和）乘以（一共有多少项），然后除二。
换到这题来，这个和就是 $\frac{1}{2}i\times \phi (i)$
{% endfold %}

它等于 $\frac{1}{2}i\times \phi (i)$。带入原式：

$=\sum _{i=1}^{n}i\times S_1(n/i)\times \frac{i\times \phi (i)}{2}$ （恭喜你又少了一个 $j$，它现在是一维了）
$=\frac{1}{2}\sum _{i=1}^n{i}^2\phi (i)\times S_1(n/i)$

$S_1(n/i)$ 显然可以整除分块，但是整除分块我们要考虑一个东西： 如何求 $i^2\phi (i)$ 的前缀和？显然，用杜教筛求一下和即可。
{% fold %}

具体怎么杜教筛 （老样子，会的跳过）

设 $f(n)=n^2\phi (n)$
关键就是要配一个好算的 $g$ ，使得 $f\times g$ 也是一个好算的函数 $h$。

由狄利克雷卷积的定义，
$h(n)=\sum _{d|n}f(d)g(\frac{n}{d})$
$=\sum _{d|n}{d}^2\phi (i)g(\frac{n}{d})$
我们发现，令 $g(n)=n^2$，就能消掉一个 $d^2$，原式继续变成：
$=\sum _{d|n}{n}^2\phi (i)$ （注意 $g(\frac{n}{d})$ 头顶上还有一个 $n$，乘出来变成 $n^2$）
=$n^2\sum _{d|n}varphi(d)=n^3$

所以说， $h(n)=f\times (g(n)=n^2)=(n^3)$

然后用杜教筛的传统艺能就行了。
{% endfold}

总复杂度非常恐怖，但是的确能通过 $n=10^{1}0$ 的数据。一定注意处处取膜

代码

{% fold %}

#include <bits/stdc++.h>
using namespace std;
namespace Flandre_Scarlet
{
	#define N   20000007ll
	#define mod 1000000007ll
	#define i6  166666668ll
	#define i2  500000004ll
	#define int long long
	#define F(i,l,r) for(int i=l;i<=r;++i)
	#define D(i,r,l) for(int i=r;i>=l;--i)
	#define Fs(i,l,r,c) for(int i=l;i<=r;c)
	#define Ds(i,r,l,c) for(int i=r;i>=l;c)
	#define MEM(x,a) memset(x,a,sizeof(x))
	#define FK(x) MEM(x,0)
	#define Tra(i,u) for(int i=G.Start(u),v=G.To(i);~i;i=G.Next(i),v=G.To(i))
	#define p_b push_back
	#define sz(a) ((int)a.size())
	#define iter(a,p) (a.begin()+p)
	int I()
	{
	    int x=0;char c=getchar();int f=1;
	    while(c<'0' or c>'9') f=(c=='-')?-1:1,c=getchar();
	    while(c>='0' and c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
	    return (x=(f==1)?x:-x);
	}
	void Rd(int cnt,...)
	{
	    va_list args; va_start(args,cnt);
	    F(i,1,cnt) {int* x=va_arg(args,int*);(*x)=I();}
	    va_end(args);
	}
	bitset<N> notp; int primes[N/5];
	int phi[N];
	void Init()
	{
		int n=2e7;
		notp[1]=1; phi[1]=1;
		int &cnt=primes[0];
		F(i,2,n)
		{
			if (!notp[i]) {primes[++cnt]=i; phi[i]=i-1;}
			for(int j=1;j<=cnt and i*primes[j]<=n;++j)
			{
				int u=primes[j];
				notp[i*u]=1;
				if (i%u==0){phi[i*u]=phi[i]*u; break;}
				else phi[i*u]=phi[i]*phi[u];
			}
		}
		F(i,1,n) phi[i]=(phi[i-1]+1ll*phi[i]*i%mod*i%mod)%mod;
	}

	int n;
	void Input()
	{
		n=I();
	}

	map<int,int> rec; 
	int S1(int n) {n%=mod; return 1ll*n%mod*(n+1)%mod*i2%mod;}
	int S2(int n) {n%=mod; return 1ll*n%mod*(n+1)%mod*(2ll*n+1)%mod*i6%mod;}
	int S3(int n) {n%=mod; return 1ll*S1(n)*S1(n)%mod;} // 我就是这里没取膜，调了两小时
	int f(int n) //杜教筛求 f(n)=n^2*phi(n) 的前缀和
	{
		if (n<=2e7) return phi[n];
		if (rec[n]) return rec[n];
		int ans=0;
		for(int l=2,r;l<=n;l=r+1)
		{
			r=n/(n/l);
			ans=(ans+f(n/l)*(S2(r)-S2(l-1)))%mod;
		} 
		n%=mod;
		return rec[n]=(S1(n)*S1(n)-ans)%mod;
	}
	void Soviet()
	{
		int ans=0;
		for(int l=2,r;l<=n;l=r+1)
		{
			r=n/(n/l);
			ans=(ans+S1(n/l)*(f(r)-f(l-1))%mod*i2%mod)%mod;
		}// 整除分块+杜教筛
		printf("%lld\n",((2ll*ans+S1(n))%mod+mod)%mod);
	}

	#define Flan void
	Flan IsMyWife()
	{
		// freopen("ans.txt","w",stdout);
		Init();
		Input();
		Soviet();
	}
	#undef int //long long
}
int main()
{
	Flandre_Scarlet::IsMyWife();
	getchar();getchar();
	return 0;
}

{% endfold %}