primecount 0.2.1

Rust wrapper for https://github.com/kimwalisch/primecount
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152

///
/// @file  pi_lmo3.cpp
/// @brief Simple demonstration implementation of the
///        Lagarias-Miller-Odlyzko prime counting algorithm.
///        This implementation uses the segmented sieve of
///        Eratosthenes to calculate S2(x).
///
///        Lagarias-Miller-Odlyzko formula:
///        pi(x) = pi(y) + S1(x, a) + S2(x, a) - 1 - P2(x, a)
///        with y = x^(1/3), a = pi(y)
///
/// Copyright (C) 2021 Kim Walisch, <kim.walisch@gmail.com>
///
/// This file is distributed under the BSD License. See the COPYING
/// file in the top level directory.
///

#include <primecount-internal.hpp>
#include <imath.hpp>
#include <generate.hpp>
#include <PhiTiny.hpp>
#include <S.hpp>

#include <stdint.h>
#include <algorithm>
#include <vector>

using std::min;
using std::max;
using std::vector;
using namespace primecount;

namespace {

/// Calculate the contribution of the special leaves.
/// This implementation uses segmentation which reduces the
/// algorithm's memory usage to O(x^(1/3) * log^2 x).
///
int64_t S2(int64_t x,
           int64_t y,
           int64_t c,
           int64_t pi_y,
           const vector<int32_t>& primes,
           const vector<int32_t>& lpf,
           const vector<int32_t>& mu)
{
  int64_t limit = x / y;
  int64_t segment_size = isqrt(limit);
  int64_t s2 = 0;

  vector<char> sieve(segment_size);
  vector<int64_t> next(primes.begin(), primes.end());
  vector<int64_t> phi(primes.size(), 0);

  // segmented sieve of Eratosthenes
  for (int64_t low = 1; low < limit; low += segment_size)
  {
    // current segment [low, high[
    int64_t high = min(low + segment_size, limit);

    std::fill(sieve.begin(), sieve.end(), 1);

    // phi(y, b) nodes with b <= c do not contribute to S2, so
    // we sieve out the multiples of the first c primes
    for (int64_t b = 1; b <= c; b++)
    {
      int64_t k = next[b];
      for (int64_t prime = primes[b]; k < high; k += prime)
        sieve[k - low] = 0;
      next[b] = k;
    }

    for (int64_t b = c + 1; b < pi_y; b++)
    {
      int64_t prime = primes[b];
      int64_t min_m = max(x / (prime * high), y / prime);
      int64_t max_m = min(x / (prime * low), y);
      int64_t i = 0;

      // Obviously if (prime >= max_m) then (prime >= lpf[max_m])
      // hence (prime < lpf[m]) will always evaluate to false
      // and no special leaves are possible.
      if (prime >= max_m)
        break;

      for (int64_t m = max_m; m > min_m; m--)
      {
        if (mu[m] != 0 && prime < lpf[m])
        {
          // We have found a special leaf. Compute it's contribution
          // phi(x / (primes[b] * m), b - 1) by counting the number
          // of unsieved elements <= x / (primes[b] * m) after having
          // removed the multiples of the first b - 1 primes.
          //
          for (int64_t xpm = x / (prime * m); i <= xpm - low; i++)
            phi[b] += sieve[i];

          s2 -= mu[m] * phi[b];
        }
      }

      // Count the remaining unsieved elements in this segment,
      // we need their count in the next segment.
      for (; i < high - low; i++)
        phi[b] += sieve[i];

      // remove the multiples of the b-th prime
      int64_t k = next[b];
      for (; k < high; k += prime * 2)
        sieve[k - low] = 0;
      next[b] = k;
    }
  }

  return s2;
}

} // namespace

namespace primecount {

/// Calculate the number of primes below x using the
/// Lagarias-Miller-Odlyzko algorithm.
/// Run time: O(x^(2/3))
/// Memory usage: O(x^(1/3) * (log x)^2)
///
int64_t pi_lmo3(int64_t x)
{
  if (x < 2)
    return 0;

  bool threads = 1;
  double alpha = get_alpha_lmo(x);
  int64_t x13 = iroot<3>(x);
  int64_t y = (int64_t) (x13 * alpha);
  int64_t c = PhiTiny::get_c(y);
  int64_t p2 = P2(x, y, threads);

  auto primes = generate_primes<int32_t>(y);
  auto lpf = generate_lpf(y);
  auto mu = generate_moebius(y);

  int64_t pi_y = primes.size() - 1;
  int64_t s1 = S1(x, y, c, threads);
  int64_t s2 = S2(x, y, c, pi_y, primes, lpf, mu);
  int64_t phi = s1 + s2;
  int64_t sum = phi + pi_y - 1 - p2;

  return sum;
}

} // namespace