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 4f1ca63c13984d58f075eac667bd2a29d0dadd89
Author: Dana Jacobsen <d...@acm.org>
Date:   Tue Jun 19 16:15:52 2012 -0600

    Add rational Chebyshev approx for Ei
---
 lib/Math/Prime/Util.pm |  24 +++++-----
 sieve.c                |   6 ++-
 t/18-functions.t       |  32 +++++++------
 util.c                 | 123 +++++++++++++++++++++++++++++--------------------
 4 files changed, 110 insertions(+), 75 deletions(-)

diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 62bc48e..c83c531 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -1,7 +1,7 @@
 package Math::Prime::Util;
 use strict;
 use warnings;
-use Carp qw/croak confess/;
+use Carp qw/croak confess carp/;
 
 BEGIN {
   $Math::Prime::Util::AUTHORITY = 'cpan:DANAJ';
@@ -694,15 +694,17 @@ others they succeed in a remarkably short time.
 
   my $Ei = ExponentialIntegral($x);
 
-Given a non-zero floating point input C<x>, this returns the exponential 
integral
-of C<x>, defined as the integral of C<e^t/t dt> from C<-infinity> to C<x>.
-Depending on the input, the integral
-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.
+Given a non-zero floating point input C<x>, this returns the real-valued
+exponential integral of C<x>, defined as the integral of C<e^t/t dt>
+from C<-infinity> to C<x>.
+Depending on the input, the integral is calculated using
+continued fractions (C<x E<lt> -1>),
+rational Chebyshev approximation (C< -1 E<lt> x E<lt> 0>),
+a convergent series (small positive C<x>),
+or an asymptotic divergent series (large positive C<x>).
+
+Accuracy should be at least 14 digits.
 
-Accuracy is typically 15 digits unless between C<-.1> and C<0>, where
-accuracy is degraded.
 
 =head2 LogarithmicIntegral
 
@@ -720,9 +722,7 @@ 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.
+Accuracy should be at least 14 digits.
 
 
 =head2 RiemannR
diff --git a/sieve.c b/sieve.c
index 1ddc013..61931b4 100644
--- a/sieve.c
+++ b/sieve.c
@@ -97,7 +97,11 @@ static void sieve_prefill(unsigned char* mem, UV startd, UV 
endd)
   UV nbytes = endd - startd + 1;
   MPUassert( (mem != 0) && (endd >= startd), "sieve_prefill bad arguments");
 
-  /* TODO: Big ranges would benefit from copy doubling */
+  /* Walk the memory, tiling in the presieve area using memcpy.
+   * This is pretty fast, but it might still benefit from using copy
+   * doubling (where we copy to the memory, then copy memory to memory
+   * doubling in size each time), as memcpy usually loves big chunks.
+   */
   while (startd <= endd) {
     UV pstartd = startd % PRESIEVE_SIZE;
     UV sieve_bytes = PRESIEVE_SIZE - pstartd;
diff --git a/t/18-functions.t b/t/18-functions.t
index 4e5bda7..04c447d 100644
--- a/t/18-functions.t
+++ b/t/18-functions.t
@@ -8,7 +8,7 @@ use Math::Prime::Util qw/prime_count ExponentialIntegral 
LogarithmicIntegral Rie
 my $use64 = Math::Prime::Util::_maxbits > 32;
 my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
 
-plan tests => 4 + 1 + 12 + 11 + 9;
+plan tests => 4 + 1 + 1 + 16 + 11 + 9;
 
 eval { ExponentialIntegral(0); };
 like($@, qr/invalid/i, "Ei(0) is invalid");
@@ -21,20 +21,26 @@ like($@, qr/invalid/i, "R(-1) is invalid");
 
 cmp_ok( LogarithmicIntegral(1), '<=', 0 - (~0 * ~0), "li(1) is -infinity" );
 
+# Example used in Math::Cephes
+cmp_closeto( ExponentialIntegral(2.2), 5.732614700, 1e-06, "Ei(2.2)");
 
 my %eivals = (
-    -10   => -.00000415696892968532438,
-     -0.1 => -1.8229239584193906660809,
-      0.693147180559945 => 1.0451637801174927848446,           # log2
-      1   =>  1.8951178163559367554665,
-      1.5 =>  3.3012854491297978379574,
-      2   =>  4.9542343560018901633795,
-      5   =>  40.185275355803177455091,
-      10  =>  2492.2289762418777591384,
-      12  =>  14959.532666397528852292,
-      20  =>  25615652.664056588820481,
-      40  =>  6039718263611241.5783592,
-      41  =>  16006649143245041.110700,
+         -10 =>  -0.00000415696892968532438,
+        -0.5 =>  -0.55977359477616,
+        -0.1 =>  -1.8229239584193906660809,
+      -0.001 =>  -6.33153936413615,
+    -0.00001 => -10.9357198000436956,
+ -0.00000001 => -17.843465089050832587,
+ 0.693147180559945 => 1.0451637801174927848446,           # log2
+         1   =>  1.8951178163559367554665,
+         1.5 =>  3.3012854491297978379574,
+         2   =>  4.9542343560018901633795,
+         5   =>  40.185275355803177455091,
+         10  =>  2492.2289762418777591384,
+         12  =>  14959.532666397528852292,
+         20  =>  25615652.664056588820481,
+         40  =>  6039718263611241.5783592,
+         41  =>  16006649143245041.110700,
 );
 
 while (my($n, $ein) = each (%eivals)) {
diff --git a/util.c b/util.c
index 8e49e1b..6b7e728 100644
--- a/util.c
+++ b/util.c
@@ -785,11 +785,13 @@ UV nth_prime(UV n)
 /*
  * See:
  *  "Multiple-Precision Exponential Integral and Related Functions"
- *      by David M. Smith, and
+ *      by David M. Smith
  *  "On the Evaluation of the Complex-Valued Exponential Integral"
  *      by Vincent Pegoraro and Philipp Slusallek
  *  "Numerical Recipes" 3rd edition
  *      by William H. Press et al.
+ *  "Rational Chevyshev Approximations for the Exponential Integral E_1(x)"
+ *      by W. J. Cody and Henry C. Thacher, Jr.
  *
  * Any mistakes here are completely my fault.  This code has not been
  * verified for anything serious.  For better results, see:
@@ -798,8 +800,6 @@ UV nth_prime(UV n)
  * undoubtedly produce more reliable results than this code does (I don't
  * know of any obvious issues with this code, but it just hasn't been used
  * by many people).
- *
- * TODO: Verify error bounds at different ranges.
  */
 
 static double const euler_mascheroni = 
0.57721566490153286060651209008240243104215933593992;
@@ -815,9 +815,8 @@ double ExponentialIntegral(double x) {
 
   if (x == 0) croak("Invalid input to ExponentialIntegral:  x must be != 0");
 
-  if (x < 0) {
-    /* continued fraction */
-    /* This loses accuracy very quickly below -0.001 or so */
+  if (x < -1) {
+    /* Continued fraction, good for x < -1 */
     double old;
     double lc = 0;
     double ld = 1 / (1 - x);
@@ -830,6 +829,32 @@ double ExponentialIntegral(double x) {
       if ( fabs(val-old) <= tol*fabs(val) )
         break;
     }
+  } else if (x < 0) {
+    /* Rational Chebyshev approximation (Cody, Thacher), good for -1 < x < 0 */
+#if 0
+    static const double C2p[3] = { -4.43668255, 4.42054938, 3.16274620 };
+    static const double C2q[3] = {  7.68641124, 5.65655216, 1.00000000 };
+    double sumn = C2p[0] - x*(C2p[1] - x*C2p[2]);
+    double sumd = C2q[0] - x*(C2q[1] - x*C2q[2]);
+#else
+    static const double C6p[7] = { -148151.02102575750838086,
+                                    150260.59476436982420737,
+                                     89904.972007457256553251,
+                                     15924.175980637303639884,
+                                      2150.0672908092918123209,
+                                       116.69552669734461083368,
+                                         5.0196785185439843791020 };
+    static const double C6q[7] = {  256664.93484897117319268,
+                                    184340.70063353677359298,
+                                     52440.529172056355429883,
+                                      8125.8035174768735759866,
+                                       750.43163907103936624165,
+                                        40.205465640027706061433,
+                                         1.0000000000000000000000 };
+    double sumn = 
C6p[0]-x*(C6p[1]-x*(C6p[2]-x*(C6p[3]-x*(C6p[4]-x*(C6p[5]-x*C6p[6])))));
+    double sumd = 
C6q[0]-x*(C6q[1]-x*(C6q[2]-x*(C6q[3]-x*(C6q[4]-x*(C6q[5]-x*C6q[6])))));
+#endif
+    val = log(-x) - sumn/sumd;
   } else if (x < -log(tol)) {
     /* Convergent series */
     fact_n = 1;
@@ -883,47 +908,47 @@ double LogarithmicIntegral(double x) {
 /*
  * Storing the first 10-20 Zeta values makes sense.  Past that it is purely
  * to avoid making the call to generate them ourselves.  We could cache the
- * calculated values. */
+ * calculated values. These all have 1 subtracted from them.  */
 static const double riemann_zeta_table[] = {
-  1.6449340668482264364724151666460251892,  /* zeta(2) */
-  1.2020569031595942853997381615114499908,
-  1.0823232337111381915160036965411679028,
-  1.0369277551433699263313654864570341681,
-  1.0173430619844491397145179297909205279,
-  1.0083492773819228268397975498497967596,
-  1.0040773561979443393786852385086524653,
-  1.0020083928260822144178527692324120605,
-  1.0009945751278180853371459589003190170,
-  1.0004941886041194645587022825264699365,
-  1.0002460865533080482986379980477396710,
-  1.0001227133475784891467518365263573957,
-  1.0000612481350587048292585451051353337,
-  1.0000305882363070204935517285106450626,
-  1.0000152822594086518717325714876367220,
-  1.0000076371976378997622736002935630292,  /* zeta(17)  Past here all we're */
-  1.0000038172932649998398564616446219397,  /* zeta(18)  getting is speed.   */
-  1.0000019082127165539389256569577951013,
-  1.0000009539620338727961131520386834493,
-  1.0000004769329867878064631167196043730,
-  1.0000002384505027277329900036481867530,
-  1.0000001192199259653110730677887188823,
-  1.0000000596081890512594796124402079358,
-  1.0000000298035035146522801860637050694,
-  1.0000000149015548283650412346585066307,
-  1.0000000074507117898354294919810041706,
-  1.0000000037253340247884570548192040184,
-  1.0000000018626597235130490064039099454,
-  1.0000000009313274324196681828717647350,
-  1.0000000004656629065033784072989233251,
-  1.0000000002328311833676505492001455976,
-  1.0000000001164155017270051977592973835,
-  1.0000000000582077208790270088924368599,
-  1.0000000000291038504449709968692942523,
-  1.0000000000145519218910419842359296322,
-  1.0000000000072759598350574810145208690,
-  1.0000000000036379795473786511902372363,
-  1.0000000000018189896503070659475848321,
-  1.0000000000009094947840263889282533118,
+  0.6449340668482264364724151666460251892,  /* zeta(2) */
+  0.2020569031595942853997381615114499908,
+  0.0823232337111381915160036965411679028,
+  0.0369277551433699263313654864570341681,
+  0.0173430619844491397145179297909205279,
+  0.0083492773819228268397975498497967596,
+  0.0040773561979443393786852385086524653,
+  0.0020083928260822144178527692324120605,
+  0.0009945751278180853371459589003190170,
+  0.0004941886041194645587022825264699365,
+  0.0002460865533080482986379980477396710,
+  0.0001227133475784891467518365263573957,
+  0.0000612481350587048292585451051353337,
+  0.0000305882363070204935517285106450626,
+  0.0000152822594086518717325714876367220,
+  0.0000076371976378997622736002935630292,  /* zeta(17)  Past here all we're */
+  0.0000038172932649998398564616446219397,  /* zeta(18)  getting is speed.   */
+  0.0000019082127165539389256569577951013,
+  0.0000009539620338727961131520386834493,
+  0.0000004769329867878064631167196043730,
+  0.0000002384505027277329900036481867530,
+  0.0000001192199259653110730677887188823,
+  0.0000000596081890512594796124402079358,
+  0.0000000298035035146522801860637050694,
+  0.0000000149015548283650412346585066307,
+  0.0000000074507117898354294919810041706,
+  0.0000000037253340247884570548192040184,
+  0.0000000018626597235130490064039099454,
+  0.0000000009313274324196681828717647350,
+  0.0000000004656629065033784072989233251,
+  0.0000000002328311833676505492001455976,
+  0.0000000001164155017270051977592973835,
+  0.0000000000582077208790270088924368599,
+  0.0000000000291038504449709968692942523,
+  0.0000000000145519218910419842359296322,
+  0.0000000000072759598350574810145208690,
+  0.0000000000036379795473786511902372363,
+  0.0000000000018189896503070659475848321,
+  0.0000000000009094947840263889282533118,
 };
 #define NPRECALC_ZETA 
(sizeof(riemann_zeta_table)/sizeof(riemann_zeta_table[0]))
 
@@ -939,7 +964,7 @@ static double evaluate_zeta(double x) {
    * improve quickly.
    */
 
-  term = 1.0;          y = term-c;  t = sum+y;  c = (t-sum)-y;  sum = t;
+  /* term = 1.0;          y = term-c;  t = sum+y;  c = (t-sum)-y;  sum = t; */
   term = 1.0/exp2(x);  y = term-c;  t = sum+y;  c = (t-sum)-y;  sum = t;
 
   for (k = 3; k <= 100000; k++) {
@@ -971,7 +996,7 @@ double RiemannR(double x) {
   for (k = 1; k <= (int)NPRECALC_ZETA; k++) {
     zeta = riemann_zeta_table[k+1-2];
     part_term *= flogx / k;
-    term = part_term / (k * zeta);
+    term = part_term / (k + k * zeta);
     y = term-c;  t = sum+y;  c = (t-sum)-y;  sum = t;
     if (fabs(term) < tol) break;
     /* printf("R  after adding %.8lf, sum = %.8lf\n", term, sum); */
@@ -981,7 +1006,7 @@ double RiemannR(double x) {
     if (fabs(term) < tol) break;
     zeta = evaluate_zeta(k+1);
     part_term *= flogx / k;
-    term = part_term / (k * zeta);
+    term = part_term / (k + k * zeta);
     y = term-c;  t = sum+y;  c = (t-sum)-y;  sum = t;
     /* printf("R  after adding %.8lf, sum = %.8lf\n", term, sum); */
   }

-- 
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

Reply via email to