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, limit = min(int(pow(N, 1.0 / 3) * pow(log(N), 4.0 / 3) * 0.75), v);
  46. 	const int B = N / K, B2 = N / limit;
  47.  
  48. 	clock_t clk = clock();
  49.  
  50. 	sieve(B2);
  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<int> roughs;
  87. 	for (int i = n_6 + 1; i <= v; ++i)
  88. 		if (s0[i] != s0[i - 1]) roughs.push_back(i);
  89. 	roughs.push_back(v + 1);
  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 k = 0, l;
  94. 		lf[i] = l1[i] - l0[i] + s0[u] * sf[u];
  95. 		for (; (l = roughs[k]) <= t; ++k)
  96. 			lf[i] -= lf[i * l] + (sf[l] - sf[l - 1]) * l0[i * l];
  97. 		for (; (l = roughs[k]) <= u; ++k)
  98. 			lf[i] -= sf[divide(m, l)] + (sf[l] - sf[l - 1]) * s0[divide(m, l)];
  99. 	}
  100.  
  101. 	fprintf(stderr, "sieve_large done %lf\n", double(clock() - clk) / CLOCKS_PER_SEC);
  102. 	clk = clock();
  103.  
  104. 	for (int id = pi[n_6]; id; --id) {
  105. 		const u32 p = primes[id], t = v / p;
  106. 		const u64 M = N / p;
  107. 		for (int i = 1; i <= t; ++i) lf[i] -= lf[i * p];
  108. 		for (int i = t + 1; i <= v; ++i) lf[i] -= sf[divide(M, i)];
  109. 		for (int i = v, j = t; j; --j)
  110. 			for (int c = j * p; i >= c; --i) sf[i] -= sf[j];
  111. 		for (int j = 1, i = p; j <= t; ++j) {
  112. 			const u32 c = p * sf[j];
  113. 			for (int e = min<u32>(v + 1, i + p); i < e; ++i) sf[i] += c;
  114. 		}
  115. 		for (int i = v; i > t; --i) lf[i] += p * sf[divide(M, i)];
  116. 		for (int i = t; i; --i) lf[i] += p * lf[i * p];
  117. 	}
  118.  
  119. 	fprintf(stderr, "sieve_small done %lf\n", double(clock() - clk) / CLOCKS_PER_SEC);
  120. 	clk = clock();
  121.  
  122. 	u32 ans = lf[1];
  123.  
  124. 	const auto ff = [&](i64 n, int k) -> double {
  125. 		double P = 1;
  126. 		i64 x = 1;
  127. 		while (k <= pi[v] && x * primes[k] <= n) {
  128. 			x *= primes[k];
  129. 			P = P * primes[k] / (primes[k] - 1);
  130. 			++k;
  131. 		}
  132. 		return P;
  133. 		};
  134.  
  135. 	vector<u32> f(v + 1);
  136. 	for (int id = 1; id <= pi[v]; ++id) f[primes[id]] += primes[id] - 1;
  137. 	for (int i = 1; i <= v; ++i) f[i] += f[i - 1];
  138.  
  139. 	int lim = 0;
  140. 	vector<vector<pair<u64, u64>>> d1(pi[v] + 1);
  141. 	function <void(u64, int, u64)> dfs4 = [&](u64 d, int beg, u64 coeff) -> void {
  142. 		u64 m = divide(N, d);
  143. 		double g = d * 1. / coeff;
  144. 		if (floor(g) == floor(g * ff(m, beg + 1))) return;
  145. 		for (int i = beg + 1; i <= pi[v]; ++i) {
  146. 			u64 pr = primes[i];
  147. 			if (pr * pr > m) break;
  148. 			double g1 = g * pr / (pr - 1);
  149. 			if (floor(g1) != floor(g))
  150. 				d1[i].push_back(make_pair(d, coeff)), d* pr <= limit && (lim = max(lim, i));
  151. 			u64 q = 1;
  152. 			for (; q * pr <= m; q *= pr)
  153. 				dfs4(d * q * pr, i, coeff * q * (pr - 1));
  154. 		}
  155. 		int left = min<i64>(v, (floor(g) + 1) * 1. / (floor(g) + 1 - g));
  156. 		left = min<i64>(left, m);
  157. 		int right = max<i64>(primes[beg], sqrt(m));
  158. 		if (left > right) ans += coeff * (f[left] - f[right]);
  159. 		};
  160. 	dfs4(1, 0, 1);
  161.  
  162. 	fprintf(stderr, "dfs done %lf\n", double(clock() - clk) / CLOCKS_PER_SEC);
  163. 	clk = clock();
  164.  
  165. 	for (int id = 1; id <= lim; ++id) {
  166. 		const u32 p = primes[id], t = v / p;
  167. 		const u64 M = N / p;
  168. 		for (auto d : d1[id])
  169. 			ans += d.second * (d.first <= v ? lf[d.first] : sf[divide(N, d.first)]);
  170. 		for (int i = 1; i <= t; ++i) lf[i] -= p * lf[i * p];
  171. 		for (int i = t + 1; i <= v; ++i) lf[i] -= p * sf[divide(M, i)];
  172. 		for (int i = v, j = t; j; --j)
  173. 			for (int c = j * p; i >= c; --i) sf[i] -= p * sf[j];
  174. 		for (int j = 1, i = p; j <= t; ++j)
  175. 			for (int e = min<int>(v + 1, i + p); i < e; ++i) sf[i] += sf[j];
  176. 		for (int i = v; i > t; --i) lf[i] += sf[divide(M, i)];
  177. 		for (int i = t; i; --i) lf[i] += lf[i * p];
  178. 		for (auto d : d1[id])
  179. 			ans -= d.second * (d.first <= v ? lf[d.first] : sf[divide(N, d.first)]);
  180. 	}
  181.  
  182. 	fprintf(stderr, "dp done %lf\n", double(clock() - clk) / CLOCKS_PER_SEC);
  183. 	clk = clock();
  184.  
  185. 	vector<u32> bit(B2 + 1);
  186. 	const auto add = [&](int x, const u32 cnt) -> void {
  187. 		while (x <= B2) bit[x] += cnt, x += x & -x;
  188. 		};
  189. 	const auto query = [&](int x) -> u32 {
  190. 		u32 ans = 0;
  191. 		while (x) ans += bit[x], x ^= x & -x;
  192. 		return ans;
  193. 		};
  194. 	add(1, 1);
  195. 	for (int id = pi[B2]; id > pi[v]; --id)
  196. 		for (auto r : h[id]) add(r, phi[r]);
  197. 	for (int id = pi[v]; id > lim; --id) {
  198. 		const u32 p = primes[id];
  199. 		const u64 m = divide(N, p);
  200. 		//while (h[i].first > id) add(h[i].second, phi[h[i].second]), --i;
  201. 		for (auto d : d1[id]) {
  202. 			u32 n = m / d.first, phi0 = p - 1;
  203. 			while (n) {
  204. 				ans += phi0 * d.second * query(n);
  205. 				n /= p;
  206. 				phi0 *= p;
  207. 			}
  208. 		}
  209. 		for (auto r : h[id]) add(r, phi[r]);
  210. 		//ans += (p - 1) * d.second * (1 + f[divide(m, d.first)] - f[p] + p);
  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 1.92s 92496KB
comments (?)
stdin
copy 
1000000000000
stdout
copy 
3875137357
stderr
copy 
sieve_large done 0.533531
sieve_small done 0.106246
dfs done 0.709759
dp done 0.532496
bit done 0.018459
https://ideone.com/fTu5Kq
language:
C++14 (gcc 8.3)
created:
2 days 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!