1 ///
2 /// @file phi_xa.cpp
3 /// @brief Test the partial sieve function phi(x, a)
4 /// which counts the numbers <= x that are not divisible
5 /// by any of the first a primes.
6 ///
7 /// Copyright (C) 2021 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 <primesieve.hpp>
16 #include <imath.hpp>
17
18 #include <stdint.h>
19 #include <iostream>
20 #include <cstdlib>
21 #include <vector>
22 #include <random>
23
24 using std::size_t;
25 using namespace primecount;
26
check(size_t x,size_t a,size_t phi_xa,size_t cnt)27 void check(size_t x, size_t a, size_t phi_xa, size_t cnt)
28 {
29 bool OK = (phi_xa == cnt);
30 std::cout << "phi(" << x << ", " << a << ") = " << phi_xa;
31 std::cout << " " << (OK ? "OK" : "ERROR") << "\n";
32 if (!OK)
33 std::exit(1);
34 }
35
check2(size_t x,size_t a,size_t phi_xa,size_t cnt)36 void check2(size_t x, size_t a, size_t phi_xa, size_t cnt)
37 {
38 if (phi_xa != cnt)
39 {
40 bool OK = (phi_xa == cnt);
41 std::cout << "phi(" << x << ", " << a << ") = " << phi_xa;
42 std::cout << " " << (OK ? "OK" : "ERROR") << "\n";
43 std::exit(1);
44 }
45 // Reduce logging because it is slow
46 else if (a % 101 == 0)
47 {
48 bool OK = (phi_xa == cnt);
49 std::cout << "phi(" << x << ", " << a << ") = " << phi_xa;
50 std::cout << " " << (OK ? "OK" : "ERROR") << "\n";
51 }
52 }
53
main()54 int main()
55 {
56 std::random_device rd;
57 std::mt19937 gen(rd());
58
59 {
60 std::uniform_int_distribution<size_t> dist(20000000, 30000000);
61 size_t size = dist(gen);
62 size_t x = size - 1;
63 size_t cnt = size - 1;
64 primesieve::iterator it;
65 std::vector<char> sieve(size, 1);
66
67 // test with small a values
68 for (size_t a = 1;; a++)
69 {
70 auto prime = it.next_prime();
71 if (prime * prime > x)
72 break;
73
74 // remove primes[a] and its multiples
75 for (auto j = prime; j <= x; j += prime)
76 {
77 cnt -= (sieve[j] == 1);
78 sieve[j] = 0;
79 }
80
81 int64_t phi_xa = phi(x, a);
82 check(x, a, phi_xa, cnt);
83 }
84 }
85
86 {
87 std::uniform_int_distribution<size_t> dist(100000, 200000);
88 size_t size = dist(gen);
89 size_t x = size - 1;
90 size_t cnt = size - 1;
91 primesieve::iterator it;
92 std::vector<char> sieve(size, 1);
93
94 // test with large a values
95 for (size_t a = 1;; a++)
96 {
97 auto prime = it.next_prime();
98 if (prime > x)
99 break;
100
101 // remove primes[a] and its multiples
102 for (auto j = prime; j <= x; j += prime)
103 {
104 cnt -= (sieve[j] == 1);
105 sieve[j] = 0;
106 }
107
108 int64_t phi_xa = phi(x, a);
109 check2(x, a, phi_xa, cnt);
110 }
111 }
112
113 {
114 std::cout << "Testing phi(x, a) multi-threading" << std::endl;
115
116 int64_t iters = 500;
117 int64_t sum1 = 0;
118 int64_t sum2 = 0;
119
120 #pragma omp parallel for reduction(+: sum1)
121 for (int64_t i = 0; i < iters; i++)
122 sum1 += pi_legendre(10000000 + i, 1);
123
124 for (int64_t i = 0; i < iters; i++)
125 sum2 += pi_legendre(10000000 + i, 1);
126
127 if (sum1 == sum2)
128 {
129 std::cout << "Multi-thread sum: " << sum1 << " == Single-thread sum: " << sum2 << " OK" << std::endl;
130 std::cout << "phi(x, a) multi-threading: no data races detected!" << std::endl;
131 }
132 else
133 {
134 std::cout << "Multi-thread sum: " << sum1 << " != Single-thread sum: " << sum2 << " ERROR" << std::endl;
135 std::exit(1);
136 }
137 }
138
139 std::cout << std::endl;
140 std::cout << "All tests passed successfully!" << std::endl;
141
142 return 0;
143 }
144