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

///
/// @file  pi_gourdon.cpp
/// @brief Implementation of Xavier Gourdon's prime counting
///        algorithm. Xavier Gourdon's algorithm is an improved
///        version of the Deleglise-Rivat algorithm, according to my
///        benchmarks Gourdon's algorithm runs up to 2x faster than
///        the Deleglise-Rivat algorithm.
///
///        Xavier Gourdon formula:
///        pi(x) = A - B + C + D + Phi0 + Sigma
///
/// 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 <gourdon.hpp>
#include <primecount.hpp>
#include <primecount-internal.hpp>
#include <imath.hpp>
#include <PhiTiny.hpp>
#include <print.hpp>

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

namespace primecount {

/// Calculate the number of primes below x using
/// Xavier Gourdon's algorithm.
/// Run time: O(x^(2/3) / (log x)^2)
/// Memory usage: O(x^(1/3) * (log x)^3)
///
int64_t pi_gourdon_64(int64_t x,
                      int threads,
                      bool is_print)
{
  if (x < 2)
    return 0;

  auto alpha = get_alpha_gourdon(x);
  double alpha_y = alpha.first;
  double alpha_z = alpha.second;
  int64_t x13 = iroot<3>(x);
  int64_t sqrtx = isqrt(x);
  int64_t y = (int64_t)(x13 * alpha_y);

  // x^(1/3) < y < x^(1/2)
  y = std::max(y, x13 + 1);
  y = std::min(y, sqrtx - 1);
  y = std::max(y, (int64_t) 1);

  int64_t k = PhiTiny::get_k(x);
  int64_t z = (int64_t)(y * alpha_z);

  // y <= z < x^(1/2)
  z = std::max(z, y);
  z = std::min(z, sqrtx - 1);
  z = std::max(z, (int64_t) 1);

  if (is_print)
  {
    print("");
    print("=== pi_gourdon_64(x) ===");
    print("pi(x) = A - B + C + D + Phi0 + Sigma");
    print_gourdon(x, y, z, k, threads);
  }

  int64_t sigma = Sigma(x, y, threads, is_print);
  int64_t phi0 = Phi0(x, y, z, k, threads, is_print);
  int64_t b = B(x, y, threads, is_print);
  int64_t ac = AC(x, y, z, k, threads, is_print);
  int64_t d_approx = D_approx(x, sigma, phi0, ac, b);
  int64_t d = D(x, y, z, k, d_approx, threads, is_print);
  int64_t sum = ac - b + d + phi0 + sigma;

  return sum;
}

#if defined(HAVE_INT128_T)

/// Calculate the number of primes below x using
/// Xavier Gourdon's algorithm.
/// Run time: O(x^(2/3) / (log x)^2)
/// Memory usage: O(x^(1/3) * (log x)^3)
///
int128_t pi_gourdon_128(int128_t x,
                        int threads,
                        bool is_print)
{
  if (x < 2)
    return 0;

  auto alpha = get_alpha_gourdon(x);
  double alpha_y = alpha.first;
  double alpha_z = alpha.second;
  maxint_t limit = get_max_x(alpha_y);

  if (x > limit)
    throw primecount_error("pi(x): x must be <= " + to_string(limit));

  int64_t x13 = iroot<3>(x);
  int64_t sqrtx = isqrt(x);
  int64_t y = (int64_t)(x13 * alpha_y);

  // x^(1/3) < y < x^(1/2)
  y = std::max(y, x13 + 1);
  y = std::min(y, sqrtx - 1);
  y = std::max(y, (int64_t) 1);

  int64_t k = PhiTiny::get_k(x);
  int64_t z = (int64_t)(y * alpha_z);

  // y <= z < x^(1/2)
  z = std::max(z, y);
  z = std::min(z, sqrtx - 1);
  z = std::max(z, (int64_t) 1);

  if (is_print)
  {
    print("");
    print("=== pi_gourdon_128(x) ===");
    print("pi(x) = A - B + C + D + Phi0 + Sigma");
    print_gourdon(x, y, z, k, threads);
  }

  int128_t sigma = Sigma(x, y, threads, is_print);
  int128_t phi0 = Phi0(x, y, z, k, threads, is_print);
  int128_t b = B(x, y, threads, is_print);
  int128_t ac = AC(x, y, z, k, threads, is_print);
  int128_t d_approx = D_approx(x, sigma, phi0, ac, b);
  int128_t d = D(x, y, z, k, d_approx, threads, is_print);
  int128_t sum = ac - b + d + phi0 + sigma;

  return sum;
}

#endif

} // namespace