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