From 29051b9210473e2b785af5fe7a04bf4b0f722b6f Mon Sep 17 00:00:00 2001 From: nipzu Date: Fri, 9 Feb 2024 19:04:43 +0200 Subject: [PATCH] Use the Gram series for calculating Ri (#144) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit RiemannR(x) = 1 + \sum_{k=1}^{∞} ln(x)^k / (zeta(k + 1) * k * k!) --- include/primesieve/forward.hpp | 1 + include/primesieve/nthPrimeApprox.hpp | 3 - src/LookupTables.cpp | 136 +++++++++++++++ src/nthPrimeApprox.cpp | 234 +++++++------------------- test/Li.cpp | 112 ------------ test/moebius.cpp | 99 ----------- 6 files changed, 198 insertions(+), 387 deletions(-) delete mode 100644 test/Li.cpp delete mode 100644 test/moebius.cpp diff --git a/include/primesieve/forward.hpp b/include/primesieve/forward.hpp index 7538e2a76..ddb25dd64 100644 --- a/include/primesieve/forward.hpp +++ b/include/primesieve/forward.hpp @@ -18,6 +18,7 @@ namespace primesieve { extern const Array bitValues; extern const Array bruijnBitValues; +extern const Array zetaInv; int get_num_threads(); int get_sieve_size(); diff --git a/include/primesieve/nthPrimeApprox.hpp b/include/primesieve/nthPrimeApprox.hpp index 94d27e132..a8488e67d 100644 --- a/include/primesieve/nthPrimeApprox.hpp +++ b/include/primesieve/nthPrimeApprox.hpp @@ -12,9 +12,6 @@ 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 primePiApprox(uint64_t x); diff --git a/src/LookupTables.cpp b/src/LookupTables.cpp index fffdd7265..ce2009ba1 100644 --- a/src/LookupTables.cpp +++ b/src/LookupTables.cpp @@ -107,4 +107,140 @@ const WheelInit wheel210Init[210] = { 1, 47 }, { 0, 47 } }; +/// Precomputed values of zetaInv[k] = 1 / zeta(k + 1). +/// Used in the calculation of the Riemann R function and its derivative. +/// These values are calculated up to a precision of 128 bits. +/// +const Array zetaInv = +{ + 0.00000000000000000000000000000000000000L, + 0.60792710185402662866327677925836583343L, + 0.83190737258070746868312627882153073442L, + 0.92393840292159016702375049404068247277L, + 0.96438734042926245912643658844498457124L, + 0.98295259226458041980489656499392413295L, + 0.99171985583844431042818593149755069165L, + 0.99593920112551514683483647280554532401L, + 0.99799563273076215686467613210509999632L, + 0.99900641306903078175222531809290878576L, + 0.99950605549762467678582298009453697739L, + 0.99975397399038468206770164303673058471L, + 0.99987730170913952450133620378723676486L, + 0.99993875561604559519730175829325041151L, + 0.99996941269930456117242340889418930870L, + 0.99998471797413523168338429153076556171L, + 0.99999236286068844254826438672217335846L, + 0.99999618272130667240681034896437598819L, + 0.99999809179092471488437703328457254051L, + 0.99999904603887616989781139055273723228L, + 0.99999952306724067715893767591651988929L, + 0.99999976154955413089570313635363076685L, + 0.99999988078008824807797958727410486625L, + 0.99999994039181450187651056224846508482L, + 0.99999997019649737359651508890102089892L, + 0.99999998509844539369129175912634580637L, + 0.99999999254928826567767626959951349247L, + 0.99999999627466598908965648972675175642L, + 0.99999999813734027995645223273129800534L, + 0.99999999906867256844770260269784066458L, + 0.99999999953433709371346353509327556591L, + 0.99999999976716881668655981073561305391L, + 0.99999999988358449828654737128147523751L, + 0.99999999994179227912436112987729664108L, + 0.99999999997089614955587603724140433426L, + 0.99999999998544807810916977419479027320L, + 0.99999999999272404016499545857700011549L, + 0.99999999999636202045263458370494986090L, + 0.99999999999818101034969624277576308610L, + 0.99999999999909050521597443825250885864L, + 0.99999999999954525262169599139249816326L, + 0.99999999999977262631541758644726975846L, + 0.99999999999988631315923199013976285333L, + 0.99999999999994315658012372737508161143L, + 0.99999999999997157829023110778923903545L, + 0.99999999999998578914517196859517856066L, + 0.99999999999999289457260478919777422111L, + 0.99999999999999644728630866289894847613L, + 0.99999999999999822364315642088282797029L, + 0.99999999999999911182157890691919855133L, + 0.99999999999999955591078968561886079545L, + 0.99999999999999977795539492019585090391L, + 0.99999999999999988897769748589339895390L, + 0.99999999999999994448884875154519064428L, + 0.99999999999999997224442437863875904455L, + 0.99999999999999998612221219027476742979L, + 0.99999999999999999306110609545584635070L, + 0.99999999999999999653055304783407738733L, + 0.99999999999999999826527652395242343096L, + 0.99999999999999999913263826198800662795L, + 0.99999999999999999956631913099793495144L, + 0.99999999999999999978315956550027802158L, + 0.99999999999999999989157978275057585938L, + 0.99999999999999999994578989137543354593L, + 0.99999999999999999997289494568776531168L, + 0.99999999999999999998644747284389883546L, + 0.99999999999999999999322373642195481090L, + 0.99999999999999999999661186821097920322L, + 0.99999999999999999999830593410549020083L, + 0.99999999999999999999915296705274530021L, + 0.99999999999999999999957648352637271666L, + 0.99999999999999999999978824176318638057L, + 0.99999999999999999999989412088159319766L, + 0.99999999999999999999994706044079660134L, + 0.99999999999999999999997353022039830147L, + 0.99999999999999999999998676511019915106L, + 0.99999999999999999999999338255509957560L, + 0.99999999999999999999999669127754978787L, + 0.99999999999999999999999834563877489392L, + 0.99999999999999999999999917281938744701L, + 0.99999999999999999999999958640969372348L, + 0.99999999999999999999999979320484686179L, + 0.99999999999999999999999989660242343087L, + 0.99999999999999999999999994830121171549L, + 0.99999999999999999999999997415060585772L, + 0.99999999999999999999999998707530292891L, + 0.99999999999999999999999999353765146443L, + 0.99999999999999999999999999676882573227L, + 0.99999999999999999999999999838441286611L, + 0.99999999999999999999999999919220643311L, + 0.99999999999999999999999999959610321653L, + 0.99999999999999999999999999979805160832L, + 0.99999999999999999999999999989902580413L, + 0.99999999999999999999999999994951290212L, + 0.99999999999999999999999999997475645103L, + 0.99999999999999999999999999998737822558L, + 0.99999999999999999999999999999368911276L, + 0.99999999999999999999999999999684455644L, + 0.99999999999999999999999999999842227819L, + 0.99999999999999999999999999999921113916L, + 0.99999999999999999999999999999960556955L, + 0.99999999999999999999999999999980278483L, + 0.99999999999999999999999999999990139239L, + 0.99999999999999999999999999999995069626L, + 0.99999999999999999999999999999997534810L, + 0.99999999999999999999999999999998767411L, + 0.99999999999999999999999999999999383702L, + 0.99999999999999999999999999999999691858L, + 0.99999999999999999999999999999999845926L, + 0.99999999999999999999999999999999922970L, + 0.99999999999999999999999999999999961481L, + 0.99999999999999999999999999999999980747L, + 0.99999999999999999999999999999999990370L, + 0.99999999999999999999999999999999995192L, + 0.99999999999999999999999999999999997593L, + 0.99999999999999999999999999999999998803L, + 0.99999999999999999999999999999999999398L, + 0.99999999999999999999999999999999999706L, + 0.99999999999999999999999999999999999850L, + 0.99999999999999999999999999999999999932L, + 0.99999999999999999999999999999999999962L, + 0.99999999999999999999999999999999999989L, + 0.99999999999999999999999999999999999991L, + 1.00000000000000000000000000000000000000L, + 1.00000000000000000000000000000000000000L, + 1.00000000000000000000000000000000000000L, + 1.00000000000000000000000000000000000000L, + 1.00000000000000000000000000000000000000L +}; + } // namespace diff --git a/src/nthPrimeApprox.cpp b/src/nthPrimeApprox.cpp index 757acb5d7..152c32810 100644 --- a/src/nthPrimeApprox.cpp +++ b/src/nthPrimeApprox.cpp @@ -22,203 +22,117 @@ #include #include +#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 +/// Calculate the Riemann R function which is a very accurate +/// approximation of the number of primes below x. +/// http://mathworld.wolfram.com/RiemannPrimeCountingFunction.html +/// The calculation is done with the Gram series: +/// RiemannR(x) = 1 + \sum_{k=1}^{∞} ln(x)^k / (zeta(k + 1) * k * k!) /// -long double li(long double x) +long double Ri(long double x) { - if (x <= 1) + if (x < 0.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 epsilon = std::numeric_limits::epsilon(); + long double old_sum = -1; + long double sum = 1; + long double term = 1; long double logx = std::log(x); - int k = 0; - for (int n = 1; true; n++) + for (int k = 1; k < 128 && std::abs(old_sum - sum) >= epsilon; k++) { - 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; + long double k_inv = 1.0L / (long double)k; + term *= logx * k_inv; + old_sum = sum; + sum += term * k_inv * primesieve::zetaInv[k]; } - 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) + // For k >= 128, approximate zeta(k + 1) by 1. + for (int k = 128; std::abs(old_sum - sum) >= epsilon; k++) { - 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; + long double k_inv = 1.0L / (long double)k; + term *= logx * k_inv; + old_sum = sum; + sum += term * k_inv; } - return t; + return sum; } -/// 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 +/// Calculate the derivative of the Riemann R function. +/// RiemannR'(x) = 1/x * \sum_{k=1}^{∞} ln(x)^(k-1) / (zeta(k + 1) * k!) /// -long double Ri(long double x) +long double Ri_prime(long double x) { - if (x <= 1) + if (x < 0.1) return 0; + long double epsilon = std::numeric_limits::epsilon(); + long double old_sum = -1; 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); + long double term = 1; + long double logx = std::log(x); - for (int n = 1; n < terms; n++) + for (int k = 1; k < 128 && std::abs(old_sum - sum) >= epsilon; k++) { - 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; + long double k_inv = 1.0L / (long double)k; + term *= logx * k_inv; + old_sum = sum; + sum += term * primesieve::zetaInv[k]; + } - sum += term; - old_term = term; - } + // For k >= 128, approximate zeta(k + 1) by 1. + for (int k = 128; std::abs(old_sum - sum) >= epsilon; k++) + { + long double k_inv = 1.0L / (long double)k; + term *= logx * k_inv; + old_sum = sum; + sum += term; } - return sum; + return sum / (x * logx); } /// 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) +/// zn+1 = zn - (Ri(zn) - x) / Ri'(zn) /// long double Ri_inverse(long double x) { if (x < 2) return 0; - long double t = Li_inverse(x); + long double logx = std::log(x); + long double loglogx = std::log(logx); + + // Calculate an initial approximation for the inverse + long double t = logx + 0.5L * loglogx; + if (x > 1600) + t += 0.5L * loglogx - 1.0L + (loglogx - 2.0L) / logx; + if (x > 1200000) + t -= (loglogx * loglogx - 6.0L * loglogx + 11.0L) / (2.0L * logx * logx); + t *= x; + long double old_term = std::numeric_limits::infinity(); while (true) { - long double term = (Ri(t) - x) * std::log(t); + long double term = (Ri(t) - x) / Ri_prime(t); // Not converging anymore if (std::abs(term) >= std::abs(old_term)) @@ -235,27 +149,11 @@ long double Ri_inverse(long double x) namespace primesieve { -uint64_t Li(uint64_t x) -{ - return (uint64_t) ::Li((long double) x); -} - uint64_t Ri(uint64_t x) { return (uint64_t) ::Ri((long double) x); } -uint64_t Li_inverse(uint64_t x) -{ - auto res = ::Li_inverse((long double) x); - - // Prevent 64-bit integer overflow - if (res > (long double) std::numeric_limits::max()) - return std::numeric_limits::max(); - - return (uint64_t) res; -} - uint64_t Ri_inverse(uint64_t x) { auto res = ::Ri_inverse((long double) x); @@ -275,12 +173,7 @@ uint64_t Ri_inverse(uint64_t x) /// uint64_t primePiApprox(uint64_t x) { - // Li(x) is faster but less accurate than Ri(x). - // For small x speed is more important than accuracy. - if (x < 1e10) - return Li(x); - else - return Ri(x); + return Ri(x); } /// nthPrimeApprox(n) is a very accurate approximation of the nth @@ -290,12 +183,7 @@ uint64_t primePiApprox(uint64_t x) /// uint64_t nthPrimeApprox(uint64_t n) { - // Li_inverse(n) is faster but less accurate than Ri_inverse(n). - // For small n speed is more important than accuracy. - if (n < 1e8) - return Li_inverse(n); - else - return Ri_inverse(n); + return Ri_inverse(n); } } // namespace diff --git a/test/Li.cpp b/test/Li.cpp deleted file mode 100644 index bccaa45d6..000000000 --- a/test/Li.cpp +++ /dev/null @@ -1,112 +0,0 @@ -/// -/// @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 -#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); - } - } - - { - uint64_t x = std::numeric_limits::max() / 10; - uint64_t res = Li_inverse(x); - if (res != std::numeric_limits::max()) - { - std::cout << "Li_inverse(" << x << ") != UINT64_MAX, failed to prevent integer overflow!" << std::endl; - std::exit(1); - } - } - - std::cout << std::endl; - std::cout << "All tests passed successfully!" << std::endl; - - return 0; -} diff --git a/test/moebius.cpp b/test/moebius.cpp deleted file mode 100644 index 579ac0d97..000000000 --- a/test/moebius.cpp +++ /dev/null @@ -1,99 +0,0 @@ -/// -/// @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; -}