51Nod 1238 最小公倍数之和 V3

出一个数N，输出小于等于N的所有数，两两之间的最小公倍数之和。

相当于计算这段程序（程序中的lcm(i,j)表示i与j的最小公倍数）：
由于结果很大，输出Mod 1000000007的结果。
```
G=0;
for(i=1;i< N;i++)
for(j=1;j<=N;j++)
{
 G = (G + lcm(i,j)) % 1000000007;
}
```
Input
输入一个数N。(2 <= N <= 10^10)
Output
输出G Mod 1000000007的结果。
Input示例
4
Output示例
72

- 推导：

$$G(x)=\sum^{n}_{i=1}\sum^{n}_{j=1}lcm(i, j)$$
$$G(x)=\sum^{n}_{i=1}\sum^{n}_{j=1}\frac{i\cdot j}{gcd(i, j)}$$
$$G(x)=\sum^{n}_{d=1}\frac{1}{d}\sum^{n}_{i=1}\sum^{n}_{j=1}{i\cdot j}\cdot [gcd(i, j)=d]$$
$$G(x)=\sum^{n}_{d=1}d\cdot \sum^{\frac{n}{d}}_{i=1}\sum^{\frac{n}{d}}_{j=1}{i\cdot j}\cdot [gcd(i, j)=1]$$

----------

- 引理1：$$\sum^{n}_{i=1}i\cdot [gcd(i,n)=1]=\frac{n\cdot\phi(n)}{2}$$
证明：由已知：\
若有 $$gcd(n,i) = 1$$
根据辗转相除：
$$gcd(n,n-i)=1$$
那么对于任意i与n互质,都有一个n-i满足与n互质，且和为n。
又因为小于i与i互质的数的个数为$\phi(i)$ 则所有小于n与n互质的 数的和为$$\frac{n\cdot \phi(n)}{2}$$
则有：$$\sum^{n}_{i=1}i\cdot [gcd(i,n)=1]=\frac{n\cdot\phi(n)}{2}$$  得证

----------

- 继续推导
$$G(x)=\sum^{n}_{d=1}d\cdot \sum^{\frac{n}{d}}_{i=1}{i\cdot }\sum^{\frac{n}{d}}_{j=1}{j}\cdot [gcd(i, j)=1]$$
因为$i=1, j=1$时被算了两次，扣去1：
$$G(x)=\sum^{n}_{d=1}d(2\sum^{\frac{n}{d}}_{i=1}{i\cdot}\sum^{i}_{j=1}{j}\cdot [gcd(i, j)=1] - 1)$$
由引理1，又因为$i=1$时有$gcd(1, 1) = 1$则有：$$G(x)=\sum^{n}_{d=1}d(2 \sum^{\frac{n}{d}}_{i=1}{i\cdot }\frac{i\cdot \phi(i)+[i=1]}{2} - 1)$$
$$G(x)=\sum^{n}_{d=1}d\cdot \sum^{\frac{n}{d}}_{i=1}{i^2\phi(i)}$$
设$$h(x)={x^2 \phi(x)}$$
可以知道，$h(x)$ 为积性函数，
考虑快速求出$h(x)$的前缀和。

----------

- 引理2：$$\sum_{d|n}\phi(d)=n$$ 
证明: 对于1~n-1中所有的数将其按照与n的最大公因数分类，相同的分为一类（即gcd(i, n)相同的分为一类）那么有：
 - 1、所有$x \in [1, n)$ 与d的最大公因数集合为n的因数集合。因为每一个数可能和n有的gcd只有可能是n的因数。而且每一个n的因数自己和n的gcd一定是它本身，那么对于n的每一个因数至少都能找到一个数使它与n的gcd为这个n的因数。
 - 2、对于每一个n的因数m，有$$\sum^{n}_{i=1}[gcd(i,n)=m]=\phi(m)$$ 考虑每一个合法的i，对于$\frac{i}{m}$和$\frac{n}{m}$，$\frac{i}{m}$一定与$\frac{n}{m}$互质。而对于每一个因数m，$\frac{n}{m}$ 一定也是n的因数，那么我们只需要统计下式即可：
$$\sum_{m|n}\phi(\frac{n}{m}) = \sum_{m|n}\phi(m)$$
得证。

----------

- 引理3：$\sum^{n}_{i=1}i^3=\frac{1}{4}n^2(n+1)^2$
证明：如下图，观察可得

$$\sum^{n}_{i=1}i^3 = S_{square} =(\sum^{n}_{i=1}i)^2=\frac{n^2(n+1)^2}{4}$$

![图片标题](https://leanote.com/api/file/getImage?fileId=5a615b4cab64413654001a69)
得证。

----------

那么有：
$$\sum^{n}_{d=1}\sum_{i|d}\phi(i)\cdot d^2=\sum^{n}_{d=1}d^2\cdot \sum_{i|d}\phi(i)=\sum^{n}_{d=1}
d^3=\frac{n^2(n+1)^2}{4}$$
又有：$$\sum^{n}_{d=1}\sum_{i|d}\phi(i)\cdot d^2=\sum^{n}_{d=1}\sum_{i|d}\frac{h(i)}{i^2}\cdot d^2= \sum^{n}_{d=1}\sum_{i|d}h(i)\cdot (\frac{d}{i})^2$$
考虑上式在枚举1~$n$ 中所有数的因数，不是很好求。我们反着计算对于每一个数，它对它倍数的贡献。那么我们这样，枚举每一个数的倍数统计答案：
$$\sum^{n}_{d=1}\sum_{i|d}h(i)\cdot (\frac{d}{i})^2=\sum^{n}_{d=1}\sum^{\frac{n}{d}}_{i=1}h(i)\cdot (\frac{d\cdot i}{i})^2 \\ =\sum^{n}_{d=1}\sum^{\frac{n}{d}}_{i=1}h(i)\cdot d^2=\sum^{n}_{d=1}d^2\cdot \sum^{\frac{n}{d}}_{i=1}h(i)$$
又因为，原式为$\sum^{n}_{d=1}d^3$，设$s(n)=\sum^{n}_{i=1}h(i)$，有：
$$\sum^{n}_{d=1}d^2\cdot \sum^{\frac{n}{d}}_{i=1}h(i)=\sum^{n}_{d=1}d^3$$
$$\sum^{n}_{d=1}d^2\cdot s(\frac{n}{d})=\sum^{n}_{d=1}d^3$$
$$\sum^{n}_{d=2}d^2\cdot s(\frac{n}{d}) + s(n)=\sum^{n}_{d=1}d^3$$
$$ s(n)=\sum^{n}_{d=1}d^3 - \sum^{n}_{d=2}d^2\cdot s(\frac{n}{d})\\ =\frac{n^2(n+1)^2}{4} - \sum^{n}_{d=2}d^2\cdot s(\frac{n}{d})$$
这样我们成功变成了杜教筛的形式，可以在$O(n^{\frac{2}{3}})$的时间内计算

AC代码：
```
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<tr1/unordered_map>
#define LL long long
#define sqr(x) ((x) * (x))
using namespace std;
using namespace tr1;
const int maxn = 6e6 + 6;
const LL mod = 1e9 + 7;
int pri[maxn], vis[maxn];
LL n, m, ans, phi[maxn], inv4, inv2, inv6, s[maxn];
unordered_map<LL, LL> hsh;
LL s3(LL x) {return (((sqr(x % mod) % mod) * (sqr((x + 1) % mod) % mod)) % mod) * inv4 % mod;}
LL s2(LL x) {return x %= mod, (((x * (x + 1)) % mod) * ((2 * x + 1) % mod) % mod) * inv6 % mod;}
LL fpow(LL a, LL b, LL c) {
	LL ans = 1;
	while(b) {
		if(b & 1) ans = ans * a % c;
		a = a * a % c, b >>= 1;
	}
	return ans;
}
void GetPri(int N) {
	vis[1] = phi[1] = s[1] = 1;
	for(register int i = 1; i <= N; ++i) {
		if(!vis[i]) pri[++pri[0]] = i, phi[i] = i - 1;
		for(register int j = 1; j <= pri[0] && i * pri[j] <= N; ++j) {
			vis[i * pri[j]] = 1;
			phi[i * pri[j]] = phi[i] * phi[pri[j]];
			if(i % pri[j] == 0) {
				phi[i * pri[j]] = phi[i] * pri[j];
				break;
			}
		}
	}
	for(register int i = 1; i <= N; ++i) {
		(s[i] = s[i - 1] + (((phi[i] * i) % mod) * i) % mod) %= mod;
	}
		
}

LL GetSum(LL x) {
	if(x <= m) return s[x];
	if(hsh[x]) return hsh[x];
	LL ret = s3(x);
	LL l = 2, r = 1, a;
	while(r < x) {
		l = r + 1, r = x / (x / l);
		a = GetSum(x / l), ret = (ret - (s2(r) - s2(l - 1)) * a % mod ) % mod;
		hsh[x / l] = a;
	}
	return (ret + mod) % mod;
}

int main() {
	scanf("%lld", &n);
	m = (int)pow(n, 2.0 / 3.0) + 100;
	inv4 = fpow(4, mod - 2, mod);
	inv6 = fpow(6, mod - 2, mod);
	inv2 = fpow(2, mod - 2, mod);
	GetPri(m);
	LL l = 1, r = 0;
	while(r < n) {
		l = r + 1, r = n / (n / l);
		(ans += (GetSum(n / r) % mod) * (((((l + r) % mod) * ((r - l + 1) % mod) % mod) * inv2) % mod) % mod) %= mod;
	}
	printf("%lld", (ans + mod) % mod);
	return 0;
} 
```

Rockdu's Blog

导航

最近发表

大佬们的友链