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

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

commit f3e7800385a667316c535f605e834e090ce5289b
Author: Dana Jacobsen <d...@acm.org>
Date:   Wed Jan 22 21:16:32 2014 -0800

    More moving functions out of Util
---
 XS.xs                       |  59 +++++++++-------
 lib/Math/Prime/Util.pm      | 168 ++++++++++++++++----------------------------
 lib/Math/Prime/Util/PP.pm   | 134 +++++++++++++++++++++++++++--------
 lib/Math/Prime/Util/PPFE.pm | 153 ++++++++++++++++++++++++++++++++++++----
 t/80-pp.t                   |   2 +-
 t/92-release-pod-coverage.t |   4 +-
 util.c                      |  27 +++----
 7 files changed, 352 insertions(+), 195 deletions(-)

diff --git a/XS.xs b/XS.xs
index 5aee2e3..a7ccee5 100644
--- a/XS.xs
+++ b/XS.xs
@@ -686,7 +686,7 @@ factor(IN SV* svn)
       switch (ix) {
         case 0:  _vcallsubn(aTHX_ gimme_v, VCALL_ROOT, "_generic_factor", 1);  
   break;
         case 1:  _vcallsubn(aTHX_ gimme_v, VCALL_ROOT, "_generic_factor_exp", 
1); break;
-        default: _vcallsubn(aTHX_ gimme_v, VCALL_ROOT, "_generic_divisors", 
1);   break;
+        default: _vcallsubn(aTHX_ gimme_v, VCALL_GMP|VCALL_PP, "divisors", 1); 
  break;
       }
       return; /* skip implicit PUTBACK */
     }
@@ -706,7 +706,7 @@ divisor_sum(IN SV* svn, ...)
       UV sigma = divisor_sum(n, k);
       if (sigma != 0)  XSRETURN_UV(sigma);   /* sigma 0 means overflow */
     }
-    _vcallsub("_generic_divisor_sum");
+    _vcallsub_with_gmp("divisor_sum");
     return; /* skip implicit PUTBACK */
 
 void
@@ -785,29 +785,22 @@ kronecker(IN SV* sva, IN SV* svb)
     _vcallsub_with_gmp("kronecker");
     return; /* skip implicit PUTBACK */
 
-double
+NV
 _XS_ExponentialIntegral(IN SV* x)
   ALIAS:
     _XS_LogarithmicIntegral = 1
     _XS_RiemannZeta = 2
     _XS_RiemannR = 3
-    _XS_chebyshev_theta = 4
-    _XS_chebyshev_psi = 5
   PREINIT:
-    double ret;
+    NV nv, ret;
   CODE:
-    if (ix < 4) {
-      NV nv = SvNV(x);
-      switch (ix) {
-        case 0: ret = (NV) _XS_ExponentialIntegral(nv); break;
-        case 1: ret = (NV) _XS_LogarithmicIntegral(nv); break;
-        case 2: ret = (NV) ld_riemann_zeta(nv); break;
-        case 3:
-        default:ret = (NV) _XS_RiemannR(nv); break;
-      }
-    } else {
-      UV uv = SvUV(x);
-      ret = (NV) chebyshev_function(uv, ix-4);
+    nv = SvNV(x);
+    switch (ix) {
+      case 0: ret = (NV) _XS_ExponentialIntegral(nv); break;
+      case 1: ret = (NV) _XS_LogarithmicIntegral(nv); break;
+      case 2: ret = (NV) ld_riemann_zeta(nv); break;
+      case 3:
+      default:ret = (NV) _XS_RiemannR(nv); break;
     }
     RETVAL = ret;
   OUTPUT:
@@ -866,12 +859,15 @@ void
 carmichael_lambda(IN SV* svn)
   ALIAS:
     mertens = 1
-    exp_mangoldt = 2
-    znprimroot = 3
+    liouville = 2
+    chebyshev_theta = 3
+    chebyshev_psi = 4
+    exp_mangoldt = 5
+    znprimroot = 6
   PREINIT:
     int status;
   PPCODE:
-    status = _validate_int(aTHX_ svn, (ix > 1) ? 1 : 0);
+    status = _validate_int(aTHX_ svn, (ix >= 5) ? 1 : 0);
     switch (ix) {
       case 0: if (status == 1) XSRETURN_UV(carmichael_lambda(my_svuv(svn)));
               _vcallsub_with_gmp("carmichael_lambda");
@@ -879,11 +875,24 @@ carmichael_lambda(IN SV* svn)
       case 1: if (status == 1) XSRETURN_IV(mertens(my_svuv(svn)));
               _vcallsub_with_pp("mertens");
               break;
-      case 2: if (status ==-1) XSRETURN_UV(1);
-              if (status == 1) XSRETURN_UV(exp_mangoldt(my_svuv(svn)));
-              _vcallsub("_generic_exp_mangoldt");
+      case 2: if (status == 1) {
+                UV factors[MPU_MAX_FACTORS+1];
+                int nfactors = factor(my_svuv(svn), factors);
+                RETURN_NPARITY( (nfactors & 1) ? -1 : 1 );
+              }
+              _vcallsub_with_gmp("liouville");
               break;
-      case 3:
+      case 3: if (status == 1) XSRETURN_NV(chebyshev_function(my_svuv(svn),0));
+              _vcallsub_with_pp("chebyshev_theta");
+              break;
+      case 4: if (status == 1) XSRETURN_NV(chebyshev_function(my_svuv(svn),1));
+              _vcallsub_with_pp("chebyshev_psi");
+              break;
+      case 5: if (status != 0)
+                XSRETURN_UV( (status == -1) ? 1 : exp_mangoldt(my_svuv(svn)) );
+              _vcallsub_with_gmp("exp_mangoldt");
+              break;
+      case 6:
       default:if (status != 0) {
                 UV r, n = my_svuv(svn);
                 if (status == -1) n = -(IV)n;
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 91ca2a8..79a64f1 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -94,15 +94,9 @@ BEGIN {
 
     *next_prime    = \&Math::Prime::Util::_generic_next_prime;
     *prev_prime    = \&Math::Prime::Util::_generic_prev_prime;
-    *exp_mangoldt  = \&Math::Prime::Util::_generic_exp_mangoldt;
     *prime_count   = \&Math::Prime::Util::_generic_prime_count;
-    *divisor_sum   = \&Math::Prime::Util::_generic_divisor_sum;
     *factor        = \&Math::Prime::Util::_generic_factor;
     *factor_exp    = \&Math::Prime::Util::_generic_factor_exp;
-    *divisors      = \&Math::Prime::Util::_generic_divisors;
-    *forprimes     = sub (&$;$) { _generic_forprimes(@_); }; ## no critic 
qw(ProhibitSubroutinePrototypes)
-    *forcomposites = sub (&$;$) { _generic_forcomposites(@_); }; ## no critic 
qw(ProhibitSubroutinePrototypes)
-    *fordivisors   = sub (&$) { _generic_fordivisors(@_); }; ## no critic 
qw(ProhibitSubroutinePrototypes)
   };
 
   # aliases for deprecated names.  Will eventually be removed.
@@ -299,6 +293,10 @@ sub primes {
   return $sref;
 }
 
+#############################################################################
+# Random primes.  These are front end functions that do input validation,
+# load the RandomPrimes module, and call its function.
+
 sub random_prime {
   my($low,$high) = @_;
   if (scalar @_ > 1) {
@@ -361,6 +359,29 @@ sub random_proven_prime_with_cert {
   return Math::Prime::Util::RandomPrimes::random_proven_prime_with_cert($bits);
 }
 
+sub miller_rabin_random {
+  my($n, $k, $seed) = @_;
+  _validate_num($n) || _validate_positive_integer($n);
+  _validate_num($k) || _validate_positive_integer($k);
+
+  return 1 if $k <= 0;
+  return (is_prob_prime($n) > 0) if $n < 100;
+  return 0 unless $n & 1;
+
+  return Math::Prime::Util::GMP::miller_rabin_random($n, $k)
+    if $_HAVE_GMP
+    && defined &Math::Prime::Util::GMP::miller_rabin_random;
+
+  require Math::Prime::Util::RandomPrimes;
+  return Math::Prime::Util::RandomPrimes::miller_rabin_random($n, $k, $seed);
+}
+
+
+#############################################################################
+# These functions almost always return bigints, so there is no XS
+# implementation.  Try to run the GMP version, and if it isn't available,
+# load PP and call it.
+
 sub primorial {
   my($n) = @_;
   _validate_num($n) || _validate_positive_integer($n);
@@ -396,19 +417,36 @@ sub consecutive_integer_lcm {
   return Math::Prime::Util::PP::consecutive_integer_lcm($n);
 }
 
-sub _generic_divisor_sum {
+# See 2011+ FLINT and Fredrik Johansson's work for state of the art.
+# Very crude timing estimates (ignores growth rates).
+#   Perl-comb   partitions(10^5)  ~ 300 seconds  ~200,000x slower than Pari
+#   GMP-comb    partitions(10^6)  ~ 120 seconds    ~1,000x slower than Pari
+#   Pari        partitions(10^8)  ~ 100 seconds
+#   Bober 0.6   partitions(10^9)  ~  20 seconds       ~50x faster than Pari
+#   Johansson   partitions(10^12) ~  10 seconds     >1000x faster than Pari
+sub partitions {
   my($n) = @_;
+  return 1 if defined $n && $n <= 0;
   _validate_num($n) || _validate_positive_integer($n);
+
+  if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::partitions) {
+    return _reftyped($_[0],Math::Prime::Util::GMP::partitions($n));
+  }
+
   require Math::Prime::Util::PP;
-  return Math::Prime::Util::PP::divisor_sum(@_);
+  return Math::Prime::Util::PP::partitions($n);
 }
 
-                                   # Need proto for the block
-sub _generic_forprimes (&$;$) {    ## no critic 
qw(ProhibitSubroutinePrototypes)
+
+#############################################################################
+# forprimes, forcomposites, fordivisors.
+# These are used when the XS code can't handle it.
+
+sub _generic_forprimes {
   my($sub, $beg, $end) = @_;
   if (!defined $end) { $end = $beg; $beg = 2; }
-  _validate_num($beg) || _validate_positive_integer($beg);
-  _validate_num($end) || _validate_positive_integer($end);
+  _validate_positive_integer($beg);
+  _validate_positive_integer($end);
   $beg = 2 if $beg < 2;
   {
     my $pp;
@@ -420,11 +458,11 @@ sub _generic_forprimes (&$;$) {    ## no critic 
qw(ProhibitSubroutinePrototypes)
   }
 }
 
-sub _generic_forcomposites(&$;$) { ## no critic 
qw(ProhibitSubroutinePrototypes)
+sub _generic_forcomposites {
   my($sub, $beg, $end) = @_;
   if (!defined $end) { $end = $beg; $beg = 4; }
-  _validate_num($beg) || _validate_positive_integer($beg);
-  _validate_num($end) || _validate_positive_integer($end);
+  _validate_positive_integer($beg);
+  _validate_positive_integer($end);
   $beg = 4 if $beg < 4;
   $end = Math::BigInt->new(''.~0) if ref($end) ne 'Math::BigInt' && $end == ~0;
   {
@@ -439,9 +477,9 @@ sub _generic_forcomposites(&$;$) { ## no critic 
qw(ProhibitSubroutinePrototypes)
   }
 }
 
-sub _generic_fordivisors (&$) {    ## no critic 
qw(ProhibitSubroutinePrototypes)
+sub _generic_fordivisors {
   my($sub, $n) = @_;
-  _validate_num($n) || _validate_positive_integer($n);
+  _validate_positive_integer($n);
   my @divisors = divisors($n);
   {
     my $pp;
@@ -453,6 +491,9 @@ sub _generic_fordivisors (&$) {    ## no critic 
qw(ProhibitSubroutinePrototypes)
   }
 }
 
+#############################################################################
+# Iterators
+
 sub prime_iterator {
   my($start) = @_;
   $start = 0 unless defined $start;
@@ -475,72 +516,6 @@ sub prime_iterator_object {
   return Math::Prime::Util::PrimeIterator->new($start);
 }
 
-# Exponential of Mangoldt function (A014963).
-# Return p if n = p^m [p prime, m >= 1], 1 otherwise.
-sub _generic_exp_mangoldt {
-  my($n) = @_;
-  return 1 if defined $n && $n <= 1;  # n <= 1
-  _validate_num($n) || _validate_positive_integer($n);
-
-  return 2 if ($n & ($n-1)) == 0;     # n power of 2
-  return 1 unless $n & 1;             # even n can't be p^m
-
-  my @pe = factor_exp($n);
-  return 1 if scalar @pe > 1;
-  return $pe[0]->[0];
-}
-
-sub liouville {
-  my($n) = @_;
-  _validate_num($n) || _validate_positive_integer($n);
-  my $l = (-1) ** scalar factor($n);
-  return $l;
-}
-
-# See 2011+ FLINT and Fredrik Johansson's work for state of the art.
-# Very crude timing estimates (ignores growth rates).
-#   Perl-comb   partitions(10^5)  ~ 300 seconds  ~200,000x slower than Pari
-#   GMP-comb    partitions(10^6)  ~ 120 seconds    ~1,000x slower than Pari
-#   Pari        partitions(10^8)  ~ 100 seconds
-#   Bober 0.6   partitions(10^9)  ~  20 seconds       ~50x faster than Pari
-#   Johansson   partitions(10^12) ~  10 seconds     >1000x faster than Pari
-sub partitions {
-  my($n) = @_;
-  return 1 if defined $n && $n <= 0;
-  _validate_num($n) || _validate_positive_integer($n);
-
-  if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::partitions) {
-    return _reftyped($_[0],Math::Prime::Util::GMP::partitions($n));
-  }
-
-  require Math::Prime::Util::PP;
-  return Math::Prime::Util::PP::partitions($n);
-}
-
-sub chebyshev_theta {
-  my($n) = @_;
-  _validate_num($n) || _validate_positive_integer($n);
-  return _XS_chebyshev_theta($n) if $n <= $_XS_MAXVAL;
-  my $sum = 0.0;
-  forprimes { $sum += log($_); } $n;
-  return $sum;
-}
-sub chebyshev_psi {
-  my($n) = @_;
-  _validate_num($n) || _validate_positive_integer($n);
-  return 0 if $n <= 1;
-  return _XS_chebyshev_psi($n) if $n <= $_XS_MAXVAL;
-  my ($sum, $logn, $sqrtn) = (0.0, log($n), int(sqrt($n)));
-  forprimes {
-    my $logp = log($_);
-    $sum += $logp * int($logn/$logp+1e-15);
-  } $sqrtn;
-  forprimes {
-    $sum += log($_);
-  } $sqrtn+1, $n;
-  return $sum;
-}
-
 #############################################################################
 # Front ends to functions.
 #
@@ -622,14 +597,6 @@ sub _generic_factor_exp {
   return (map { [$_, $exponents{$_}] } @factors);
 }
 
-sub _generic_divisors {
-  my($n) = @_;
-  _validate_num($n) || _validate_positive_integer($n);
-
-  require Math::Prime::Util::PP;
-  return Math::Prime::Util::PP::divisors($n);
-}
-
 
 sub lucas_sequence {
   my($n, $P, $Q, $k) = @_;
@@ -639,9 +606,11 @@ sub lucas_sequence {
     _validate_num($testP) || _validate_positive_integer($testP); }
   { my $testQ = (!defined $Q || $Q >= 0) ? $Q : -$Q;
     _validate_num($testQ) || _validate_positive_integer($testQ); }
+
   return _XS_lucas_sequence($n, $P, $Q, $k)
     if ref($_[0]) ne 'Math::BigInt' && $n <= $_XS_MAXVAL
     && ref($_[3]) ne 'Math::BigInt' && $k <= $_XS_MAXVAL;
+
   if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::lucas_sequence) {
     return map { ($_ > ''.~0) ? Math::BigInt->new(''.$_) : $_ }
            Math::Prime::Util::GMP::lucas_sequence($n, $P, $Q, $k);
@@ -651,23 +620,6 @@ sub lucas_sequence {
          Math::Prime::Util::PP::lucas_sequence($n, $P, $Q, $k);
 }
 
-sub miller_rabin_random {
-  my($n, $k, $seed) = @_;
-  _validate_num($n) || _validate_positive_integer($n);
-  _validate_num($k) || _validate_positive_integer($k);
-
-  return 1 if $k <= 0;
-  return (is_prob_prime($n) > 0) if $n < 100;
-  return 0 unless $n & 1;
-
-  return Math::Prime::Util::GMP::miller_rabin_random($n, $k)
-    if $_HAVE_GMP
-    && defined &Math::Prime::Util::GMP::miller_rabin_random;
-
-  require Math::Prime::Util::RandomPrimes;
-  return Math::Prime::Util::RandomPrimes::miller_rabin_random($n, $k, $seed);
-}
-
 
 #############################################################################
 
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 4317ad8..461e767 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -40,6 +40,7 @@ BEGIN {
   use constant BTWO            => Math::BigInt->new(2);
   use constant B_PRIM759       => Math::BigInt->new("64092011671807087969");
   use constant B_PRIM235       => Math::BigInt->new("30");
+  use constant PI_TIMES_8      => 25.13274122871834590770114707;
 }
 
 {
@@ -82,6 +83,29 @@ sub _upgrade_to_float {
   return Math::BigFloat->new($_[0]);
 }
 
+# Get the accuracy of variable x, or the max default from BigInt/BigFloat
+# One might think to use ref($x)->accuracy() but numbers get upgraded and
+# downgraded willy-nilly, and it will do the wrong thing from the user's
+# perspective.
+sub _find_big_acc {
+  my($x) = @_;
+
+  $b = $x->accuracy();
+  return $b if defined $b;
+
+  my ($i,$f) = (Math::BigInt->accuracy(), Math::BigFloat->accuracy());
+  return (($i > $f) ? $i : $f)   if defined $i && defined $f;
+  return $i if defined $i;
+  return $f if defined $f;
+
+  ($i,$f) = (Math::BigInt->div_scale(), Math::BigFloat->div_scale());
+  return (($i > $f) ? $i : $f)   if defined $i && defined $f;
+  return $i if defined $i;
+  return $f if defined $f;
+  return 18;
+}
+
+
 sub _validate_num {
   my($n, $min, $max) = @_;
   croak "Parameter must be defined" if !defined $n;
@@ -556,7 +580,7 @@ sub consecutive_integer_lcm {
 
   my $max = (MPU_32BIT) ? 22 : (OLD_PERL_VERSION) ? 37 : 46;
   my $pn = ref($n) ? ref($n)->new(1) : ($n >= $max) ? Math::BigInt->bone() : 1;
-  foreach my $p (@{primes($n)}) {
+  for (my $p = 2; $p <= $n; $p = next_prime($p)) {
     my($p_power, $pmin) = ($p, int($n/$p));
     $p_power *= $p while $p_power <= $pmin;
     $pn *= $p_power;
@@ -715,6 +739,25 @@ sub mertens {
   return $sum;
 }
 
+sub liouville {
+  my($n) = @_;
+  my $l = (-1) ** scalar factor($n);
+  return $l;
+}
+
+# Exponential of Mangoldt function (A014963).
+# Return p if n = p^m [p prime, m >= 1], 1 otherwise.
+sub exp_mangoldt {
+  my($n) = @_;
+  return 1 if defined $n && $n <= 1;  # n <= 1
+  return 2 if ($n & ($n-1)) == 0;     # n power of 2
+  return 1 unless $n & 1;             # even n can't be p^m
+
+  my @pe = Math::Prime::Util::factor_exp($n);
+  return 1 if scalar @pe > 1;
+  return $pe[0]->[0];
+}
+
 sub carmichael_lambda {
   my($n) = @_;
   return euler_phi($n) if $n < 8;                # = phi(n) for n < 8
@@ -739,19 +782,13 @@ my @_ds_overflow =  # We'll use BigInt math if the input 
is larger than this.
    : ( 50,           845404560,      52560,    1548,   252,   84);
 sub divisor_sum {
   my($n, $k) = @_;
-  return 1 if defined $n && $n == 1;
+  return 1 if $n == 1;
 
   if (defined $k && ref($k) eq 'CODE') {
     my $sum = $n-$n;
-    if (ref($n) eq 'Math::BigInt') {
-      # If the original number was a bigint, make sure all divisors are.
-      foreach my $d (Math::Prime::Util::divisors($n)) {
-        $sum += $k->(Math::BigInt->new("$d"));
-      }
-    } else {
-      foreach my $d (Math::Prime::Util::divisors($n)) {
-        $sum += $k->($d);
-      }
+    my $refn = ref($n);
+    foreach my $d (Math::Prime::Util::divisors($n)) {
+      $sum += $k->( $refn ? $refn->new("$d") : $d );
     }
     return $sum;
   }
@@ -1185,10 +1222,15 @@ sub prime_count_lower {
 
   my $result;
   if ($x > 1000_000_000_000 && 
Math::Prime::Util::prime_get_config()->{'assume_rh'}) {
+    # Schoenfeld bound
     my $lix = LogarithmicIntegral($x);
     my $sqx = sqrt($x);
-    # Schoenfeld bound:    (constant is 8 * Pi)
-    $result = $lix - (($sqx*$flogx) / 25.13274122871834590770114707);
+    if (ref($x) eq 'Math::BigFloat') {
+      my $xdigits = _find_big_acc($x);
+      $result = $lix - (($sqx*$flogx) / (Math::BigFloat->bpi($xdigits)*8));
+    } else {
+      $result = $lix - (($sqx*$flogx) / PI_TIMES_8);
+    }
   } elsif ($x < 599) {
     $result = $x / ($flogx - 0.7);   # For smaller numbers this works out well.
   } else {
@@ -1237,10 +1279,15 @@ sub prime_count_upper {
 
   my $result;
   if ($x > 10000_000_000_000 && 
Math::Prime::Util::prime_get_config()->{'assume_rh'}) {
+    # Schoenfeld bound
     my $lix = LogarithmicIntegral($x);
     my $sqx = sqrt($x);
-    # Schoenfeld bound:    (constant is 8 * Pi)
-    $result = $lix + (($sqx*$flogx) / 25.13274122871834590770114707);
+    if (ref($x) eq 'Math::BigFloat') {
+      my $xdigits = _find_big_acc($x);
+      $result = $lix + (($sqx*$flogx) / (Math::BigFloat->bpi($xdigits)*8));
+    } else {
+      $result = $lix + (($sqx*$flogx) / PI_TIMES_8);
+    }
   } elsif ($x <  1621) { $result = ($x / ($flogx - 1.048)) + 1.0; }
     elsif ($x <  5000) { $result = ($x / ($flogx - 1.071)) + 1.0; }
     elsif ($x < 15900) { $result = ($x / ($flogx - 1.098)) + 1.0; }
@@ -2883,7 +2930,7 @@ sub ecm_factor {
 
 sub divisors {
   my($n) = @_;
-  _validate_num($n) || _validate_positive_integer($n);
+  _validate_positive_integer($n);
 
   # In scalar context, returns sigma_0(n).  Very fast.
   return Math::Prime::Util::divisor_sum($n,0) unless wantarray;
@@ -2918,6 +2965,33 @@ sub divisors {
 }
 
 
+sub chebyshev_theta {
+  my($n) = @_;
+  my $sum = 0.0;
+  for (my $p = 2; $p <= $n; $p = next_prime($p)) {
+    $sum += log($p);
+  }
+  return $sum;
+}
+
+sub chebyshev_psi {
+  my($n) = @_;
+  return 0 if $n <= 1;
+
+  my ($sum, $p, $logn, $sqrtn) = (0.0, 2, log($n), int(sqrt($n)));
+
+  for ( ; $p <= $sqrtn; $p = next_prime($p)) {
+    my $logp = log($p);
+    $sum += $logp * int($logn/$logp+1e-15);
+  }
+
+  for ( ; $p <= $n; $p = next_prime($p)) {
+    $sum += log($p);
+  }
+
+  return $sum;
+}
+
 sub ExponentialIntegral {
   my($x) = @_;
   return - MPU_INFINITY if $x == 0;
@@ -2932,8 +3006,8 @@ sub ExponentialIntegral {
       do { require Math::BigFloat; Math::BigFloat->import(); }
         if !defined $Math::BigFloat::VERSION;
       $x = Math::BigFloat->new("$x") if ref($x) ne 'Math::BigFloat';
-      $wantbf = 1;
-      $xdigits = $x->accuracy || Math::BigFloat->accuracy() || 
Math::BigFloat->div_scale();
+      $wantbf = _find_big_acc($x);
+      $xdigits = $wantbf;
     }
     my $rnd = 0;  # MPFR_RNDN;
     my $bit_precision = int($xdigits * 3.322) + 4;
@@ -2944,7 +3018,7 @@ sub ExponentialIntegral {
     Math::MPFR::Rmpfr_set_prec($eix, $bit_precision);
     Math::MPFR::Rmpfr_eint($eix, $rx, $rnd);
     my $strval = Math::MPFR::Rmpfr_get_str($eix, 10, 0, $rnd);
-    return ($wantbf)  ?  Math::BigFloat->new($strval)  :  0.0 + $strval;
+    return ($wantbf)  ?  Math::BigFloat->new($strval,$wantbf)  :  0.0 + 
$strval;
   }
 
   $x = Math::BigFloat->new("$x") if defined $bignum::VERSION && ref($x) ne 
'Math::BigFloat';
@@ -3031,8 +3105,8 @@ sub LogarithmicIntegral {
     my $wantbf = 0;
     my $xdigits = 18;
     if (defined $bignum::VERSION || ref($x) =~ /^Math::Big/) {
-      $wantbf = 1;
-      $xdigits = $x->accuracy || Math::BigInt->accuracy() || 
Math::BigInt->div_scale();
+      $wantbf = _find_big_acc($x);
+      $xdigits = $wantbf;
     }
     $xdigits += length(int(log(0.0+"$x"))) + 1;
     my $rnd = 0;  # MPFR_RNDN;
@@ -3045,9 +3119,7 @@ sub LogarithmicIntegral {
     Math::MPFR::Rmpfr_set_prec($lix, $bit_precision);
     Math::MPFR::Rmpfr_eint($lix, $rx, $rnd);
     my $strval = Math::MPFR::Rmpfr_get_str($lix, 10, 0, $rnd);
-    return Math::BigFloat->new($strval, ($x->accuracy || 
Math::BigInt->accuracy() || Math::BigInt->div_scale()))
-      if $wantbf;
-    return 0.0 + $strval;
+    return ($wantbf)  ?  Math::BigFloat->new($strval,$wantbf)  :  0.0 + 
$strval;
   }
 
   if ($x == 2) {
@@ -3066,7 +3138,7 @@ sub LogarithmicIntegral {
   my $xdigits = 0;
   my $finalacc = 0;
   if (ref($x) =~ /^Math::Big/) {
-    $xdigits = $x->accuracy || Math::BigFloat->accuracy() || 
Math::BigFloat->div_scale();
+    $xdigits = _find_big_acc($x);
     my $xlen = length($x->bfloor->bstr());
     $xdigits = $xlen if $xdigits < $xlen;
     $finalacc = $xdigits;
@@ -3179,8 +3251,8 @@ sub RiemannZeta {
         $x->accuracy($xacc) if $xacc;
       }
       $x = Math::BigFloat->new("$x") if ref($x) ne 'Math::BigFloat';
-      $wantbf = 1;
-      $xdigits = $x->accuracy || Math::BigFloat->accuracy() || 
Math::BigFloat->div_scale();
+      $wantbf = _find_big_acc($x);
+      $xdigits = $wantbf;
     }
     my $rnd = 0;  # MPFR_RNDN;
     my $bit_precision = int($xdigits * 3.322) + 7;
@@ -3194,7 +3266,7 @@ sub RiemannZeta {
     Math::MPFR::Rmpfr_zeta($zetax, $rx, $rnd);
     Math::MPFR::Rmpfr_sub_ui($zetax, $zetax, 1, $rnd);
     my $strval = Math::MPFR::Rmpfr_get_str($zetax, 10, $xdigits, $rnd);
-    return ($wantbf)  ?  Math::BigFloat->new($strval)  :  0.0 + $strval;
+    return ($wantbf)  ?  Math::BigFloat->new($strval,$wantbf)  :  0.0 + 
$strval;
   }
 
   if (defined $bignum::VERSION || ref($x) =~ /^Math::Big/) {
@@ -3268,8 +3340,8 @@ sub RiemannR {
         $x->accuracy($xacc) if $xacc;
       }
       $x = Math::BigFloat->new("$x") if ref($x) ne 'Math::BigFloat';
-      $wantbf = 1;
-      $xdigits = $x->accuracy || Math::BigFloat->accuracy() || 
Math::BigFloat->div_scale();
+      $wantbf = _find_big_acc($x);
+      $xdigits = $wantbf;
     }
     my $rnd = 0;  # MPFR_RNDN;
     my $bit_precision = int($xdigits * 3.322) + 8;  # Add some extra
@@ -3310,7 +3382,7 @@ sub RiemannR {
       Math::MPFR::Rmpfr_add($rsum, $rsum, $rterm, $rnd);
     }
     my $strval = Math::MPFR::Rmpfr_get_str($rsum, 10, $xdigits, $rnd);
-    return ($wantbf)  ?  Math::BigFloat->new($strval)  :  0.0 + $strval;
+    return ($wantbf)  ?  Math::BigFloat->new($strval,$wantbf)  :  0.0 + 
$strval;
   }
 
   if (defined $bignum::VERSION || ref($x) =~ /^Math::Big/) {
diff --git a/lib/Math/Prime/Util/PPFE.pm b/lib/Math/Prime/Util/PPFE.pm
index 98ed907..5787af5 100644
--- a/lib/Math/Prime/Util/PPFE.pm
+++ b/lib/Math/Prime/Util/PPFE.pm
@@ -1,9 +1,13 @@
-# The PP front end, only loaded if XS is not used.  It is intended to load
-# directly into the package namespace.
-
+package Math::Prime::Util::PPFE;
 use strict;
 use warnings;
 use Math::Prime::Util::PP;
+use Carp qw/carp croak confess/;
+
+# The PP front end, only loaded if XS is not used.
+# It is intended to load directly into the MPU namespace.
+
+package Math::Prime::Util;
 
 *_validate_num = \&Math::Prime::Util::PP::_validate_num;
 *_prime_memfreeall = \&Math::Prime::Util::PP::_prime_memfreeall;
@@ -11,12 +15,6 @@ use Math::Prime::Util::PP;
 *prime_precalc  = \&Math::Prime::Util::PP::prime_precalc;
 
 
-sub mertens {
-  my($n) = @_;
-  _validate_positive_integer($n);
-  return Math::Prime::Util::PP::mertens(@_);
-}
-
 sub moebius {
   if (scalar @_ <= 1) {
     my($n) = @_;
@@ -47,18 +45,35 @@ sub jordan_totient {
   _validate_positive_integer($k);
   return 0 if defined $n && $n < 0;
   _validate_positive_integer($n);
-  return Math::Prime::Util::PP::jordan_totient(@_);
+  return Math::Prime::Util::PP::jordan_totient($k, $n);
 }
 sub carmichael_lambda {
   my($n) = @_;
   _validate_positive_integer($n);
-  return Math::Prime::Util::PP::carmichael_lambda(@_);
+  return Math::Prime::Util::PP::carmichael_lambda($n);
+}
+sub mertens {
+  my($n) = @_;
+  _validate_positive_integer($n);
+  return Math::Prime::Util::PP::mertens($n);
+}
+sub liouville {
+  my($n) = @_;
+  _validate_positive_integer($n);
+  return Math::Prime::Util::PP::liouville($n);
+}
+sub exp_mangoldt {
+  my($n) = @_;
+  return 1 if defined $n && $n <= 1;
+  _validate_positive_integer($n);
+  return Math::Prime::Util::PP::exp_mangoldt($n);
 }
 
+
 sub nth_prime {
   my($n) = @_;
   _validate_positive_integer($n);
-  return Math::Prime::Util::PP::nth_prime(@_);
+  return Math::Prime::Util::PP::nth_prime($n);
 }
 sub nth_prime_lower {
   my($n) = @_;
@@ -90,7 +105,7 @@ sub prime_count_approx {
   _validate_positive_integer($n);
   return Math::Prime::Util::PP::prime_count_approx($n);
 }
-  
+
 
 sub is_prime {
   my($n) = @_;
@@ -261,6 +276,19 @@ sub ecm_factor {
   Math::Prime::Util::PP::ecm_factor($n, $B1, $B2, $ncurves);
 }
 
+sub divisors {
+  my($n) = @_;
+  _validate_positive_integer($n);
+  return Math::Prime::Util::PP::divisors($n);
+}
+
+sub divisor_sum {
+  my($n, $k) = @_;
+  _validate_positive_integer($n);
+  _validate_positive_integer($k) if defined $k && ref($k) ne 'CODE';
+  return Math::Prime::Util::PP::divisor_sum($n, $k);
+}
+
 sub gcd {
   return Math::Prime::Util::PP::gcd(@_);
 }
@@ -275,4 +303,103 @@ sub legendre_phi {
   return Math::Prime::Util::PP::legendre_phi($x, $a);
 }
 
+sub chebyshev_theta {
+  my($n) = @_;
+  _validate_positive_integer($n);
+  return Math::Prime::Util::PP::chebyshev_theta($n);
+}
+sub chebyshev_psi {
+  my($n) = @_;
+  _validate_positive_integer($n);
+  return Math::Prime::Util::PP::chebyshev_psi($n);
+}
+
+#############################################################################
+
+sub forprimes (&$;$) {    ## no critic qw(ProhibitSubroutinePrototypes)
+  my($sub, $beg, $end) = @_;
+  if (!defined $end) { $end = $beg; $beg = 2; }
+  _validate_num($beg) || _validate_positive_integer($beg);
+  _validate_num($end) || _validate_positive_integer($end);
+  $beg = 2 if $beg < 2;
+  {
+    my $pp;
+    local *_ = \$pp;
+    for (my $p = next_prime($beg-1);  $p <= $end;  $p = next_prime($p)) {
+      $pp = $p;
+      $sub->();
+    }
+  }
+}
+
+sub forcomposites(&$;$) { ## no critic qw(ProhibitSubroutinePrototypes)
+  my($sub, $beg, $end) = @_;
+  if (!defined $end) { $end = $beg; $beg = 4; }
+  _validate_num($beg) || _validate_positive_integer($beg);
+  _validate_num($end) || _validate_positive_integer($end);
+  $beg = 4 if $beg < 4;
+  $end = Math::BigInt->new(''.~0) if ref($end) ne 'Math::BigInt' && $end == ~0;
+  {
+    my $pp;
+    local *_ = \$pp;
+    for ( ; $beg <= $end ; $beg++ ) {
+      if (!is_prime($beg)) {
+        $pp = $beg;
+        $sub->();
+      }
+    }
+  }
+}
+
+sub fordivisors (&$) {    ## no critic qw(ProhibitSubroutinePrototypes)
+  my($sub, $n) = @_;
+  _validate_num($n) || _validate_positive_integer($n);
+  my @divisors = divisors($n);
+  {
+    my $pp;
+    local *_ = \$pp;
+    foreach my $d (@divisors) {
+      $pp = $d;
+      $sub->();
+    }
+  }
+}
+
 1;
+
+__END__
+
+=pod
+
+=head1 NAME
+
+Math::Prime::Util::PPFE - PP front end for Math::Prime::Util
+
+=head1 SYNOPSIS
+
+This loads the PP code and adds input validation front ends.  It is only
+meant to be used when XS is not used.
+
+=head1 DESCRIPTION
+
+Loads PP module and implements PP front-end functions for all XS code.
+This is used only if the XS code is not loaded.
+
+=head1 SEE ALSO
+
+L<Math::Prime::Util>
+
+L<Math::Prime::Util::PP>
+
+=head1 AUTHORS
+
+Dana Jacobsen E<lt>d...@acm.orge<gt>
+
+
+=head1 COPYRIGHT
+
+Copyright 2014 by Dana Jacobsen E<lt>d...@acm.orge<gt>
+
+This program is free software; you can redistribute it and/or modify it under 
the same terms as Perl itself.
+
+=cut
diff --git a/t/80-pp.t b/t/80-pp.t
index 627fdae..b6033ec 100644
--- a/t/80-pp.t
+++ b/t/80-pp.t
@@ -301,7 +301,7 @@ require_ok 'Math::Prime::Util::PrimalityProving';
     *moebius        = \&Math::Prime::Util::PP::moebius;
     *euler_phi      = \&Math::Prime::Util::PP::euler_phi;
     *mertens        = \&Math::Prime::Util::PP::mertens;
-    *exp_mangoldt   = \&Math::Prime::Util::_generic_exp_mangoldt;
+    *exp_mangoldt   = \&Math::Prime::Util::PP::exp_mangoldt;
 
     *RiemannR            = \&Math::Prime::Util::PP::RiemannR;
     *RiemannZeta         = \&Math::Prime::Util::PP::RiemannZeta;
diff --git a/t/92-release-pod-coverage.t b/t/92-release-pod-coverage.t
index 1eff345..ba0a023 100644
--- a/t/92-release-pod-coverage.t
+++ b/t/92-release-pod-coverage.t
@@ -17,7 +17,9 @@ eval "use Test::Pod::Coverage 1.08";
 plan skip_all => "Test::Pod::Coverage 1.08 required for testing POD coverage"
   if $@;
 
-my @modules = Test::Pod::Coverage::all_modules();
+my @modules = grep { $_ ne 'Math::Prime::Util::PPFE' }
+              Test::Pod::Coverage::all_modules();
+
 plan tests => scalar @modules;
 
 #my $ppsubclass = { trustme => [mpu_public_regex()] };
diff --git a/util.c b/util.c
index e8e8e31..796a8c6 100644
--- a/util.c
+++ b/util.c
@@ -222,7 +222,7 @@ int _XS_is_prime(UV n)
 
 UV next_prime(UV n)
 {
-  UV d, m, sieve_size, next;
+  UV m, sieve_size, next;
   const unsigned char* sieve;
 
   if (n < 30*NPRIME_SIEVE30) {
@@ -249,7 +249,7 @@ UV next_prime(UV n)
 UV prev_prime(UV n)
 {
   const unsigned char* sieve;
-  UV d, m, prev;
+  UV m, prev;
 
   if (n < 30*NPRIME_SIEVE30)
     return prev_prime_in_sieve(prime_sieve30, n);
@@ -646,12 +646,12 @@ UV prime_count_upper(UV n)
   fn     = (long double) n;
   flogn  = logl(n);
 
-  for (i = 0; i < NUPPER_THRESH; i++)
+  for (i = 0; i < (int)NUPPER_THRESH; i++)
     if (n < _upper_thresh[i].thresh)
       break;
 
-  if (i < NUPPER_THRESH)   a = _upper_thresh[i].aval;
-  else                     a = 2.334;   /* Dusart 2010, page 2 */
+  if (i < (int)NUPPER_THRESH) a = _upper_thresh[i].aval;
+  else                        a = 2.334;   /* Dusart 2010, page 2 */
 
   upper = fn/flogn * (1.0 + 1.0/flogn + a/(flogn*flogn));
   return (UV) ceill(upper);
@@ -1516,12 +1516,8 @@ long double ld_riemann_zeta(long double x) {
     return sum;
   }
 
-  if (x > 2000.0) {
-    /* 1) zeta(2000)-1 is about 8.7E-603, which is far less than a IEEE-754
-     *    64-bit double can represent.  A 128-bit quad could go to ~16000.
-     * 2) pow / powl start getting obnoxiously slow with values like -7500. */
+  if (x > 17000.0)
     return 0.0;
-  }
 
 #if 0
   {
@@ -1564,8 +1560,8 @@ long double ld_riemann_zeta(long double x) {
     for (i = 2; i < 11; i++) {
       b = powl( i, -x );
       s += b;
-      if (fabsl(b/s) < LDBL_EPSILON)
-         return s;
+      if (fabsl(b) < fabsl(LDBL_EPSILON * s))
+        return s;
     }
     s = s + b*w/(x-1.0) - 0.5 * b;
     a = 1.0;
@@ -1575,8 +1571,7 @@ long double ld_riemann_zeta(long double x) {
       b /= w;
       t = a*b/A[i];
       s = s + t;
-      t = fabsl(t/s);
-      if (t < LDBL_EPSILON)
+      if (fabsl(t) < fabsl(LDBL_EPSILON * s))
         break;
       a *= x + k + 1.0;
       b /= w;
@@ -1601,8 +1596,8 @@ long double _XS_RiemannR(long double x) {
     part_term *= flogx / k;
     term = part_term / (k + k * ld_riemann_zeta(k+1));
     KAHAN_SUM(sum, term);
-    /* printf("R  after adding %.15lg, sum = %.15lg\n", term, sum); */
-    if (fabsl(term/sum) < LDBL_EPSILON) break;
+    /* printf("R %5d after adding %.18Lg, sum = %.19Lg\n", k, term, sum); */
+    if (fabsl(term) < fabsl(LDBL_EPSILON*sum)) break;
   }
 
   return 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