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
143
|
///
/// @file phi.cpp
/// @brief Test the partial sieve function phi(x, a)
/// which counts the numbers <= x that are not divisible
/// by any of the first a primes.
///
/// 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.hpp>
#include <primecount-internal.hpp>
#include <primesieve.hpp>
#include <imath.hpp>
#include <stdint.h>
#include <iostream>
#include <cstdlib>
#include <vector>
#include <random>
using std::size_t;
using namespace primecount;
void check(size_t x, size_t a, size_t phi_xa, size_t cnt)
{
bool OK = (phi_xa == cnt);
std::cout << "phi(" << x << ", " << a << ") = " << phi_xa;
std::cout << " " << (OK ? "OK" : "ERROR") << "\n";
if (!OK)
std::exit(1);
}
void check2(size_t x, size_t a, size_t phi_xa, size_t cnt)
{
if (phi_xa != cnt)
{
bool OK = (phi_xa == cnt);
std::cout << "phi(" << x << ", " << a << ") = " << phi_xa;
std::cout << " " << (OK ? "OK" : "ERROR") << "\n";
std::exit(1);
}
// Reduce logging because it is slow
else if (a % 101 == 0)
{
bool OK = (phi_xa == cnt);
std::cout << "phi(" << x << ", " << a << ") = " << phi_xa;
std::cout << " " << (OK ? "OK" : "ERROR") << "\n";
}
}
int main()
{
std::random_device rd;
std::mt19937 gen(rd());
{
std::uniform_int_distribution<size_t> dist(20000000, 30000000);
size_t size = dist(gen);
size_t x = size - 1;
size_t cnt = size - 1;
primesieve::iterator it;
std::vector<char> sieve(size, 1);
// test with small a values
for (size_t a = 1;; a++)
{
auto prime = it.next_prime();
if (prime * prime > x)
break;
// remove primes[a] and its multiples
for (auto j = prime; j <= x; j += prime)
{
cnt -= (sieve[j] == 1);
sieve[j] = 0;
}
int64_t phi_xa = phi(x, a);
check(x, a, phi_xa, cnt);
}
}
{
std::uniform_int_distribution<size_t> dist(100000, 200000);
size_t size = dist(gen);
size_t x = size - 1;
size_t cnt = size - 1;
primesieve::iterator it;
std::vector<char> sieve(size, 1);
// test with large a values
for (size_t a = 1;; a++)
{
auto prime = it.next_prime();
if (prime > x)
break;
// remove primes[a] and its multiples
for (auto j = prime; j <= x; j += prime)
{
cnt -= (sieve[j] == 1);
sieve[j] = 0;
}
int64_t phi_xa = phi(x, a);
check2(x, a, phi_xa, cnt);
}
}
{
std::cout << "Testing phi(x, a) multi-threading" << std::endl;
int64_t iters = 500;
int64_t sum1 = 0;
int64_t sum2 = 0;
#pragma omp parallel for reduction(+: sum1)
for (int64_t i = 0; i < iters; i++)
sum1 += pi_legendre(10000000 + i, 1);
for (int64_t i = 0; i < iters; i++)
sum2 += pi_legendre(10000000 + i, 1);
if (sum1 == sum2)
{
std::cout << "Multi-thread sum: " << sum1 << " == Single-thread sum: " << sum2 << " OK" << std::endl;
std::cout << "phi(x, a) multi-threading: no data races detected!" << std::endl;
}
else
{
std::cout << "Multi-thread sum: " << sum1 << " != Single-thread sum: " << sum2 << " ERROR" << std::endl;
std::exit(1);
}
}
std::cout << std::endl;
std::cout << "All tests passed successfully!" << std::endl;
return 0;
}
|
