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