[libmath-prime-util-perl] 18/33: More moving functions out of Util
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:51:42 UTC 2015
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 <dana at 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>dana at acm.orgE<gt>
+
+
+=head1 COPYRIGHT
+
+Copyright 2014 by Dana Jacobsen E<lt>dana at 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
More information about the Pkg-perl-cvs-commits
mailing list