diff --git a/lib/primesieve/CMakeLists.txt b/lib/primesieve/CMakeLists.txt index fe639d82..d92fe472 100644 --- a/lib/primesieve/CMakeLists.txt +++ b/lib/primesieve/CMakeLists.txt @@ -70,6 +70,7 @@ set(LIB_SRC src/api-c.cpp src/MemoryPool.cpp src/PrimeGenerator.cpp src/nthPrime.cpp + src/nthPrimeApprox.cpp src/ParallelSieve.cpp src/popcount.cpp src/PreSieve.cpp diff --git a/lib/primesieve/ChangeLog b/lib/primesieve/ChangeLog index 2fb20de5..ca8b49e0 100644 --- a/lib/primesieve/ChangeLog +++ b/lib/primesieve/ChangeLog @@ -1,6 +1,8 @@ -Changes in version 11.2, 02/01/2024 +Changes in version 11.2, 08/01/2024 =================================== +* nthPrime.cpp: Rewritten using faster nth prime approximation. +* nthPrimeApprox.cpp: Logarithmic integral and Riemann R implementations. * cmake/libatomic.cmake: Fix failed to find libatomic #141. * .github/workflows/ci.yml: Port AppVeyor CI tests to GitHub Actions. * doc/C_API.md: Fix off by 1 error in OpenMP example #137. diff --git a/lib/primesieve/include/primesieve/PrimeSieve.hpp b/lib/primesieve/include/primesieve/PrimeSieve.hpp index 59b58d91..ba5fb529 100644 --- a/lib/primesieve/include/primesieve/PrimeSieve.hpp +++ b/lib/primesieve/include/primesieve/PrimeSieve.hpp @@ -78,6 +78,7 @@ class PrimeSieve // nth prime uint64_t nthPrime(uint64_t); uint64_t nthPrime(int64_t, uint64_t); + uint64_t negativeNthPrime(int64_t, uint64_t); // Count counts_t& getCounts(); uint64_t getCount(int) const; diff --git a/lib/primesieve/include/primesieve/nthPrimeApprox.hpp b/lib/primesieve/include/primesieve/nthPrimeApprox.hpp new file mode 100644 index 00000000..58e68891 --- /dev/null +++ b/lib/primesieve/include/primesieve/nthPrimeApprox.hpp @@ -0,0 +1,23 @@ +/// +/// @file nthPrimeApprox.hpp +/// +/// Copyright (C) 2024 Kim Walisch, +/// +/// This file is distributed under the BSD License. See the COPYING +/// file in the top level directory. +/// + +#include +#include + +namespace primesieve { + +Vector generate_moebius(int64_t max); +uint64_t Li(uint64_t x); +uint64_t Li_inverse(uint64_t x); +uint64_t Ri(uint64_t x); +uint64_t Ri_inverse(uint64_t x); +uint64_t primesApprox(uint64_t x); +uint64_t nthPrimeApprox(uint64_t n); + +} // namespace diff --git a/lib/primesieve/src/nthPrime.cpp b/lib/primesieve/src/nthPrime.cpp index 631ced9c..e524220d 100644 --- a/lib/primesieve/src/nthPrime.cpp +++ b/lib/primesieve/src/nthPrime.cpp @@ -1,7 +1,7 @@ /// /// @file nthPrime.cpp /// -/// Copyright (C) 2022 Kim Walisch, +/// Copyright (C) 2024 Kim Walisch, /// /// This file is distributed under the BSD License. See the COPYING /// file in the top level directory. @@ -13,84 +13,31 @@ #include #include #include +#include #include #include #include #include - -using namespace primesieve; +#include namespace { -void checkLimit(uint64_t start) -{ - if (start >= get_max_stop()) - throw primesieve_error("nth prime > 2^64"); -} - -void checkLowerLimit(uint64_t stop) -{ - if (stop == 0) - throw primesieve_error("nth prime < 2 is impossible"); -} - -bool sieveBackwards(int64_t n, int64_t count, uint64_t stop) -{ - return (count >= n) && - !(count == n && stop < 2); -} +/// PrimePi(2^64) +const uint64_t max_n = 425656284035217743ull; -// Prime count approximation -int64_t pix(int64_t n) +/// Average prime gap near n +inline uint64_t avgPrimeGap(uint64_t n) { double x = (double) n; - x = std::max(4.0, x); - double pix = x / std::log(x); - return (int64_t) pix; -} - -uint64_t nthPrimeDist(int64_t n, int64_t count, uint64_t start) -{ - double x = (double) (n - count); - - x = std::abs(x); - x = std::max(x, 4.0); - - // rough pi(x) approximation + x = std::max(8.0, x); double logx = std::log(x); - double loglogx = std::log(logx); - double pix = x * (logx + loglogx - 1); - - // correct start if sieving backwards to - // get more accurate approximation - if (count >= n) - { - double st = start - pix; - st = std::max(0.0, st); - start = (uint64_t) st; - } - - // approximate the nth prime using: - // start + n * log(start + pi(n) / loglog(n)) - double startPix = start + pix / loglogx; - startPix = std::max(4.0, startPix); - double logStartPix = std::log(startPix); - double dist = std::max(pix, x * logStartPix); - - // ensure start + dist <= nth prime - if (count < n) - dist -= std::sqrt(dist) * std::log(logStartPix) * 2; - // ensure start + dist >= nth prime - if (count > n) - dist += std::sqrt(dist) * std::log(logStartPix) * 2; - - // if n is very small: - // ensure start + dist >= nth prime - double primeGap = maxPrimeGap(startPix); - dist = std::max(dist, primeGap); - - return (uint64_t) dist; + // When we buffer primes using primesieve::iterator, we + // want to make sure we buffer primes up to the nth + // prime. Therefore we use +2 here, better to buffer + // slightly too many primes than not enough primes. + double primeGap = logx + 2; + return (uint64_t) primeGap; } } // namespace @@ -104,60 +51,124 @@ uint64_t PrimeSieve::nthPrime(uint64_t n) uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start) { + if (n < 0) + return negativeNthPrime(n, start); + else if (n == 0) + n = 1; // like Mathematica + else if ((uint64_t) n > max_n) + throw primesieve_error("nth_prime(n): n must be <= " + std::to_string(max_n)); + setStart(start); auto t1 = std::chrono::system_clock::now(); + uint64_t nApprox = checkedAdd(primesApprox(start), n); + nApprox = std::min(nApprox, max_n); + uint64_t primeApprox = nthPrimeApprox(nApprox); + primeApprox = std::max(primeApprox, start); + int64_t countApprox = 0; + uint64_t prime = 0; - if (n == 0) - n = 1; // like Mathematica - else if (n > 0) + // Only use multi-threading if the sieving distance is sufficiently + // large. For small n this if statement also avoids calling + // countPrimes() and hence the initialization overhead of + // O(x^0.5 log log x^0.5) occurs only once (instead of twice) when + // using primesieve::iterator further down. + if (primeApprox - start > isqrt(primeApprox) / 10) + { + // Count primes > start start = checkedAdd(start, 1); - else if (n < 0) - start = checkedSub(start, 1); - - uint64_t stop = start; - uint64_t dist = nthPrimeDist(n, 0, start); - uint64_t nthPrimeGuess = checkedAdd(start, dist); - - int64_t count = 0; - int64_t tinyN = 100000; - tinyN = std::max(tinyN, pix(isqrt(nthPrimeGuess))); + primeApprox = std::max(start, primeApprox); + countApprox = countPrimes(start, primeApprox); + start = primeApprox; + } - while ((n - count) > tinyN || - sieveBackwards(n, count, stop)) + // Here we are very close to the nth prime < sqrt(nth_prime), + // we simply iterate over the primes until we find it. + if (countApprox < n) { - if (count < n) - { - checkLimit(start); - dist = nthPrimeDist(n, count, start); - stop = checkedAdd(start, dist); - count += countPrimes(start, stop); - start = checkedAdd(stop, 1); - } - if (sieveBackwards(n, count, stop)) + start = checkedAdd(start, 1); + uint64_t dist = (n - countApprox) * avgPrimeGap(primeApprox); + uint64_t stop = checkedAdd(start, dist); + primesieve::iterator iter(start, stop); + for (int64_t i = countApprox; i < n; i++) + prime = iter.next_prime(); + } + else // if (countApprox >= n) + { + uint64_t dist = (countApprox - n) * avgPrimeGap(primeApprox); + uint64_t stop = checkedSub(start, dist); + primesieve::iterator iter(start, stop); + for (int64_t i = countApprox; i >= n; i--) { - checkLowerLimit(stop); - dist = nthPrimeDist(n, count, stop); - start = checkedSub(start, dist); - count -= (int64_t) countPrimes(start, stop); - stop = checkedSub(start, 1); + prime = iter.prev_prime(); + if_unlikely(prime == 0) + throw primesieve_error("nth_prime(n): invalid n, nth prime < 2 is impossible!"); } } - if (n < 0) - count -= 1; + auto t2 = std::chrono::system_clock::now(); + std::chrono::duration seconds = t2 - t1; + seconds_ = seconds.count(); - // here start < nth prime, - // hence we can sieve forward the remaining - // distance and find the nth prime - ASSERT(count < n); + return prime; +} + +/// Used for n < 0 +uint64_t PrimeSieve::negativeNthPrime(int64_t n, uint64_t start) +{ + ASSERT(n < 0); + n = -n; - checkLimit(start); - dist = nthPrimeDist(n, count, start) * 2; - stop = checkedAdd(start, dist); + if ((uint64_t) n >= start) + throw primesieve_error("nth_prime(n): abs(n) must be < start"); + else if ((uint64_t) n > max_n) + throw primesieve_error("nth_prime(n): abs(n) must be <= " + std::to_string(max_n)); + + setStart(start); + auto t1 = std::chrono::system_clock::now(); + uint64_t nApprox = checkedSub(primesApprox(start), n); + nApprox = std::min(nApprox, max_n); + uint64_t primeApprox = nthPrimeApprox(nApprox); + primeApprox = std::min(primeApprox, start); uint64_t prime = 0; + int64_t countApprox = 0; + + // Only use multi-threading if the sieving distance is sufficiently + // large. For small n this if statement also avoids calling + // countPrimes() and hence the initialization overhead of + // O(x^0.5 log log x^0.5) occurs only once (instead of twice) when + // using primesieve::iterator further down. + if (start - primeApprox > isqrt(start) / 10) + { + // Count primes < start + start = checkedSub(start, 1); + primeApprox = std::min(primeApprox, start); + countApprox = countPrimes(primeApprox, start); + start = primeApprox; + } - for (primesieve::iterator it(start, stop); count < n; count++) - prime = it.next_prime(); + // Here we are very close to the nth prime < sqrt(nth_prime), + // we simply iterate over the primes until we find it. + if (countApprox >= n) + { + uint64_t dist = (countApprox - n) * avgPrimeGap(start); + uint64_t stop = checkedAdd(start, dist); + primesieve::iterator iter(start, stop); + for (int64_t i = countApprox; i >= n; i--) + prime = iter.next_prime(); + } + else // if (countApprox < n) + { + start = checkedSub(start, 1); + uint64_t dist = (n - countApprox) * avgPrimeGap(start); + uint64_t stop = checkedSub(start, dist); + primesieve::iterator iter(start, stop); + for (int64_t i = countApprox; i < n; i++) + { + prime = iter.prev_prime(); + if_unlikely(prime == 0) + throw primesieve_error("nth_prime(n): invalid n, nth prime < 2 is impossible!"); + } + } auto t2 = std::chrono::system_clock::now(); std::chrono::duration seconds = t2 - t1; diff --git a/lib/primesieve/src/nthPrimeApprox.cpp b/lib/primesieve/src/nthPrimeApprox.cpp new file mode 100644 index 00000000..b55f214b --- /dev/null +++ b/lib/primesieve/src/nthPrimeApprox.cpp @@ -0,0 +1,277 @@ +/// +/// @file nthPrimeApprox.cpp +/// This file contains implementations of the logarithmic +/// integral and the Riemann R function which are very +/// accurate approximations of PrimePi(x). Please note that +/// most of this code has been copied from the primecount +/// project: https://github.com/kimwalisch/primecount +/// +/// Note that while the Riemann R function is extremely +/// accurate it is much slower than other simpler PrimePi(x) +/// approximations. When speed matters, e.g. for allocating +/// a vector of primes, we avoid using the functions defined +/// in this file. Currently, the functions defined in this +/// file are only used in nthPrime.cpp where accuracy is of +/// utmost importance. +/// +/// Copyright (C) 2024 Kim Walisch, +/// +/// This file is distributed under the BSD License. See the COPYING +/// file in the top level directory. +/// + +#include +#include + +#include +#include +#include + +namespace primesieve { + +/// Generate a vector with Möbius function values. +/// This implementation is based on code by Rick Sladkey: +/// https://mathoverflow.net/q/99545 +/// +Vector generate_moebius(int64_t max) +{ + int64_t sqrt = isqrt(max); + int64_t size = max + 1; + Vector mu(size); + std::fill(mu.begin(), mu.end(), 1); + + for (int64_t i = 2; i <= sqrt; i++) + { + if (mu[i] == 1) + { + for (int64_t j = i; j < size; j += i) + mu[j] *= (int32_t) -i; + for (int64_t j = i * i; j < size; j += i * i) + mu[j] = 0; + } + } + + for (int64_t i = 2; i < size; i++) + { + if (mu[i] == i) + mu[i] = 1; + else if (mu[i] == -i) + mu[i] = -1; + else if (mu[i] < 0) + mu[i] = 1; + else if (mu[i] > 0) + mu[i] = -1; + } + + return mu; +} + +} // namespace + +namespace { + +/// Calculate the logarithmic integral using +/// Ramanujan's formula: +/// https://en.wikipedia.org/wiki/Logarithmic_integral_function#Series_representation +/// +long double li(long double x) +{ + if (x <= 1) + return 0; + + long double gamma = 0.577215664901532860606512090082402431L; + long double sum = 0; + long double inner_sum = 0; + long double factorial = 1; + long double p = -1; + long double q = 0; + long double power2 = 1; + long double logx = std::log(x); + int k = 0; + + for (int n = 1; true; n++) + { + p *= -logx; + factorial *= n; + q = factorial * power2; + power2 *= 2; + + for (; k <= (n - 1) / 2; k++) + inner_sum += 1.0L / (2 * k + 1); + + auto old_sum = sum; + sum += (p / q) * inner_sum; + + // Not converging anymore + if (std::abs(sum - old_sum) < std::numeric_limits::epsilon()) + break; + } + + return gamma + std::log(logx) + std::sqrt(x) * sum; +} + +/// Calculate the offset logarithmic integral which is a very +/// accurate approximation of the number of primes <= x. +/// Li(x) > pi(x) for 24 <= x <= ~ 10^316 +/// +long double Li(long double x) +{ + long double li2 = 1.045163780117492784844588889194613136L; + + if (x <= li2) + return 0; + else + return li(x) - li2; +} + +/// Calculate the inverse offset logarithmic integral which +/// is a very accurate approximation of the nth prime. +/// Li^-1(x) < nth_prime(x) for 7 <= x <= 10^316 +/// +/// This implementation computes Li^-1(x) as the zero of the +/// function f(z) = Li(z) - x using the Newton–Raphson method. +/// Note that Li'(z) = 1 / log(z). +/// https://math.stackexchange.com/a/853192 +/// +/// Newton–Raphson method: +/// zn+1 = zn - (f(zn) / f'(zn)). +/// zn+1 = zn - (Li(zn) - x) / (1 / log(zn)) +/// zn+1 = zn - (Li(zn) - x) * log(zn) +/// +long double Li_inverse(long double x) +{ + if (x < 2) + return 0; + + long double t = x * std::log(x); + long double old_term = std::numeric_limits::infinity(); + + while (true) + { + long double term = (Li(t) - x) * std::log(t); + + // Not converging anymore + if (std::abs(term) >= std::abs(old_term)) + break; + + t -= term; + old_term = term; + } + + return t; +} + +/// Calculate the Riemann R function which is a very accurate +/// approximation of the number of primes below x. +/// RiemannR(x) = \sum_{n=1}^{∞} μ(n)/n * li(x^(1/n)) +/// http://mathworld.wolfram.com/RiemannPrimeCountingFunction.html +/// +long double Ri(long double x) +{ + if (x <= 1) + return 0; + + long double sum = 0; + long double old_term = std::numeric_limits::infinity(); + auto terms = (int) (std::log2(x) * 2 + 10); + auto mu = primesieve::generate_moebius(terms); + + for (int n = 1; n < terms; n++) + { + if (mu[n]) + { + long double root = std::pow(x, 1.0L / n); + long double term = (li(root) * mu[n]) / n; + + // Not converging anymore + if (std::abs(term) >= std::abs(old_term)) + break; + + sum += term; + old_term = term; + } + } + + return sum; +} + +/// Calculate the inverse Riemann R function which is a very +/// accurate approximation of the nth prime. +/// This implementation computes Ri^-1(x) as the zero of the +/// function f(z) = Ri(z) - x using the Newton–Raphson method. +/// Note that Ri'(z) = 1 / log(z). +/// https://math.stackexchange.com/a/853192 +/// +/// Newton–Raphson method: +/// zn+1 = zn - (f(zn) / f'(zn)). +/// zn+1 = zn - (Ri(zn) - x) / (1 / log(zn)) +/// zn+1 = zn - (Ri(zn) - x) * log(zn) +/// +long double Ri_inverse(long double x) +{ + if (x < 2) + return 0; + + long double t = Li_inverse(x); + long double old_term = std::numeric_limits::infinity(); + + while (true) + { + long double term = (Ri(t) - x) * std::log(t); + + // Not converging anymore + if (std::abs(term) >= std::abs(old_term)) + break; + + t -= term; + old_term = term; + } + + return t; +} + +} // namespace + +namespace primesieve { + +uint64_t Li(uint64_t x) +{ + return (uint64_t) ::Li((long double) x); +} + +uint64_t Li_inverse(uint64_t x) +{ + return (uint64_t) ::Li_inverse((long double) x); +} + +uint64_t Ri(uint64_t x) +{ + return (uint64_t) ::Ri((long double) x); +} + +uint64_t Ri_inverse(uint64_t x) +{ + return (uint64_t) ::Ri_inverse((long double) x); +} + +uint64_t primesApprox(uint64_t x) +{ + // Li(x) is faster but less accurate than Ri(x). + // For small n speed is more important than accuracy. + if (x < 1e10) + return Li(x); + else + return Ri(x); +} + +uint64_t nthPrimeApprox(uint64_t n) +{ + // Li_inverse(x) is faster but less accurate than Ri_inverse(x). + // For small n speed is more important than accuracy. + if (n < 1e8) + return Li_inverse(n); + else + return Ri_inverse(n); +} + +} // namespace diff --git a/lib/primesieve/test/Li.cpp b/lib/primesieve/test/Li.cpp new file mode 100644 index 00000000..f1f674a6 --- /dev/null +++ b/lib/primesieve/test/Li.cpp @@ -0,0 +1,101 @@ +/// +/// @file Li.cpp +/// @brief Test the offset logarithmic integral function. +/// Li(x) = li(x) - li(2) +/// +/// Copyright (C) 2024 Kim Walisch, +/// +/// This file is distributed under the BSD License. See the COPYING +/// file in the top level directory. +/// + +#include + +#include +#include +#include +#include +#include + +using std::max; +using std::size_t; +using namespace primesieve; + +std::vector Li_table = +{ + 5, // Li(10^1) + 29, // Li(10^2) + 176, // Li(10^3) + 1245, // Li(10^4) + 9628, // Li(10^5) + 78626, // Li(10^6) + 664917, // Li(10^7) + 5762208, // Li(10^8) + 50849233, // Li(10^9) + 455055613, // Li(10^10) + 4118066399ll, // Li(10^11) + 37607950279ll, // Li(10^12) + 346065645809ll, // Li(10^13) + 3204942065690ll, // Li(10^14) + 29844571475286ll // Li(10^15) +}; + +void check(bool OK) +{ + std::cout << " " << (OK ? "OK" : "ERROR") << "\n"; + if (!OK) + std::exit(1); +} + +int main() +{ + uint64_t x = 1; + for (size_t i = 0; i < Li_table.size(); i++) + { + x *= 10; + std::cout << "Li(" << x << ") = " << Li(x); + check(Li(x) == Li_table[i]); + } + + x = 1; + for (size_t i = 0; i < Li_table.size(); i++) + { + x *= 10; + std::cout << "Li_inverse(" << Li_table[i] << ") = " << Li_inverse(Li_table[i]); + check(Li_inverse(Li_table[i]) <= x && + Li_inverse(Li_table[i] + 1) > x); + } + + // Sanity checks for small values of Li(x) + for (x = 0; x < 300000; x++) + { + uint64_t lix = Li(x); + double logx = std::log(max((double) x, 2.0)); + + if ((x >= 11 && lix < x / logx) || + (x >= 2 && lix > x * logx)) + { + std::cout << "Li(" << x << ") = " << lix << " ERROR" << std::endl; + std::exit(1); + } + } + + // Sanity checks for small values of Li_inverse(x) + for (x = 2; x < 30000; x++) + { + uint64_t res = Li_inverse(x); + double logx = std::log((double) x); + + if (res < x || + (x >= 4 && res > x * logx * logx)) + { + std::cout << "Li_inverse(" << x << ") = " << res << " ERROR" << std::endl; + std::exit(1); + } + } + + std::cout << std::endl; + std::cout << "All tests passed successfully!" << std::endl; + + return 0; +} diff --git a/lib/primesieve/test/Ri.cpp b/lib/primesieve/test/Ri.cpp new file mode 100644 index 00000000..ad71c321 --- /dev/null +++ b/lib/primesieve/test/Ri.cpp @@ -0,0 +1,127 @@ +/// +/// @file Ri.cpp +/// @brief Test the Riemann R function. +/// +/// Copyright (C) 2024 Kim Walisch, +/// +/// This file is distributed under the BSD License. See the COPYING +/// file in the top level directory. +/// + +#include + +#include +#include +#include +#include +#include + +using std::max; +using std::size_t; +using namespace primesieve; + +std::vector Ri_table = +{ + 4, // Ri(10^1) + 25, // Ri(10^2) + 168, // Ri(10^3) + 1226, // Ri(10^4) + 9587, // Ri(10^5) + 78527, // Ri(10^6) + 664667, // Ri(10^7) + 5761551, // Ri(10^8) + 50847455, // Ri(10^9) + 455050683, // Ri(10^10) + 4118052494ll, // Ri(10^11) + 37607910542ll, // Ri(10^12) + 346065531065ll, // Ri(10^13) + 3204941731601ll // Ri(10^14) +}; + +void check(bool OK) +{ + std::cout << " " << (OK ? "OK" : "ERROR") << "\n"; + if (!OK) + std::exit(1); +} + +int main() +{ + uint64_t x = 1; + for (size_t i = 0; i < Ri_table.size(); i++) + { + x *= 10; + std::cout << "Ri(" << x << ") = " << Ri(x); + check(Ri(x) == Ri_table[i]); + } + + x = 1; + for (size_t i = 0; i < Ri_table.size(); i++) + { + x *= 10; + std::cout << "Ri_inverse(" << Ri_table[i] << ") = " << Ri_inverse(Ri_table[i]); + check(Ri_inverse(Ri_table[i]) < x && + Ri_inverse(Ri_table[i] + 1) >= x); + } + + // Sanity checks for tiny values of Ri(x) + for (x = 0; x < 10000; x++) + { + uint64_t rix = Ri(x); + double logx = std::log(max((double) x, 2.0)); + + if ((x >= 20 && rix < x / logx) || + (x >= 2 && rix > x * logx)) + { + std::cout << "Ri(" << x << ") = " << rix << " ERROR" << std::endl; + std::exit(1); + } + } + + // Sanity checks for small values of Ri(x) + for (; x < 100000; x += 101) + { + uint64_t rix = Ri(x); + double logx = std::log(max((double) x, 2.0)); + + if ((x >= 20 && rix < x / logx) || + (x >= 2 && rix > x * logx)) + { + std::cout << "Ri(" << x << ") = " << rix << " ERROR" << std::endl; + std::exit(1); + } + } + + // Sanity checks for tiny values of Ri_inverse(x) + for (x = 2; x < 1000; x++) + { + uint64_t res = Ri_inverse(x); + double logx = std::log((double) x); + + if (res < x || + (x >= 5 && res > x * logx * logx)) + { + std::cout << "Ri_inverse(" << x << ") = " << res << " ERROR" << std::endl; + std::exit(1); + } + } + + // Sanity checks for small values of Ri_inverse(x) + for (; x < 100000; x += 101) + { + uint64_t res = Ri_inverse(x); + double logx = std::log((double) x); + + if (res < x || + (x >= 5 && res > x * logx * logx)) + { + std::cout << "Ri_inverse(" << x << ") = " << res << " ERROR" << std::endl; + std::exit(1); + } + } + + std::cout << std::endl; + std::cout << "All tests passed successfully!" << std::endl; + + return 0; +} diff --git a/lib/primesieve/test/moebius.cpp b/lib/primesieve/test/moebius.cpp new file mode 100644 index 00000000..579ac0d9 --- /dev/null +++ b/lib/primesieve/test/moebius.cpp @@ -0,0 +1,99 @@ +/// +/// @file moebius.cpp +/// @brief Test the generate_moebius(n) function +/// @link https://en.wikipedia.org/wiki/Moebius_function +/// +/// Copyright (C) 2024 Kim Walisch, +/// +/// This file is distributed under the BSD License. See the COPYING +/// file in the top level directory. +/// + +#include + +#include +#include +#include +#include + +using std::size_t; +using namespace primesieve; + +/// Values of mu(n) for the first 1000 integers +/// https://oeis.org/A008683/b008683.txt +std::vector moebius = +{ + 0, 1, -1, -1, 0, -1, 1, -1, 0, 0, 1, -1, 0, -1, 1, 1, 0, -1, 0, -1, + 0, 1, 1, -1, 0, 0, 1, 0, 0, -1, -1, -1, 0, 1, 1, 1, 0, -1, 1, 1, + 0, -1, -1, -1, 0, 0, 1, -1, 0, 0, 0, 1, 0, -1, 0, 1, 0, 1, 1, -1, + 0, -1, 1, 0, 0, 1, -1, -1, 0, 1, -1, -1, 0, -1, 1, 0, 0, 1, -1, -1, + 0, 0, 1, -1, 0, 1, 1, 1, 0, -1, 0, 1, 0, 1, 1, 1, 0, -1, 0, 0, + 0, -1, -1, -1, 0, -1, 1, -1, 0, -1, -1, 1, 0, -1, -1, 1, 0, 0, 1, 1, + 0, 0, 1, 1, 0, 0, 0, -1, 0, 1, -1, -1, 0, 1, 1, 0, 0, -1, -1, -1, + 0, 1, 1, 1, 0, 1, 1, 0, 0, -1, 0, -1, 0, 0, -1, 1, 0, -1, 1, 1, + 0, 1, 0, -1, 0, -1, 1, -1, 0, 0, -1, 0, 0, -1, -1, 0, 0, 1, 1, -1, + 0, -1, -1, 1, 0, 1, -1, 1, 0, 0, -1, -1, 0, -1, 1, -1, 0, -1, 0, -1, + 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, -1, 0, 1, 1, 1, 0, 1, 1, 1, + 0, 1, -1, -1, 0, 0, 1, -1, 0, -1, -1, -1, 0, -1, 0, 1, 0, 1, -1, -1, + 0, -1, 0, 0, 0, 0, -1, 1, 0, 1, 0, -1, 0, 1, 1, -1, 0, -1, -1, 1, + 0, 0, 1, -1, 0, 1, -1, 1, 0, -1, 0, -1, 0, -1, 1, 0, 0, -1, 1, 0, + 0, -1, -1, -1, 0, -1, -1, 1, 0, 0, -1, 1, 0, -1, 0, 1, 0, 0, 1, 1, + 0, 1, 1, 1, 0, 1, 0, -1, 0, 1, -1, -1, 0, -1, 1, 0, 0, -1, -1, 1, + 0, 1, -1, 1, 0, 0, 1, 1, 0, 1, 1, -1, 0, 0, 1, 1, 0, -1, 0, 1, + 0, 1, 0, 0, 0, -1, 1, -1, 0, -1, 0, 0, 0, -1, -1, 1, 0, -1, 1, -1, + 0, 0, 1, 0, 0, 1, -1, -1, 0, 0, -1, 1, 0, -1, -1, 0, 0, 1, 0, -1, + 0, 1, 1, -1, 0, -1, 1, 0, 0, -1, 1, 1, 0, 1, 1, 1, 0, -1, 1, -1, + 0, -1, -1, 1, 0, 0, -1, 1, 0, -1, -1, 1, 0, 1, 0, 1, 0, 1, -1, -1, + 0, -1, 1, 0, 0, 0, -1, 1, 0, -1, -1, -1, 0, -1, -1, -1, 0, 1, -1, -1, + 0, 0, -1, -1, 0, 1, 1, 1, 0, -1, 0, 1, 0, 1, 1, -1, 0, -1, 1, 0, + 0, -1, 1, -1, 0, -1, 1, -1, 0, 1, -1, 1, 0, 1, -1, 0, 0, 0, 1, -1, + 0, 1, 1, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, -1, 0, 0, 1, -1, -1, + 0, 1, 1, -1, 0, 1, -1, 0, 0, -1, 1, 1, 0, 0, 1, 1, 0, 1, -1, 1, + 0, -1, 0, -1, 0, 0, 1, 1, 0, 0, -1, 0, 0, 1, -1, 1, 0, 1, 1, 0, + 0, -1, 1, 1, 0, 1, 1, -1, 0, 0, 0, 1, 0, 1, 1, -1, 0, -1, 0, 1, + 0, -1, 1, -1, 0, 1, 1, 0, 0, -1, 1, -1, 0, 1, -1, 0, 0, -1, 0, 1, + 0, 1, -1, 1, 0, 0, 1, -1, 0, 1, -1, 1, 0, -1, 0, -1, 0, 1, -1, -1, + 0, -1, -1, 0, 0, 0, -1, -1, 0, -1, -1, 1, 0, -1, 1, -1, 0, -1, -1, -1, + 0, 0, 1, 1, 0, 0, 1, -1, 0, 1, 0, -1, 0, 1, 1, 1, 0, 0, -1, 0, + 0, -1, -1, -1, 0, -1, -1, -1, 0, 1, 0, -1, 0, -1, -1, 1, 0, 0, -1, -1, + 0, -1, 1, -1, 0, -1, 0, 1, 0, 1, -1, 1, 0, -1, 1, 0, 0, -1, -1, 1, + 0, 1, -1, -1, 0, 1, 0, 1, 0, 1, 1, -1, 0, 0, 1, 1, 0, 1, 1, 1, + 0, -1, 0, 1, 0, -1, 1, 1, 0, -1, -1, 0, 0, 1, 1, -1, 0, 1, 1, -1, + 0, 1, 0, 1, 0, 0, 0, -1, 0, 0, -1, 1, 0, -1, 1, 0, 0, 1, 0, -1, + 0, -1, -1, -1, 0, 1, 1, 0, 0, 1, 0, -1, 0, 1, -1, 1, 0, -1, 1, -1, + 0, -1, -1, 1, 0, 0, 1, 1, 0, -1, 1, 1, 0, -1, 0, 0, 0, -1, 1, 1, + 0, 1, -1, 0, 0, 1, -1, -1, 0, 1, -1, 1, 0, 1, 1, -1, 0, -1, 1, 1, + 0, 0, 1, 1, 0, -1, -1, 1, 0, -1, 0, -1, 0, 1, -1, 1, 0, 1, 1, 0, + 0, -1, -1, -1, 0, 0, -1, -1, 0, -1, -1, 1, 0, 0, -1, 1, 0, 0, 1, -1, + 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, -1, -1, 0, 0, -1, 1, -1, + 0, -1, 1, -1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, -1, 0, 0, -1, 1, 1, + 0, -1, 0, -1, 0, -1, 1, -1, 0, 1, -1, 0, 0, 1, -1, 1, 0, -1, 1, 1, + 0, 1, -1, -1, 0, 1, -1, -1, 0, 0, 1, -1, 0, 1, 1, -1, 0, 1, 0, -1, + 0, 1, 1, 1, 0, 0, 1, 0, 0, -1, 1, 0, 0, 1, 1, -1, 0, -1, -1, 1, + 0, -1, -1, 1, 0, 0, -1, -1, 0, 1, 0, 1, 0, -1, 0, 1, 0, -1, 1, 1, + 0, 0, -1, 0, 0, 1, 1, -1, 0, -1, -1, -1, 0, 1, 1, 0, 0, -1, -1, 1, + 0, 0, 1, -1, 0, 1, -1, -1, 0, 1, 0, -1, 0, 1, -1, 1, 0, -1, 1, 0 +}; + +void check(bool OK) +{ + std::cout << " " << (OK ? "OK" : "ERROR") << "\n"; + if (!OK) + std::exit(1); +} + +int main() +{ + auto mu = generate_moebius(moebius.size() - 1); + + for (size_t i = 1; i < mu.size(); i++) + { + std::cout << "mu(" << i << ") = " << mu[i]; + check(mu[i] == moebius[i]); + } + + std::cout << std::endl; + std::cout << "All tests passed successfully!" << std::endl; + + return 0; +}