【51Nod 1238】最小公倍数之和 V3

http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1238
设\(A(n)=\sum\limits_{i=1}^n\frac{in}{(i,n)}\)，则\(ans=\sum\limits_{i=1}^n\left(2A(i)-i\right)\)

\[\begin{aligned} A(n)=&n\sum_{d|n}\sum_{i=1}^{\frac nd}i\left[\left(i,\frac nd\right)=1\right]\\ =&\frac n2\sum_{d|n}\left(\varphi(d)d+[d=1]\right)\\ =&\frac n2\sum_{d|n}\varphi(d)d+\frac n2\\ ans=&\sum_{i=1}^n\left(2A(i)-i\right)\\ =&\sum_{i=1}^ni\sum_{d|i}\varphi(d)d\\ =&\sum_{i=1}^ni\sum_{d=1}^{\left\lfloor\frac ni\right\rfloor}\varphi(d)d^2\\ \end{aligned} \]

设\(S(n)=\sum\limits_{i=1}^n\varphi(i)i^2\)，求出\(O\left(\sqrt n\right)\)不同下取整取值的\(S\)就可以算出答案了，所以现在重点是杜教筛\(S(n)\)。
先让\(f(n)=\varphi(n)n^2\)卷上\(f(n)=n^2\)

\[\sum_{i=1}^{n}\sum_{d|i}\varphi(d)d^2\left(\frac id\right)^2=\sum_{i=1}^ni^2\sum_{d|i}\varphi(d)=\sum_{i=1}^ni^3 \]

\[\begin{aligned} \sum_{i=1}^ni^3=&\sum_{i=1}^{n}\sum_{d|i}\varphi(d)d^2\left(\frac id\right)^2\\ =&\sum_{i=1}^ni^2\sum_{d=1}^{\left\lfloor\frac ni\right\rfloor}\varphi(d)d^2\\ =&\sum_{i=1}^ni^2S\left(\left\lfloor\frac ni\right\rfloor\right)\\ \end{aligned} \]

如果要求\(S(n)\)，\(S(n)=\sum\limits_{i=1}^ni^3-\sum\limits_{i=2}^ni^2S\left(\left\lfloor\frac ni\right\rfloor\right)\)，时间复杂度\(O\left(n^{\frac 23}\right)\)。
PS：\(\sum\limits_{i=1}^ni^3=\left(\frac{n(n+1)}{2}\right)^2\)

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;

const int N = 100003;
const int B = 4641589;
const int p = 1000000007;
const int ni2 = 500000004;
const int ni6 = 166666668;

ll n;
bool notp[B + 1];
int phi[B + 1], sum[B + 1], num = 0, prime[B + 1];

void Euler_shai() {
	sum[1] = phi[1] = 1;
	for (int i = 2; i <= B; ++i) {
		if (!notp[i]) prime[++num] = i, phi[i] = i - 1;
		for (int j = 1; j <= num && prime[j] * i <= B; ++j) {
			notp[prime[j] * i] = true;
			if (i % prime[j] == 0) {
				phi[prime[j] * i] = prime[j] * phi[i];
				break;
			} else
				phi[prime[j] * i] = (prime[j] - 1) * phi[i];
		}
		sum[i] = (1ll * phi[i] * i % p * i % p + sum[i - 1]) % p;
	}
}

int ps[B];
ll ndx;
int S(ll x) {return (ndx = n / x) <= B ? sum[ndx] : ps[x];}

void DJ_shai() {
	for (ll i = n, y; i >= 1; i = n / (y + 1)) {
		y = n / i;
		if (y <= B) continue;
		int &s = ps[i];
		s = y % p * (y % p) % p * ((y + 1) % p) % p * ((y + 1) % p) % p * ni2 % p * ni2 % p;
		for (ll j = 2, pre = 1, spre = 1, sj; j <= y; spre = sj, pre = j, ++j) {
			j = y / (y / j);
			sj = j % p * ((j + 1) % p) % p * ((j * 2 + 1) % p) % p * ni6 % p;
			((s -= (sj - spre + p) % p * S(i * j) % p) += p) %= p;
		}
	}
}

main() {
	scanf("%lld", &n);
	Euler_shai();
	
	int ans = 0;
	DJ_shai();
	
	for (ll i = 1, pre = 0; i <= n; pre = i, ++i) {
		i = n / (n / i);
		(ans += (i + pre + 1) % p * ((i - pre) % p) % p * ni2 % p * S(i) % p) %= p;
	}
	
	printf("%d\n", ans);
	return 0;
}