1 ///
2 /// @file pi_deleglise_rivat.cpp
3 /// @brief 64-bit and 128-bit parallel implementations of the
4 /// Deleglise-Rivat prime counting algorithm.
5 ///
6 /// Deleglise-Rivat formula:
7 /// pi(x) = pi(y) + S1(x, a) + S2(x, a) - 1 - P2(x, a)
8 /// S2(x, a) = S2_trivial(x, a) + S2_easy(x, a) + S2_hard(x, a)
9 /// with y = alpha * x^(1/3), a = pi(y)
10 ///
11 /// This implementation is based on the paper:
12 /// Tomás Oliveira e Silva, Computing pi(x): the combinatorial
13 /// method, Revista do DETUA, vol. 4, no. 6, March 2006,
14 /// pp. 759-768.
15 ///
16 /// Copyright (C) 2021 Kim Walisch, <kim.walisch@gmail.com>
17 ///
18 /// This file is distributed under the BSD License. See the COPYING
19 /// file in the top level directory.
20 ///
21
22 #include <primecount.hpp>
23 #include <primecount-internal.hpp>
24 #include <imath.hpp>
25 #include <PhiTiny.hpp>
26 #include <int128_t.hpp>
27 #include <print.hpp>
28 #include <S.hpp>
29
30 #include <stdint.h>
31 #include <string>
32
33 using namespace primecount;
34
35 namespace {
36
37 /// Calculate the contribution of the special leaves
38 template <typename T>
S2_trivial(T x,int64_t y,int64_t z,int64_t c,int threads)39 T S2(T x,
40 int64_t y,
41 int64_t z,
42 int64_t c,
43 T s2_approx,
44 int threads,
45 bool is_print)
46 {
47 T s2_trivial = S2_trivial(x, y, z, c, threads, is_print);
48 T s2_easy = S2_easy(x, y, z, c, threads, is_print);
49 T s2_hard_approx = s2_approx - (s2_trivial + s2_easy);
50 T s2_hard = S2_hard(x, y, z, c, s2_hard_approx, threads, is_print);
51 T s2 = s2_trivial + s2_easy + s2_hard;
52
53 return s2;
54 }
55
56 } // namespace
57
58 namespace primecount {
59
60 /// Calculate the number of primes below x using the
61 /// Deleglise-Rivat algorithm.
62 /// Run time: O(x^(2/3) / (log x)^2)
63 /// Memory usage: O(x^(1/3) * (log x)^3)
64 ///
65 int64_t pi_deleglise_rivat_64(int64_t x,
66 int threads,
67 bool is_print)
68 {
69 if (x < 2)
70 return 0;
71
72 double alpha = get_alpha_deleglise_rivat(x);
73 int64_t x13 = iroot<3>(x);
74 int64_t y = (int64_t) (x13 * alpha);
75 int64_t z = x / y;
76 int64_t pi_y = pi_noprint(y, threads);
77 int64_t c = PhiTiny::get_c(y);
78
79 if (is_print)
80 {
81 print("");
82 print("=== pi_deleglise_rivat_64(x) ===");
83 print("pi(x) = S1 + S2 + pi(y) - 1 - P2");
84 print(x, y, z, c, threads);
85 }
86
87 int64_t p2 = P2(x, y, threads, is_print);
88 int64_t s1 = S1(x, y, c, threads, is_print);
89 int64_t s2_approx = S2_approx(x, pi_y, p2, s1);
90 int64_t s2 = S2(x, y, z, c, s2_approx, threads, is_print);
91 int64_t phi = s1 + s2;
S2_trivial(int64_t x,int64_t y,int64_t z,int64_t c,int threads,bool is_print)92 int64_t sum = phi + pi_y - 1 - p2;
93
94 return sum;
95 }
96
97 #if defined(HAVE_INT128_T)
98
99 /// Calculate the number of primes below x using the
100 /// Deleglise-Rivat algorithm.
101 /// Run time: O(x^(2/3) / (log x)^2)
102 /// Memory usage: O(x^(1/3) * (log x)^3)
103 ///
104 int128_t pi_deleglise_rivat_128(int128_t x,
105 int threads,
106 bool is_print)
107 {
108 if (x < 2)
109 return 0;
110
111 double alpha = get_alpha_deleglise_rivat(x);
112 maxint_t limit = get_max_x(alpha);
113
114 if (x > limit)
115 throw primecount_error("pi(x): x must be <= " + to_string(limit));
116
S2_trivial(int128_t x,int64_t y,int64_t z,int64_t c,int threads,bool is_print)117 int64_t y = (int64_t) (iroot<3>(x) * alpha);
118 int64_t z = (int64_t) (x / y);
119 int64_t pi_y = pi_noprint(y, threads);
120 int64_t c = PhiTiny::get_c(y);
121
122 if (is_print)
123 {
124 print("");
125 print("=== pi_deleglise_rivat_128(x) ===");
126 print("pi(x) = S1 + S2 + pi(y) - 1 - P2");
127 print(x, y, z, c, threads);
128 }
129
130 int128_t p2 = P2(x, y, threads, is_print);
131 int128_t s1 = S1(x, y, c, threads, is_print);
132 int128_t s2_approx = S2_approx(x, pi_y, p2, s1);
133 int128_t s2 = S2(x, y, z, c, s2_approx, threads, is_print);
134 int128_t phi = s1 + s2;
135 int128_t sum = phi + pi_y - 1 - p2;
136
137 return sum;
138 }
139
140 #endif
141
142 } // namespace
143