fork download 
  1. // This can be achieved in Time complexity O(N^(2/3)) by presieving values of phi.
  2. // Space complexity is O(N^(2/3)) ie 4D + I.
  3. // We use that phi_sum(d) floor(N/d)^2 over 1 <= d <= N is the same as
  4. // what we need except it counts each pair where i != j twice (but we can easily
  5. // correct for this).
  6.  
  7. // Note This is not Deleglise Rivat. Deleglise Rivat achives similar time complexity but space compleity would be O(N^(1/3))
  8.  
  9. // O(N^(2/3)) precompute O(1) per query but query can be only among N,floor(N/2),floor(N/3),floor(N/4),...,floor(N/D)
  10. #include <bits/stdc++.h>
  11. using namespace std;
  12. using ll = long long;
  13.  
  14. // Prime modulus for this problem
  15. #define MOD 998244353
  16. // Inverse of 2 modulo MOD
  17. #define INV2 499122177
  18. #define N 100000000000 // 10^11
  19. // We must have D * I <= N < (D + 1) * I.
  20. // Optimal performance is achieved with I approximately N^(1/3).
  21. #define I 4641
  22. #define D 21547080
  23.  
  24. // Least prime factor of each integer. Necessary for linear sieve.
  25. int32_t lpf[D + 1];
  26. // Prime numbers below D. Will be much smaller than D + 1 in actuality.
  27. int32_t primes[D + 1];
  28. // Euler totient function below D.
  29. int32_t phi[D + 1];
  30. // Summatory totient function.
  31. int32_t phi_sum[D + 1];
  32. // large_phi_sum_compressed[i] will be the sum of the euler totient function up to N / i.
  33. int32_t large_phi_sum_compressed[I + 1];
  34.  
  35. // Multiply 64-bits integers modulo MOD.
  36. ll mul_mod(ll a, ll b) {
  37.   return ((ll)a * (ll)b) % MOD;
  38. }
  39.  
  40. // nth triangular number modulo MOD.
  41. int32_t T_mod(ll n) {
  42.   return mul_mod(mul_mod(n % MOD, (n + 1) % MOD), INV2);
  43. }
  44.  
  45. int main() {
  46.   phi[1] = 1;
  47.   int32_t nprimes = 0;
  48.   for (int32_t k = 2; k <= D; ++k) {
  49.     if (lpf[k] == 0) {
  50.       phi[k] = k - 1;
  51.       lpf[k] = k;
  52.       primes[nprimes++] = k;
  53.     }
  54.  
  55.     for (ll j = 0;
  56.          primes[j] < lpf[k] && (ll)primes[j] * (ll)k <= (ll)D;
  57.          ++j) {
  58.       lpf[primes[j] * k] = primes[j];
  59.       phi[primes[j] * k] = (primes[j] - 1) * phi[k];
  60.     }
  61.  
  62.     if ((ll)lpf[k] * (ll)k <= (ll)D) {
  63.       lpf[lpf[k] * k] = lpf[k];
  64.       phi[lpf[k] * k] = lpf[k] * phi[k];
  65.     }
  66.   }
  67.   for (int32_t d = 1; d <= D; ++d) {
  68.     phi_sum[d] = (phi_sum[d - 1] + phi[d]) % MOD;
  69.     // cout << "etf_sum(" << d << ")=" << phi_sum[d] << endl;
  70.   }
  71.  
  72.   large_phi_sum_compressed[I] = phi_sum[D];
  73.  
  74.   // Use recurrence phi_sum(N) = T(N) - sum_d >=2 phi_sum(N/d) to compute large_phi_sum_compressed
  75.   // values.
  76.   cout << "Now computing all larges, amortised O(N^(2/3)) for all query or O(1) per query" << endl;
  77.   for (ll i = I - 1; i > 0; --i) {
  78.     const ll x_i = N / i;
  79.     large_phi_sum_compressed[i] = T_mod(x_i);
  80.  
  81.     ll d = 2;
  82.     for (; d * i < I; ++d) {
  83.       large_phi_sum_compressed[i] = (large_phi_sum_compressed[i] - large_phi_sum_compressed[d * i] + MOD) % MOD;
  84.     }
  85.  
  86.     for (; d * d <= x_i; ++d) {
  87.       large_phi_sum_compressed[i] = (large_phi_sum_compressed[i] - phi_sum[x_i / d] + MOD) % MOD;
  88.     }
  89.     for (ll k = x_i / d; k > 0; --k) {
  90.       ll m = mul_mod((x_i / k - x_i / (k + 1)) % MOD, phi_sum[k]);
  91.       large_phi_sum_compressed[i] = (large_phi_sum_compressed[i] - m + MOD) % MOD;
  92.     }
  93.     // cout << "etf_sum(" << N << "/" << i << "=" << N/i << ")=" << large_phi_sum_compressed[i] <<
  94.     // endl;  // == etf_sum(N/i)
  95.   }
  96.  
  97.   ll ii = 3456;
  98.   if (ii < I) cout << "etf_sum(" << N << "/" << ii << "=" << N / ii << ")=" << large_phi_sum_compressed[ii]  << endl;
  99.   else cout << "etf_sum(" << N << "/" << ii << "=" << N / ii << ")=" << phi_sum[N / ii] << endl;
  100.   cout << clock() / (double)CLOCKS_PER_SEC << endl;
  101. }
  102.  
Success #stdin #stdout 2.36s 261660KB
comments (?)
stdin 
Standard input is empty
stdout 
Now computing all larges, amortised O(N^(2/3)) for all query or O(1) per query
etf_sum(100000000000/3456=28935185)=534740299
2.34478
https://ideone.com/XQ3xc6
language:
C++14 (gcc 8.3)
created:
2 years ago
visibility:
secret

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!