驱魔错误球球了
查看原帖
驱魔错误球球了
519384
Link_Cut_Y楼主2023/8/2 15:46

取模,打错字了。

写的 PN 筛。

#include <unordered_map>
#include <iostream>
#include <cstring>
#include <cstdio>
#define int long long
		
using namespace std;
		
const int N = 6000010;
const int mod = 1e9 + 7;
		
unordered_map<int, int> Gsum;
int phi[N], p[N], s[N];
int lim, ans, n, cnt, inv2, inv6;
bool is_prime[N];

int power(int a, int b = mod - 2) {
	int ans = 1;
	for (; b; b >>= 1, a = a * a % mod)
		if (b & 1) ans = ans * a % mod;
	return ans;
}
void init(int n) {
	phi[1] = 1;
	for (int i = 2; i <= n; i ++ ) {
		if (!is_prime[i]) p[ ++ cnt] = i, phi[i] = i - 1;
		for (int j = 1; j <= cnt and i * p[j] <= n; j ++ ) {
			is_prime[i * p[j]] = true;
			if (i % p[j] == 0) { phi[i * p[j]] = phi[i] * p[j]; break; }
			phi[i * p[j]] = phi[i] * phi[p[j]];
		}
	}
	for (int i = 1; i <= n; i ++ )
		s[i] = (s[i - 1] + i * phi[i] % mod) % mod;
}
int get_G(int n) {
	if (n <= lim) return s[n];
	if (Gsum[n]) return Gsum[n];
	int ans = n % mod * (n + 1) % mod * (n * 2 % mod + 1) % mod * inv6 % mod;
	for (int l = 2, r; l <= n; l = r + 1) {
		r = n / (n / l);
		ans -= (l % mod + r % mod) % mod * (r % mod - l % mod + 1 + mod) % mod * inv2 % mod * get_G(n / l) % mod;
		if (ans < 0) ans += mod;
	}
	return Gsum[n] = ans;
}
void dfs(int now, int num, int val) {
//	cout << now << ' ' << num << ' ' << val << endl;
//	system("pause");
	if (now > cnt or num * p[now] > n) {
		if (n / num <= lim) ans = (ans + val % mod * s[n / num] % mod) % mod;
		else ans = (ans + val % mod * Gsum[n / num] % mod) % mod;
		return;
	}
	dfs(now + 1, num, val); 
	int u = (p[now] - 1) % mod * p[now] % mod * p[now] % mod, tmp = p[now] * p[now];
	for (int i = 2; num * tmp <= n; i ++ , tmp = tmp * p[now]) {
		dfs(now + 1, num * tmp, val % mod * (i - 1) % mod * u % mod);
		u = u % mod * p[now] % mod;
	}
}
		
signed main() {
//	freopen("debug.txt", "w", stdout);
	inv2 = power(2);
	inv6 = power(6);
	scanf("%lld", &n);
	lim = 6000000;
	init(lim); get_G(n);
	dfs(1, 1, 1);
	printf("%lld\n", ans);
//	cout << "OK" << endl;
	return 0;
}
2023/8/2 15:46
加载中...