This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.08
in repository libmath-prime-util-perl.
commit fd61c4c7f6f06c4d37ac263952e2f32da5c5364a
Author: Dana Jacobsen <d...@acm.org>
Date: Mon Jun 18 16:54:19 2012 -0600
Accuracy for math functions
---
Changes | 1 +
TODO | 2 --
lib/Math/Prime/Util.pm | 10 ++++++++++
util.c | 21 +++++++++++----------
4 files changed, 22 insertions(+), 12 deletions(-)
diff --git a/Changes b/Changes
index aa4f9da..2eeaf46 100644
--- a/Changes
+++ b/Changes
@@ -2,6 +2,7 @@ Revision history for Perl extension Math::Prime::Util.
0.08 20 June 2012
- Added thread scaffolding.
+ - Accuracy improvement and measurements for math functions.
0.07 17 June 2012
- Fixed a bug in next_prime found by Lou Godio (thank you VERY much!).
diff --git a/TODO b/TODO
index 303f8c6..93c3a6a 100644
--- a/TODO
+++ b/TODO
@@ -37,5 +37,3 @@
- Iterator or tie?
- Get rid of erat_simple and bitarray.h. They're only there for comparison.
-
-- Need to be more aware of precision for the math operations.
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index e225335..fa0c2c6 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -672,6 +672,9 @@ is calculated using continued fractions (negative C<x>), a
convergent series
(small positive C<x>), or an asymptotic divergent series (large positive C<x>).
This function only considers the real part.
+Accuracy is typically 15 digits unless between C<-.1> and C<0>, where
+accuracy is degraded.
+
=head2 LogarithmicIntegral
my $li = LogarithmicIntegral($x)
@@ -688,6 +691,11 @@ may be defined as C<Li(x) = li(x) - li(2)>.
This function is implemented as C<li(x) = Ei(ln x)> after handling special
values.
+Accuracy is typically 15 digits unless between C<0.9> and C<1>, where
+accuracy is degraded. These are not typically values one would use for
+prime number functionality, so should have little impact.
+
+
=head2 RiemannR
my $r = RiemannR($x);
@@ -696,6 +704,8 @@ Given a positive non-zero floating point input, returns the
floating
point value of Riemann's R function. Riemann's R function gives a very close
approximation to the prime counting function.
+Accuracy should be at least 14 digits.
+
=head1 LIMITATIONS
diff --git a/util.c b/util.c
index 30c4c8d..8e49e1b 100644
--- a/util.c
+++ b/util.c
@@ -806,7 +806,7 @@ static double const euler_mascheroni =
0.577215664901532860606512090082402431042
static double const li2 = 1.045163780117492784844588889194613136522615578151;
double ExponentialIntegral(double x) {
- double const tol = 0.0000000001;
+ double const tol = 1e-16;
double val, term, fact_n;
double y, t;
double sum = 0.0;
@@ -817,11 +817,12 @@ double ExponentialIntegral(double x) {
if (x < 0) {
/* continued fraction */
+ /* This loses accuracy very quickly below -0.001 or so */
double old;
double lc = 0;
double ld = 1 / (1 - x);
val = ld * (-exp(x));
- for (n = 1; n <= 2000; n++) {
+ for (n = 1; n <= 100000; n++) {
lc = 1 / (2*n + 1 - x - n*n*lc);
ld = 1 / (2*n + 1 - x - n*n*ld);
old = val;
@@ -836,7 +837,7 @@ double ExponentialIntegral(double x) {
y = euler_mascheroni-c; t = sum+y; c = (t-sum)-y; sum = t;
y = log(x)-c; t = sum+y; c = (t-sum)-y; sum = t;
- for (n = 1; n <= 100; n++) {
+ for (n = 1; n <= 200; n++) {
fact_n *= x/n;
term = fact_n/n;
y = term-c; t = sum+y; c = (t-sum)-y; sum = t;
@@ -852,10 +853,10 @@ double ExponentialIntegral(double x) {
term = 1.0;
y = 1.0-c; t = sum+y; c = (t-sum)-y; sum = t;
- for (n = 1; n <= 100; n++) {
+ for (n = 1; n <= 200; n++) {
last_term = term;
term *= n/x;
- if (term < (tol/val)) break;
+ if (term < tol) break;
if (term < last_term) {
y = term-c; t = sum+y; c = (t-sum)-y; sum = t;
/* printf("A after adding %.8lf, sum = %.8lf\n", term, sum); */
@@ -927,7 +928,7 @@ static const double riemann_zeta_table[] = {
#define NPRECALC_ZETA
(sizeof(riemann_zeta_table)/sizeof(riemann_zeta_table[0]))
static double evaluate_zeta(double x) {
- double const tol = 1e-14;
+ double const tol = 1e-16;
double y, t;
double sum = 0.0;
double c = 0.0;
@@ -942,7 +943,7 @@ static double evaluate_zeta(double x) {
term = 1.0/exp2(x); y = term-c; t = sum+y; c = (t-sum)-y; sum = t;
for (k = 3; k <= 100000; k++) {
- if (term < tol) break;
+ if (fabs(term) < tol) break;
if (k == 4) term = 1.0 / exp2(2*x);
else if (k == 8) term = 1.0 / exp2(4*x);
else if (k == 16) term = 1.0 / exp2(8*x);
@@ -953,7 +954,7 @@ static double evaluate_zeta(double x) {
}
double RiemannR(double x) {
- double const tol = 1e-14;
+ double const tol = 1e-16;
double y, t, part_term, term, flogx, zeta;
double sum = 0.0;
double c = 0.0;
@@ -972,12 +973,12 @@ double RiemannR(double x) {
part_term *= flogx / k;
term = part_term / (k * zeta);
y = term-c; t = sum+y; c = (t-sum)-y; sum = t;
- if (term < tol) break;
+ if (fabs(term) < tol) break;
/* printf("R after adding %.8lf, sum = %.8lf\n", term, sum); */
}
/* Finish with function */
for (; k <= 10000; k++) {
- if (term < tol) break;
+ if (fabs(term) < tol) break;
zeta = evaluate_zeta(k+1);
part_term *= flogx / k;
term = part_term / (k * zeta);
--
Alioth's /usr/local/bin/git-commit-notice on
/srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git
_______________________________________________
Pkg-perl-cvs-commits mailing list
Pkg-perl-cvs-commits@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits
