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

http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1238
设(A(n)=sumlimits_{i=1}^nfrac{in}{(i,n)})，则(ans=sumlimits_{i=1}^nleft(2A(i)-i ight))

[egin{aligned} A(n)=&nsum_{d|n}sum_{i=1}^{frac nd}ileft[left(i,frac nd ight)=1 ight]\ =&frac n2sum_{d|n}left(varphi(d)d+[d=1] ight)\ =&frac n2sum_{d|n}varphi(d)d+frac n2\ ans=&sum_{i=1}^nleft(2A(i)-i ight)\ =&sum_{i=1}^nisum_{d|i}varphi(d)d\ =&sum_{i=1}^nisum_{d=1}^{leftlfloorfrac ni ight floor}varphi(d)d^2\ end{aligned} ]

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

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

[egin{aligned} sum_{i=1}^ni^3=&sum_{i=1}^{n}sum_{d|i}varphi(d)d^2left(frac id ight)^2\ =&sum_{i=1}^ni^2sum_{d=1}^{leftlfloorfrac ni ight floor}varphi(d)d^2\ =&sum_{i=1}^ni^2Sleft(leftlfloorfrac ni ight floor ight)\ end{aligned} ]

如果要求(S(n))，(S(n)=sumlimits_{i=1}^ni^3-sumlimits_{i=2}^ni^2Sleft(leftlfloorfrac ni ight floor ight))，时间复杂度(Oleft(n^{frac 23} ight))。
PS：(sumlimits_{i=1}^ni^3=left(frac{n(n+1)}{2} ight)^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
", ans);
	return 0;
}