[libmath-prime-util-perl] 17/59: Many changes for bignum support
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:54 UTC 2015
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 3e23f1f0bb4510b3e2b2d38e98e466310031e568
Author: Dana Jacobsen <dana at acm.org>
Date: Sun Jul 1 02:50:01 2012 -0600
Many changes for bignum support
---
Changes | 4 +
TODO | 8 +-
XS.xs | 28 +--
lib/Math/Prime/Util.pm | 466 +++++++++++++++++++++++++++++++++++++++++-----
lib/Math/Prime/Util/PP.pm | 328 +++++++++++---------------------
t/03-init.t | 36 ++--
t/10-isprime.t | 2 +-
t/11-primes.t | 2 +-
t/12-nextprime.t | 2 +-
t/13-primecount.t | 2 +-
t/14-nthprime.t | 2 +-
t/16-randomprime.t | 2 +-
t/17-pseudoprime.t | 2 +-
t/18-functions.t | 2 +-
t/30-relations.t | 2 +-
t/50-factoring.t | 2 +-
t/80-pp.t | 43 ++---
util.c | 26 +--
util.h | 21 ++-
19 files changed, 607 insertions(+), 373 deletions(-)
diff --git a/Changes b/Changes
index b13078c..def2989 100644
--- a/Changes
+++ b/Changes
@@ -4,6 +4,10 @@ Revision history for Perl extension Math::Prime::Util.
- Make miller_rabin return 0 if given even number.
- Add tests for PP.
- Cleanup of some tests.
+ - The XS miller_rabin code now works with large base > n.
+ - Moved prime_count_* and nth_prime_* into Util.pm. This lets them
+ work with big numbers.
+ - factor and all_factor should correctly work with bigints.
0.09 25 June 2012
- Pure Perl code added. Passes all tests. Used only if the XSLoader
diff --git a/TODO b/TODO
index f04f361..265ad6f 100644
--- a/TODO
+++ b/TODO
@@ -25,12 +25,8 @@
- Move .c / .h files into separate directory.
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.
+
+- redo _XS_factor to return a sorted array
diff --git a/XS.xs b/XS.xs
index 842e2b3..b828e94 100644
--- a/XS.xs
+++ b/XS.xs
@@ -40,26 +40,8 @@ prime_count(IN UV low, IN UV high = 0)
RETVAL
UV
-prime_count_lower(IN UV n)
-
-UV
-prime_count_upper(IN UV n)
-
-UV
-prime_count_approx(IN UV n)
-
-UV
nth_prime(IN UV n)
-UV
-nth_prime_lower(IN UV n)
-
-UV
-nth_prime_upper(IN UV n)
-
-UV
-nth_prime_approx(IN UV n)
-
int
is_prime(IN UV n)
@@ -78,7 +60,7 @@ _get_prime_cache_size()
RETVAL
int
-_maxbits()
+_XS_prime_maxbits()
CODE:
RETVAL = BITS_PER_WORD;
OUTPUT:
@@ -212,7 +194,7 @@ erat_primes(IN UV low, IN UV high)
void
-XS_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 +353,10 @@ int
is_prob_prime(IN UV n)
double
-XS_ExponentialIntegral(double x)
+_XS_ExponentialIntegral(double x)
double
-XS_LogarithmicIntegral(double x)
+_XS_LogarithmicIntegral(double x)
double
-XS_RiemannR(double x)
+_XS_RiemannR(double x)
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index f7511d6..a70bc6c 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -31,43 +31,37 @@ BEGIN {
# Load PP code. Nothing exported.
require Math::Prime::Util::PP;
+ # There is no GMP module yet
+ $_Config{'gmp'} = 0;
eval {
require XSLoader;
XSLoader::load(__PACKAGE__, $Math::Prime::Util::VERSION);
prime_precalc(0);
$_Config{'xs'} = 1;
- $_Config{'maxbits'} = _maxbits();
+ $_Config{'maxbits'} = _XS_prime_maxbits();
1;
} or do {
-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;
+ $_Config{'maxbits'} = Math::Prime::Util::PP::_PP_prime_maxbits();
+ carp "Using Pure Perl implementation: $@";
*_prime_memfreeall = \&Math::Prime::Util::PP::_prime_memfreeall;
*prime_memfree = \&Math::Prime::Util::PP::prime_memfree;
*prime_precalc = \&Math::Prime::Util::PP::prime_precalc;
*prime_count = \&Math::Prime::Util::PP::prime_count;
- *prime_count_upper = \&Math::Prime::Util::PP::prime_count_upper;
- *prime_count_lower = \&Math::Prime::Util::PP::prime_count_lower;
- *prime_count_approx = \&Math::Prime::Util::PP::prime_count_approx;
*nth_prime = \&Math::Prime::Util::PP::nth_prime;
- *nth_prime_upper = \&Math::Prime::Util::PP::nth_prime_upper;
- *nth_prime_lower = \&Math::Prime::Util::PP::nth_prime_lower;
- *nth_prime_approx = \&Math::Prime::Util::PP::nth_prime_approx;
*is_prime = \&Math::Prime::Util::PP::is_prime;
*next_prime = \&Math::Prime::Util::PP::next_prime;
*prev_prime = \&Math::Prime::Util::PP::prev_prime;
+ # Perhaps these can be moved here?
*miller_rabin = \&Math::Prime::Util::PP::miller_rabin;
*is_prob_prime = \&Math::Prime::Util::PP::is_prob_prime;
- *factor = \&Math::Prime::Util::PP::factor;
+ # These probably shouldn't even be exported
*trial_factor = \&Math::Prime::Util::PP::trial_factor;
*fermat_factor = \&Math::Prime::Util::PP::fermat_factor;
*holf_factor = \&Math::Prime::Util::PP::holf_factor;
@@ -75,10 +69,6 @@ die $@;
*pbrent_factor = \&Math::Prime::Util::PP::pbrent_factor;
*prho_factor = \&Math::Prime::Util::PP::prho_factor;
*pminus1_factor = \&Math::Prime::Util::PP::pminus1_factor;
-
- *RiemannR = \&Math::Prime::Util::PP::RiemannR;
- *LogarithmicIntegral = \&Math::Prime::Util::PP::LogarithmicIntegral;
- *ExponentialIntegral = \&Math::Prime::Util::PP::ExponentialIntegral;
}
}
END {
@@ -99,7 +89,13 @@ if ($_Config{'maxbits'} == 32) {
sub prime_get_config {
my %config = %_Config;
+
+ $config{'precalc_to'} = ($_Config{'xs'})
+ ? _get_prime_cache_size
+ : Math::Prime::Util::PP::_get_prime_cache_size;
+
return \%config;
+
}
sub _validate_positive_integer {
@@ -113,6 +109,27 @@ sub _validate_positive_integer {
1;
}
+my @_primes_small = (
+ 0,2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,
+ 101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,
+ 193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,
+ 293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,
+ 409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499);
+my @_prime_count_small = (
+ 0,0,1,2,2,3,3,4,4,4,4,5,5,6,6,6,6,7,7,8,8,8,8,9,9,9,9,9,9,10,10,
+ 11,11,11,11,11,11,12,12,12,12,13,13,14,14,14,14,15,15,15,15,15,15,
+ 16,16,16,16,16,16,17,17,18,18,18,18,18,18,19);
+#my @_prime_next_small = (
+# 2,2,3,5,5,7,7,11,11,11,11,13,13,17,17,17,17,19,19,23,23,23,23,
+# 29,29,29,29,29,29,31,31,37,37,37,37,37,37,41,41,41,41,43,43,47,
+# 47,47,47,53,53,53,53,53,53,59,59,59,59,59,59,61,61,67,67,67,67,67,67,71);
+
+
+
+
+
+#############################################################################
+
sub primes {
my $optref = (ref $_[0] eq 'HASH') ? shift : {};
croak "no parameters to primes" unless scalar @_ > 0;
@@ -247,14 +264,24 @@ sub random_ndigit_prime {
sub all_factors {
my $n = shift;
+ my $use_bigint = (ref($n) =~ /^Math::Big/);
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;
+ if (!$use_bigint) {
+ foreach my $f2 (@all) {
+ $all_factors{$f1*$f2} = 1 if ($f1*$f2) < $n;
+ }
+ } else {
+ # Many of the factors will be numified after coming back, so we need
+ # to make sure we're using bigints when we calculate the product.
+ foreach my $f2 (@all) {
+ my $product = Math::BigInt->new("$f1") * Math::BigInt->new("$f2");
+ $all_factors{$product} = 1 if $product < $n;
+ }
}
$all_factors{$f1} = 1;
}
@@ -262,6 +289,7 @@ sub all_factors {
return @factors;
}
+
# A008683 Moebius function mu(n)
# A030059, A013929, A030229, A002321, A005117, A013929 all relate.
@@ -271,7 +299,6 @@ sub all_factors {
sub moebius {
my($n) = @_;
_validate_positive_integer($n, 1);
- return if $n <= 0;
return 1 if $n == 1;
# Quick check for small replicated factors
@@ -285,18 +312,35 @@ sub moebius {
return (((scalar @factors) % 2) == 0) ? 1 : -1;
}
+
+# Euler Phi, aka Euler Totient. A000010
+
sub euler_phi {
my($n) = @_;
+ # SAGE defines this to be 0 for all n <= 0. Others choose differently.
+ return 0 if defined $n && $n <= 0; # Following SAGE's logic here.
_validate_positive_integer($n);
- return 1 if $n <= 1; # phi(0) is disputed
+ return 1 if $n <= 1;
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);
}
+
+ # Alternate way if you want to avoid divisions.
+ #my $totient = 1;
+ #foreach my $factor (@factors) {
+ # $totient *= ($factor - 1);
+ # while ($factor_mult{$factor} > 1) {
+ # $totient *= $factor;
+ # $factor_mult{$factor}--;
+ # }
+ #}
+
$totient;
}
@@ -308,10 +352,20 @@ sub euler_phi {
# based on the input (XS, GMP, PP).
#############################################################################
+# Doing a sub here like:
#
-# 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 foo { my($n) = @_; _validate_positive_integer($n);
+# return _XS_... if $_Config{'xs'} && $n <= $_Config{'maxparam'}; }
#
+# takes about 0.7uS on my machine. Operations like is_prime and factor run
+# on small input (under 100_000) typically take a lot less time than this. So
+# The overhead for these is significantly more than just the XS call itself.
+#
+# The plan for some of these functions will be to invert the operation. That
+# is, the XS functions will look at the input and make a call here if the input
+# is large.
+
+
#sub is_prime {
# my($n) = @_;
# _validate_positive_integer($n);
@@ -323,27 +377,267 @@ sub euler_phi {
sub factor {
my($n) = @_;
_validate_positive_integer($n);
-
- return XS_factor($n) if $_Config{'xs'} && $n <= $_Config{'maxparam'};
+
+ #return _XS_factor($n) if $_Config{'xs'} && $n <= $_Config{'maxparam'};
+ if ($_Config{'xs'} && $n <= $_Config{'maxparam'}) {
+ my @factors = sort {$a<=>$b} _XS_factor($n);
+ return @factors;
+ }
$n = $n->as_number() if ref($n) =~ /^Math::BigFloat/;
$n = $n->numify if $n <= ~0;
Math::Prime::Util::PP::factor($n);
}
+#############################################################################
+
+sub prime_count_approx {
+ my($x) = @_;
+ _validate_positive_integer($x);
+
+ return $_prime_count_small[$x] if $x <= $#_prime_count_small;
+
+ # 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%
+
+ # return int( (prime_count_upper($x) + prime_count_lower($x)) / 2);
+
+ # return int( LogarithmicIntegral($x) );
+
+ # return int( LogarithmicIntegral($x) - LogarithmicIntegral(sqrt($x))/2 );
+
+ return int(RiemannR($x)+0.5);
+}
+
+sub prime_count_lower {
+ my($x) = @_;
+ _validate_positive_integer($x);
+
+ return $_prime_count_small[$x] if $x <= $#_prime_count_small;
+
+ if (ref($x) =~ /^Math::BigInt/ && !defined $Math::BigFloat::VERSION) {
+ require Math::BigFloat; Math::BigFloat->import;
+ $x = new Math::BigFloat "$x";
+ }
+
+ my $flogx = log($x);
+
+ # Chebyshev: 1*x/logx x >= 17
+ # Rosser & Schoenfeld: x/(logx-1/2) x >= 67
+ # Dusart 1999: x/logx*(1+1/logx+1.8/logxlogx) x >= 32299
+
+ # For smaller numbers this works out well.
+ return int( $x / ($flogx - 0.7) ) if $x < 599;
+
+ my $a;
+ # Hand tuned for small numbers (< 60_000M)
+ if ($x < 2700) { $a = 0.30; }
+ elsif ($x < 5500) { $a = 0.90; }
+ elsif ($x < 19400) { $a = 1.30; }
+ elsif ($x < 32299) { $a = 1.60; }
+ elsif ($x < 176000) { $a = 1.80; }
+ elsif ($x < 315000) { $a = 2.10; }
+ elsif ($x < 1100000) { $a = 2.20; }
+ elsif ($x < 4500000) { $a = 2.31; }
+ elsif ($x < 233000000) { $a = 2.36; }
+ elsif ($x < 5433800000) { $a = 2.32; }
+ elsif ($x <60000000000) { $a = 2.15; }
+ else { $a = 1.80; } # Dusart 1999, page 14
+
+ return int( ($x/$flogx) * (1.0 + 1.0/$flogx + $a/($flogx*$flogx)) );
+}
+
+sub prime_count_upper {
+ my($x) = @_;
+ _validate_positive_integer($x);
+
+ return $_prime_count_small[$x] if $x <= $#_prime_count_small;
+
+ if (ref($x) =~ /^Math::BigInt/ && !defined $Math::BigFloat::VERSION) {
+ require Math::BigFloat;
+ $x = new Math::BigFloat "$x";
+ }
+
+ # Chebyshev: 1.25506*x/logx x >= 17
+ # Rosser & Schoenfeld: x/(logx-3/2) x >= 67
+ # Dusart 1999: x/logx*(1+1/logx+2.51/logxlogx) x >= 355991
+
+ my $flogx = log($x);
+
+ # These work out well for small values
+ return int( ($x / ($flogx - 1.048)) + 1.0 ) if $x < 1621;
+ return int( ($x / ($flogx - 1.071)) + 1.0 ) if $x < 5000;
+ return int( ($x / ($flogx - 1.098)) + 1.0 ) if $x < 15900;
+
+ my $a;
+ # Hand tuned for small numbers (< 60_000M)
+ if ($x < 24000) { $a = 2.30; }
+ elsif ($x < 59000) { $a = 2.48; }
+ elsif ($x < 350000) { $a = 2.52; }
+ elsif ($x < 355991) { $a = 2.54; }
+ elsif ($x < 356000) { $a = 2.51; }
+ elsif ($x < 3550000) { $a = 2.50; }
+ elsif ($x < 3560000) { $a = 2.49; }
+ elsif ($x < 5000000) { $a = 2.48; }
+ elsif ($x < 8000000) { $a = 2.47; }
+ elsif ($x < 13000000) { $a = 2.46; }
+ elsif ($x < 18000000) { $a = 2.45; }
+ elsif ($x < 31000000) { $a = 2.44; }
+ elsif ($x < 41000000) { $a = 2.43; }
+ elsif ($x < 48000000) { $a = 2.42; }
+ elsif ($x < 119000000) { $a = 2.41; }
+ elsif ($x < 182000000) { $a = 2.40; }
+ elsif ($x < 192000000) { $a = 2.395; }
+ elsif ($x < 213000000) { $a = 2.390; }
+ elsif ($x < 271000000) { $a = 2.385; }
+ elsif ($x < 322000000) { $a = 2.380; }
+ elsif ($x < 400000000) { $a = 2.375; }
+ elsif ($x < 510000000) { $a = 2.370; }
+ elsif ($x < 682000000) { $a = 2.367; }
+ elsif ($x <60000000000) { $a = 2.362; }
+ else { $a = 2.51; }
+
+ return int( ($x/$flogx) * (1.0 + 1.0/$flogx + $a/($flogx*$flogx)) + 1.0 );
+}
+
+#############################################################################
+
+sub nth_prime_approx {
+ my($n) = @_;
+ _validate_positive_integer($n);
+
+ return $_primes_small[$n] if $n <= $#_primes_small;
+
+ if (ref($n) =~ /^Math::BigInt/ && !defined $Math::BigFloat::VERSION) {
+ require Math::BigFloat;
+ $n = new Math::BigFloat "$n";
+ }
+
+ my $flogn = log($n);
+ my $flog2n = log($flogn);
+
+ # Cipolla 1902:
+ # m=0 fn * ( flogn + flog2n - 1 );
+ # m=1 + ((flog2n - 2)/flogn) );
+ # m=2 - (((flog2n*flog2n) - 6*flog2n + 11) / (2*flogn*flogn))
+ # + O((flog2n/flogn)^3)
+ #
+ # Shown in Dusart 1999 page 12, as well as other sources such as:
+ # http://www.emis.de/journals/JIPAM/images/153_02_JIPAM/153_02.pdf
+ # where the main issue you run into is that you're doing polynomial
+ # interpolation, so it oscillates like crazy with many high-order terms.
+ # Hence I'm leaving it at m=2.
+ #
+
+ my $approx = $n * ( $flogn + $flog2n - 1
+ + (($flog2n - 2)/$flogn)
+ - ((($flog2n*$flog2n) - 6*$flog2n + 11) / (2*$flogn*$flogn))
+ );
+
+ # Apply a correction to help keep values close.
+ my $order = $flog2n/$flogn;
+ $order = $order*$order*$order * $n;
+
+ if ($n < 259) { $approx += 10.4 * $order; }
+ elsif ($n < 775) { $approx += 7.52* $order; }
+ elsif ($n < 1271) { $approx += 5.6 * $order; }
+ elsif ($n < 2000) { $approx += 5.2 * $order; }
+ elsif ($n < 4000) { $approx += 4.3 * $order; }
+ elsif ($n < 12000) { $approx += 3.0 * $order; }
+ elsif ($n < 150000) { $approx += 2.1 * $order; }
+ elsif ($n < 200000000) { $approx += 0.0 * $order; }
+ else { $approx += -0.010 * $order; }
+
+ if ( ($approx >= ~0) && (ref($n) !~ /^Math::Big/) ) {
+ return $_Config{'maxprime'} if $n <= $_Config{'maxprimeidx'};
+ croak "nth_prime_approx($n) overflow";
+ }
+
+ return int($approx + 0.5);
+}
+
+# The nth prime will be greater than or equal to this number
+sub nth_prime_lower {
+ my($n) = @_;
+ _validate_positive_integer($n);
+
+ return $_primes_small[$n] if $n <= $#_primes_small;
+
+ if (ref($n) =~ /^Math::BigInt/ && !defined $Math::BigFloat::VERSION) {
+ require Math::BigFloat;
+ $n = new Math::BigFloat "$n";
+ }
+
+ my $flogn = log($n);
+ my $flog2n = log($flogn); # Note distinction between log_2(n) and log^2(n)
+
+ # Dusart 1999 page 14, for all n >= 2
+ my $lower = $n * ($flogn + $flog2n - 1.0 + (($flog2n-2.25)/$flogn));
+
+ if ( ($lower >= ~0) && (ref($n) !~ /^Math::Big/) ) {
+ return $_Config{'maxprime'} if $n <= $_Config{'maxprimeidx'};
+ croak "nth_prime_lower($n) overflow";
+ }
+
+ return int($lower);
+}
+
+# The nth prime will be less or equal to this number
+sub nth_prime_upper {
+ my($n) = @_;
+ _validate_positive_integer($n);
+
+ return $_primes_small[$n] if $n <= $#_primes_small;
+
+ if (ref($n) =~ /^Math::BigInt/ && !defined $Math::BigFloat::VERSION) {
+ require Math::BigFloat;
+ $n = new Math::BigFloat "$n";
+ }
+
+ my $flogn = log($n);
+ my $flog2n = log($flogn); # Note distinction between log_2(n) and log^2(n)
+
+ my $upper;
+ if ($n >= 39017) { # Dusart 1999 page 14
+ $upper = $n * ( $flogn + $flog2n - 0.9484 );
+ } elsif ($n >= 27076) { # Dusart 1999 page 14
+ $upper = $n * ( $flogn + $flog2n - 1.0 + (($flog2n-1.80)/$flogn) );
+ } elsif ($n >= 7022) { # Robin 1983
+ $upper = $n * ( $flogn + 0.9385 * $flog2n );
+ } else {
+ $upper = $n * ( $flogn + $flog2n );
+ }
+
+ if ( ($upper >= ~0) && (ref($n) !~ /^Math::Big/) ) {
+ return $_Config{'maxprime'} if $n <= $_Config{'maxprimeidx'};
+ croak "nth_prime_upper($n) overflow";
+ }
+
+ return int($upper + 1.0);
+}
+
+
+#############################################################################
+
+
+#############################################################################
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) {
+ if (ref($n) =~ /^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);
+ return _XS_RiemannR($n);
# We could make a new object, like:
# require Math::BigFloat;
@@ -356,14 +650,14 @@ 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) {
+ if (ref($n) =~ /^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);
+ return _XS_ExponentialIntegral($n);
}
sub LogarithmicIntegral {
@@ -381,6 +675,7 @@ sub LogarithmicIntegral {
ExponentialIntegral(log($n));
}
+#############################################################################
use Math::Prime::Util::MemFree;
@@ -479,6 +774,10 @@ Version 0.09
my $rand_prime = random_prime(100, 10000); # random prime within a range
my $rand_prime = random_ndigit_prime(6); # random 6-digit prime
+ # Euler phi on large number
+ use bigint; say euler_phi( 801294088771394680000412 );
+ # returns 391329671260448564651280
+
=head1 DESCRIPTION
@@ -487,19 +786,58 @@ methods, is_prime, prime_count, nth_prime, approximations and bounds for
the prime_count and nth prime, next_prime and prev_prime, factoring utilities,
and more.
-All routines currently work in native integers (32-bit or 64-bit). Bignum
-support may be added later. If you need bignum support for these types of
-functions inside Perl now, I recommend L<Math::Pari>.
-
The default sieving and factoring are intended to be (and currently are)
the fastest on CPAN, including L<Math::Prime::XS>, L<Math::Prime::FastSieve>,
-and L<Math::Factor::XS>. It seems to be faster than L<Math::Pari> for
-everything except factoring certain 16-20 digit numbers.
+L<Math::Factor::XS>, and L<Math::Prime::TiedArray>. For numbers in the 10-20
+digit range, it is often orders of magnitude faster. Typically it is faster
+than L<Math::Pari> for 64-bit operations, with the exception of factoring
+certain 16-20 digit numbers.
+
+The main development of the module has been for working with Perl UVs, so
+32-bit or 64-bit. Bignum support is limited. If you need full bignum
+support for these types of functions inside Perl now, I recommend L<Math::Pari>.
The module is thread-safe and allows concurrency between Perl threads while
still sharing a prime cache. It is not itself multithreaded. The one caveat
-is on Win32 where you must use C<precalc> if the function will use primes
-(C<primes>, C<prime_count> greater than 900M, C<nth_prime> greater than 45M).
+is on Win32 (non-Cygwin) where you must use C<precalc> if the function will
+use primes (C<primes>, C<prime_count> greater than 900M,
+C<nth_prime> greater than 45M).
+
+
+=head1 BIGNUM SUPPORT
+
+A number of the functions support big numbers, but currently not all. The
+ones that do:
+
+ prime_count_lower
+ prime_count_upper
+ prime_count_approx
+ nth_prime_lower
+ nth_prime_upper
+ nth_prime_approx
+ factor
+ all_factors
+ moebius
+ euler_phi
+ ExponentialIntegral
+ LogarithmicIntegral
+ RiemannR
+
+These still do not:
+
+ is_prime
+ is_prob_prime
+ miller_rabin
+ primes
+ next_prime
+ prev_prime
+ prime_count
+ nth_prime
+ random_prime
+ random_ndigit_prime
+
+It is possible to call the L<Math::Prime::Util::PP> versions directly, though
+performance may be very suboptimal.
=head1 FUNCTIONS
@@ -758,6 +1096,20 @@ and C<-1^t> if C<n> is a product of C<t> distinct primes. This is an
important function in prime number theory.
+=head2 euler_phi
+
+ say "The Euler totient of $n is ", euler_phi($n);
+
+Returns the Euler totient function (also called Euler's phi or phi function)
+for an integer value. This is an arithmetic function that counts the number
+of positive integers less than or equal to C<n> that are relatively prime to
+C<n>. Given the definition used, C<euler_phi> will return 0 for all
+C<n E<lt> 1>. This follows the logic used by SAGE. Mathematic/WolframAlpha
+returns 0 for input 0, but returns C<euler_phi(-n)> for C<n E<lt> 0>.
+
+
+
+
=head1 UTILITY FUNCTIONS
=head2 prime_precalc
@@ -791,6 +1143,22 @@ your routines were inside an eval that died, things will still get cleaned up.
If you call another function that uses a MemFree object, the cache will stay
in place because you still have an object.
+=head2 prime_get_config
+
+ my $cached_up_to = prime_get_config->{'precalc_to'};
+
+Returns a reference to a hash of the current settings. The hash is copy of
+the configuration, so changing it has no effect. The settings include:
+
+ precalc_to primes up to this number are calculated
+ maxbits the maximum number of bits for native operations
+ xs 0 or 1, indicating the XS code is running
+ gmp 0 or 1, indicating GMP code is available
+ maxparam the largest value for most functions, without bigint
+ maxdigits the max digits in a number, without bigint
+ maxprime the largest representable prime, without bigint
+ maxprimeidx the index of maxprime, without bigint
+
=head1 FACTORING FUNCTIONS
@@ -800,17 +1168,21 @@ in place because you still have an object.
my @factors = factor(3_369_738_766_071_892_021);
# returns (204518747,16476429743)
-Produces the prime factors of a positive number input. They may not be in
-numerical order. The special cases of C<n = 0> and C<n = 1> will
-return C<n>, which guarantees multiplying the factors together will
-always result in the input value, though those are the only cases where
-the returned factors are not prime.
+Produces the prime factors of a positive number input, in numerical order.
+The special cases of C<n = 0> and C<n = 1> will return C<n>, which
+guarantees multiplying the factors together will always result in the
+input value, though those are the only cases where the returned factors
+are not prime.
The current algorithm is to use trial division for small numbers, while large
numbers go through a sequence of small trials, SQUFOF, Pollard's Rho, Hart's
one line factorization, and finally trial division for any survivors. This
process is repeated for each non-prime factor.
+While factoring works on bigints, the algorithms are currently set up for
+smaller numbers, and bignum support is all in pure Perl. Hence, it will be
+somewhat slow for "easy" numbers and too slow for "hard" numbers.
+
=head2 all_factors
@@ -945,11 +1317,6 @@ Print pseudoprimes base 17:
I have not completed testing all the functions near the word size limit
(e.g. C<2^32> for 32-bit machines). Please report any problems you find.
-The extra factoring algorithms are mildly interesting but really limited by
-not being big number aware. Assuming a desktop PC, every 32-bit number
-should be factored by the main routine in a few microseconds, and 64-bit
-numbers should be a few milliseconds at worst.
-
Perl versions earlier than 5.8.0 have issues with 64-bit that show up in the
factoring tests. The test suite will try to determine if your Perl is broken.
If you use later versions of Perl, or Perl 5.6.2 32-bit, or Perl 5.6.2 64-bit
@@ -1023,7 +1390,7 @@ The differences are in the implementations:
- 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.
+ 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
@@ -1052,7 +1419,8 @@ the extreme end). Pari's underlying algorithms and code are very
sophisticated, and will always be more so than this module, and of course
supports bignums which is a huge advantage. Small numbers factor much, much
faster with Math::Prime::Util. Pari passes M::P::U in speed somewhere in the
-16 digit range and rapidly increases its lead.
+16 digit range and rapidly increases its lead. For bignums, there is no
+question that Math::Pari is far superior at this point.
The presentation here:
L<http://math.boisestate.edu/~liljanab/BOISECRYPTFall09/Jacobsen.pdf>
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 74047d9..5043faa 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -28,7 +28,7 @@ BEGIN {
: 32;
no Config;
}
-sub _maxbits { $_uv_size }
+sub _PP_prime_maxbits { $_uv_size }
# If $n < $_half_word, then $n*$n will be exact.
my $_half_word = (~0 == 18446744073709551615) ? 4294967296 : # 64-bit
@@ -276,7 +276,7 @@ sub primes {
sub next_prime {
my($n) = @_;
croak "Input must be a positive integer" unless _is_positive_int($n);
- return 0 if ($n >= ((_maxbits == 32) ? 4294967291 : 18446744073709551557))
+ return 0 if ($n >= ((_PP_prime_maxbits == 32) ? 4294967291 : 18446744073709551557))
&& (!defined $bigint::VERSION);
return $_prime_next_small[$n] if $n <= $#_prime_next_small;
@@ -354,177 +354,6 @@ sub prime_count {
$count;
}
-sub prime_count_lower {
- my($x) = @_;
- croak "Input must be a positive integer" unless _is_positive_int($x);
-
- return $_prime_count_small[$x] if $x <= $#_prime_count_small;
-
- my $flogx = log($x);
-
- return int( $x / ($flogx - 0.7) ) if $x < 599;
-
- my $a;
- if ($x < 2700) { $a = 0.30; }
- elsif ($x < 5500) { $a = 0.90; }
- elsif ($x < 19400) { $a = 1.30; }
- elsif ($x < 32299) { $a = 1.60; }
- elsif ($x < 176000) { $a = 1.80; }
- elsif ($x < 315000) { $a = 2.10; }
- elsif ($x < 1100000) { $a = 2.20; }
- elsif ($x < 4500000) { $a = 2.31; }
- elsif ($x < 233000000) { $a = 2.36; }
- elsif ($x < 5433800000) { $a = 2.32; }
- elsif ($x <60000000000) { $a = 2.15; }
- else { $a = 1.80; } # Dusart 1999, page 14
-
- return int( ($x/$flogx) * (1.0 + 1.0/$flogx + $a/($flogx*$flogx)) );
-}
-
-sub prime_count_upper {
- my($x) = @_;
- croak "Input must be a positive integer" unless _is_positive_int($x);
-
- return $_prime_count_small[$x] if $x <= $#_prime_count_small;
-
- my $flogx = log($x);
-
- return int( ($x / ($flogx - 1.048)) + 1.0 ) if $x < 1621;
- return int( ($x / ($flogx - 1.071)) + 1.0 ) if $x < 5000;
- return int( ($x / ($flogx - 1.098)) + 1.0 ) if $x < 15900;
-
- my $a;
- if ($x < 24000) { $a = 2.30; }
- elsif ($x < 59000) { $a = 2.48; }
- elsif ($x < 350000) { $a = 2.52; }
- elsif ($x < 355991) { $a = 2.54; }
- elsif ($x < 356000) { $a = 2.51; }
- elsif ($x < 3550000) { $a = 2.50; }
- elsif ($x < 3560000) { $a = 2.49; }
- elsif ($x < 5000000) { $a = 2.48; }
- elsif ($x < 8000000) { $a = 2.47; }
- elsif ($x < 13000000) { $a = 2.46; }
- elsif ($x < 18000000) { $a = 2.45; }
- elsif ($x < 31000000) { $a = 2.44; }
- elsif ($x < 41000000) { $a = 2.43; }
- elsif ($x < 48000000) { $a = 2.42; }
- elsif ($x < 119000000) { $a = 2.41; }
- elsif ($x < 182000000) { $a = 2.40; }
- elsif ($x < 192000000) { $a = 2.395; }
- elsif ($x < 213000000) { $a = 2.390; }
- elsif ($x < 271000000) { $a = 2.385; }
- elsif ($x < 322000000) { $a = 2.380; }
- elsif ($x < 400000000) { $a = 2.375; }
- elsif ($x < 510000000) { $a = 2.370; }
- elsif ($x < 682000000) { $a = 2.367; }
- elsif ($x <60000000000) { $a = 2.362; }
- else { $a = 2.51; }
-
- return int( ($x/$flogx) * (1.0 + 1.0/$flogx + $a/($flogx*$flogx)) + 1.0 );
-}
-
-sub prime_count_approx {
- my($x) = @_;
- croak "Input must be a positive integer" unless _is_positive_int($x);
-
- return $_prime_count_small[$x] if $x <= $#_prime_count_small;
-
- #return int( (prime_count_upper($x) + prime_count_lower($x)) / 2);
- return int(RiemannR($x)+0.5);
-}
-
-sub nth_prime_lower {
- my($n) = @_;
- croak "Input must be a positive integer" unless _is_positive_int($n);
-
- return $_primes_small[$n] if $n <= $#_primes_small;
-
- my $flogn = log($n);
- my $flog2n = log($flogn);
-
- # Dusart 1999 page 14, for all n >= 2
- my $lower = $n * ($flogn + $flog2n - 1.0 + (($flog2n-2.25)/$flogn));
-
- if ( ($lower >= ~0) && (!defined $bignum::VERSION) ) {
- if (_maxbits == 32) {
- return 4294967291 if $n <= 203280221;
- } else {
- return 18446744073709551557 if $n <= 425656284035217743;
- }
- croak "nth_prime_lower($n) overflow";
- }
- return int($lower);
-}
-
-sub nth_prime_upper {
- my($n) = @_;
- croak "Input must be a positive integer" unless _is_positive_int($n);
-
- return $_primes_small[$n] if $n <= $#_primes_small;
-
- my $flogn = log($n);
- my $flog2n = log($flogn);
-
- my $upper;
- if ($n >= 39017) {
- $upper = $n * ( $flogn + $flog2n - 0.9484 ); # Dusart 1999 page 14
- } elsif ($n >= 27076) {
- $upper = $n * ( $flogn + $flog2n - 1.0 + (($flog2n-1.80)/$flogn) );
- } elsif ($n >= 7022) {
- $upper = $n * ( $flogn + 0.9385 * $flog2n ); # Robin 1983
- } else {
- $upper = $n * ( $flogn + $flog2n );
- }
-
- if ( ($upper >= ~0) && (!defined $bignum::VERSION) ) {
- if (_maxbits == 32) {
- return 4294967291 if $n <= 203280221;
- } else {
- return 18446744073709551557 if $n <= 425656284035217743;
- }
- croak "nth_prime_upper($n) overflow";
- }
- return int($upper + 1.0);
-}
-
-sub nth_prime_approx {
- my($n) = @_;
- croak "Input must be a positive integer" unless _is_positive_int($n);
-
- return $_primes_small[$n] if $n <= $#_primes_small;
-
- my $flogn = log($n);
- my $flog2n = log($flogn);
-
- my $approx = $n * ( $flogn + $flog2n - 1
- + (($flog2n - 2)/$flogn)
- - ((($flog2n*$flog2n) - 6*$flog2n + 11) / (2*$flogn*$flogn))
- );
-
- my $order = $flog2n/$flogn;
- $order = $order*$order*$order * $n;
-
- if ($n < 259) { $approx += 10.4 * $order; }
- elsif ($n < 775) { $approx += 7.52* $order; }
- elsif ($n < 1271) { $approx += 5.6 * $order; }
- elsif ($n < 2000) { $approx += 5.2 * $order; }
- elsif ($n < 4000) { $approx += 4.3 * $order; }
- elsif ($n < 12000) { $approx += 3.0 * $order; }
- elsif ($n < 150000) { $approx += 2.1 * $order; }
- elsif ($n < 200000000) { $approx += 0.0 * $order; }
- else { $approx += -0.010 * $order; }
-
- if ( ($approx >= ~0) && (!defined $bignum::VERSION) ) {
- if (_maxbits == 32) {
- return 4294967291 if $n <= 203280221;
- } else {
- return 18446744073709551557 if $n <= 425656284035217743;
- }
- croak "nth_prime_approx($n) overflow";
- }
-
- return int($approx + 0.5);
-}
sub nth_prime {
my($n) = @_;
@@ -533,7 +362,7 @@ sub nth_prime {
return $_primes_small[$n] if $n <= $#_primes_small;
if (!defined $bigint::VERSION) {
- if (_maxbits == 32) {
+ if (_PP_prime_maxbits == 32) {
croak "nth_prime($n) overflow" if $n > 203280221;
} else {
croak "nth_prime($n) overflow" if $n > 425656284035217743;
@@ -698,12 +527,16 @@ sub is_prob_prime {
return 2 if $n < 121;
my @bases;
+ # We can do this, if MR accepts bases larger than the number
+ #if ($n < 316349281) { @bases = (11000544, 31481107); }
if ($n < 9080191) { @bases = (31, 73); }
elsif ($n < 4759123141) { @bases = (2, 7, 61); }
elsif ($n < 105936894253) { @bases = (2, 1005905886, 1340600841); }
elsif ($n < 31858317218647) { @bases = (2, 642735, 553174392, 3046413974); }
elsif ($n < 3071837692357849) { @bases = (2, 75088, 642735, 203659041, 3613982119); }
- else { @bases = (2, 325, 9375, 28178, 450775, 9780504, 1795265022); }
+ elsif ($n < 18446744073709551615) { @bases = (2, 325, 9375, 28178, 450775, 9780504, 1795265022); }
+ # Seriously we need BPSW here.
+ else { @bases = (2, 325, 9375, 28178, 450775, 9780504, 1795265022, 15, 52, 311); }
# Run our selected set of Miller-Rabin strong probability tests
my $prob_prime = miller_rabin($n, @bases);
@@ -763,6 +596,7 @@ sub factor {
my @factors = _basic_factor($n);
return @factors if $n < 4;
+ # Use 'n = int($n/7)' instead of 'n/=7' to not "upgrade" n to a Math::BigFloat.
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); }
@@ -770,23 +604,46 @@ sub factor {
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;
+ if ($n < (31*31)) {
+ push @factors, $n if $n != 1;
+ return @factors;
+ }
my @nstack = ($n);
while (@nstack) {
$n = pop @nstack;
+ # Don't use bignum on $n if it has gotten small enough.
+ $n = $n->numify if ref($n) =~ /^Math::Big/ && $n <= ~0;
#print "Looking at $n with stack ", join(",", at nstack), "\n";
while ( ($n >= (31*31)) && !_is_prime7($n) ) {
my @ftry;
- @ftry = prho_factor($n, 64*1024);
- @ftry = holf_factor($n, 64*1024) if scalar @ftry == 1;
+ my $holf_rounds = 0;
+ if ($n <= ~0) {
+ $holf_rounds = 64*1024;
+ #warn "trying holf 64k on $n\n";
+ @ftry = holf_factor($n, $holf_rounds);
+ }
+ if (scalar @ftry < 2) {
+ foreach my $add (3, 5, 7, 11, 13) {
+ #warn "trying prho 64k {$add} on $n\n" if scalar @ftry < 2;
+ @ftry = prho_factor($n, 64*1024, $add) if scalar @ftry < 2;
+ }
+ }
+ if (scalar @ftry < 2) {
+ #warn "trying holf 128k on $n\n";
+ @ftry = holf_factor($n, 128*1024, $holf_rounds);
+ $holf_rounds += 128*1024;
+ }
+ if (scalar @ftry < 2) {
+ #warn "trying prho 128k {17} on $n\n";
+ @ftry = prho_factor($n, 128*1024, 17);
+ }
if (scalar @ftry > 1) {
#print " split into ", join(",", at ftry), "\n";
$n = shift @ftry;
push @nstack, @ftry;
} else {
+ #warn "trial factor $n\n";
push @factors, trial_factor($n);
#print " trial into ", join(",", at factors), "\n";
$n = 1;
@@ -795,7 +652,7 @@ sub factor {
}
push @factors, $n if $n != 1;
}
- @factors;
+ sort {$a<=>$b} @factors;
}
# TODO:
@@ -803,37 +660,58 @@ sub fermat_factor { trial_factor(@_) }
sub squfof_factor { trial_factor(@_) }
sub prho_factor {
- my($n, $rounds) = @_;
+ my($n, $rounds, $a) = @_;
croak "Input must be a positive integer" unless _is_positive_int($n);
$rounds = 4*1024*1024 unless defined $rounds;
+ $a = 3 unless defined $a;
my @factors = _basic_factor($n);
return @factors if $n < 4;
my $inloop = 0;
- my $a = 3;
my $U = 7;
my $V = 7;
- for my $i (1 .. $rounds) {
- # U^2+a % n
- $U = _mulmod($U, $U, $n);
- $U = (($n-$U) > $a) ? $U+$a : $U+$a-$n;
- # V^2+a % n
- $V = _mulmod($V, $V, $n);
- $V = (($n-$V) > $a) ? $V+$a : $V+$a-$n;
- # V^2+a % n
- $V = _mulmod($V, $V, $n);
- $V = (($n-$V) > $a) ? $V+$a : $V+$a-$n;
- my $f = _gcd_ui( ($U > $V) ? $U-$V : $V-$U, $n );
- if ($f == $n) {
- last if $inloop++; # We've been here before
- } elsif ($f != 1) {
- push @factors, $f;
- push @factors, int($n/$f);
- croak "internal error in prho" unless ($f * int($n/$f)) == $n;
- return @factors;
+ if ( ($n < $_half_word) || (ref($n) =~ /^Math::Big/) ) {
+
+ for my $i (1 .. $rounds) {
+ $U = ($U * $U + $a) % $n;
+ $V = ($V * $V + $a) % $n;
+ $V = ($V * $V + $a) % $n;
+ my $f = _gcd_ui( ($U > $V) ? $U-$V : $V-$U, $n );
+ if ($f == $n) {
+ last if $inloop++; # We've been here before
+ } elsif ($f != 1) {
+ push @factors, $f;
+ push @factors, int($n/$f);
+ croak "internal error in prho" unless ($f * int($n/$f)) == $n;
+ return @factors;
+ }
}
+
+ } else {
+
+ for my $i (1 .. $rounds) {
+ # U^2+a % n
+ $U = _mulmod($U, $U, $n);
+ $U = (($n-$U) > $a) ? $U+$a : $U+$a-$n;
+ # V^2+a % n
+ $V = _mulmod($V, $V, $n);
+ $V = (($n-$V) > $a) ? $V+$a : $V+$a-$n;
+ # V^2+a % n
+ $V = _mulmod($V, $V, $n);
+ $V = (($n-$V) > $a) ? $V+$a : $V+$a-$n;
+ my $f = _gcd_ui( ($U > $V) ? $U-$V : $V-$U, $n );
+ if ($f == $n) {
+ last if $inloop++; # We've been here before
+ } elsif ($f != 1) {
+ push @factors, $f;
+ push @factors, int($n/$f);
+ croak "internal error in prho" unless ($f * int($n/$f)) == $n;
+ return @factors;
+ }
+ }
+
}
push @factors, $n;
@factors;
@@ -847,22 +725,40 @@ sub pbrent_factor {
my @factors = _basic_factor($n);
return @factors if $n < 4;
- my $a = 11;
+ my $a = 1;
my $Xi = 2;
my $Xm = 2;
- for my $i (1 .. $rounds) {
- # Xi^2+a % n
- $Xi = _mulmod($Xi, $Xi, $n);
- $Xi = (($n-$Xi) > $a) ? $Xi+$a : $Xi+$a-$n;
- my $f = _gcd_ui( ($Xi > $Xm) ? $Xi-$Xm : $Xm-$Xi, $n );
- if ( ($f != 1) && ($f != $n) ) {
- push @factors, $f;
- push @factors, int($n/$f);
- croak "internal error in pbrent" unless ($f * int($n/$f)) == $n;
- return @factors;
+ if ( ($n < $_half_word) || (ref($n) =~ /^Math::Big/) ) {
+
+ for my $i (1 .. $rounds) {
+ $Xi = ($Xi * $Xi + $a) % $n;
+ my $f = _gcd_ui( ($Xi > $Xm) ? $Xi-$Xm : $Xm-$Xi, $n );
+ if ( ($f != 1) && ($f != $n) ) {
+ push @factors, $f;
+ push @factors, int($n/$f);
+ croak "internal error in pbrent" unless ($f * int($n/$f)) == $n;
+ return @factors;
+ }
+ $Xm = $Xi if ($i & ($i-1)) == 0; # i is a power of 2
+ }
+
+ } else {
+
+ for my $i (1 .. $rounds) {
+ # Xi^2+a % n
+ $Xi = _mulmod($Xi, $Xi, $n);
+ $Xi = (($n-$Xi) > $a) ? $Xi+$a : $Xi+$a-$n;
+ my $f = _gcd_ui( ($Xi > $Xm) ? $Xi-$Xm : $Xm-$Xi, $n );
+ if ( ($f != 1) && ($f != $n) ) {
+ push @factors, $f;
+ push @factors, int($n/$f);
+ croak "internal error in pbrent" unless ($f * int($n/$f)) == $n;
+ return @factors;
+ }
+ $Xm = $Xi if ($i & ($i-1)) == 0; # i is a power of 2
}
- $Xm = $Xi if ($i & ($i-1)) == 0; # i is a power of 2
+
}
push @factors, $n;
@factors;
@@ -894,20 +790,22 @@ sub pminus1_factor {
}
sub holf_factor {
- my($n, $rounds) = @_;
+ my($n, $rounds, $startrounds) = @_;
croak "Input must be a positive integer" unless _is_positive_int($n);
$rounds = 64*1024*1024 unless defined $rounds;
+ $startrounds = 1 unless defined $startrounds;
my @factors = _basic_factor($n);
return @factors if $n < 4;
- for my $i (1 .. $rounds) {
+ for my $i ($startrounds .. $rounds) {
my $s = int(sqrt($n * $i));
$s++ if ($s * $s) != ($n * $i);
my $m = ($s < $_half_word) ? ($s*$s) % $n : _mulmod($s, $s, $n);
# Check for perfect square
my $mcheck = $m & 127;
next if (($mcheck*0x8bc40d7d) & ($mcheck*0xa1e2f5d1) & 0x14020a);
+ # ... 82% of non-squares were rejected by the bloom filter
my $f = int(sqrt($m));
next unless $f*$f == $m;
$f = _gcd_ui( ($s > $f) ? $s - $f : $f - $s, $n);
diff --git a/t/03-init.t b/t/03-init.t
index 079b07a..3dbf12e 100644
--- a/t/03-init.t
+++ b/t/03-init.t
@@ -1,27 +1,27 @@
#!/usr/bin/env perl
use strict;
use warnings;
-use Math::Prime::Util qw/prime_precalc prime_memfree/;
+use Math::Prime::Util qw/prime_precalc prime_memfree prime_get_config/;
use Test::More tests => 3 + 3 + 3 + 6;
my $bigsize = 10_000_000;
-# How do we test these? Use some private functions and assume things about
-# internal behavior. Not the best idea.
+# This is still a slightly dubious assumption, that the precalc size _must_
+# go up when we request it.
-can_ok( 'Math::Prime::Util', '_get_prime_cache_size' );
+can_ok( 'Math::Prime::Util', 'prime_get_config' );
-my $init_size = Math::Prime::Util::_get_prime_cache_size;
+my $init_size = prime_get_config->{'precalc_to'};
prime_precalc($bigsize);
-cmp_ok( Math::Prime::Util::_get_prime_cache_size, '>', $init_size, "Internal space grew after large precalc" );
+cmp_ok( prime_get_config->{'precalc_to'}, '>', $init_size, "Internal space grew after large precalc" );
prime_memfree;
-is( Math::Prime::Util::_get_prime_cache_size, $init_size, "Internal space went back to original size after memfree" );
+is( prime_get_config->{'precalc_to'}, $init_size, "Internal space went back to original size after memfree" );
# Now do the object way.
@@ -30,30 +30,30 @@ is( Math::Prime::Util::_get_prime_cache_size, $init_size, "Internal space went b
my $mf = Math::Prime::Util::MemFree->new;
isa_ok $mf, 'Math::Prime::Util::MemFree';
prime_precalc($bigsize);
- cmp_ok( Math::Prime::Util::_get_prime_cache_size, '>', $init_size, "Internal space grew after large precalc" );
+ cmp_ok( prime_get_config->{'precalc_to'}, '>', $init_size, "Internal space grew after large precalc" );
}
-is( Math::Prime::Util::_get_prime_cache_size, $init_size, "Memory released after MemFree object goes out of scope");
+is( prime_get_config->{'precalc_to'}, $init_size, "Memory released after MemFree object goes out of scope");
# Wrap multiple calls, make sure we wait until the last one is done.
{
my $mf = Math::Prime::Util::MemFree->new;
prime_precalc($bigsize);
- cmp_ok( Math::Prime::Util::_get_prime_cache_size, '>', $init_size, "Internal space grew after large precalc" );
+ cmp_ok( prime_get_config->{'precalc_to'}, '>', $init_size, "Internal space grew after large precalc" );
{
my $mf2 = Math::Prime::Util::MemFree->new;
prime_precalc( 2 * $bigsize );
}
- cmp_ok( Math::Prime::Util::_get_prime_cache_size, '>', $init_size, "Memory not freed yet because a MemFree object still live." );
+ cmp_ok( prime_get_config->{'precalc_to'}, '>', $init_size, "Memory not freed yet because a MemFree object still live." );
}
-is( Math::Prime::Util::_get_prime_cache_size, $init_size, "Memory released after last MemFree object goes out of scope");
+is( prime_get_config->{'precalc_to'}, $init_size, "Memory released after last MemFree object goes out of scope");
# Show how an eval death can leak
-eval { prime_precalc($bigsize); cmp_ok( Math::Prime::Util::_get_prime_cache_size, '>', $init_size, "Internal space grew after large precalc" ); prime_memfree; };
-is( Math::Prime::Util::_get_prime_cache_size, $init_size, "Memory freed after successful eval");
+eval { prime_precalc($bigsize); cmp_ok( prime_get_config->{'precalc_to'}, '>', $init_size, "Internal space grew after large precalc" ); prime_memfree; };
+is( prime_get_config->{'precalc_to'}, $init_size, "Memory freed after successful eval");
-eval { prime_precalc($bigsize); cmp_ok( Math::Prime::Util::_get_prime_cache_size, '>', $init_size, "Internal space grew after large precalc" ); die; prime_memfree; };
-isnt( Math::Prime::Util::_get_prime_cache_size, $init_size, "Memory normally not freed after eval die");
+eval { prime_precalc($bigsize); cmp_ok( prime_get_config->{'precalc_to'}, '>', $init_size, "Internal space grew after large precalc" ); die; prime_memfree; };
+isnt( prime_get_config->{'precalc_to'}, $init_size, "Memory normally not freed after eval die");
prime_memfree;
-eval { my $mf = Math::Prime::Util::MemFree->new; prime_precalc($bigsize); cmp_ok( Math::Prime::Util::_get_prime_cache_size, '>', $init_size, "Internal space grew after large precalc" ); die; };
-is( Math::Prime::Util::_get_prime_cache_size, $init_size, "Memory is freed after eval die using object scoper");
+eval { my $mf = Math::Prime::Util::MemFree->new; prime_precalc($bigsize); cmp_ok( prime_get_config->{'precalc_to'}, '>', $init_size, "Internal space grew after large precalc" ); die; };
+is( prime_get_config->{'precalc_to'}, $init_size, "Memory is freed after eval die using object scoper");
diff --git a/t/10-isprime.t b/t/10-isprime.t
index 3e96b1d..4ca9e60 100644
--- a/t/10-isprime.t
+++ b/t/10-isprime.t
@@ -5,7 +5,7 @@ use warnings;
use Test::More;
use Math::Prime::Util qw/is_prime/;
-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/11-primes.t b/t/11-primes.t
index a63b213..23a3434 100644
--- a/t/11-primes.t
+++ b/t/11-primes.t
@@ -5,7 +5,7 @@ use warnings;
use Test::More;
use Math::Prime::Util qw/primes prime_count/;
-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 => 12+3 + 12 + 1 + 19 + ($use64 ? 1 : 0) + 1 + 13*5;
diff --git a/t/12-nextprime.t b/t/12-nextprime.t
index dd8a4eb..5f15a4b 100644
--- a/t/12-nextprime.t
+++ b/t/12-nextprime.t
@@ -5,7 +5,7 @@ use warnings;
use Test::More;
use Math::Prime::Util qw/next_prime prev_prime/;
-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 => 3573*2 + 3*2 + 6 + 2 + 148 + 148 + 1;
diff --git a/t/13-primecount.t b/t/13-primecount.t
index f711515..95682c3 100644
--- a/t/13-primecount.t
+++ b/t/13-primecount.t
@@ -5,7 +5,7 @@ use warnings;
use Test::More;
use Math::Prime::Util qw/prime_count prime_count_lower prime_count_upper prime_count_approx/;
-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};
# Powers of 2: http://oeis.org/A007053/b007053.txt
diff --git a/t/14-nthprime.t b/t/14-nthprime.t
index 19f660d..a5c1840 100644
--- a/t/14-nthprime.t
+++ b/t/14-nthprime.t
@@ -5,7 +5,7 @@ use warnings;
use Test::More;
use Math::Prime::Util qw/primes nth_prime nth_prime_lower nth_prime_upper nth_prime_approx/;
-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};
my $broken64 = (18446744073709550592 == ~0);
diff --git a/t/16-randomprime.t b/t/16-randomprime.t
index 277b09e..eff57e0 100644
--- a/t/16-randomprime.t
+++ b/t/16-randomprime.t
@@ -8,7 +8,7 @@ use Test::More;
#sub rand { return 0.5; }
use Math::Prime::Util qw/random_prime random_ndigit_prime is_prime/;
-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 => 13+4 + 6*1 + 11*3 + 10*2 + 12*2 + 2*$use64 + 10*2;
diff --git a/t/17-pseudoprime.t b/t/17-pseudoprime.t
index ae07392..f2f04ed 100644
--- a/t/17-pseudoprime.t
+++ b/t/17-pseudoprime.t
@@ -5,7 +5,7 @@ use warnings;
use Test::More;
use Math::Prime::Util qw/is_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};
# small primes
diff --git a/t/18-functions.t b/t/18-functions.t
index 04c447d..a6415eb 100644
--- a/t/18-functions.t
+++ b/t/18-functions.t
@@ -5,7 +5,7 @@ use warnings;
use Test::More;
use Math::Prime::Util qw/prime_count ExponentialIntegral LogarithmicIntegral RiemannR/;
-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 => 4 + 1 + 1 + 16 + 11 + 9;
diff --git a/t/30-relations.t b/t/30-relations.t
index d1b907a..282bc3f 100644
--- a/t/30-relations.t
+++ b/t/30-relations.t
@@ -12,7 +12,7 @@ use Math::Prime::Util qw/primes
my @trials = qw/1 2 3 4 5 6 7 17 57 89 102 1337 8573 84763 784357 1000001 2573622/;
-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 => 5 * scalar @trials;
diff --git a/t/50-factoring.t b/t/50-factoring.t
index 7590155..c826416 100644
--- a/t/50-factoring.t
+++ b/t/50-factoring.t
@@ -5,7 +5,7 @@ use warnings;
use Test::More;
use Math::Prime::Util qw/factor all_factors is_prime/;
-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};
diff --git a/t/80-pp.t b/t/80-pp.t
index 87c1a07..a19dc19 100644
--- a/t/80-pp.t
+++ b/t/80-pp.t
@@ -4,6 +4,7 @@ use warnings;
# This is a subset of our tests. You really should run the whole test suite
# on the PP code. What this will do is basic regression testing.
+my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
use Test::More;
my @small_primes = qw/
@@ -128,13 +129,10 @@ my %pivals32 = (
4294967295 => 203280221,
);
my %pivals_small = map { $_ => $pivals32{$_} }
- grep {$_ <= 2000000}
+ grep {$_ <= 200000}
keys %pivals32;
my %pi_intervals = (
- "868396 to 9478505" => 563275,
- "1118105 to 9961674" => 575195,
- "24689 to 7973249" => 535368,
"1e10 +2**16" => 2821,
"17 to 13" => 0,
"3 to 17" => 6,
@@ -145,6 +143,14 @@ my %pi_intervals = (
"191912783 +247" => 1,
"191912784 +246" => 0,
);
+my %extra_pi_intervals = (
+ "868396 to 9478505" => 563275,
+ "1118105 to 9961674" => 575195,
+ "24689 to 7973249" => 535368,
+);
+# Add extra intervals to pi_intervals if we're doing release testing
+ at pi_intervals{keys %extra_pi_intervals} = values %extra_pi_intervals if $extra;
+
# Remove any entries where the high value is too large for us
# ikegami++ for the delete from a hash slice idea
delete @pi_intervals{ grep { (parse_range($_))[1] > ~0 } keys %pi_intervals };
@@ -161,7 +167,7 @@ my %nthprimes32 = (
100000000 => 2038074743,
);
my %nthprimes_small = map { $_ => $nthprimes32{$_} }
- grep {$_ <= 2000000}
+ grep {$_ <= 200000}
keys %nthprimes32;
my %eivals = (
@@ -212,27 +218,21 @@ plan tests => 1 +
2*(1087 + @primes + @composites) +
3 + scalar(keys %small_single) + scalar(keys %small_range) +
2*scalar(keys %primegaps) + 8 + 148 + 148 + 1 +
- 3*scalar(keys %pivals32) + scalar(keys %pivals_small) + scalar(keys %pi_intervals) +
- 2*scalar(keys %pivals_small) + 3*scalar(keys %nthprimes32) + scalar(keys %nthprimes_small) +
+ scalar(keys %pivals_small) + scalar(keys %pi_intervals) +
+ 2*scalar(keys %pivals_small) + scalar(keys %nthprimes_small) +
4 + $num_pseudoprimes +
scalar(keys %eivals) + scalar(keys %livals) + scalar(keys %rvals) +
scalar @primes + 3*scalar @composites +
0;
-use Math::Prime::Util qw/primes/;
+use Math::Prime::Util qw/primes prime_count_approx prime_count_lower/;
require_ok 'Math::Prime::Util::PP';
# This function skips some setup
undef *primes;
*primes = \&Math::Prime::Util::PP::primes;
*prime_count = \&Math::Prime::Util::PP::prime_count;
- *prime_count_upper = \&Math::Prime::Util::PP::prime_count_upper;
- *prime_count_lower = \&Math::Prime::Util::PP::prime_count_lower;
- *prime_count_approx = \&Math::Prime::Util::PP::prime_count_approx;
*nth_prime = \&Math::Prime::Util::PP::nth_prime;
- *nth_prime_upper = \&Math::Prime::Util::PP::nth_prime_upper;
- *nth_prime_lower = \&Math::Prime::Util::PP::nth_prime_lower;
- *nth_prime_approx = \&Math::Prime::Util::PP::nth_prime_approx;
*is_prime = \&Math::Prime::Util::PP::is_prime;
*next_prime = \&Math::Prime::Util::PP::next_prime;
@@ -312,13 +312,6 @@ is(next_prime(1234567890), 1234567891, "next_prime(1234567890) == 1234567891)");
###############################################################################
-while (my($n, $pin) = each (%pivals32)) {
- cmp_ok( prime_count_upper($n), '>=', $pin, "Pi($n) <= upper estimate" );
- cmp_ok( prime_count_lower($n), '<=', $pin, "Pi($n) >= lower estimate" );
- my $approx_range = abs($pin - prime_count_approx($n));
- my $range_limit = ($n <= 100000000) ? 100 : 500;
- cmp_ok( $approx_range, '<=', $range_limit, "prime_count_approx($n) within $range_limit");
-}
while (my($n, $pin) = each (%pivals_small)) {
is( prime_count($n), $pin, "Pi($n) = $pin" );
}
@@ -336,14 +329,6 @@ while (my($n, $pin) = each (%pivals_small)) {
cmp_ok( nth_prime($next), '>=', $n, "nth_prime($next) >= $n");
}
-while (my($n, $nth) = each (%nthprimes32)) {
- cmp_ok( nth_prime_upper($n), '>=', $nth, "nth_prime($n) <= upper estimate" );
- cmp_ok( nth_prime_lower($n), '<=', $nth, "nth_prime($n) >= lower estimate" );
-
- my $approx = nth_prime_approx($n);
- my $percent_limit = ($n >= 775) ? 1 : 2;
- cmp_ok( abs($nth - $approx) / $nth, '<=', $percent_limit/100.0, "nth_prime_approx($n) = $approx within $percent_limit\% of $nth");
-}
while (my($n, $nth) = each (%nthprimes_small)) {
is( nth_prime($n), $nth, "nth_prime($n) = $nth" );
}
diff --git a/util.c b/util.c
index 174d68f..f57c3d0 100644
--- a/util.c
+++ b/util.c
@@ -393,7 +393,7 @@ static const unsigned char prime_count_small[] =
#define NPRIME_COUNT_SMALL (sizeof(prime_count_small)/sizeof(prime_count_small[0]))
static const double F1 = 1.0;
-UV prime_count_lower(UV x)
+UV _XS_prime_count_lower(UV x)
{
double fx, flogx;
double a = 1.80; /* Dusart 1999, page 14 */
@@ -427,7 +427,7 @@ UV prime_count_lower(UV x)
}
-UV prime_count_upper(UV x)
+UV _XS_prime_count_upper(UV x)
{
double fx, flogx;
double a = 2.51; /* Dusart 1999, page 14 */
@@ -495,7 +495,7 @@ UV prime_count_upper(UV x)
}
-UV prime_count_approx(UV x)
+UV _XS_prime_count_approx(UV x)
{
/*
* A simple way:
@@ -522,8 +522,8 @@ UV prime_count_approx(UV x)
if (x < NPRIME_COUNT_SMALL)
return prime_count_small[x];
- //R = XS_LogarithmicIntegral(x) - XS_LogarithmicIntegral(sqrt(x))/2;
- R = XS_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)))
@@ -615,7 +615,7 @@ static const unsigned short primes_small[] =
#define NPRIMES_SMALL (sizeof(primes_small)/sizeof(primes_small[0]))
/* The nth prime will be greater than or equal to this number */
-UV nth_prime_lower(UV n)
+UV _XS_nth_prime_lower(UV n)
{
double fn = (double) n;
double flogn, flog2n, lower;
@@ -644,7 +644,7 @@ UV nth_prime_lower(UV n)
/* The nth prime will be less or equal to this number */
-UV nth_prime_upper(UV n)
+UV _XS_nth_prime_upper(UV n)
{
double fn = (double) n;
double flogn, flog2n, upper;
@@ -678,7 +678,7 @@ UV nth_prime_upper(UV n)
}
-UV nth_prime_approx(UV n)
+UV _XS_nth_prime_approx(UV n)
{
double fn, flogn, flog2n, order, approx;
@@ -761,7 +761,7 @@ UV nth_prime(UV n)
return primes_small[n];
/* Determine a bound on the nth prime. We know it comes before this. */
- upper_limit = nth_prime_upper(n);
+ upper_limit = _XS_nth_prime_upper(n);
MPUassert(upper_limit > 0, "nth_prime got an upper limit of 0");
/* Get the primary cache, and ensure it is at least this large. If the
@@ -835,7 +835,7 @@ UV nth_prime(UV n)
static double const euler_mascheroni = 0.57721566490153286060651209008240243104215933593992;
static double const li2 = 1.045163780117492784844588889194613136522615578151;
-double XS_ExponentialIntegral(double x) {
+double _XS_ExponentialIntegral(double x) {
double const tol = 1e-16;
double val, term, fact_n;
double y, t;
@@ -927,12 +927,12 @@ double XS_ExponentialIntegral(double x) {
return val;
}
-double XS_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 XS_ExponentialIntegral(log(x));
+ return _XS_ExponentialIntegral(log(x));
}
/*
@@ -1008,7 +1008,7 @@ static double evaluate_zeta(double x) {
return sum;
}
-double XS_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 b93c075..efab28c 100644
--- a/util.h
+++ b/util.h
@@ -9,19 +9,20 @@ extern UV next_trial_prime(UV x);
extern UV next_prime(UV x);
extern UV prev_prime(UV x);
-extern UV prime_count_lower(UV x);
-extern UV prime_count_upper(UV x);
-extern UV prime_count_approx(UV x);
extern UV prime_count(UV low, UV high);
-
-extern UV nth_prime_lower(UV n);
-extern UV nth_prime_upper(UV x);
-extern UV nth_prime_approx(UV x);
extern UV nth_prime(UV x);
-extern double XS_ExponentialIntegral(double x);
-extern double XS_LogarithmicIntegral(double x);
-extern double XS_RiemannR(double x);
+/* These have been moved into the main Util.pm */
+extern UV _XS_prime_count_lower(UV x);
+extern UV _XS_prime_count_upper(UV x);
+extern UV _XS_prime_count_approx(UV x);
+extern UV _XS_nth_prime_lower(UV n);
+extern UV _XS_nth_prime_upper(UV x);
+extern UV _XS_nth_prime_approx(UV 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
More information about the Pkg-perl-cvs-commits
mailing list