1 ///
2 /// @file S2_xy.cpp
3 /// @brief Test the computation of the special leaves
4 /// S2(x, y) used in the Lagarias-Miller-Odlyzko
5 /// and Deleglise-Rivat prime counting algorithms.
6 ///
7 /// Copyright (C) 2020 Kim Walisch, <kim.walisch@gmail.com>
8 ///
9 /// This file is distributed under the BSD License. See the COPYING
10 /// file in the top level directory.
11 ///
12
13 #include <primecount.hpp>
14 #include <primecount-internal.hpp>
15 #include <PhiTiny.hpp>
16 #include <generate.hpp>
17 #include <imath.hpp>
18 #include <S.hpp>
19
20 #include <stdint.h>
21 #include <iostream>
22 #include <cstdlib>
23 #include <vector>
24 #include <random>
25
26 using namespace primecount;
27
check(bool OK)28 void check(bool OK)
29 {
30 std::cout << " " << (OK ? "OK" : "ERROR") << "\n";
31 if (!OK)
32 std::exit(1);
33 }
34
main()35 int main()
36 {
37 std::random_device rd;
38 std::mt19937 gen(rd());
39 std::uniform_int_distribution<int> dist(1, 10000000);
40 int threads = 1;
41
42 // test small x
43 for (int i = 1; i < 30000; i++)
44 {
45 int64_t x = i;
46 double alpha = get_alpha_deleglise_rivat(x);
47 int64_t y = (int64_t) (alpha * iroot<3>(x));
48 int64_t pi_y = pi_noprint(y, threads);
49 int64_t z = x / y;
50 int64_t c = PhiTiny::get_c(y);
51 int64_t s2 = 0;
52
53 auto primes = generate_primes<int32_t>(y);
54 auto lpf = generate_lpf(y);
55 auto mu = generate_moebius(y);
56
57 // special leaves
58 for (int64_t b = c + 1; b < pi_y; b++)
59 for (int64_t m = (y / primes[b]) + 1; m <= y; m++)
60 if (lpf[m] > primes[b])
61 s2 -= mu[m] * phi(x / (primes[b] * m), b - 1);
62
63 std::cout << "S2(" << x << ", " << y << ") = " << s2;
64
65 check(s2 == S2_trivial(x, y, z, c, threads) +
66 S2_easy(x, y, z, c, threads) +
67 S2_hard(x, y, z, c, Ri(x), threads));
68 }
69
70 // test random x
71 for (int i = 0; i < 500; i++)
72 {
73 int64_t x = dist(gen);
74 double alpha = get_alpha_deleglise_rivat(x);
75 int64_t y = (int64_t) (alpha * iroot<3>(x));
76 int64_t pi_y = pi_noprint(y, threads);
77 int64_t z = x / y;
78 int64_t c = PhiTiny::get_c(y);
79 int64_t s2 = 0;
80
81 auto primes = generate_primes<int32_t>(y);
82 auto lpf = generate_lpf(y);
83 auto mu = generate_moebius(y);
84
85 // special leaves
86 for (int64_t b = c + 1; b < pi_y; b++)
87 for (int64_t m = (y / primes[b]) + 1; m <= y; m++)
88 if (lpf[m] > primes[b])
89 s2 -= mu[m] * phi(x / (primes[b] * m), b - 1);
90
91 std::cout << "S2(" << x << ", " << y << ") = " << s2;
92
93 check(s2 == S2_trivial(x, y, z, c, threads) +
94 S2_easy(x, y, z, c, threads) +
95 S2_hard(x, y, z, c, Ri(x), threads));
96 }
97
98 std::cout << std::endl;
99 std::cout << "All tests passed successfully!" << std::endl;
100
101 return 0;
102 }
103