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

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

commit cc71cd27829abb53b6523bd3b7811858ebe9f796
Author: Dana Jacobsen <d...@acm.org>
Date:   Sat Jun 23 11:47:20 2012 -0600

    Pure Perl working on all tests
---
 Changes                   |   2 +-
 lib/Math/Prime/Util.pm    |  27 +++--
 lib/Math/Prime/Util/PP.pm | 267 +++++++++++++++++++++++++++++++++++++++++++---
 util.c                    |   2 +-
 4 files changed, 273 insertions(+), 25 deletions(-)

diff --git a/Changes b/Changes
index 99cde38..0a6c18f 100644
--- a/Changes
+++ b/Changes
@@ -1,7 +1,7 @@
 Revision history for Perl extension Math::Prime::Util.
 
 0.09  23 June 2012
-    - Started work on pure perl code.
+    - Pure Perl code.  Passes all tests, but 1 to 13000x slower.
 
 0.08  22 June 2012
     - Added thread safety and tested good concurrency.
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index d2d5e83..565396b 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -69,7 +69,9 @@ BEGIN {
     *prho_factor    = \&Math::Prime::Util::PP::prho_factor;
     *pminus1_factor = \&Math::Prime::Util::PP::pminus1_factor;
 
-    # TODO: Ei(x), li(x), R(x)
+    *RiemannR            = \&Math::Prime::Util::PP::RiemannR;
+    *LogarithmicIntegral = \&Math::Prime::Util::PP::LogarithmicIntegral;
+    *ExponentialIntegral = \&Math::Prime::Util::PP::ExponentialIntegral;
     # TODO: prime_count is horribly, horribly slow
     # TODO: We should have some tests to verify XS vs. PP.
   }
@@ -85,7 +87,7 @@ my $_maxprimeidx=(_maxbits == 32) ?  203280221 :   
425656284035217743;
 
 
 sub primes {
-  my $optref = {};  $optref = shift if ref $_[0] eq 'HASH';
+  my $optref = (ref $_[0] eq 'HASH')  ?  shift  :  {};
   croak "no parameters to primes" unless scalar @_ > 0;
   croak "too many parameters to primes" unless scalar @_ <= 2;
   my $low = (@_ == 2)  ?  shift  :  2;
@@ -96,8 +98,8 @@ sub primes {
   if ( (!defined $low) || (!defined $high) ||
        ($low =~ tr/0123456789//c) || ($high =~ tr/0123456789//c)
      ) {
-    $low  = 'undef' unless defined $low;
-    $high = 'undef' unless defined $high;
+    $low  = '<undef>' unless defined $low;
+    $high = '<undef>' unless defined $high;
     croak "Parameters [ $low $high ] must be positive integers";
   }
 
@@ -155,6 +157,7 @@ sub primes {
   return $sref;
 }
 
+
 sub random_prime {
   my $low = (@_ == 2)  ?  shift  :  2;
   my $high = shift;
@@ -170,8 +173,8 @@ sub random_prime {
 
   # Make sure we have a valid range.
   # TODO: this is is killing performance with large numbers
-  $low = next_prime($low-1);
-  $high = ($high < ~0) ? prev_prime($high+1) : prev_prime($high);
+  $low = next_prime($low - 1);
+  $high = ($high < ~0)  ?  prev_prime($high + 1)  :  prev_prime($high);
   return $low if ($low == $high) && is_prime($low);
   return if $low >= $high;
 
@@ -198,18 +201,19 @@ sub random_prime {
       do {
         $prime = $low + $irandf->($range);
         croak "Random function broken?" if $loop_limit-- < 0;
-      }  while ( !($prime%2) || !($prime%3) || !is_prime($prime) );
+      }  while ( !($prime % 2 ) || !($prime % 3) || !is_prime($prime) );
     } else {
       do {
         my $rand = ( ($irandf->(4294967295) << 32) + $irandf->(4294967295) ) % 
$range;
         $prime = $low + $rand;
         croak "Random function broken?" if $loop_limit-- < 0;
-      }  while ( !($prime%2) || !($prime%3) || !is_prime($prime) );
+      }  while ( !($prime % 2) || !($prime % 3) || !is_prime($prime) );
     }
   }
   return $prime;
 }
 
+
 sub random_ndigit_prime {
   my $digits = shift;
   if ((!defined $digits) || ($digits > $_maxdigits) || ($digits < 1)) {
@@ -223,12 +227,14 @@ sub random_ndigit_prime {
 
 # Perhaps a random_nbit_prime ?   Definition?
 
+
 sub all_factors {
   my $n = shift;
   my @factors = factor($n);
   my %all_factors;
   foreach my $f1 (@factors) {
     next if $f1 >= $n;
+    # We're adding to %all_factors in the loop, so grab the keys now.
     my @all = keys %all_factors;;
     foreach my $f2 (@all) {
       $all_factors{$f1*$f2} = 1 if ($f1*$f2) < $n;
@@ -835,21 +841,26 @@ C<is_prime>: my impressions:
    Math::Primality           Very slow      Very slow
 
 The differences are in the implementations:
+
    - L<Math::Prime::FastSieve> only works in a sieved range, which is really
      fast if you can do it (M::P::U will do the same if you call
      C<prime_precalc>).  Larger inputs just need too much time and memory
      for the sieve.
+
    - L<Math::Primality> uses a GMP BPSW test which is overkill for our 64-bit
      range.  It's generally an order of magnitude or two slower than any
      of the others.  
+
    - L<Math::Pari> has some very effective code, but it has some overhead to 
get
      to it from Perl.  That means for small numbers it is relatively slow: an
      order of magnitude slower than M::P::XS and M::P::Util (though arguably
      this is only important for benchmarking since "slow" is ~2 microseconds).
      Large numbers transition over to smarter tests so don't slow down much.
+
    - L<Math::Prime::XS> does trial divisions, which is wonderful if the input
      has a small factor (or is small itself).  But it can take 1000x longer
      if given a large prime.
+
    - L<Math::Prime::Util> looks in the sieve for a fast bit lookup if that
      exists (default up to 30,000 but it can be expanded, e.g.
      C<prime_precalc>), uses trial division for numbers higher than this but
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 8e27c3f..93b604d 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -99,6 +99,44 @@ sub is_prime {
   2*_is_prime7($n);
 }
 
+sub _sieve_erat {
+  my($end) = @_;
+  my $last = int(($end+1)/2);
+
+  my $sieve = '';
+  my $n = 3;
+  while ( ($n*$n) <= $end ) {
+    my $s = $n*$n;
+    while ($s <= $end) {
+      vec($sieve, $s >> 1, 1) = 1;
+      $s += 2*$n;
+    }
+    do { $n += 2 } while vec($sieve, $n >> 1, 1) != 0;
+  }
+  # Mark 1 as composite
+  vec($sieve, 0, 1) = 1;
+  $sieve;
+}
+# Uses 8x more memory, but about 50% faster
+sub _sieve_erat_array {
+  my($end) = @_;
+  my $last = int(($end+1)/2);
+
+  my @sieve;
+  my $n = 3;
+  while ( ($n*$n) <= $end ) {
+    my $s = $n*$n;
+    while ($s <= $end) {
+      $sieve[$s>>1] = 1;
+      $s += 2*$n;
+    }
+    do { $n += 2 } while $sieve[$n>>1];
+  }
+  # Mark 1 as composite
+  $sieve[0] = 1;
+  \@sieve;
+}
+
 sub primes {
   my $optref = (ref $_[0] eq 'HASH')  ?  shift  :  {};
   croak "no parameters to primes" unless scalar @_ > 0;
@@ -119,14 +157,23 @@ sub primes {
   push @$sref, 5  if ($low <= 5) && ($high >= 5);
   $low = 7 if $low < 7;
 
-  my $n = $low;
-  my $base = 30 * int($n/30);
-  my $in = 0;  $in++ while ($n - $base) > $_prime_indices[$in];
-  $n = $base + $_prime_indices[$in];
-  while ($n <= $high) {
-    push @$sref, $n  if _is_prime7($n);
-    if (++$in == 8) {  $base += 30; $in = 0;  }
+  if ($low == 7) {
+    my $sieve = _sieve_erat_array($high);
+    my $n = 7;
+    while ($n <= $high) {
+      push @$sref, $n  if !$sieve->[$n>>1];
+      $n += 2;
+    }
+  } else {
+    my $n = $low;
+    my $base = 30 * int($n/30);
+    my $in = 0;  $in++ while ($n - $base) > $_prime_indices[$in];
     $n = $base + $_prime_indices[$in];
+    while ($n <= $high) {
+      push @$sref, $n  if _is_prime7($n);
+      if (++$in == 8) {  $base += 30; $in = 0;  }
+      $n = $base + $_prime_indices[$in];
+    }
   }
   $sref;
 }
@@ -203,14 +250,27 @@ sub prime_count {
   $count++ if ($low <= 5) && ($high >= 5);
   $low = 7 if $low < 7;
 
-  my $n = $low;
-  my $base = 30 * int($n/30);
-  my $in = 0;  $in++ while ($n - $base) > $_prime_indices[$in];
-  $n = $base + $_prime_indices[$in];
-  while ($n <= $high) {
-    $count++ if _is_prime7($n);
-    if (++$in == 8) {  $base += 30; $in = 0;  }
+  if ($low == 7) {
+    my $sieve = _sieve_erat_array($high);
+    my $n = 7;
+    while ($n <= $high) {
+      $count++ if !$sieve->[$n>>1];
+      $n += 4;
+      if ($n <= $high) {
+        $count++ if !$sieve->[$n>>1];
+        $n += 2;
+      }
+    }
+  } else {
+    my $n = $low;
+    my $base = 30 * int($n/30);
+    my $in = 0;  $in++ while ($n - $base) > $_prime_indices[$in];
     $n = $base + $_prime_indices[$in];
+    while ($n <= $high) {
+      $count++ if _is_prime7($n);
+      if (++$in == 8) {  $base += 30; $in = 0;  }
+      $n = $base + $_prime_indices[$in];
+    }
   }
   $count;
 }
@@ -290,7 +350,8 @@ sub prime_count_approx {
 
   return $_prime_count_small[$x] if $x <= $#_prime_count_small;
 
-  return int( (prime_count_upper($x) + prime_count_lower($x)) / 2);
+  #return int( (prime_count_upper($x) + prime_count_lower($x)) / 2);
+  return int(RiemannR($x)+0.5);
 }
 
 sub nth_prime_lower {
@@ -542,6 +603,182 @@ sub pbrent_factor { trial_factor(@_) }
 sub prho_factor { trial_factor(@_) }
 sub pminus1_factor { trial_factor(@_) }
 
+
+my $_const_euler = 0.57721566490153286060651209008240243104215933593992;
+my $_const_li2 = 1.045163780117492784844588889194613136522615578151;
+
+sub ExponentialIntegral {
+  my($x) = @_;
+  my $tol = 1e-16;
+  my $sum = 0.0;
+  my($y, $t);
+  my $c = 0.0;
+
+  croak "Invalid input to ExponentialIntegral:  x must be != 0" if $x == 0;
+
+  my $val; # The result from one of the four methods
+
+  if ($x < -1) {
+    # Continued fraction
+    my $lc = 0;
+    my $ld = 1 / (1 - $x);
+    $val = $ld * (-exp($x));
+    for my $n (1 .. 100000) {
+      $lc = 1 / (2*$n + 1 - $x - $n*$n*$lc);
+      $ld = 1 / (2*$n + 1 - $x - $n*$n*$ld);
+      my $old = $val;
+      $val *= $ld/$lc;
+      last if abs($val - $old) <= ($tol * abs($val));
+    }
+  } elsif ($x < 0) {
+    # Rational Chebyshev approximation
+    my @C6p = ( -148151.02102575750838086,
+                 150260.59476436982420737,
+                  89904.972007457256553251,
+                  15924.175980637303639884,
+                   2150.0672908092918123209,
+                    116.69552669734461083368,
+                      5.0196785185439843791020);
+    my @C6q = (  256664.93484897117319268,
+                 184340.70063353677359298,
+                  52440.529172056355429883,
+                   8125.8035174768735759866,
+                    750.43163907103936624165,
+                     40.205465640027706061433,
+                      1.0000000000000000000000);
+    my $sumn = 
$C6p[0]-$x*($C6p[1]-$x*($C6p[2]-$x*($C6p[3]-$x*($C6p[4]-$x*($C6p[5]-$x*$C6p[6])))));
+    my $sumd = 
$C6q[0]-$x*($C6q[1]-$x*($C6q[2]-$x*($C6q[3]-$x*($C6q[4]-$x*($C6q[5]-$x*$C6q[6])))));
+    $val = log(-$x) - ($sumn / $sumd);
+  } elsif ($x < -log($tol)) {
+    # Convergent series
+    my $fact_n = 1;
+    $y = $_const_euler-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+    $y = log($x)-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+    for my $n (1 .. 200) {
+      $fact_n *= $x/$n;
+      my $term = $fact_n / $n;
+      $y = $term-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+      last if $term < $tol;
+    }
+    $val = $sum;
+  } else {
+    # Asymptotic divergent series
+    $val = exp($x) / $x;
+    $y = 1.0-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+    my $term = 1.0;
+    for my $n (1 .. 200) {
+      my $last_term = $term;
+      $term *= $n/$x;
+      last if $term < $tol;
+      if ($term < $last_term) {
+        $y = $term-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+      } else {
+        $y = (-$last_term/3)-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+        last;
+      }
+    }
+    $val *= $sum;
+  }
+  $val;
+}
+
+sub LogarithmicIntegral {
+  my($x) = @_;
+  return 0 if $x == 0;
+  return 0+(-Infinity) if $x == 1;
+  return $_const_li2 if $x == 2;
+  croak "Invalid input to LogarithmicIntegral:  x must be > 0" if $x <= 0;
+  ExponentialIntegral(log($x));
+}
+
+# Riemann Zeta function for integers, used for computing Riemann R
+my @_Riemann_Zeta_Table = (
+  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,
+  0.0000038172932649998398564616446219397,
+  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,
+);
+
+# Compute Riemann Zeta function.  Slow and inaccurate near x = 1, but improves
+# very rapidly (x = 5 is quite reasonable).
+sub _evaluate_zeta {
+  my($x) = @_;
+  my $tol = 1e-16;
+  my $sum = 0.0;
+  my($y, $t);
+  my $c = 0.0;
+
+  $y = (1.0/(2.0**$x))-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+
+  for my $k (3 .. 100000) {
+    my $term = 1.0 / ($k ** $x);
+    $y = $term-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+    last if abs($term) < $tol;
+  }
+  $sum;
+}
+
+# Riemann R function
+sub RiemannR {
+  my($x) = @_;
+  my $tol = 1e-16;
+  my $sum = 0.0;
+  my($y, $t);
+  my $c = 0.0;
+
+  croak "Invalid input to ReimannR:  x must be > 0" if $x <= 0;
+
+  $y = 1.0-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+  my $flogx = log($x);
+  my $part_term = 1.0;
+
+  for my $k (1 .. 10000) {
+    # Small k from table, larger k from function
+    my $zeta = ($k <= $#_Riemann_Zeta_Table) ? $_Riemann_Zeta_Table[$k+1-2]
+                                             : _evaluate_zeta($k+1);
+    $part_term *= $flogx / $k;
+    my $term = $part_term / ($k + $k * $zeta);
+    $y = $term-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+    last if abs($term) < $tol;
+  }
+  $sum;
+}
+
 1;
 
 __END__
diff --git a/util.c b/util.c
index a318f53..d3d040a 100644
--- a/util.c
+++ b/util.c
@@ -919,7 +919,7 @@ 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");
+  if (x <= 0) croak("Invalid input to LogarithmicIntegral:  x must be > 0");
   return ExponentialIntegral(log(x));
 }
 

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