fork(1) download
copy 
  1. #include <cstdio>
  2. #include <cassert>
  3. #include <cmath>
  4.  
  5. #include <algorithm>
  6. #include <vector>
  7.  
  8. #include <tuple>
  9.  
  10. using namespace std;
  11.  
  12. using i64 = long long;
  13. using i128 = __int128_t;
  14. using u8 = unsigned char;
  15. using u32 = unsigned;
  16. using u64 = unsigned long long;
  17. using u128 = __uint128_t;
  18.  
  19. inline int isqrt(i64 n) {
  20.   return sqrtl(n);
  21. }
  22.  
  23. inline int icbrt(i64 n) {
  24.   return cbrtl(n);
  25. }
  26.  
  27. template <typename F>
  28. struct InplaceFenwickTree {
  29.   using SmallSum = typename F::SmallSum;
  30.   using LargeSum = typename F::LargeSum;
  31.   using Integer = typename F::Integer;
  32.  
  33.   InplaceFenwickTree(i64 N, int S, int L, vector<SmallSum>* smalls, vector<LargeSum>* larges)
  34.       : N(N), S(S), L(L), smalls(*smalls), larges(*larges), small_sum(SmallSum(0)) {}
  35.  
  36.   void merge_from_cumulative_form() {
  37.     for (int i = (S - L) + 1; i <= S - 1; ++i) larges[i] -= larges[i + 1];
  38.     if (L > 0) larges[S] -= smalls[S];
  39.     for (int i = S; i > 0; --i) smalls[i] -= smalls[i - 1];
  40.     merge();
  41.   }
  42.  
  43.   void merge() {
  44.     for (int i = 1; i <= S; ++i) {
  45.       int parent = i + (i & -i);
  46.       if (parent <= S) smalls[parent] += smalls[i];
  47.     }
  48.     small_sum = sum_lo(S);
  49.     for (int i = 1; i <= L; ++i) {
  50.       int parent = i + (i & -i);
  51.       if (parent <= L) larges[S + 1 - parent] += larges[S + 1 - i];
  52.     }
  53.   }
  54.  
  55.   void split_lo() {
  56.     for (int i = S; i > 0; --i) {
  57.       int parent = i + (i & -i);
  58.       if (parent <= S) smalls[parent] -= smalls[i];
  59.     }
  60.     for (int i = 1; i <= S; ++i) smalls[i] += smalls[i - 1];
  61.   }
  62.  
  63.   void split_hi() {
  64.     for (int i = L; i > 0; --i) {
  65.       int parent = i + (i & -i);
  66.       if (parent <= L) larges[S + 1 - parent] -= larges[S + 1 - i];
  67.     }
  68.     if (L > 0) larges[S] += smalls[S];
  69.     for (int i = S - 1; i > S - L; --i) larges[i] += larges[i + 1];
  70.   }
  71.  
  72.   void split() {
  73.     split_lo();
  74.     split_hi();
  75.   }
  76.  
  77.   SmallSum sum_lo(int n) const {
  78.     SmallSum ret = SmallSum(0);
  79.     for (; n > 0; n &= n - 1) ret += smalls[n];
  80.     return ret;
  81.   }
  82.  
  83.   LargeSum sum_hi(int n) const {
  84.     LargeSum ret = small_sum;
  85.     for (n = S + 1 - n; n > 0; n &= n - 1) ret += larges[S + 1 - n];
  86.     return ret;
  87.   }
  88.  
  89.   void add_lo(int n, Integer coeff) {
  90.     for (; n <= S; n += n & -n) smalls[n] += coeff;
  91.     small_sum += coeff;
  92.   }
  93.  
  94.   void add_hi(int n, Integer coeff) {
  95.     for (n = S + 1 - n; n <= L; n += n & -n) larges[S + 1 - n] += coeff;
  96.   }
  97.  
  98.   const i64 N;
  99.   const int S, L;
  100.   vector<SmallSum>& smalls;
  101.   vector<LargeSum>& larges;
  102.   SmallSum small_sum;
  103. };
  104.  
  105. template <typename F>
  106. class PrimeSummation {
  107.   using SmallSum = typename F::SmallSum;
  108.   using LargeSum = typename F::LargeSum;
  109.   using Integer = typename F::Integer;
  110.  
  111. public:
  112.   PrimeSummation(i64 N, const vector<int>& primes)
  113.       : N(N),
  114.         sqrt_N(isqrt(N)), cbrt_N(icbrt(N)), sixth_root_N(icbrt(sqrt_N)),
  115.         threshold(icbrt(N)), L(N / (threshold + 1)),
  116.         primes(primes),
  117.         smalls(sqrt_N + 1, SmallSum(0)), larges(sqrt_N + 1, LargeSum(0)),
  118.         fenwick(N, sqrt_N, sqrt_N - threshold, &this->smalls, &this->larges) {
  119.  
  120.     assert(N < (i64(1) << 62));
  121.     assert(N >= 1); assert(i64(cbrt_N + 1) * i64(cbrt_N + 1) > L);
  122.     for (int i = 1; i <= sqrt_N; ++i) smalls[i] = F::sum(i) - F::one();
  123.     for (int i = 1; i <= sqrt_N; ++i) larges[i] = F::sum(N / i) - F::one();
  124.     const int pi3 = upper_bound(primes.begin(), primes.end(), cbrt_N) - primes.begin();
  125.     const int pi4 = upper_bound(primes.begin(), primes.end(), isqrt(sqrt_N)) - primes.begin();
  126.     const int pi6 = upper_bound(primes.begin(), primes.end(), sixth_root_N) - primes.begin();
  127.     clock_t t0 = clock();
  128.     calc(0, pi6);
  129.     clock_t t1 = clock();
  130.     fprintf(stderr, "S: %.4f seconds\n", double(t1 - t0) / CLOCKS_PER_SEC);
  131.     fenwick.merge_from_cumulative_form();
  132.     calc_intermediate(pi6, pi4, [&] (int n) { return fenwick.sum_lo(n); });
  133.     fenwick.split_lo();
  134.     clock_t t21 = clock();
  135.     fprintf(stderr, "M: %.4f seconds\n", double(t21 - t1) / CLOCKS_PER_SEC);
  136.     calc_intermediate(pi4, pi3, [&] (int n) { return smalls[n]; });
  137.     fenwick.split_hi();
  138.     clock_t t22 = clock();
  139.     fprintf(stderr, "M: %.4f seconds\n", double(t22 - t21) / CLOCKS_PER_SEC);
  140.     calc(pi3, primes.size());
  141.     clock_t t3 = clock();
  142.     fprintf(stderr, "L: %.4f seconds\n", double(t3 - t22) / CLOCKS_PER_SEC);
  143.   }
  144.  
  145. private:
  146.   void calc(const int pi_lo, const int pi_hi) {
  147.     for (int pi = pi_lo; pi < pi_hi; ++pi) {
  148.       const int p = primes[pi];
  149.       const i64 M = N / p, q = i64(p) * p;
  150.       const int w = sqrt_N / p, v = min<i64>(sqrt_N, N / q);
  151.       const SmallSum psum = smalls[p - 1];
  152.       const Integer c = F()(p);
  153.       for (int i = 1; i <= w; ++i) larges[i] -= (larges[i * p] - psum) * c;
  154.       const int t = min(isqrt(M), v);
  155.       for (int i = w + 1; i <= t; ++i) larges[i] -= (smalls[M / i] - psum) * c;
  156.       for (int i = t + 1, q = M / i, r = M % i; i <= v; ++i, r -= q) {
  157.         if (r < 0) r += i, --q;
  158.         larges[i] -= (smalls[q] - psum) * c;
  159.       }
  160.       for (int j = sqrt_N / p; j >= p; --j) {
  161.         SmallSum coeff = (smalls[j] - psum) * c;
  162.         for (int i = j * p, e = min(sqrt_N + 1, (j + 1) * p); i < e; ++i) smalls[i] -= coeff;
  163.       }
  164.     }
  165.   }
  166.  
  167.   void rec(const i64 n, const i64 M, size_t pi_beg, const Integer coeff) {
  168.     if (n <= sqrt_N) fenwick.add_lo(n, coeff);
  169.     else fenwick.add_hi(M, coeff);
  170.     for (size_t pi = pi_beg; pi < primes.size(); ++pi) {
  171.       int p = primes[pi];
  172.       if (i128(n) * p > L) break;
  173.       i64 nn = n, nM = M; const Integer c = F()(p);
  174.       Integer ncoeff = coeff;
  175.       while (i128(nn) * p <= L) rec(nn *= p, nM /= p, pi + 1, ncoeff *= c);
  176.     }
  177.   }
  178.  
  179.   template <typename Func>
  180.   void calc_intermediate(const int pi_lo, const int pi_hi, Func func) {
  181.     for (int pi = pi_lo; pi < pi_hi; ++pi) {
  182.       const int p = primes[pi];
  183.       const SmallSum psum = func(p - 1);
  184.       const Integer c = F()(p);
  185.       const int w = threshold / p;
  186.       for (int i = 1; i <= w; ++i) larges[i] -= (larges[i * p] - psum) * c;
  187.       const i64 M = N / p;
  188.       const int u = min<i64>(threshold, M / p);
  189.       const int v = min<int>(u, sqrt_N / p);
  190.       for (int i = w + 1; i <= v; ++i) larges[i] -= (fenwick.sum_hi(i * p) - psum) * c;
  191.       for (int i = v + 1; i <= u; ++i) larges[i] -= (func(M / i) - psum) * c;
  192.       for (size_t pj = pi; pj < primes.size(); ++pj) {
  193.         const int q = primes[pj];
  194.         if (i64(p) * q > L) break;
  195.         rec(i64(p) * q, M / q, pj, -Integer(c) * F()(q));
  196.       }
  197.     }
  198.   }
  199.  
  200. private:
  201.   const i64 N; const int sqrt_N, cbrt_N, sixth_root_N;
  202.   const int threshold; const i64 L;
  203.   const vector<int>& primes;
  204.  
  205. public:
  206.   vector<SmallSum> smalls;
  207.   vector<LargeSum> larges;
  208.  
  209. private:
  210.   InplaceFenwickTree<F> fenwick;
  211. };
  212.  
  213. template <typename U, typename DU>
  214. struct FastDiv {
  215.   static constexpr int S = sizeof(U) * 8;
  216.   FastDiv() {}
  217.   FastDiv(U n) : m(n) {
  218.     s = __lg(n - 1);
  219.     x = ((DU(1) << (s + S)) + n - 1) / n;
  220.   }
  221.   friend U operator / (U n, FastDiv d) {
  222.     return (DU(n) * d.x >> d.s) >> S;
  223.   }
  224.   friend U operator % (U n, FastDiv d) { return n - n / d * d.m; }
  225.   U m, x; int s;
  226. };
  227.  
  228. template <typename F>
  229. class MultiplicativeSummation {
  230.   using SmallSum = typename F::SmallSum;
  231.   using LargeSum = typename F::LargeSum;
  232.   using Integer = typename F::Integer;
  233.  
  234. public:
  235.   MultiplicativeSummation(i64 N,
  236.         const vector<int>& primes,
  237.         vector<SmallSum>* smalls,
  238.         vector<LargeSum>* larges)
  239.       : N(N),
  240.         sqrt_N(isqrt(N)), cbrt_N(icbrt(N)), sixth_root_N(icbrt(sqrt_N)),
  241.         threshold(icbrt(N)), L(N / (threshold + 1)),
  242.         primes(primes),
  243.         smalls(*smalls), larges(*larges),
  244.         fenwick(N, sqrt_N, sqrt_N - threshold, &(*smalls), &(*larges)) {
  245.  
  246.     assert(N < (i64(1) << 62));
  247.     assert(N >= 1); assert(i64(cbrt_N + 1) * i64(cbrt_N + 1) > L);
  248.     const int pi3 = upper_bound(primes.begin(), primes.end(), cbrt_N) - primes.begin();
  249.     const int pi6 = upper_bound(primes.begin(), primes.end(), sixth_root_N) - primes.begin();
  250.  
  251.     clock_t t0 = clock();
  252.     calc_larges(pi3);
  253.     clock_t t1 = clock();
  254.     fprintf(stderr, "L: %.4f seconds\n", double(t1 - t0) / CLOCKS_PER_SEC);
  255.     fenwick.merge();
  256.     calc_intermediate(pi6, pi3);
  257.     fenwick.split();
  258.     clock_t t2 = clock();
  259.     fprintf(stderr, "M: %.4f seconds\n", double(t2 - t1) / CLOCKS_PER_SEC);
  260.     calc_smalls(pi6);
  261.     clock_t t3 = clock();
  262.     fprintf(stderr, "S: %.4f seconds\n", double(t3 - t2) / CLOCKS_PER_SEC);
  263.   }
  264.  
  265. private:
  266.   void calc_larges(size_t pi3) {
  267.     for (int i = 1; i <= threshold; ++i) {
  268.       const i64 M = N / i;
  269.       LargeSum s = F::one() + (larges[i] - smalls[cbrt_N]);
  270.       for (size_t pi = pi3; pi < primes.size(); ++pi) {
  271.         int p = primes[pi];
  272.         if (i64(p) * p > M) break;
  273.         i64 d = i64(i) * p;
  274.         s += F()(p, 1) * ((d <= sqrt_N ? larges[d] : smalls[M / p]) - smalls[p]);
  275.         s += F()(p, 2);
  276.       }
  277.       larges[i] = s;
  278.     }
  279.     for (int i = threshold + 1; i <= sqrt_N - 1; ++i) larges[i] -= larges[i + 1];
  280.     if (sqrt_N > threshold) larges[sqrt_N] -= smalls[sqrt_N];
  281.     for (int i = sqrt_N; i > cbrt_N; --i) smalls[i] -= smalls[i - 1];
  282.     for (int i = cbrt_N; i > 1; --i) smalls[i] = 0;
  283.     smalls[1] = F::one();
  284.   }
  285.  
  286.   void rec(const i64 n, const i64 M, size_t pi_beg, const Integer coeff) {
  287.     if (n <= sqrt_N) fenwick.add_lo(n, coeff);
  288.     else fenwick.add_hi(M, coeff);
  289.     for (size_t pi = pi_beg; pi < primes.size(); ++pi) {
  290.       int p = primes[pi];
  291.       if (i128(n) * p > L) break;
  292.       i64 nn = n, nM = M;
  293.       for (int e = 1; i128(nn) * p <= L; ++e) {
  294.         rec(nn *= p, nM /= p, pi + 1, coeff * F()(p, e));
  295.       }
  296.     }
  297.   }
  298.  
  299.   void calc_intermediate(const int pi_lo, const int pi3) {
  300.     std::pair<int, SmallSum> memo[65] = {}; // should be initialized with 0.
  301.     for (int pi = pi3 - 1; pi >= pi_lo; --pi) {
  302.       const int p = primes[pi]; auto fp = FastDiv<u32, u64>(p);
  303.       for (int i = 1; i <= threshold; ++i) {
  304.         LargeSum s = larges[i];
  305.         i64 d = i64(i) * p; int e = 1;
  306.         for (; d <= threshold; d *= p) s += F()(p, e++) * larges[d];
  307.         for (; d <= sqrt_N; d *= p) s += F()(p, e++) * fenwick.sum_hi(d);
  308.         int se = e;
  309.         for (int M = N / d; memo[e].first != M; ++e, M = M / fp) memo[e] = {M, F()(p, e) * fenwick.sum_lo(M)};
  310.         for (; se < e; --e) memo[e - 1].second += memo[e].second;
  311.         larges[i] = s + memo[e].second;
  312.       }
  313.       i64 q = 1; i64 M = N;
  314.       for (int e = 1; i128(q) * p <= L; ++e) rec(q *= p, M /= p, pi + 1, F()(p, e));
  315.     }
  316.   }
  317.  
  318.   void calc_smalls(const int pi_hi) {
  319.     std::pair<int, SmallSum> memo[65] = {}; // should be initialized with 0.
  320.     for (int pi = pi_hi - 1; pi >= 0; --pi) {
  321.       const int p = primes[pi]; auto fp = FastDiv<u32, u64>(p);
  322.       for (int i = 1; i <= sqrt_N; ++i) {
  323.         LargeSum s = larges[i];
  324.         i64 d = i64(i) * p; int e = 1;
  325.         for (; d <= sqrt_N; d *= p) s += F()(p, e++) * larges[d];
  326.         int se = e;
  327.         for (int M = N / d; memo[e].first != M; ++e, M = M / fp) memo[e] = {M, F()(p, e) * smalls[M]};
  328.         for (; se < e; --e) memo[e - 1].second += memo[e].second;
  329.         larges[i] = s + memo[e].second;
  330.       }
  331.       for (int j = sqrt_N / p; j >= 1; --j) {
  332.         const Integer coeff = smalls[j];
  333.         i64 q = p;
  334.         for (int e = 1; j * q <= sqrt_N; q *= p, ++e) {
  335.           const Integer c = F()(p, e) * coeff;
  336.           const int end = min<i64>(sqrt_N + 1, i64(j + 1) * q);
  337.           for (int i = j * q; i < end; ++i) smalls[i] += c;
  338.         }
  339.       }
  340.     }
  341.   }
  342.  
  343. private:
  344.   const i64 N; const int sqrt_N, cbrt_N, sixth_root_N;
  345.   const int threshold; const i64 L;
  346.   const vector<int>& primes;
  347.  
  348. public:
  349.   vector<SmallSum>& smalls;
  350.   vector<LargeSum>& larges;
  351.  
  352. private:
  353.   InplaceFenwickTree<F> fenwick;
  354. };
  355.  
  356. struct Ring {
  357.   using LargeSum = i64;
  358.   using SmallSum = int;
  359.   using Integer = int;
  360.   SmallSum operator () (int, int e) const { return e == 1 ? -SmallSum(1) : 0; }
  361.   static Integer one() { return SmallSum(1); }
  362. };
  363.  
  364. struct Op {
  365.   using LargeSum = i64;
  366.   using SmallSum = int;
  367.   using Integer = int;
  368.   Integer operator () (int p) const { return Integer(1); }
  369.   static LargeSum sum(i64 N) { return LargeSum(N); }
  370.   static Integer one() { return Integer(1); }
  371. };
  372.  
  373. vector<int> prime_sieve(int N) {
  374.   if (N <= 1) return vector<int>();
  375.   const int v = isqrt(N);
  376.   vector<bool> isp(N + 1, 1);
  377.   for (int p = 2; p <= v; ++p) if (isp[p]) {
  378.     for (int i = p * p; i <= N; i += p) isp[i] = false;
  379.   }
  380.   const int rsize = N >= 60194 ? N / (log(N) - 1.1)
  381.                                : max(1., N / (log(N) - 1.11)) + 1;
  382.   vector<int> primes(rsize); int ps = 0;
  383.   for (int i = 2; i <= N; ++i) if (isp[i]) primes[ps++] = i;
  384.   assert(ps <= rsize);
  385.   primes.resize(ps);
  386.   return primes;
  387. }
  388.  
  389. void verify() {
  390.   for (i64 N = 1; N <= 1e6; ++N) {
  391.     const int v = isqrt(N);
  392.     auto primes = prime_sieve(v);
  393.     auto ps = PrimeSummation<Op>(N, primes);
  394.     auto ms = MultiplicativeSummation<Ring>(N, primes, &ps.smalls, &ps.larges);
  395.     assert(ms.larges[1] == Op::LargeSum(N));
  396.     if (N % 100 == 0) printf("OK: %llu\n", N);
  397.   }
  398. }
  399.  
  400. /*
  401. S: 0.0580 seconds
  402. M: 0.1119 seconds
  403. M: 0.1056 seconds
  404. L: 0.0501 seconds
  405. 37607912018
  406. L: 0.0517 seconds
  407. M: 0.4052 seconds
  408. S: 0.1782 seconds
  409. -522626
  410. 0.97678
  411. 0.980 sec
  412.  
  413. S: 0.2800 seconds
  414. M: 0.6059 seconds
  415. M: 0.5038 seconds
  416. L: 0.2281 seconds
  417. 346065536839
  418. L: 0.2337 seconds
  419. M: 2.0540 seconds
  420. S: 0.9004 seconds
  421. 10000000000000
  422. 4.84525
  423. 4.853 sec
  424.  
  425. S: 1.1435 seconds
  426. M: 2.5596 seconds
  427. M: 2.4708 seconds
  428. L: 1.0569 seconds
  429. 3204941750802
  430. L: 1.0805 seconds
  431. M: 9.8681 seconds
  432. S: 3.8567 seconds
  433. 100000000000000
  434. 22.15263
  435. 22.171 sec
  436.  
  437. S: 4.7559 seconds
  438. M: 10.7700 seconds
  439. M: 12.5694 seconds
  440. L: 4.9361 seconds
  441. 29844570422669
  442. L: 5.1237 seconds
  443. M: 45.7336 seconds
  444. S: 16.3929 seconds
  445. 1000000000000000
  446. 100.65623
  447. 100.703 sec
  448.  
  449. */
  450.  
  451. void solve() {
  452.   const i64 N = 1e13;
  453.   const int v = isqrt(N);
  454.   auto primes = prime_sieve(v);
  455.   clock_t beg = clock();
  456.   auto ps = PrimeSummation<Op>(N, primes);
  457.   printf("%lld\n", ps.larges[1]);
  458.   for (int i = 1; i <= v; ++i) ps.smalls[i] *= -Op::Integer(1);
  459.   for (int i = 1; i <= v; ++i) ps.larges[i] *= -Op::Integer(1);
  460.   auto ms = MultiplicativeSummation<Ring>(N, primes, &ps.smalls, &ps.larges);
  461.   printf("%lld\n", ms.larges[1]);
  462.   clock_t end = clock();
  463.   printf("%.5f\n", double(end - beg) / CLOCKS_PER_SEC);
  464. }
  465.  
  466. int main() {
  467.   clock_t beg = clock();
  468.   solve();
  469.   clock_t end = clock();
  470.   fprintf(stderr, "%.3f sec\n", double(end - beg) / CLOCKS_PER_SEC);
  471.   return 0;
  472. }
  473.  
Success #stdin #stdout #stderr 8.2s 15264KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
346065536839
599582
8.22094
stderr
copy 
S: 0.3948 seconds
M: 0.9185 seconds
M: 0.9635 seconds
L: 0.4750 seconds
L: 0.4604 seconds
M: 3.1667 seconds
S: 1.8016 seconds
8.237 sec
https://ideone.com/9F9amv
Multiplicative Sum [2]
language:
C++14 (gcc 8.3)
created:
5 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!