fork download
copy 
  1. #pragma GCC optimize(2)
  2. #include <iostream>
  3. #include <vector>
  4. #include <cmath>
  5. #include <utility>
  6. #include <functional>
  7. #include <algorithm> 
  8. #include <ctime>
  9. #include <bitset>
  10. using namespace std;
  11. using u32 = unsigned int;
  12. using u64 = unsigned long long;
  13. using i64 = long long;
  14. vector<int> pi, pos;
  15. vector<u32> phi, primes;
  16. const int MAX = 8000000;
  17. static bitset<MAX> is_prime;
  18. void sieve(int v) {
  19. 	pi.resize(v + 1);
  20. 	phi.resize(v + 1);
  21. 	pos.resize(v + 1);
  22.  
  23. 	primes.push_back(1);
  24. 	phi[1] = 1;
  25. 	for (int i = 2, t; i <= v; ++i) {
  26. 		if (!phi[i]) pos[i] = primes.size(), primes.push_back(i), phi[i] = i - 1;
  27. 		for (int j = 1; j <= pos[i] && (t = i * primes[j]) <= v; ++j) {
  28. 			pos[t] = j;
  29. 			if (j == pos[i]) {
  30. 				phi[t] = phi[i] * primes[j];
  31. 				break;
  32. 			}
  33. 			else
  34. 				phi[t] = phi[i] * (primes[j] - 1);
  35. 		}
  36. 	}
  37. 	for (int id = 1; id < primes.size(); ++id) pi[primes[id]] = 1;
  38. 	for (int i = 1; i <= v; ++i) pi[i] += pi[i - 1];
  39. }
  40. u32 solve(const u64 N) {
  41. 	const int v = sqrt(N + 0.5);
  42. 	const int n_6 = pow(N + 0.5, 1.0 / 6) * 0.4;
  43. 	const int n_4 = sqrt(v + 0.5);
  44. 	const int K = n_6 * v;
  45. 	const int B = N / K;
  46. 	const int limit = min<i64>(cbrt(N) * log(N), v);
  47. 	const int B2 = N / limit;
  48.  
  49. 	clock_t clk = clock();
  50.  
  51. 	sieve(v);
  52.  
  53. 	vector<u32> s0(v + 1), s1(v + 1), l0(v + 1), l1(v + 1);
  54. 	const auto divide = [](u64 n, u64 d) -> u64 {return double(n) / d; };
  55. 	for (int i = 1; i <= v; ++i) s0[i] = i, s1[i] = u64(i) * (i + 1) / 2;
  56. 	for (int i = 1; i <= v; ++i) l0[i] = N / i, l1[i] = (N / i) * (N / i + 1) / 2;
  57. 	for (int id = 1; id <= pi[n_6]; ++id) {
  58. 		const u32 p = primes[id], t = v / p;
  59. 		const i64 M = N / p;
  60. 		for (int i = 1; i <= t; ++i)
  61. 			l0[i] -= l0[i * p], l1[i] -= l1[i * p] * p;
  62. 		for (int i = t + 1; i <= v; ++i)
  63. 			l0[i] -= s0[divide(M, i)], l1[i] -= s1[divide(M, i)] * p;
  64. 		for (int i = v, j = t; j; --j)
  65. 			for (int e = j * p; i >= e; --i)
  66. 				s0[i] -= s0[j], s1[i] -= s1[j] * p;
  67. 	}
  68.  
  69. 	vector<u32> sf(v + 1), lf(v + 1);
  70. 	function <void(u64, int, u32)> dfs = [&](u64 n, int beg, u32 coeff) -> void {
  71. 		if (n <= v) sf[n] += coeff - n + 1;
  72. 		else lf[divide(N, n)] += coeff - n + 1;
  73. 		int m = K / n;
  74. 		for (int i = beg + 1; i <= pi[v]; ++i) {
  75. 			if (primes[i] > m) break;
  76. 			u64 q = 1;
  77. 			for (; q * primes[i] <= m; q *= primes[i])
  78. 				dfs(n * q * primes[i], i, coeff * q * (primes[i] - 1));
  79. 		}
  80. 		};
  81. 	dfs(1, pi[n_6], 1);
  82. 	for (int i = 1; i <= v; ++i) sf[i] += sf[i - 1];
  83. 	lf[v] += sf[v];
  84. 	for (int i = v - 1; i > B; --i) lf[i] += lf[i + 1];
  85. 	for (int i = 1; i <= v; ++i) sf[i] += s1[i] - s0[i], lf[i] += l1[i] - l0[i];
  86.  
  87. 	vector<int> roughs;
  88. 	for (int i = n_6 + 1; i <= v; ++i)
  89. 		if (s0[i] != s0[i - 1]) roughs.push_back(i);
  90. 	roughs.push_back(v + 1);
  91. 	for (int i = B; i >= 1; --i) {
  92. 		const u64 m = divide(N, i);
  93. 		const int u = sqrt(m + 0.5), t = divide(v, i);
  94. 		int k = 0, l;
  95. 		lf[i] = l1[i] - l0[i] + s0[u] * sf[u];
  96. 		for (; (l = roughs[k]) <= t; ++k)
  97. 			lf[i] -= lf[i * l] + (sf[l] - sf[l - 1]) * l0[i * l];
  98. 		for (; (l = roughs[k]) <= u; ++k)
  99. 			lf[i] -= sf[divide(m, l)] + (sf[l] - sf[l - 1]) * s0[divide(m, l)];
  100. 	}
  101.  
  102. 	fprintf(stderr, "sieve done %lf\n", double(clock() - clk) / CLOCKS_PER_SEC);
  103. 	clk = clock();
  104.  
  105. 	u32 ans = 0;
  106.  
  107. 	const auto ff = [&](i64 n, int k) -> double {
  108. 		double P = 1;
  109. 		i64 x = 1;
  110. 		while (k <= pi[v] && x * primes[k] <= n) {
  111. 			x *= primes[k];
  112. 			P = P * primes[k] / (primes[k] - 1);
  113. 			++k;
  114. 		}
  115. 		return P;
  116. 		};
  117.  
  118. 	vector<u32> f(pi[v] + 1);
  119. 	for (int id = 1; id <= pi[v]; ++id) f[id] += f[id - 1] + primes[id] - 1;
  120.  
  121. 	int lim = 0;
  122. 	vector<vector<pair<u64, u32>>> d1(pi[v] + 1);
  123. 	function <void(u64, int, u64)> dfs1 = [&](u64 d, int beg, u64 coeff) -> void {
  124. 		u64 m = divide(N, d);
  125. 		double g = d * 1. / coeff;
  126. 		if (int(g) == int(g * ff(m, beg + 1))) return;
  127. 		for (int i = beg + 1; i <= pi[v]; ++i) {
  128. 			u64 pr = primes[i];
  129. 			if (pr * pr > m) break;
  130. 			double g1 = g * pr / (pr - 1);
  131. 			if (int(g1) != int(g))
  132. 				d1[i].push_back(make_pair(d, coeff)), d * pr <= limit && (lim = max(lim, i));
  133. 			u64 q = 1;
  134. 			for (; q * pr <= m; q *= pr)
  135. 				dfs1(d * q * pr, i, coeff * q * (pr - 1));
  136. 		}
  137. 		int left = pi[min<i64>(min<i64>(v, (int(g) + 1) * 1. / (int(g) + 1 - g)), m)];
  138. 		int right = max(beg, pi[sqrt(m)]);
  139. 		if (left > right) ans += coeff * (f[left] - f[right]);
  140. 		};
  141. 	dfs1(1, 0, 1);
  142.  
  143. 	fprintf(stderr, "dfs done %lf\n", double(clock() - clk) / CLOCKS_PER_SEC);
  144. 	clk = clock();
  145.  
  146. 	for (int i = 1; i <= v; ++i) l0[i] = lf[i], s0[i] = sf[i];
  147. 	for (int id = pi[n_6]; id; --id) {
  148. 		const u32 p = primes[id], t = v / p;
  149. 		const u64 M = N / p;
  150. 		if (id <= lim) {
  151. 			for (auto d : d1[id])
  152. 				ans -= d.second * (d.first <= v ? l0[d.first] : s0[divide(N, d.first)]);
  153. 		}
  154. 		for (int i = 1; i <= t; ++i) l0[i] -= l0[i * p];
  155. 		for (int i = t + 1; i <= v; ++i) l0[i] -= s0[divide(M, i)];
  156. 		for (int i = v, j = t; j; --j)
  157. 			for (int c = j * p; i >= c; --i) s0[i] -= s0[j];
  158. 		for (int j = 1, i = p; j <= t; ++j) {
  159. 			const u32 c = p * s0[j];
  160. 			for (int e = min<u32>(v + 1, i + p); i < e; ++i) s0[i] += c;
  161. 		}
  162. 		for (int i = v; i > t; --i) l0[i] += p * s0[divide(M, i)];
  163. 		for (int i = t; i; --i) l0[i] += p * l0[i * p];
  164. 		if (id <= lim) {
  165. 			for (auto d : d1[id])
  166. 				ans += d.second * (d.first <= v ? l0[d.first] : s0[divide(N, d.first)]);
  167. 		}
  168. 	}
  169. 	ans += l0[1];
  170. 	for (int id = pi[n_6] + 1; id <= lim; ++id) {
  171. 		const u32 p = primes[id], t = v / p;
  172. 		const u64 M = N / p;
  173. 		for (auto d : d1[id])
  174. 			ans += d.second * (d.first <= v ? lf[d.first] : sf[divide(N, d.first)]);
  175. 		for (int i = 1; i <= t; ++i) lf[i] -= p * lf[i * p];
  176. 		for (int i = t + 1; i <= v; ++i) lf[i] -= p * sf[divide(M, i)];
  177. 		for (int i = v, j = t; j; --j)
  178. 			for (int c = j * p; i >= c; --i) sf[i] -= p * sf[j];
  179. 		for (int j = 1, i = p; j <= t; ++j)
  180. 			for (int e = min<int>(v + 1, i + p); i < e; ++i) sf[i] += sf[j];
  181. 		for (int i = v; i > t; --i) lf[i] += sf[divide(M, i)];
  182. 		for (int i = t; i; --i) lf[i] += lf[i * p];
  183. 		for (auto d : d1[id])
  184. 			ans -= d.second * (d.first <= v ? lf[d.first] : sf[divide(N, d.first)]);
  185. 	}
  186.  
  187. 	fprintf(stderr, "dp done %lf\n", double(clock() - clk) / CLOCKS_PER_SEC);
  188. 	clk = clock();
  189.  
  190. 	vector<u32> bit(B2 + 1);
  191. 	const auto add = [&](int x, const u32 cnt) -> void {
  192. 		while (x <= B2) bit[x] += cnt, x += x & -x;
  193. 		};
  194. 	const auto query = [&](int x) -> u32 {
  195. 		u32 ans = 0;
  196. 		while (x) ans += bit[x], x ^= x & -x;
  197. 		return ans;
  198. 		};
  199.  
  200. 	add(1, 1);
  201.  
  202. 	int v1 = sqrt(B2), v0 = (B2 - 1) / 2;
  203. 	for (int id = 2, p = 3; p <= v1; p = primes[++id])
  204. 			for (int j = p * p >> 1; j <= v0; j += p)
  205. 				is_prime[j] = true;
  206. 	for (int i = (v + 1) / 2; i <= v0; ++i)
  207. 		if (!is_prime[i])
  208. 			add(2 * i + 1, 2 * i), primes.push_back(2 * i + 1);
  209.  
  210. 	function <void(u64, int, u32)> dfs2 = [&](u64 n, int id, u32 fn) {
  211. 		add(n, fn);
  212. 		const int m = B2 / n;
  213. 		for (int i = id + 1; i < primes.size(); ++i) {
  214. 			const int pr = primes[i];
  215. 			if (pr > m) break;
  216. 			u64 q = 1;
  217. 			for (; q * pr <= m; q *= pr)
  218. 				dfs2(n * q * pr, i, fn * q * (pr - 1));
  219. 		}
  220. 	};
  221.  
  222. 	for (int id = pi[v]; id > lim; --id) {
  223. 		const u32 p = primes[id];
  224. 		const u64 m = divide(N, p);
  225. 		for (auto d : d1[id]) {
  226. 			u32 n = m / d.first, phi0 = p - 1;
  227. 			while (n) {
  228. 				ans += phi0 * d.second * query(n);
  229. 				n /= p;
  230. 				phi0 *= p;
  231. 			}
  232. 		}
  233. 		u64 q = 1;
  234. 		while (p * q <= B2) {
  235. 			dfs2(p * q, id, q * (p - 1));
  236. 			q *= p;
  237. 		}
  238. 	}
  239. 	fprintf(stderr, "bit done %lf\n", double(clock() - clk) / CLOCKS_PER_SEC);
  240. 	return N * (N + 1) / 2 - ans;
  241. }
  242. signed main() {
  243. 	u64 n;
  244. 	cin >> n;
  245. 	cout << solve(n);
  246. 	return 0;
  247. }
Success #stdin #stdout #stderr 4.07s 184932KB
comments (?)
stdin
copy 
5000000000000
stdout
copy 
1132365761
stderr
copy 
sieve done 1.500418
dfs done 1.204779
dp done 1.221135
bit done 0.128766
https://ideone.com/SJcFiy
language:
C++14 (gcc 8.3)
created:
5 hours ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!