mirror of
https://github.com/kimwalisch/primecount.git
synced 2026-06-06 18:59:47 +09:00
276 lines
6.7 KiB
C++
276 lines
6.7 KiB
C++
///
|
|
/// @file Li.cpp
|
|
/// @brief Test the Eulerian logarithmic integral function.
|
|
/// Li(x) = li(x) - li(2)
|
|
///
|
|
/// Copyright (C) 2025 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-internal.hpp>
|
|
#include <calculator.hpp>
|
|
#include <int128_t.hpp>
|
|
#include <imath.hpp>
|
|
|
|
#include <stdint.h>
|
|
#include <iostream>
|
|
#include <cmath>
|
|
#include <cstdlib>
|
|
#include <string>
|
|
#include <vector>
|
|
#include <array>
|
|
|
|
using std::size_t;
|
|
using namespace primecount;
|
|
|
|
/// Generated using Mathematica:
|
|
/// Table[IntegerPart[li[k]-li[2]], {k, 2, 99}]
|
|
std::array<int64_t, 100> Li_tiny =
|
|
{
|
|
0, 0, 0, 1, 1, 2, 3, 3, 4, 4,
|
|
5, 5, 5, 6, 6, 7, 7, 7, 8, 8,
|
|
8, 9, 9, 9, 10, 10, 10, 11, 11, 11,
|
|
11, 12, 12, 12, 13, 13, 13, 13, 14, 14,
|
|
14, 15, 15, 15, 15, 16, 16, 16, 16, 17,
|
|
17, 17, 17, 18, 18, 18, 18, 19, 19, 19,
|
|
19, 20, 20, 20, 20, 21, 21, 21, 21, 22,
|
|
22, 22, 22, 23, 23, 23, 23, 23, 24, 24,
|
|
24, 24, 25, 25, 25, 25, 25, 26, 26, 26,
|
|
26, 27, 27, 27, 27, 27, 28, 28, 28, 28
|
|
};
|
|
|
|
std::vector<int64_t> 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)
|
|
};
|
|
|
|
#if defined(HAVE_FLOAT128)
|
|
|
|
std::vector<std::string> Li_f128 =
|
|
{
|
|
"29844571475286", // Li(10^15)
|
|
"279238344248555", // Li(10^16)
|
|
"2623557165610820", // Li(10^17)
|
|
"24739954309690413", // Li(10^18)
|
|
"234057667376222381", // Li(10^19)
|
|
"2220819602783663482", // Li(20^20)
|
|
"21127269486616126181", // Li(20^21)
|
|
"201467286691248261497", // Li(20^22)
|
|
"1925320391614054155137", // Li(20^23)
|
|
"18435599767366347775143", // Li(20^24)
|
|
"176846309399198930392618", // Li(20^25)
|
|
"1699246750872593033005722", // Li(20^26)
|
|
"16352460426842189113085404", // Li(20^27)
|
|
"157589269275974838158399970", // Li(20^28)
|
|
"1520698109714276717287880526", // Li(20^29)
|
|
"14692398897720447639079087668" // Li(20^30)
|
|
};
|
|
|
|
#endif
|
|
|
|
void check(bool OK)
|
|
{
|
|
std::cout << " " << (OK ? "OK" : "ERROR") << "\n";
|
|
if (!OK)
|
|
std::exit(1);
|
|
}
|
|
|
|
int main()
|
|
{
|
|
for (int64_t x = 0; x < (int64_t) Li_tiny.size(); x++)
|
|
{
|
|
std::cout << "Li(" << x << ") = " << Li(x);
|
|
check(Li(x) == Li_tiny[x]);
|
|
}
|
|
|
|
{
|
|
int64_t x = 10;
|
|
for (size_t i = 0; i < Li_table.size(); i++)
|
|
{
|
|
std::cout << "Li(" << x << ") = " << Li(x);
|
|
check(Li(x) == Li_table[i]);
|
|
x *= 10;
|
|
}
|
|
}
|
|
|
|
#if defined(HAVE_FLOAT128) && \
|
|
defined(HAVE_INT128_T)
|
|
|
|
{
|
|
int128_t x = ipow<15>((int128_t) 10);
|
|
for (size_t i = 0; i < Li_f128.size(); i++)
|
|
{
|
|
std::string str = to_string(Li(x));
|
|
std::cout << "Li(" << x << ") = " << str;
|
|
check(str == Li_f128[i]);
|
|
x *= 10;
|
|
}
|
|
}
|
|
|
|
#endif
|
|
|
|
for (int64_t x = 1; x < (int64_t) Li_tiny.size(); x++)
|
|
{
|
|
int64_t y = Li_tiny[x];
|
|
std::cout << "Li_inverse(" << y << ") = " << Li_inverse(y);
|
|
check(Li_inverse(y) < x &&
|
|
Li_inverse(y + 1) >= x);
|
|
}
|
|
|
|
{
|
|
int64_t x = 10;
|
|
for (size_t i = 0; i < Li_table.size(); i++)
|
|
{
|
|
int64_t y = Li_table[i];
|
|
std::cout << "Li_inverse(" << y << ") = " << Li_inverse(y);
|
|
check(Li_inverse(y) < x &&
|
|
Li_inverse(y + 1) >= x);
|
|
x *= 10;
|
|
}
|
|
}
|
|
|
|
#if defined(HAVE_FLOAT128) && \
|
|
defined(HAVE_INT128_T)
|
|
|
|
{
|
|
int128_t x = ipow<15>((int128_t) 10);
|
|
for (size_t i = 0; i < Li_f128.size(); i++)
|
|
{
|
|
int128_t y = calculator::eval<int128_t>(Li_f128[i]);
|
|
std::cout << "Li_inverse(" << y << ") = " << Li_inverse(y);
|
|
check(Li_inverse(y) < x &&
|
|
Li_inverse(y + 1) >= x);
|
|
x *= 10;
|
|
}
|
|
}
|
|
|
|
#endif
|
|
|
|
{
|
|
// Li(9760) = 1219.000098
|
|
int64_t x = 9760;
|
|
int64_t y = 1219;
|
|
std::cout << "Li(" << x << ") = " << Li(x);
|
|
check(Li(x) == y);
|
|
|
|
std::cout << "Li_inverse(" << y << ") = " << Li_inverse(y);
|
|
check(Li_inverse(y) < x &&
|
|
Li_inverse(y + 1) >= x);
|
|
}
|
|
|
|
{
|
|
// Li(9494) = 1189.9997
|
|
int64_t x = 9494;
|
|
int64_t y = 1189;
|
|
std::cout << "Li(" << x << ") = " << Li(x);
|
|
check(Li(x) == y);
|
|
|
|
std::cout << "Li_inverse(" << y << ") = " << Li_inverse(y);
|
|
check(Li_inverse(y) < x &&
|
|
Li_inverse(y + 1) >= x);
|
|
}
|
|
|
|
// Sanity checks for tiny values of Li(x)
|
|
int64_t x;
|
|
for (x = 0; x < 10000; x++)
|
|
{
|
|
int64_t Lix = Li(x);
|
|
double logx = std::log(std::max((double) x, 2.0));
|
|
|
|
if (Lix < 0 ||
|
|
(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(x)
|
|
for (; x < 100000; x += 101)
|
|
{
|
|
int64_t Lix = Li(x);
|
|
double logx = std::log(std::max((double) x, 2.0));
|
|
|
|
if (Lix < 0 ||
|
|
(x >= 11 && Lix < x / logx) ||
|
|
(x >= 2 && Lix > x * logx))
|
|
{
|
|
std::cout << "Li(" << x << ") = " << Lix << " ERROR" << std::endl;
|
|
std::exit(1);
|
|
}
|
|
}
|
|
|
|
// Sanity checks for tiny values of Li_inverse(x)
|
|
for (x = 2; x < 1000; x++)
|
|
{
|
|
int64_t res = Li_inverse(x);
|
|
double logx = std::log((double) x);
|
|
|
|
if (res < 0 ||
|
|
res < x ||
|
|
(x >= 4 && res > x * logx * logx))
|
|
{
|
|
std::cout << "Li_inverse(" << x << ") = " << res << " ERROR" << std::endl;
|
|
std::exit(1);
|
|
}
|
|
}
|
|
|
|
// Sanity checks for small values of Li_inverse(x)
|
|
for (; x < 100000; x += 101)
|
|
{
|
|
int64_t res = Li_inverse(x);
|
|
double logx = std::log((double) x);
|
|
|
|
if (res < 0 ||
|
|
res < x ||
|
|
(x >= 4 && res > x * logx * logx))
|
|
{
|
|
std::cout << "Li_inverse(" << x << ") = " << res << " ERROR" << std::endl;
|
|
std::exit(1);
|
|
}
|
|
}
|
|
|
|
{
|
|
int64_t x = pstd::numeric_limits<int64_t>::max() / 10;
|
|
int64_t res = Li_inverse(x);
|
|
if (res != pstd::numeric_limits<int64_t>::max())
|
|
{
|
|
std::cout << "Li_inverse(" << x << ") != INT64_MAX, failed to prevent integer overflow!" << std::endl;
|
|
std::exit(1);
|
|
}
|
|
}
|
|
|
|
#if defined(HAVE_INT128_T)
|
|
{
|
|
int128_t x = pstd::numeric_limits<int128_t>::max() / 10;
|
|
int128_t res = Li_inverse(x);
|
|
if (res != pstd::numeric_limits<int128_t>::max())
|
|
{
|
|
std::cout << "Li_inverse(" << x << ") != INT128_MAX, failed to prevent integer overflow!" << std::endl;
|
|
std::exit(1);
|
|
}
|
|
}
|
|
#endif
|
|
|
|
std::cout << std::endl;
|
|
std::cout << "All tests passed successfully!" << std::endl;
|
|
|
|
return 0;
|
|
}
|