Skip to content

Commit

Permalink
Improve accuracy
Browse files Browse the repository at this point in the history
  • Loading branch information
@kimwalisch
kimwalisch committed Feb 17, 2024
1 parent 3e60719 commit 42dd380
Show file tree
Hide file tree
Showing 2 changed files with 52 additions and 48 deletions.
50 changes: 26 additions & 24 deletions src/LogarithmicIntegral.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -226,67 +226,69 @@ namespace primecount {
int64_t Li(int64_t x)
{
#if defined(HAVE_FLOAT128)
if (x > 1 && std::log10(x) >= std::numeric_limits<long double>::digits10)
if (x > 1e10 && std::log10(x) >= std::numeric_limits<long double>::digits10)
return (int64_t) ::Li((__float128) x);
#endif

return (int64_t) ::Li((long double) x);
}

int64_t Li_inverse(int64_t x)
{
#if defined(HAVE_FLOAT128)
double safetyThreshold = 4;
if (x > 1 && std::log10(x) >= std::numeric_limits<long double>::digits10 - safetyThreshold)
if (x > 1e10)
{
__float128 res = ::Li_inverse((__float128) x);
if (res > (__float128) std::numeric_limits<int64_t>::max())
return std::numeric_limits<int64_t>::max();
return (int64_t) res;
double logx = std::log(x);
if (std::log10((double) x * (logx * logx)) >= std::numeric_limits<long double>::digits10)
{
__float128 res = ::Li_inverse((__float128) x);
if (res > (__float128) std::numeric_limits<int64_t>::max())
return std::numeric_limits<int64_t>::max();
else
return (int64_t) res;
}
}
#endif

long double res = ::Li_inverse((long double) x);

// Prevent integer overflow
if (res > (long double) std::numeric_limits<int64_t>::max())
return std::numeric_limits<int64_t>::max();

return (int64_t) res;
else
return (int64_t) res;
}

#ifdef HAVE_INT128_T

int128_t Li(int128_t x)
{
#if defined(HAVE_FLOAT128)
if (x > 1 && std::log10(x) >= std::numeric_limits<long double>::digits10)
if (x > 1e10 && std::log10(x) >= std::numeric_limits<long double>::digits10)
return (int128_t) ::Li((__float128) x);
#endif

return (int128_t) ::Li((long double) x);
}

int128_t Li_inverse(int128_t x)
{
#if defined(HAVE_FLOAT128)
double safetyThreshold = 4;
if (x > 1 && std::log10(x) >= std::numeric_limits<long double>::digits10 - safetyThreshold)
if (x > 1e10)
{
__float128 res = ::Li_inverse((__float128) x);
if (res > (__float128) std::numeric_limits<int128_t>::max())
return std::numeric_limits<int128_t>::max();
return (int128_t) res;
double logx = std::log(x);
if (std::log10((double) x * (logx * logx)) >= std::numeric_limits<long double>::digits10)
{
__float128 res = ::Li_inverse((__float128) x);
if (res > (__float128) std::numeric_limits<int128_t>::max())
return std::numeric_limits<int128_t>::max();
else
return (int128_t) res;
}
}
#endif

long double res = ::Li_inverse((long double) x);

// Prevent integer overflow
if (res > (long double) std::numeric_limits<int128_t>::max())
return std::numeric_limits<int128_t>::max();

return (int128_t) res;
else
return (int128_t) res;
}

#endif
Expand Down
50 changes: 26 additions & 24 deletions src/RiemannR.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -550,67 +550,69 @@ namespace primecount {
int64_t RiemannR(int64_t x)
{
#if defined(HAVE_FLOAT128)
if (x > 1 && std::log10(x) >= std::numeric_limits<long double>::digits10)
if (x > 1e10 && std::log10(x) >= std::numeric_limits<long double>::digits10)
return (int64_t) ::RiemannR((__float128) x);
#endif

return (int64_t) ::RiemannR((long double) x);
}

int64_t RiemannR_inverse(int64_t x)
{
#if defined(HAVE_FLOAT128)
double safetyThreshold = 4;
if (x > 1 && std::log10(x) >= std::numeric_limits<long double>::digits10 - safetyThreshold)
if (x > 1e10)
{
__float128 res = ::RiemannR_inverse((__float128) x);
if (res > (__float128) std::numeric_limits<int64_t>::max())
return std::numeric_limits<int64_t>::max();
return (int64_t) res;
double logx = std::log(x);
if (std::log10((double) x * (logx * logx)) >= std::numeric_limits<long double>::digits10)
{
__float128 res = ::RiemannR_inverse((__float128) x);
if (res > (__float128) std::numeric_limits<int64_t>::max())
return std::numeric_limits<int64_t>::max();
else
return (int64_t) res;
}
}
#endif

long double res = ::RiemannR_inverse((long double) x);

// Prevent integer overflow
if (res > (long double) std::numeric_limits<int64_t>::max())
return std::numeric_limits<int64_t>::max();

return (int64_t) res;
else
return (int64_t) res;
}

#ifdef HAVE_INT128_T

int128_t RiemannR(int128_t x)
{
#if defined(HAVE_FLOAT128)
if (x > 1 && std::log10(x) >= std::numeric_limits<long double>::digits10)
if (x > 1e10 && std::log10(x) >= std::numeric_limits<long double>::digits10)
return (int128_t) ::RiemannR((__float128) x);
#endif

return (int128_t) ::RiemannR((long double) x);
}

int128_t RiemannR_inverse(int128_t x)
{
#if defined(HAVE_FLOAT128)
double safetyThreshold = 4;
if (x > 1 && std::log10(x) >= std::numeric_limits<long double>::digits10 - safetyThreshold)
if (x > 1e10)
{
__float128 res = ::RiemannR_inverse((__float128) x);
if (res > (__float128) std::numeric_limits<int128_t>::max())
return std::numeric_limits<int128_t>::max();
return (int128_t) res;
double logx = std::log(x);
if (std::log10((double) x * (logx * logx)) >= std::numeric_limits<long double>::digits10)
{
__float128 res = ::RiemannR_inverse((__float128) x);
if (res > (__float128) std::numeric_limits<int128_t>::max())
return std::numeric_limits<int128_t>::max();
else
return (int128_t) res;
}
}
#endif

long double res = ::RiemannR_inverse((long double) x);

// Prevent integer overflow
if (res > (long double) std::numeric_limits<int128_t>::max())
return std::numeric_limits<int128_t>::max();

return (int128_t) res;
else
return (int128_t) res;
}

#endif
Expand Down

0 comments on commit 42dd380

Please sign in to comment.