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

///
/// @file  pi_lmo2.cpp
/// @brief Simple demonstration implementation of the
///        Lagarias-Miller-Odlyzko prime counting algorithm.
///        This implementation uses the sieve of Eratosthenes
///        (without segmentation) 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 <vector>

using std::vector;
using namespace primecount;

namespace {

/// Calculate the contribution of the special leaves.
/// This implementation uses the sieve of Eratosthenes
/// (without segmentation).
/// Memory usage: O(x^(2/3) * (log x)^2)
///
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 s2 = 0;
  int64_t b = 1;
  vector<char> sieve(limit, 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 (; b <= c; b++)
  {
    int64_t prime = primes[b];
    for (int64_t k = prime; k < limit; k += prime)
      sieve[k] = 0;
  }

  for (; b < pi_y; b++)
  {
    int64_t prime = primes[b];
    int64_t i = 1;
    int64_t phi = 0;

    for (int64_t m = y; m > y / prime; 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; i++)
          phi += sieve[i];

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

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

  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^(2/3) / (log x)^2)
///
int64_t pi_lmo2(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