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