This is an automated email from the git hooks/post-receive script.

ppm-guest pushed a commit to annotated tag v0.06
in repository libmath-prime-util-perl.

commit 92269f671fa20a3e7a0329e4c24da822a858d25b
Author: Dana Jacobsen <d...@acm.org>
Date:   Sun Jun 10 04:31:44 2012 -0600

    Add math functions
---
 Changes                    |   4 +-
 README                     |   2 +-
 XS.xs                      |   9 ++
 examples/bench-pcapprox.pl |  37 +++++++
 examples/test-pcapprox.pl  |  72 +++++++++++++
 lib/Math/Prime/Util.pm     |  51 ++++++++-
 t/18-functions.t           |  88 +++++++++++++++
 util.c                     | 259 ++++++++++++++++++++++++++++++++++++++++++++-
 util.h                     |   4 +
 9 files changed, 514 insertions(+), 12 deletions(-)

diff --git a/Changes b/Changes
index eaf0d16..1990528 100644
--- a/Changes
+++ b/Changes
@@ -3,8 +3,10 @@ Revision history for Perl extension Math::Prime::Util.
 0.05  8 June 2012
     - Use assembler for mulmod if 64-bit and gcc and x86
     - nth_prime routines should now all croak on overflow in the same way.
-    - Segmented prime_count, things like this are efficient:
+    - Segmented prime_count, things like this are reasonably efficient:
             say prime_count( 10**16,  10**16 + 2**20 )
+    - Add Ei(x), li(x), and R(x) functions.
+    - prime_count_approx uses R(x), making it vastly more accurate.
 
 0.04  7 June 2012
     - Didn't do tests on 32-bit machine before release.  Test suite caught
diff --git a/README b/README
index bffb345..36ed54e 100644
--- a/README
+++ b/README
@@ -1,4 +1,4 @@
-Math::Prime::Util version 0.04
+Math::Prime::Util version 0.05
 
 A set of utilities related to prime numbers.  These include multiple sieving
 methods, is_prime, prime_count, nth_prime, approximations and bounds for
diff --git a/XS.xs b/XS.xs
index a73f3c2..e12bbe5 100644
--- a/XS.xs
+++ b/XS.xs
@@ -388,3 +388,12 @@ miller_rabin(IN UV n, ...)
 
 int
 is_prob_prime(IN UV n)
+
+double
+ExponentialIntegral(double x)
+
+double
+LogarithmicIntegral(double x)
+
+double
+RiemannR(double x)
diff --git a/examples/bench-pcapprox.pl b/examples/bench-pcapprox.pl
new file mode 100755
index 0000000..576e275
--- /dev/null
+++ b/examples/bench-pcapprox.pl
@@ -0,0 +1,37 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+use Math::Prime::Util ":all";
+use Benchmark qw/:all/;
+use List::Util qw/min max/;
+my $maxdigits = (~0 <= 4294967295) ? 10 : 20;
+
+my $count = shift || -5;
+
+srand(29);
+test_at_digits($_) for (5 .. $maxdigits);
+
+
+sub test_at_digits {
+  my $digits = shift;
+  die "Digits must be > 0" unless $digits > 0;
+
+  my $base = ($digits == 1) ? 0 : int(10 ** ($digits-1));
+  my $max = int(10 ** $digits);
+  $max = ~0 if $max > ~0;
+  my @nums = map { $base+int(rand($max-$base)) } (1 .. 1000);
+  my $min_num = min @nums;
+  my $max_num = max @nums;
+
+  #print "miller_rabin for 1000 random $digits-digit numbers ($min_num - 
$max_num)\n";
+
+  my $sum;
+  cmpthese($count,{
+    'lower' => sub { $sum += prime_count_lower($_) for @nums },
+    'luapprox' => sub { $sum += 
(prime_count_lower($_)+prime_count_upper($_))/2 for @nums },
+    'approx' => sub { $sum += prime_count_approx($_) for @nums },
+    'li' => sub { $sum += LogarithmicIntegral($_) for @nums },
+    'R' => sub { $sum += RiemannR($_) for @nums },
+  });
+  print "\n";
+}
diff --git a/examples/test-pcapprox.pl b/examples/test-pcapprox.pl
new file mode 100755
index 0000000..70b1041
--- /dev/null
+++ b/examples/test-pcapprox.pl
@@ -0,0 +1,72 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+use Math::Prime::Util qw/prime_count prime_count_approx prime_count_lower 
prime_count_upper LogarithmicIntegral RiemannR/;
+$| = 1;  # fast pipes
+my $use64 = Math::Prime::Util::_maxbits > 32;
+
+
+my %pivals = (
+                  10 => 4,
+                 100 => 25,
+                1000 => 168,
+               10000 => 1229,
+              100000 => 9592,
+             1000000 => 78498,
+            10000000 => 664579,
+           100000000 => 5761455,
+          1000000000 => 50847534,
+         10000000000 => 455052511,
+        100000000000 => 4118054813,
+       1000000000000 => 37607912018,
+      10000000000000 => 346065536839,
+     100000000000000 => 3204941750802,
+    1000000000000000 => 29844570422669,
+   10000000000000000 => 279238341033925,
+  100000000000000000 => 2623557157654233,
+ 1000000000000000000 => 24739954287740860,
+10000000000000000000 => 234057667276344607,
+);
+
+printf("  N    %12s  %12s  %12s\n", "pc_approx", "Li", "R");
+printf("-----  %12s  %12s  %12s\n", '-'x12,'-'x12,'-'x12);
+foreach my $n (sort {$a<=>$b} keys %pivals) {
+  my $pin  = $pivals{$n};
+  my $pca  = prime_count_approx($n);
+
+  my $pcli = ($n < 2) ? 0 : 
int(LogarithmicIntegral($n)-LogarithmicIntegral(2)+0.5);
+
+  my $r = int(RiemannR($n)+0.5);
+
+  printf "10^%2d  %12d  %12d  %12d\n", length($n)-1,
+         abs($pca-$pin), abs($pcli-$pin), abs($r-$pin);
+}
+
+# Also see http://empslocal.ex.ac.uk/people/staff/mrwatkin/zeta/encoding1.htm
+# for some ideas one how this could be made even more accurate.
+
+print "\n";
+print "Lower / Upper bounds.  Percentages.\n";
+print "\n";
+
+printf("  N    %12s  %12s  %12s  %12s\n", "lower", "upper", "SchoenfeldL", 
"SchoenfeldU");
+printf("-----  %12s  %12s  %12s  %12s\n", '-'x12,'-'x12,'-'x12,'-'x12);
+foreach my $n (sort {$a<=>$b} keys %pivals) {
+  my $pin  = $pivals{$n};
+  my $pcl  = prime_count_lower($n);
+  my $pcu  = prime_count_upper($n);
+  my ($scl,$scu) = schoenfeld($n);
+
+  #printf "10^%2d  %12d  %12d\n", length($n)-1, $pin-$pcl, $pcu-$pin;
+  printf "10^%2d  %12.7f  %12.7f  %12.7f  %12.7f\n",
+         length($n)-1, 100*($pin-$pcl)/$pin, 100*($pcu-$pin)/$pin,
+         100*($pin-$scl)/$pin, 100*($scu-$pin)/$pin;
+}
+
+sub schoenfeld {
+  my $x = shift;
+  my $lix = LogarithmicIntegral($x);
+  my $bound = ((sqrt($x)*log($x)) / 8*3.1415926535);
+  ($lix-$bound,$lix+$bound);
+}
+
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 217f8a4..0ca736a 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -5,7 +5,7 @@ use Carp qw/croak confess/;
 
 BEGIN {
   $Math::Prime::Util::AUTHORITY = 'cpan:DANAJ';
-  $Math::Prime::Util::VERSION = '0.04';
+  $Math::Prime::Util::VERSION = '0.05';
 }
 
 # parent is cleaner, and in the Perl 5.10.1 / 5.12.0 core, but not earlier.
@@ -20,6 +20,7 @@ our @EXPORT_OK = qw(
                      nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
                      random_prime random_ndigit_prime
                      factor
+                     ExponentialIntegral LogarithmicIntegral RiemannR
                    );
 our %EXPORT_TAGS = (all => [ @EXPORT_OK ]);
 
@@ -387,10 +388,9 @@ above that range.
         " primes below one quintillion.\n";
 
 Returns an approximation to the C<prime_count> function, without having to
-generate any primes.  The results are very close for small numbers, but less
-so with large ranges.  The current implementation is 0.00033% too small
-for the example, but takes under a microsecond and no memory to get the
-result.
+generate any primes.  The current implementation uses the Riemann R function
+which is quite accurate.  A slightly faster answer, but less accurate,
+can be obtained by averaging the upper and lower bounds.
 
 
 =head2 nth_prime
@@ -618,6 +618,46 @@ input value.  On some inputs they will take a very long 
time, while on
 others they succeed in a remarkably short time.
 
 
+
+=head1 MATHEMATICAL FUNCTIONS
+
+=head2 ExponentialIntegral
+
+  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.
+
+=head2 LogarithmicIntegral
+
+  my $li = LogarithmicIntegral($x)
+
+Given a positive floating point input, returns the floating point logarithmic
+integral of C<x>, defined as the integral of C<dt/ln t> from C<0> to C<x>.
+If given a negative input, the function will croak.  The function returns
+0 at C<x = 0>, and C<-infinity> at C<x = 1>.
+
+This is often known as C<li(x)>.  A related function is the offset logarithmic
+integral, sometimes known as C<Li(x)> which avoids the singularity at 1.  It
+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.
+
+=head2 RiemannR
+
+  my $r = RiemannR($x);
+
+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.
+
+
+
 =head1 LIMITATIONS
 
 I have not completed testing all the functions near the word size limit
@@ -631,6 +671,7 @@ try to determine if your Perl is broken.  This will show up 
in factoring tests.
 Perl 5.6.2 32-bit works fine, as do later versions with 32-bit and 64-bit.
 
 
+
 =head1 PERFORMANCE
 
 Counting the primes to C<10^10> (10 billion), with time in seconds.
diff --git a/t/18-functions.t b/t/18-functions.t
new file mode 100644
index 0000000..1ed634a
--- /dev/null
+++ b/t/18-functions.t
@@ -0,0 +1,88 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+
+use Test::More;
+use Math::Prime::Util qw/prime_count ExponentialIntegral LogarithmicIntegral 
RiemannR/;
+
+my $use64 = Math::Prime::Util::_maxbits > 32;
+my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
+
+#plan tests => 14*3 + ($extra ? 14 : 8) + ($use64 ? 2*18 : 0);
+
+eval { ExponentialIntegral(0); };
+like($@, qr/invalid/i, "Ei(0) is invalid");
+eval { LogarithmicIntegral(-1); };
+like($@, qr/invalid/i, "li(-1) is invalid");
+eval { RiemannR(0); };
+like($@, qr/invalid/i, "R(0) is invalid");
+eval { RiemannR(-1); };
+like($@, qr/invalid/i, "R(-1) is invalid");
+
+cmp_ok( LogarithmicIntegral(1), '<=', 0 - (~0 * ~0), "li(1) is -infinity" );
+
+
+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,
+);
+
+while (my($n, $ein) = each (%eivals)) {
+  cmp_closeto( ExponentialIntegral($n), $ein, 0.00000001 * abs($ein), "Ei($n) 
~= $ein");
+}
+
+# In pari these are:  -eint1(-log($n))
+my %livals = (
+              0 =>  0,
+           1.01 => -4.0229586739299358695031,
+              2 =>  1.0451637801174927848446,
+             10 =>  6.1655995047872979375230,
+             24 =>  11.200315795232698830550,
+           1000 =>  177.60965799015222668764,
+         100000 =>  9629.8090010507982050343,
+      100000000 =>  5762209.3754480314675691,
+     4294967295 =>  203284081.95454158906409,
+    10000000000 =>  455055614.58662307560953,
+   100000000000 =>  4118066400.6216115150394,
+);
+
+while (my($n, $lin) = each (%livals)) {
+  cmp_closeto( LogarithmicIntegral($n), $lin, 0.00000001 * abs($lin), "li($n) 
~= $lin");
+}
+
+# Values from T. R. Nicely for comparison
+my %rvals = (
+           1.01 =>  1.0060697180622924796117,
+              2 =>  1.5410090161871318832885,
+             10 =>  4.5645831410050902398658,
+           1000 =>  168.35944628116734806491,
+        1000000 =>  78527.399429127704858870,
+       10000000 =>  664667.44756474776798535,
+     4294967295 =>  203280697.51326064541983,
+    10000000000 =>  455050683.30684692446315,
+18446744073709551615 => 4.25656284014012122706963685602e17,
+);
+while (my($n, $rin) = each (%rvals)) {
+  cmp_closeto( RiemannR($n), $rin, 0.00000001 * abs($rin), "R($n) ~= $rin");
+}
+
+
+done_testing();
+
+sub cmp_closeto {
+  my $got = shift;
+  my $expect = shift;
+  my $tolerance = shift;
+  my $message = shift;
+  cmp_ok( abs($got - $expect), '<=', $tolerance, $message );
+}
diff --git a/util.c b/util.c
index c924322..c469f05 100644
--- a/util.c
+++ b/util.c
@@ -377,13 +377,30 @@ UV prime_count_upper(UV x)
 
 UV prime_count_approx(UV x)
 {
-  /* Placeholder for fancy algorithms, like Tomás Oliveira e Silva's:
+  /*
+   * A simple way:
+   *     return ((prime_count_lower(x) + prime_count_upper(x)) / 2);
+   * With the current bounds, this is ~131k at 10^10 and 436B at 10^19.
+   *
+   * The logarithmic integral works quite well, with absolute errors of
+   * ~3100 at 10^10 and ~100M at 10^19.
+   *
+   * Riemann's R function works astoundingly well, with errors of ~1828
+   * at 10^10 and 24M at 10^19.
+   *
+   * Getting fancier, one try using Riemann's pi formula:
    *     http://trac.sagemath.org/sage_trac/ticket/8135
-   * or an approximation to Li(x) plus a delta.
    */
-  UV lower = prime_count_lower(x);
-  UV upper = prime_count_upper(x);
-  return ((lower + upper) / 2);
+  double R;
+  if (x < NPRIME_COUNT_SMALL)
+    return prime_count_small[x];
+
+  R = RiemannR(x);
+  /* We could add the additional factor:
+   *   R = R - (1.0 / log(x)) + (M_1_PI * atan(M_PI/log(x)))
+   * but it's extraordinarily small, so not worth calculating here.
+   */
+  return (UV)(R+0.5);
 }
 
 
@@ -762,3 +779,235 @@ UV prime_count_seg(UV low, UV high)
 
   return count;
 }
+
+/*
+ * See:
+ *  "Multiple-Precision Exponential Integral and Related Functions"
+ *      by David M. Smith, and
+ *  "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.
+ *
+ * Any mistakes here are completely my fault.  This code has not been
+ * verified for anything serious.  For better reulsts, see:
+ *    http://www.trnicely.net/pi/pix_0000.htm
+ * which although the author claims are demonstration programs, will
+ * produce more usable results than this code does.
+ */
+
+static double const euler_mascheroni = 
0.57721566490153286060651209008240243104215933593992;
+static double const li2 = 1.045163780117492784844588889194613136522615578151;
+
+double ExponentialIntegral(double x) {
+  double const tol = 0.0000000001;
+  double val, term, x_to_the_n, fact_n;
+  double y, t;
+  double sum = 0.0;
+  double c = 0.0;
+  int n;
+
+  if (x == 0) croak("Invalid input to ExponentialIntegral:  x must be != 0");
+
+  if (x < 0) {
+    /* continued fraction */
+#if 0
+    val = 0;
+    for (n = 50; n >= 1; n--) {
+      val = -n * n / (2 * n + 1 - x + val);
+    }
+    val = -exp(x) / (1 - x + val);
+#else
+    double old;
+    double lc = 0;
+    double ld = 1 / (1 - x);
+    val = ld * (-exp(x));
+    for (n = 1; n <= 2000; n++) {
+      lc = 1 / (2*n + 1 - x - n*n*lc);
+      ld = 1 / (2*n + 1 - x - n*n*ld);
+      old = val;
+      val *= ld/lc;
+      if ( fabs(val-old) <= tol*fabs(val) )
+        break;
+    }
+#endif
+  } else if (x < -log(tol)) {
+    /* Convergent series */
+    fact_n = 1;
+
+    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++) {
+      fact_n *= x/n;
+      term = fact_n/n;
+      y = term-c;  t = sum+y;  c = (t-sum)-y;  sum = t;
+      /* printf("C  after adding %.8lf, val = %.8lf\n", term, sum); */
+      if (term < tol) break;
+    }
+    val = sum;
+    if (term > tol) printf("Tolerance not met after %d rounds for Ei(%lf)\n", 
n, x);
+  } else {
+    /* Asymptotic divergent series */
+#if 0
+    double last_term, sum;
+    last_term = 1.0;
+    x_to_the_n = 1;
+    fact_n = 1;
+
+    y = 1.0-c;  t = sum+y;  c = (t-sum)-y;  sum = t;
+
+    for (n = 1; n <= 100; n++) {
+      x_to_the_n *= x;
+      fact_n *= n;
+      term = fact_n / x_to_the_n;
+      if (term > last_term) break;
+      last_term = term;
+      y = term-c;  t = sum+y;  c = (t-sum)-y;  sum = t;
+      /* printf("A  after adding %.8lf, sum = %.8lf\n", term, sum); */
+    }
+    val = (exp(x) / x) * sum - euler_mascheroni;
+#else
+    double last_term, sum;
+
+    val = exp(x) / x;
+    term = 1.0;
+    y = 1.0-c;  t = sum+y;  c = (t-sum)-y;  sum = t;
+
+    for (n = 1; n <= 100; n++) {
+      last_term = term;
+      term *= n/x;
+      if (term < (tol/val)) 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); */
+      } else {
+        y = (-last_term/3)-c;  t = sum+y;  c = (t-sum)-y;  sum = t;
+        /* printf("A  after adding %.8lf, sum = %.8lf\n", -last_term/3, sum); 
*/
+        break;
+      }
+    }
+    val *= sum;
+#endif
+  }
+
+  return val;
+}
+
+double LogarithmicIntegral(double x) {
+  if (x == 0) return 0;
+  if (x == 1) return -INFINITY;
+  if (x == 2) return li2;
+  if (x <= 0) croak("Invalid input to ExponentialIntegral:  x must be > 0");
+  return ExponentialIntegral(log(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. */
+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,
+};
+#define NPRECALC_ZETA 
(sizeof(riemann_zeta_table)/sizeof(riemann_zeta_table[0]))
+
+static double evaluate_zeta(double x) {
+  double const tol = 1e-14;
+  double y, t;
+  double sum = 0.0;
+  double c = 0.0;
+  double term;
+  int k;
+
+  /* Simple method.  Slow and inaccurate near x=1, but iterations and accuracy
+   * improve quickly.
+   */
+
+  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++) {
+    if (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);
+    else              term = 1.0 / pow(k, x);
+    y = term-c;  t = sum+y;  c = (t-sum)-y;  sum = t;
+  }
+  return sum;
+}
+
+double RiemannR(double x) {
+  double const tol = 1e-14;
+  double y, t, part_term, term, flogx, zeta;
+  double sum = 0.0;
+  double c = 0.0;
+  int k;
+
+  if (x <= 0) croak("Invalid input to ReimannR:  x must be > 0");
+
+  y = 1.0-c;  t = sum+y;  c = (t-sum)-y;  sum = t;
+
+  flogx = log(x);
+  part_term = 1;
+
+  /* Do small k with zetas from table */
+  for (k = 1; k <= NPRECALC_ZETA; k++) {
+    zeta = riemann_zeta_table[k+1-2];
+    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;
+    /* printf("R  after adding %.8lf, sum = %.8lf\n", term, sum); */
+  }
+  /* Finish with function */
+  for (; k <= 10000; k++) {
+    if (term < tol) break;
+    zeta = evaluate_zeta(k+1);
+    part_term *= flogx / k;
+    term = part_term / (k * zeta);
+    y = term-c;  t = sum+y;  c = (t-sum)-y;  sum = t;
+    /* printf("R  after adding %.8lf, sum = %.8lf\n", term, sum); */
+  }
+
+
+  return sum;
+}
diff --git a/util.h b/util.h
index 6a09d98..cb7ff43 100644
--- a/util.h
+++ b/util.h
@@ -20,6 +20,10 @@ extern UV  nth_prime_upper(UV x);
 extern UV  nth_prime_approx(UV x);
 extern UV  nth_prime(UV x);
 
+extern double ExponentialIntegral(double x);
+extern double LogarithmicIntegral(double x);
+extern double RiemannR(double x);
+
 #define SEGMENT_CHUNK_SIZE  UVCONST(262144)
 extern unsigned char* get_prime_segment(void);
 extern void           free_prime_segment(void);

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