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

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

commit d94fc1102872ebb07e1e321aafc8c296ad2d5e25
Author: Dana Jacobsen <d...@acm.org>
Date:   Fri Jun 29 21:33:30 2012 -0600

    Lots of bignum support
---
 Makefile.PL               |   3 +
 TODO                      |  10 ++-
 XS.xs                     |   8 +-
 lib/Math/Prime/Util.pm    | 225 +++++++++++++++++++++++++++++++++++++---------
 lib/Math/Prime/Util/PP.pm |  93 +++++++++++++------
 t/15-probprime.t          |   2 +-
 t/19-moebius.t            |  55 ++++++++++++
 util.c                    |  22 ++---
 util.h                    |   6 +-
 9 files changed, 335 insertions(+), 89 deletions(-)

diff --git a/Makefile.PL b/Makefile.PL
index a9afe32..b2f5988 100644
--- a/Makefile.PL
+++ b/Makefile.PL
@@ -23,6 +23,9 @@ WriteMakefile1(
                       'Carp'             => 0,
                       'Tie::Array'       => 0,
                       'base'             => 0,
+                      'bignum'           => 0,
+                      'Math::BigInt'     => 0,
+                      'Math::BigFloat'   => 0,
                     },
 
     MIN_PERL_VERSION => 5.006002,
diff --git a/TODO b/TODO
index 828b213..f04f361 100644
--- a/TODO
+++ b/TODO
@@ -14,7 +14,7 @@
 - Is the current PP.pm setup the way we want to do it?
 
 - input validation (in XS, or do we need to make Perl wrappers for everything?)
-  We can do inpuut val in XS by looking at the NV.  But I think long term we'll
+  We can do input val in XS by looking at the NV.  But I think long term we'll
   have a little Perl front end for everything to route to bignum routines or
   regular routines.
 
@@ -26,3 +26,11 @@
   version does it in a painful way.  Something simpler to be had?
 
 - In PP nth_prime* and prime_count*, turn on bignum if bigint is on.
+
+- Add euler_phi documentation
+
+- More testing for bignum, including a test suite file
+
+- Move the approx / bounds functions into main, to support bignum properly.
+
+- need next_prime and prev_prime bignum support.  random_ndigit_prime needs it.
diff --git a/XS.xs b/XS.xs
index d7f2f32..842e2b3 100644
--- a/XS.xs
+++ b/XS.xs
@@ -212,7 +212,7 @@ erat_primes(IN UV low, IN UV high)
 
 
 void
-factor(IN UV n)
+XS_factor(IN UV n)
   PPCODE:
     if (n < 4) {
       XPUSHs(sv_2mortal(newSVuv( n ))); /* If n is 0-3, we're done. */
@@ -371,10 +371,10 @@ int
 is_prob_prime(IN UV n)
 
 double
-ExponentialIntegral(double x)
+XS_ExponentialIntegral(double x)
 
 double
-LogarithmicIntegral(double x)
+XS_LogarithmicIntegral(double x)
 
 double
-RiemannR(double x)
+XS_RiemannR(double x)
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index d81d3b5..f7511d6 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -12,6 +12,7 @@ BEGIN {
 # use parent qw( Exporter );
 use base qw( Exporter );
 our @EXPORT_OK = qw(
+                     prime_get_config
                      prime_precalc prime_memfree
                      is_prime is_prob_prime miller_rabin
                      primes
@@ -19,23 +20,29 @@ our @EXPORT_OK = qw(
                      prime_count prime_count_lower prime_count_upper 
prime_count_approx
                      nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
                      random_prime random_ndigit_prime
-                     factor all_factors
+                     factor all_factors moebius euler_phi
                      ExponentialIntegral LogarithmicIntegral RiemannR
                    );
 our %EXPORT_TAGS = (all => [ @EXPORT_OK ]);
 
-my $_pure_perl;
+my %_Config;
 
 BEGIN {
+
+  # Load PP code.  Nothing exported.
+  require Math::Prime::Util::PP;
+
   eval {
     require XSLoader;
     XSLoader::load(__PACKAGE__, $Math::Prime::Util::VERSION);
     prime_precalc(0);
-    $_pure_perl = 0;
+    $_Config{'xs'} = 1;
+    $_Config{'maxbits'} = _maxbits();
     1;
   } or do {
-    require Math::Prime::Util::PP;
-    $_pure_perl = 1;
+die $@;
+    $_Config{'xs'} = 0;
+    $_Config{'maxbits'} = Math::Prime::Util::PP::_maxbits();
     #carp "Using Pure Perl implementation";
     *_get_prime_cache_size = \&Math::Prime::Util::PP::_get_prime_cache_size;
     *_maxbits = \&Math::Prime::Util::PP::_maxbits;
@@ -72,18 +79,39 @@ BEGIN {
     *RiemannR            = \&Math::Prime::Util::PP::RiemannR;
     *LogarithmicIntegral = \&Math::Prime::Util::PP::LogarithmicIntegral;
     *ExponentialIntegral = \&Math::Prime::Util::PP::ExponentialIntegral;
-    # TODO: We should have some tests to verify XS vs. PP.
   }
 }
 END {
   _prime_memfreeall;
 }
 
-my $_maxparam  = (_maxbits == 32) ? 4294967295 : 18446744073709551615;
-my $_maxdigits = (_maxbits == 32) ? 10 : 20;
-my $_maxprime  = (_maxbits == 32) ? 4294967291 : 18446744073709551557;
-my $_maxprimeidx=(_maxbits == 32) ?  203280221 :   425656284035217743;
+if ($_Config{'maxbits'} == 32) {
+  $_Config{'maxparam'}    = 4294967295;
+  $_Config{'maxdigits'}   = 10;
+  $_Config{'maxprime'}    = 4294967291;
+  $_Config{'maxprimeidx'} = 203280221;
+} else {
+  $_Config{'maxparam'}    = 18446744073709551615;
+  $_Config{'maxdigits'}   = 20;
+  $_Config{'maxprime'}    = 18446744073709551557;
+  $_Config{'maxprimeidx'} = 425656284035217743;
+}
 
+sub prime_get_config {
+  my %config = %_Config;
+  return \%config;
+}
+
+sub _validate_positive_integer {
+  my($n, $min, $max) = @_;
+  croak "Parameter must be defined" if !defined $n;
+  croak "Parameter '$n' must be a positive integer" if $n =~ tr/0123456789//c;
+  croak "Parameter '$n' must be >= $min" if defined $min && $n < $min;
+  croak "Parameter '$n' must be <= $max" if defined $max && $n > $max;
+  croak "Parameter '$n' outside of integer range" if $n > $_Config{'maxparam'}
+                                                  && ref($n) !~ /^Math::Big/;
+  1;
+}
 
 sub primes {
   my $optref = (ref $_[0] eq 'HASH')  ?  shift  :  {};
@@ -91,25 +119,16 @@ sub primes {
   croak "too many parameters to primes" unless scalar @_ <= 2;
   my $low = (@_ == 2)  ?  shift  :  2;
   my $high = shift;
-  my $sref = [];
 
-  # Validate parameters
-  if ( (!defined $low) || (!defined $high) ||
-       ($low =~ tr/0123456789//c) || ($high =~ tr/0123456789//c)
-     ) {
-    $low  = '<undef>' unless defined $low;
-    $high = '<undef>' unless defined $high;
-    croak "Parameters [ $low $high ] must be positive integers";
-  }
-
-  # Verify the parameters are in range.
-  if (!$_pure_perl || !defined $bigint::VERSION) {
-    croak "Parameters [ $low $high ] not in range 0-$_maxparam" unless $low <= 
$_maxparam && $high <= $_maxparam;
-  }
+  _validate_positive_integer($low);
+  _validate_positive_integer($high);
 
+  my $sref = [];
   return $sref if ($low > $high) || ($high < 2);
 
-  return Math::Prime::Util::PP::primes($low,$high) if $_pure_perl;
+  if ( (!$_Config{'xs'}) || ($high > $_Config{'maxparam'}) ) {
+    return Math::Prime::Util::PP::primes($low,$high);
+  }
 
   my $method = $optref->{'method'};
   $method = 'Dynamic' unless defined $method;
@@ -162,19 +181,11 @@ sub primes {
 sub random_prime {
   my $low = (@_ == 2)  ?  shift  :  2;
   my $high = shift;
-  if ( (!defined $low) || (!defined $high) ||
-       ($low =~ tr/0123456789//c) || ($high =~ tr/0123456789//c)
-     ) {
-    $low  = 'undef' unless defined $low;
-    $high = 'undef' unless defined $high;
-    croak "Parameters [ $low $high ] must be positive integers";
-  }
-  if (!$_pure_perl || !defined $bigint::VERSION) {
-    croak "Parameters [ $low $high ] not in range 0-$_maxparam" unless $low <= 
$_maxparam && $high <= $_maxparam;
-  }
-  $low = 2 if $low < 2;
+  _validate_positive_integer($low);
+  _validate_positive_integer($high);
 
-  # Make sure we have a valid range.
+  # Tighten the range to the nearest prime.
+  $low = 2 if $low < 2;
   $low = next_prime($low - 1);
   $high = ($high < ~0)  ?  prev_prime($high + 1)  :  prev_prime($high);
   return $low if ($low == $high) && is_prime($low);
@@ -220,15 +231,12 @@ sub random_prime {
 my @_random_ndigit_ranges;
 
 sub random_ndigit_prime {
-  my $digits = shift;
-  # TODO: bigint with many digits
-  if ((!defined $digits) || ($digits > $_maxdigits) || ($digits < 1)) {
-    croak "Digits must be between 1 and $_maxdigits";
-  }
+  my($digits) = @_;
+  _validate_positive_integer($digits, 1, (defined $bigint::VERSION) ? 10000 : 
$_Config{'maxdigits'} );
   if (!defined $_random_ndigit_ranges[$digits]) {
     my $low = ($digits == 1) ? 0 : int(10 ** ($digits-1));
     my $high = int(10 ** $digits);
-    $high = ~0 if $high > ~0;
+    $high = ~0 if $high > ~0 && !defined $bigint::VERSION;
     $_random_ndigit_ranges[$digits] = [next_prime($low), prev_prime($high)];
   }
   my ($low, $high) = @{$_random_ndigit_ranges[$digits]};
@@ -237,7 +245,6 @@ sub random_ndigit_prime {
 
 # Perhaps a random_nbit_prime ?   Definition?
 
-
 sub all_factors {
   my $n = shift;
   my @factors = factor($n);
@@ -255,6 +262,126 @@ sub all_factors {
   return @factors;
 }
 
+# A008683 Moebius function mu(n)
+# A030059, A013929, A030229, A002321, A005117, A013929 all relate.
+
+# One can argue for the Omega function (A001221), Euler Phi (A000010), and
+# Merten's functions also.
+
+sub moebius {
+  my($n) = @_;
+  _validate_positive_integer($n, 1);
+  return if $n <= 0;
+  return 1 if $n == 1;
+
+  # Quick check for small replicated factors
+  return 0 if ($n >= 25) && (($n % 4) == 0 || ($n % 9) == 0 || ($n % 25) == 0);
+
+  my @factors = factor($n);
+  my %all_factors;
+  foreach my $factor (@factors) {
+    return 0 if $all_factors{$factor}++;
+  }
+  return (((scalar @factors) % 2) == 0) ? 1 : -1;
+}
+
+sub euler_phi {
+  my($n) = @_;
+  _validate_positive_integer($n);
+  return 1 if $n <= 1;  # phi(0) is disputed
+
+  my %factor_mult;
+  my @factors = grep { !$factor_mult{$_}++ } factor($n);
+  my $totient = $n;
+  foreach my $factor (@factors) {
+    # These divisions should always be exact
+    $totient = int($totient/$factor) * ($factor-1);
+  }
+  $totient;
+}
+
+
+#############################################################################
+# Front ends to functions.
+#
+# These will do input validation, then call the appropriate internal function
+# based on the input (XS, GMP, PP).
+#############################################################################
+
+#
+# This works, but sadly has a lot of overhead -- 7x more overhead than the
+# entire is_prime C function for n under 10M or so.
+#
+#sub is_prime {
+#  my($n) = @_;
+#  _validate_positive_integer($n);
+#
+#  return XS_is_prime($n) if $_Config{'xs'} && $n <= $_Config{'maxparam'};
+#  return Math::Prime::Util::PP::is_prime($n);
+#}
+
+sub factor {
+  my($n) = @_;
+  _validate_positive_integer($n);
+  
+  return XS_factor($n) if $_Config{'xs'} && $n <= $_Config{'maxparam'};
+
+  $n = $n->as_number() if ref($n) =~ /^Math::BigFloat/;
+  $n = $n->numify if $n <= ~0;
+  Math::Prime::Util::PP::factor($n);
+}
+
+
+sub RiemannR {
+  my($n) = @_;
+  croak("Invalid input to ReimannR:  x must be > 0") if $n <= 0;
+
+  if (ref($n) eq 'Math::BigInt' && !defined $Math::BigFloat::VERSION) {
+    require Math::BigFloat;
+    $n = new Math::BigFloat "$n";
+  }
+
+  return Math::Prime::Util::PP::RiemannR($n, 1e-30) if ref($n) =~ /^Math::Big/;
+  return Math::Prime::Util::PP::RiemannR($n) if !$_Config{'xs'};
+  return XS_RiemannR($n);
+
+  # We could make a new object, like:
+  #    require Math::BigFloat;
+  #    my $bign = new Math::BigFloat "$n";
+  #    my $result = Math::Prime::Util::PP::RiemannR($bign);
+  #    return $result;
+}
+
+sub ExponentialIntegral {
+  my($n) = @_;
+  croak "Invalid input to ExponentialIntegral:  x must be != 0" if $n == 0;
+
+  if (ref($n) eq 'Math::BigInt' && !defined $Math::BigFloat::VERSION) {
+    require Math::BigFloat;
+    $n = new Math::BigFloat "$n";
+  }
+
+  return Math::Prime::Util::PP::ExponentialIntegral($n, 1e-30) if ref($n) =~ 
/^Math::Big/;
+  return Math::Prime::Util::PP::ExponentialIntegral($n) if !$_Config{'xs'};
+  return XS_ExponentialIntegral($n);
+}
+
+sub LogarithmicIntegral {
+  my($n) = @_;
+  return 0 if $n == 0;
+  croak("Invalid input to LogarithmicIntegral:  x must be >= 0") if $n <= 0;
+
+  if (defined $Math::BigFloat::VERSION) {
+    return Math::BigFloat->binf('-') if $n == 1;
+    return 
Math::BigFloat->new('1.045163780117492784844588889194613136522615578151201575832909144075013205210359530172717405626383356306')
 if $n == 2;
+  } else {
+    return 0+'-inf' if $n == 1;
+    return 1.045163780117492784844588889194613136522615578151 if $n == 2;
+  }
+  ExponentialIntegral(log($n));
+}
+
+
 use Math::Prime::Util::MemFree;
 
 1;
@@ -619,6 +746,16 @@ digits between 1 and the maximum native type (10 for 
32-bit, 20 for 64-bit).
 One of the primes within that range (e.g. 1000 - 9999 for 4-digits) will be
 uniformly selected using the L<rand> function.
 
+=head2 moebius
+
+  say "$n is square free" if moebius($n) != 0;
+  $sum += moebius($_) for (1..200); say "Mertens(200) = $sum";
+
+Returns the Möbius function (also called the Moebius, Mobius, or MoebiusMu
+function) for a positive non-zero integer input.  This function is 1 if
+C<n = 1>, 0 if C<n> is not square free (i.e. C<n> has a repeated factor),
+and C<-1^t> if C<n> is a product of C<t> distinct primes.  This is an
+important function in prime number theory.
 
 
 =head1 UTILITY FUNCTIONS
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 84613d6..74047d9 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -584,6 +584,18 @@ sub _mulmod {
   $r;
 }
 
+sub _native_powmod {
+  my($n, $power, $m) = @_;
+  my $t = 1;
+  $n = $n % $m;
+  while ($power) {
+    $t = ($t * $n) % $m if ($power & 1) != 0;
+    $n = ($n * $n) % $m;
+    $power >>= 1;
+  }
+  $t;
+}
+
 sub _powmod {
   my($n, $power, $m) = @_;
   my $t = 1;
@@ -639,20 +651,41 @@ sub miller_rabin {
     $d >>= 1;
   }
 
-  foreach my $a (@bases) {
-    croak "Base $a is invalid" if $a < 2;
-    my $x = _powmod($a, $d, $n);
-    next if ($x == 1) || ($x == ($n-1));
+  if ( ($n < $_half_word) || (ref($n) =~ /^Math::Big/) ) {
+
+    foreach my $a (@bases) {
+      croak "Base $a is invalid" if $a < 2;
+      my $x = _native_powmod($a, $d, $n);
+      next if ($x == 1) || ($x == ($n-1));
+      foreach my $r (1 .. $s) {
+        $x = ($x*$x) % $n;
+        return 0 if $x == 1;
+        if ($x == ($n-1)) {
+          $a = 0;
+          last;
+        }
+      }
+      return 0 if $a != 0;
+    }
 
-    foreach my $r (1 .. $s) {
-      $x = ($x < $_half_word) ? ($x*$x) % $n : _mulmod($x, $x, $n);
-      return 0 if $x == 1;
-      if ($x == ($n-1)) {
-        $a = 0;
-        last;
+  } else {
+
+    foreach my $a (@bases) {
+      croak "Base $a is invalid" if $a < 2;
+      my $x = _powmod($a, $d, $n);
+      next if ($x == 1) || ($x == ($n-1));
+
+      foreach my $r (1 .. $s) {
+        $x = ($x < $_half_word) ? ($x*$x) % $n : _mulmod($x, $x, $n);
+        return 0 if $x == 1;
+        if ($x == ($n-1)) {
+          $a = 0;
+          last;
+        }
       }
+      return 0 if $a != 0;
     }
-    return 0 if $a != 0;
+
   }
   1;
 }
@@ -682,10 +715,15 @@ sub is_prob_prime {
 sub _basic_factor {
   # MODIFIES INPUT SCALAR
   return ($_[0]) if $_[0] < 4;
+
   my @factors;
-  while ( ($_[0] % 2) == 0 ) { push @factors, 2;  $_[0] /= 2; }
-  while ( ($_[0] % 3) == 0 ) { push @factors, 3;  $_[0] /= 3; }
-  while ( ($_[0] % 5) == 0 ) { push @factors, 5;  $_[0] /= 5; }
+  while ( ($_[0] % 2) == 0 ) { push @factors, 2;  $_[0] = int($_[0] / 2); }
+  while ( ($_[0] % 3) == 0 ) { push @factors, 3;  $_[0] = int($_[0] / 3); }
+  while ( ($_[0] % 5) == 0 ) { push @factors, 5;  $_[0] = int($_[0] / 5); }
+
+  # Stop using bignum if we can
+  $_[0] = $_[0]->numify if ref($_[0]) =~ /^Math::Big/ && $_[0] <= ~0;
+
   if ( ($_[0] > 1) && _is_prime7($_[0]) ) {
     push @factors, $_[0];
     $_[0] = 1;
@@ -706,7 +744,7 @@ sub trial_factor {
     if ( ($n % $f) == 0) {
       while ( ($n % $f) == 0) {
         push @factors, $f;
-        $n /= $f;
+        $n = int($n/$f);
       }
       $limit = int( sqrt($n) + 0.001);
     }
@@ -725,13 +763,16 @@ sub factor {
   my @factors = _basic_factor($n);
   return @factors if $n < 4;
 
-  while (($n %  7) == 0) { push @factors,  7;  $n /=  7; }
-  while (($n % 11) == 0) { push @factors, 11;  $n /= 11; }
-  while (($n % 13) == 0) { push @factors, 13;  $n /= 13; }
-  while (($n % 17) == 0) { push @factors, 17;  $n /= 17; }
-  while (($n % 19) == 0) { push @factors, 19;  $n /= 19; }
-  while (($n % 23) == 0) { push @factors, 23;  $n /= 23; }
-  while (($n % 29) == 0) { push @factors, 29;  $n /= 29; }
+  while (($n %  7) == 0) { push @factors,  7;  $n = int($n /  7); }
+  while (($n % 11) == 0) { push @factors, 11;  $n = int($n / 11); }
+  while (($n % 13) == 0) { push @factors, 13;  $n = int($n / 13); }
+  while (($n % 17) == 0) { push @factors, 17;  $n = int($n / 17); }
+  while (($n % 19) == 0) { push @factors, 19;  $n = int($n / 19); }
+  while (($n % 23) == 0) { push @factors, 23;  $n = int($n / 23); }
+  while (($n % 29) == 0) { push @factors, 29;  $n = int($n / 29); }
+
+  # Stop using bignum if possible
+  $n = $n->numify if ref($n) =~ /^Math::Big/ && $n <= ~0;
 
   my @nstack = ($n);
   while (@nstack) {
@@ -885,8 +926,8 @@ my $_const_euler = 
0.57721566490153286060651209008240243104215933593992;
 my $_const_li2 = 1.045163780117492784844588889194613136522615578151;
 
 sub ExponentialIntegral {
-  my($x) = @_;
-  my $tol = 1e-16;
+  my($x, $tol) = @_;
+  $tol = 1e-16 unless defined $tol;
   my $sum = 0.0;
   my($y, $t);
   my $c = 0.0;
@@ -1032,8 +1073,8 @@ sub _evaluate_zeta {
 
 # Riemann R function
 sub RiemannR {
-  my($x) = @_;
-  my $tol = 1e-16;
+  my($x, $tol) = @_;
+  $tol = 1e-16 unless defined $tol;
   my $sum = 0.0;
   my($y, $t);
   my $c = 0.0;
diff --git a/t/15-probprime.t b/t/15-probprime.t
index 79a8779..67c14ee 100644
--- a/t/15-probprime.t
+++ b/t/15-probprime.t
@@ -5,7 +5,7 @@ use warnings;
 use Test::More;
 use Math::Prime::Util qw/is_prime is_prob_prime miller_rabin/;
 
-my $use64 = Math::Prime::Util::_maxbits > 32;
+my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
 my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
 
 plan tests => 6 + 19 + 3573 + (5 + 29 + 22 + 23 + 16) + 15 + 27
diff --git a/t/19-moebius.t b/t/19-moebius.t
new file mode 100644
index 0000000..9ab422c
--- /dev/null
+++ b/t/19-moebius.t
@@ -0,0 +1,55 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+
+use Test::More;
+use Math::Prime::Util qw/moebius euler_phi/;
+
+my @moeb_vals = (qw/ 1 -1 -1 0 -1 1 -1 0 0 1 -1 0 -1 1 1 0 -1 0 -1 0 /);
+my %mertens = (
+        1 =>    1,
+        2 =>    0,
+        3 =>   -1,
+        4 =>   -1,
+        5 =>   -2,
+       10 =>   -1,
+      100 =>    1,
+     1000 =>    2,
+    10000 =>  -23,
+   100000 =>  -48,
+#  1000000 =>  212,
+# 10000000 => 1037,
+);
+
+my %totients = (
+          1 => 1,
+          2 => 1,
+         20 => 8,
+         36 => 12,
+     123456 => 41088, 
+     123457 => 123456,
+  123456789 => 82260072,
+);
+
+plan tests => 0 + 1
+                + scalar @moeb_vals
+                + scalar(keys %mertens)
+                + scalar(keys %totients);
+
+ok(!eval { moebius(0); }, "moebius(0)");
+
+my $i = 1;
+foreach my $m (@moeb_vals) {
+  is( moebius($i), $m, "moebius($i) == $m" );
+  $i++;
+}
+
+while (my($n, $mertens) = each (%mertens)) {
+  my $M = 0;
+  $M += moebius($_) for (1 .. $n);
+  is( $M, $mertens, "Mertens($n) == $mertens" );
+}
+
+while (my($n, $phi) = each (%totients)) {
+  is( euler_phi($n), $phi, "euler_phi($n) == $phi" );
+}
diff --git a/util.c b/util.c
index 749ddba..174d68f 100644
--- a/util.c
+++ b/util.c
@@ -508,11 +508,12 @@ UV prime_count_approx(UV x)
    * Riemann's R function works even better, with errors of ~1828 at 10^10
    * and 24M at 10^19.
    *
-   *    Method           10^10 %error  10^19 %error
-   *    ---------------  ------------  ------------
-   *    average bounds    .01%          .0002%
-   *    li(n)             .0007%        .00000004%
-   *    R(n)              .0004%        .00000001%
+   *    Method             10^10 %error  10^19 %error
+   *    -----------------  ------------  ------------
+   *    average bounds      .01%          .0002%
+   *    li(n)               .0007%        .00000004%
+   *    li(n)-li(n^.5)/2    .0004%        .00000001%
+   *    R(n)                .0004%        .00000001%
    *
    * Getting fancier, one try using Riemann's pi formula:
    *     http://trac.sagemath.org/sage_trac/ticket/8135
@@ -521,7 +522,8 @@ UV prime_count_approx(UV x)
   if (x < NPRIME_COUNT_SMALL)
     return prime_count_small[x];
 
-  R = RiemannR(x);
+  //R = XS_LogarithmicIntegral(x) - XS_LogarithmicIntegral(sqrt(x))/2;
+  R = XS_RiemannR(x);
 
   /* We could add the additional factor:
    *   R = R - (1.0 / log(x)) + (M_1_PI * atan(M_PI/log(x)))
@@ -833,7 +835,7 @@ UV nth_prime(UV n)
 static double const euler_mascheroni = 
0.57721566490153286060651209008240243104215933593992;
 static double const li2 = 1.045163780117492784844588889194613136522615578151;
 
-double ExponentialIntegral(double x) {
+double XS_ExponentialIntegral(double x) {
   double const tol = 1e-16;
   double val, term, fact_n;
   double y, t;
@@ -925,12 +927,12 @@ double ExponentialIntegral(double x) {
   return val;
 }
 
-double LogarithmicIntegral(double x) {
+double XS_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 LogarithmicIntegral:  x must be > 0");
-  return ExponentialIntegral(log(x));
+  return XS_ExponentialIntegral(log(x));
 }
 
 /*
@@ -1006,7 +1008,7 @@ static double evaluate_zeta(double x) {
   return sum;
 }
 
-double RiemannR(double x) {
+double XS_RiemannR(double x) {
   double const tol = 1e-16;
   double y, t, part_term, term, flogx, zeta;
   double sum = 0.0;
diff --git a/util.h b/util.h
index a506ec2..b93c075 100644
--- a/util.h
+++ b/util.h
@@ -19,9 +19,9 @@ 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);
+extern double XS_ExponentialIntegral(double x);
+extern double XS_LogarithmicIntegral(double x);
+extern double XS_RiemannR(double x);
 
 /* Above this value, is_prime will do deterministic Miller-Rabin */
 /* With 64-bit math, we can do much faster mulmods from 2^16-2^32 */

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