[libmath-prime-util-perl] 20/59: Lots of bignum changes, new tests, update version number
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:55 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 2188da0a629fbc2acdf2f1ccb1ddc5a2f334a2a3
Author: Dana Jacobsen <dana at acm.org>
Date: Sun Jul 1 20:33:35 2012 -0600
Lots of bignum changes, new tests, update version number
---
Changes | 2 +
MANIFEST | 1 +
README | 10 ++-
TODO | 14 ++--
XS.xs | 47 +++++++----
examples/bench-pp-sieve.pl | 2 +-
examples/test-nthapprox.pl | 6 +-
examples/test-pcapprox.pl | 13 +--
lib/Math/Prime/Util.pm | 113 ++++++++++++--------------
lib/Math/Prime/Util/MemFree.pm | 4 +-
lib/Math/Prime/Util/PP.pm | 93 ++++++++++++++++++++--
lib/Math/Prime/Util/PrimeArray.pm | 4 +-
t/02-can.t | 3 +-
t/81-bignum.t | 163 ++++++++++++++++++++++++++++++++++++++
14 files changed, 368 insertions(+), 107 deletions(-)
diff --git a/Changes b/Changes
index def2989..741dec2 100644
--- a/Changes
+++ b/Changes
@@ -8,6 +8,8 @@ Revision history for Perl extension Math::Prime::Util.
- 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.
+ - many more bigint/bignum changes
+ - factor always returns sorted results
0.09 25 June 2012
- Pure Perl code added. Passes all tests. Used only if the XSLoader
diff --git a/MANIFEST b/MANIFEST
index c633711..9003f63 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -54,6 +54,7 @@ t/31-threading.t
t/50-factoring.t
t/51-primearray.t
t/80-pp.t
+t/81-bignum.t
t/90-release-perlcritic.t
t/91-release-pod-syntax.t
t/92-release-pod-coverage.t
diff --git a/README b/README
index 9b6ebf1..6d6e1de 100644
--- a/README
+++ b/README
@@ -1,12 +1,14 @@
-Math::Prime::Util version 0.09
+Math::Prime::Util version 0.10
A set of utilities related to prime numbers. These include multiple sieving
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.
+the prime_count and nth prime, next_prime and prev_prime, moebius and totient
+functions, integer factoring, and more.
The default sieving and factoring are intended to be the fastest on CPAN,
-including Math::Prime::XS, Math::Prime::FastSieve, and Math::Factor::XS.
+including Math::Prime::XS, Math::Prime::FastSieve, Math::Factor::XS,
+Math::Primality, and Math::Prime::TiedArray. For non-bignums, it is typically
+faster than Math::Pari (and doesn't require Pari to be installed).
SYNOPSIS
diff --git a/TODO b/TODO
index 8e652c2..9736f7f 100644
--- a/TODO
+++ b/TODO
@@ -20,21 +20,25 @@
- Faster SQUFOF
-- better prime count upper/lower bounds
+- better prime count upper/lower bounds (under RH for 64-bit).
- Move .c / .h files into separate directory.
version does it in a painful way. Something simpler to be had?
-- More testing for bignum, including a test suite file
+- finish test suite for bignum. Work on making it faster.
- need next_prime and prev_prime bignum support. random_ndigit_prime needs it.
-- redo _XS_factor to return a sorted array
-
- Need to add more bignum factoring support:
- GMP prho
- GMP pminus1
+ - GMP SQUFOF
- GMP tinyqs
- That should at least make us usable for 80-90 bit numbers.
+ The first three should give us reasonable support up to ~30 digits. Adding
+ a tinyQS (e.g. yafu, msieve, flint) would both be faster and extend to ~40
+ digits (in a reasonable time). Going beyong that is going to need full
+ MPQS or SIQS. Another possibility is to see if the GMP-ECM library is
+ installed and call that, but I'm not sure how much it will help if we get
+ the above running.
- Must add BPSW if we want is_prime to be meaningful for >64-bit.
diff --git a/XS.xs b/XS.xs
index b828e94..c9f93ec 100644
--- a/XS.xs
+++ b/XS.xs
@@ -199,9 +199,11 @@ _XS_factor(IN UV n)
if (n < 4) {
XPUSHs(sv_2mortal(newSVuv( n ))); /* If n is 0-3, we're done. */
} else {
- UV tlim = 19; /* Below this we've checked */
- UV factor_stack[MPU_MAX_FACTORS+1];
- int nstack = 0;
+ UV tlim = 53; /* Below this we've checked with trial division */
+ UV tofac_stack[MPU_MAX_FACTORS+1];
+ UV factored_stack[MPU_MAX_FACTORS+1];
+ int ntofac = 0;
+ int nfactored = 0;
/* Quick trial divisions. Crude use of GCD to hopefully go faster. */
while ( (n% 2) == 0 ) { n /= 2; XPUSHs(sv_2mortal(newSVuv( 2 ))); }
if ( (n >= UVCONST(3*3)) && (gcd_ui(n, UVCONST(3234846615) != 1)) ) {
@@ -214,7 +216,6 @@ _XS_factor(IN UV n)
while ( (n%19) == 0 ) { n /= 19; XPUSHs(sv_2mortal(newSVuv( 19 ))); }
while ( (n%23) == 0 ) { n /= 23; XPUSHs(sv_2mortal(newSVuv( 23 ))); }
while ( (n%29) == 0 ) { n /= 29; XPUSHs(sv_2mortal(newSVuv( 29 ))); }
- tlim = 31;
}
if ( (n >= UVCONST(31*31)) && (gcd_ui(n, UVCONST(95041567) != 1)) ) {
while ( (n%31) == 0 ) { n /= 31; XPUSHs(sv_2mortal(newSVuv( 31 ))); }
@@ -222,7 +223,6 @@ _XS_factor(IN UV n)
while ( (n%41) == 0 ) { n /= 41; XPUSHs(sv_2mortal(newSVuv( 41 ))); }
while ( (n%43) == 0 ) { n /= 43; XPUSHs(sv_2mortal(newSVuv( 43 ))); }
while ( (n%47) == 0 ) { n /= 47; XPUSHs(sv_2mortal(newSVuv( 47 ))); }
- tlim = 53;
}
do { /* loop over each remaining factor */
while ( (n >= (tlim*tlim)) && (!is_definitely_prime(n)) ) {
@@ -230,21 +230,21 @@ _XS_factor(IN UV n)
if (n > UVCONST(10000000) ) { /* tune this */
/* For sufficiently large n, try more complex methods. */
/* SQUFOF (succeeds 98-99.9%) */
- split_success = squfof_factor(n, factor_stack+nstack, 256*1024)-1;
+ split_success = squfof_factor(n, tofac_stack+ntofac, 256*1024)-1;
/* A few rounds of Pollard rho (succeeds most of the rest) */
if (!split_success) {
- split_success = prho_factor(n, factor_stack+nstack, 400)-1;
+ split_success = prho_factor(n, tofac_stack+ntofac, 400)-1;
}
/* Some rounds of HOLF, good for close to perfect squares */
if (!split_success) {
- split_success = holf_factor(n, factor_stack+nstack, 2000)-1;
+ split_success = holf_factor(n, tofac_stack+ntofac, 2000)-1;
}
/* Less than 0.00003% of numbers make it past here. */
}
if (split_success) {
MPUassert( split_success == 1, "split factor returned more than 2 factors");
- nstack++;
- n = factor_stack[nstack];
+ ntofac++; /* Leave one on the to-be-factored stack */
+ n = tofac_stack[ntofac]; /* Set n to the other one */
} else {
/* trial divisions */
UV f = tlim;
@@ -253,7 +253,8 @@ _XS_factor(IN UV n)
while (f <= limit) {
if ( (n%f) == 0 ) {
do {
- n /= f; XPUSHs(sv_2mortal(newSVuv( f )));
+ n /= f;
+ factored_stack[nfactored++] = f;
} while ( (n%f) == 0 );
limit = (UV) (sqrt(n)+0.1);
}
@@ -263,9 +264,27 @@ _XS_factor(IN UV n)
break; /* We just factored n via trial division. Exit loop. */
}
}
- if (n != 1) XPUSHs(sv_2mortal(newSVuv( n )));
- if (nstack > 0) n = factor_stack[nstack-1];
- } while (nstack-- > 0);
+ /* n is now prime (or 1), so add to already-factored stack */
+ if (n != 1) factored_stack[nfactored++] = n;
+ /* Pop the next number off the to-factor stack */
+ if (ntofac > 0) n = tofac_stack[ntofac-1];
+ } while (ntofac-- > 0);
+ /* Now push all the factored results in sorted order */
+ {
+ int i, j;
+ for (i = 0; i < nfactored; i++) {
+ int mini = i;
+ for (j = i+1; j < nfactored; j++)
+ if (factored_stack[j] < factored_stack[mini])
+ mini = j;
+ if (mini != i) {
+ UV t = factored_stack[mini];
+ factored_stack[mini] = factored_stack[i];
+ factored_stack[i] = t;
+ }
+ XPUSHs(sv_2mortal(newSVuv( factored_stack[i] )));
+ }
+ }
}
#define SIMPLE_FACTOR(func, n, rounds) \
diff --git a/examples/bench-pp-sieve.pl b/examples/bench-pp-sieve.pl
index 681c568..a8620fa 100644
--- a/examples/bench-pp-sieve.pl
+++ b/examples/bench-pp-sieve.pl
@@ -3,7 +3,7 @@ use strict;
use warnings;
use Benchmark qw/:all/;
-use Devel::Size qw/total_size/;
+#use Devel::Size qw/total_size/;
use Math::Prime::Util;
use Math::Prime::FastSieve;
*mpu_erat = \&Math::Prime::Util::erat_primes;
diff --git a/examples/test-nthapprox.pl b/examples/test-nthapprox.pl
index 523a7b4..20b6490 100755
--- a/examples/test-nthapprox.pl
+++ b/examples/test-nthapprox.pl
@@ -3,8 +3,6 @@ use strict;
use warnings;
use Math::Prime::Util ":all";
$| = 1; # fast pipes
-my $use64 = Math::Prime::Util::_maxbits > 32;
-
my %nthprimes = (
1 => 2,
@@ -33,7 +31,7 @@ foreach my $n (sort {$a<=>$b} keys %nthprimes) {
my $nth = $nthprimes{$n};
my $ntha = nth_prime_approx($n);
- printf "10^%2d %13llu %12.7f\n", length($n)-1, abs($nth-$ntha), 100*($ntha-$nth)/$nth;
+ printf "10^%2d %13lu %12.7f\n", length($n)-1, abs($nth-$ntha), 100*($ntha-$nth)/$nth;
}
print "\n";
@@ -48,5 +46,5 @@ foreach my $n (sort {$a<=>$b} keys %nthprimes) {
my $nthu = nth_prime_upper($n);
printf "10^%2d %12.7f %12.7f\n",
- length($n)-1, 100*($nth-$nthl)/$nth, 100*($nthu-$nth)/$nth;
+ length($n)-1, 100.0*($nth-$nthl)/$nth, 100.0*($nthu-$nth)/$nth;
}
diff --git a/examples/test-pcapprox.pl b/examples/test-pcapprox.pl
index cae95d2..0f10cef 100755
--- a/examples/test-pcapprox.pl
+++ b/examples/test-pcapprox.pl
@@ -3,7 +3,6 @@ use strict;
use warnings;
use Math::Prime::Util qw/prime_count prime_count_approx prime_count_lower prime_count_upper LogarithmicIntegral RiemannR/;
$| = 1; # fast pipes
-my $use64 = Math::Prime::Util::_maxbits > 32;
my %pivals = (
@@ -28,18 +27,20 @@ my %pivals = (
10000000000000000000 => 234057667276344607,
);
-printf(" N %12s %12s %12s\n", "pc_approx", "Li", "R");
-printf("----- %12s %12s %12s\n", '-'x12,'-'x12,'-'x12);
+printf(" N %12s %12s %12s %12s\n", "pc_approx", "Li", "LiCor", "R");
+printf("----- %12s %12s %12s %12s\n", '-'x12,'-'x12,'-'x12,'-'x12);
foreach my $n (sort {$a<=>$b} keys %pivals) {
my $pin = $pivals{$n};
my $pca = prime_count_approx($n);
- my $pcli = ($n < 2) ? 0 : int(LogarithmicIntegral($n)-LogarithmicIntegral(2)+0.5);
+ my $Lisub = sub { my $x = shift; return ($x < 2) ? 0 : (LogarithmicIntegral($x)-LogarithmicIntegral(2)+0.5); };
+ my $pcli = int($Lisub->($n));
+ my $pclicor = int( $Lisub->($n) - ($Lisub->(sqrt($n)) / 2) );
my $r = int(RiemannR($n)+0.5);
- printf "10^%2d %12d %12d %12d\n", length($n)-1,
- abs($pca-$pin), abs($pcli-$pin), abs($r-$pin);
+ printf "10^%2d %12d %12d %12d %12d\n", length($n)-1,
+ abs($pca-$pin), abs($pcli-$pin), abs($pclicor-$pin), abs($r-$pin);
}
# Also see http://empslocal.ex.ac.uk/people/staff/mrwatkin/zeta/encoding1.htm
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index a70bc6c..633370f 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -5,7 +5,7 @@ use Carp qw/croak confess carp/;
BEGIN {
$Math::Prime::Util::AUTHORITY = 'cpan:DANAJ';
- $Math::Prime::Util::VERSION = '0.09';
+ $Math::Prime::Util::VERSION = '0.10';
}
# parent is cleaner, and in the Perl 5.10.1 / 5.12.0 core, but not earlier.
@@ -104,11 +104,29 @@ sub _validate_positive_integer {
croak "Parameter '$n' must be a positive integer" if $n =~ tr/0123456789//c;
croak "Parameter '$n' must be >= $min" if defined $min && $n < $min;
croak "Parameter '$n' must be <= $max" if defined $max && $n > $max;
- croak "Parameter '$n' outside of integer range" if $n > $_Config{'maxparam'}
- && ref($n) !~ /^Math::Big/;
+ if ($n > $_Config{'maxparam'}) {
+ croak "Parameter '$n' outside of integer range" if !defined $Math::BigInt::VERSION;
+ # Make $n a proper object if it isn't one already (or convert from float)
+ $_[0] = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
+ }
1;
}
+# It you use bigint then call one of the approx/bounds/math functions, you'll
+# end up with full bignum turned on. This seems non-optimal. However, if I
+# don't do this, then you'll get wrong results and end up with it turned on
+# _anyway_. As soon as anyone does something like log($n) where $n is a
+# Math::BigInt, it auto-upgrade and loads up Math::BigFloat.
+#
+# Ideally we'd notice we were causing this, and turn off Math::BigFloat after
+# we were done.
+sub _upgrade_to_float {
+ my($n) = @_;
+ return $n unless defined $Math::BigInt::VERSION || defined $Math::BigFloat::VERSION;
+ do { require Math::BigFloat; Math::BigFloat->import; } if defined $Math::BigInt::VERSION && !defined $Math::BigFloat::VERSION;
+ return Math::BigFloat->new($n);
+}
+
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,
@@ -264,14 +282,13 @@ 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;;
- if (!$use_bigint) {
+ if (!defined $bigint::VERSION) {
foreach my $f2 (@all) {
$all_factors{$f1*$f2} = 1 if ($f1*$f2) < $n;
}
@@ -280,6 +297,7 @@ sub all_factors {
# 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");
+ $product = $product->numify if $product <= ~0;
$all_factors{$product} = 1 if $product < $n;
}
}
@@ -325,22 +343,19 @@ sub euler_phi {
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;
+ # Direct from Euler's product formula. Note division will be exact.
+ #my $totient = $n;
#foreach my $factor (@factors) {
- # $totient *= ($factor - 1);
- # while ($factor_mult{$factor} > 1) {
- # $totient *= $factor;
- # $factor_mult{$factor}--;
- # }
+ # $totient = int($totient/$factor) * ($factor-1);
#}
+ # Alternate way doing multiplications only.
+ my $totient = 1;
+ foreach my $factor (@factors) {
+ $totient *= ($factor - 1);
+ $totient *= $factor for (2 .. $factor_mult{$factor});
+ }
+
$totient;
}
@@ -359,7 +374,7 @@ sub euler_phi {
#
# 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 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
@@ -378,13 +393,9 @@ sub factor {
my($n) = @_;
_validate_positive_integer($n);
- #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;
- }
+ return _XS_factor($n) if $_Config{'xs'} && $n <= $_Config{'maxparam'};
- $n = $n->as_number() if ref($n) =~ /^Math::BigFloat/;
+ $n = $n->as_number() if ref($n) eq 'Math::BigFloat';
$n = $n->numify if $n <= ~0;
Math::Prime::Util::PP::factor($n);
}
@@ -397,6 +408,9 @@ sub prime_count_approx {
return $_prime_count_small[$x] if $x <= $#_prime_count_small;
+ # Turn on high precision FP if they gave us a big number.
+ #do { require Math::BigFloat; Math::BigFloat->import; } if defined $Math::BigInt::VERSION && !defined $Math::BigFloat::VERSION;
+
# Method 10^10 %error 10^19 %error
# ----------------- ------------ ------------
# average bounds .01% .0002%
@@ -419,10 +433,7 @@ sub prime_count_lower {
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";
- }
+ $x = _upgrade_to_float($x) if defined $Math::BigInt::VERSION && ref($x) ne 'Math::BigFloat';
my $flogx = log($x);
@@ -457,10 +468,7 @@ sub prime_count_upper {
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";
- }
+ $x = _upgrade_to_float($x) if defined $Math::BigInt::VERSION && ref($x) ne 'Math::BigFloat';
# Chebyshev: 1.25506*x/logx x >= 17
# Rosser & Schoenfeld: x/(logx-3/2) x >= 67
@@ -512,10 +520,7 @@ sub nth_prime_approx {
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";
- }
+ $n = _upgrade_to_float($n) if defined $Math::BigInt::VERSION && ref($n) ne 'Math::BigFloat';
my $flogn = log($n);
my $flog2n = log($flogn);
@@ -552,7 +557,7 @@ sub nth_prime_approx {
elsif ($n < 200000000) { $approx += 0.0 * $order; }
else { $approx += -0.010 * $order; }
- if ( ($approx >= ~0) && (ref($n) !~ /^Math::Big/) ) {
+ if ( ($approx >= ~0) && (ref($approx) ne 'Math::BigFloat') ) {
return $_Config{'maxprime'} if $n <= $_Config{'maxprimeidx'};
croak "nth_prime_approx($n) overflow";
}
@@ -567,10 +572,7 @@ sub nth_prime_lower {
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";
- }
+ $n = _upgrade_to_float($n) if defined $Math::BigInt::VERSION && ref($n) ne 'Math::BigFloat';
my $flogn = log($n);
my $flog2n = log($flogn); # Note distinction between log_2(n) and log^2(n)
@@ -578,7 +580,7 @@ sub nth_prime_lower {
# 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/) ) {
+ if ( ($lower >= ~0) && (ref($lower) ne 'Math::BigFloat') ) {
return $_Config{'maxprime'} if $n <= $_Config{'maxprimeidx'};
croak "nth_prime_lower($n) overflow";
}
@@ -593,10 +595,7 @@ sub nth_prime_upper {
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";
- }
+ $n = _upgrade_to_float($n) if defined $Math::BigInt::VERSION && ref($n) ne 'Math::BigFloat';
my $flogn = log($n);
my $flog2n = log($flogn); # Note distinction between log_2(n) and log^2(n)
@@ -612,7 +611,7 @@ sub nth_prime_upper {
$upper = $n * ( $flogn + $flog2n );
}
- if ( ($upper >= ~0) && (ref($n) !~ /^Math::Big/) ) {
+ if ( ($upper >= ~0) && (ref($upper) ne 'Math::BigFloat') ) {
return $_Config{'maxprime'} if $n <= $_Config{'maxprimeidx'};
croak "nth_prime_upper($n) overflow";
}
@@ -630,12 +629,7 @@ sub RiemannR {
my($n) = @_;
croak("Invalid input to ReimannR: x must be > 0") if $n <= 0;
- 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, 1e-30) if defined $Math::BigFloat::VERSION;
return Math::Prime::Util::PP::RiemannR($n) if !$_Config{'xs'};
return _XS_RiemannR($n);
@@ -650,12 +644,7 @@ sub ExponentialIntegral {
my($n) = @_;
croak "Invalid input to ExponentialIntegral: x must be != 0" if $n == 0;
- 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, 1e-30) if defined $Math::BigFloat::VERSION;
return Math::Prime::Util::PP::ExponentialIntegral($n) if !$_Config{'xs'};
return _XS_ExponentialIntegral($n);
}
@@ -698,7 +687,7 @@ Math::Prime::Util - Utilities related to prime numbers, including fast sieves an
=head1 VERSION
-Version 0.09
+Version 0.10
=head1 SYNOPSIS
@@ -1105,7 +1094,7 @@ 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>.
+also returns 0 for input 0, but returns C<euler_phi(-n)> for C<n E<lt> 0>.
diff --git a/lib/Math/Prime/Util/MemFree.pm b/lib/Math/Prime/Util/MemFree.pm
index 7575960..37511bf 100644
--- a/lib/Math/Prime/Util/MemFree.pm
+++ b/lib/Math/Prime/Util/MemFree.pm
@@ -4,7 +4,7 @@ use warnings;
BEGIN {
$Math::Prime::Util::MemFree::AUTHORITY = 'cpan:DANAJ';
- $Math::Prime::Util::MemFree::VERSION = '0.08';
+ $Math::Prime::Util::MemFree::VERSION = '0.10';
}
use base qw( Exporter );
@@ -44,7 +44,7 @@ Math::Prime::Util::MemFree - An auto-free object for Math::Prime::Util
=head1 VERSION
-Version 0.08
+Version 0.10
=head1 SYNOPSIS
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 5043faa..e208a1c 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -5,7 +5,7 @@ use Carp qw/carp croak confess/;
BEGIN {
$Math::Prime::Util::PP::AUTHORITY = 'cpan:DANAJ';
- $Math::Prime::Util::PP::VERSION = '0.09';
+ $Math::Prime::Util::PP::VERSION = '0.10';
}
# The Pure Perl versions of all the Math::Prime::Util routines.
@@ -76,8 +76,7 @@ my @_prevwheel30 = (29,29,1,1,1,1,1,1,7,7,7,7,11,11,13,13,13,13,17,17,19,19,19,1
sub _is_prime7 { # n must not be divisible by 2, 3, or 5
my($n) = @_;
- # TODO: bignum on 32-bit
- return is_prob_prime($n) if (~0 == 18446744073709551615) && ($n > 10_000_000);
+ return is_prob_prime($n) if $n > 10_000_000;
foreach my $i (qw/7 11 13 17 19 23 29/) {
return 2 if $i*$i > $n;
@@ -111,11 +110,17 @@ sub _is_prime7 { # n must not be divisible by 2, 3, or 5
sub is_prime {
my($n) = @_;
croak "Input must be an integer" unless defined $_[0];
+ if ($n > ~0) {
+ croak "Input out of range" if !defined $Math::BigInt::VERSION;
+ $n = Math::BigInt->new($n) unless ref($n) =~ /^Math::Big/;
+ } else {
+ $n = $n->numify if $n < ~0 && ref($n) =~ /^Math::Big/;
+ }
+
return 2 if ($n == 2) || ($n == 3) || ($n == 5); # 2, 3, 5 are prime
return 0 if $n < 7; # everything else below 7 is composite
# multiples of 2,3,5 are composite
return 0 if (($n % 2) == 0) || (($n % 3) == 0) || (($n % 5) == 0);
-
return _is_prime7($n);
}
@@ -145,6 +150,7 @@ sub _sieve_erat_string {
# Prefill with 3 and 5 already marked.
my $whole = int( ($end>>1) / 15);
+ croak "Sieve too large" if $whole > 1_145_324_612; # ~32 GB string
my $sieve = "100010010010110" . "011010010010110" x $whole;
# Make exactly the number of entries requested, never more.
substr($sieve, ($end>>1)+1) = '';
@@ -232,6 +238,27 @@ sub _sieve_segment {
\$sieve;
}
+sub trial_primes {
+ my($low,$high) = @_;
+ if (!defined $high) {
+ $high = $low;
+ $low = 2;
+ }
+ croak "Input must be a positive integer" unless _is_positive_int($low)
+ && _is_positive_int($high);
+
+ return if $low > $high;
+
+ my @primes;
+ $low-- if $low >= 2;
+ my $curprime = next_prime($low);
+ while ($curprime <= $high) {
+ push @primes, $curprime;
+ $curprime = next_prime($curprime);
+ }
+ return \@primes;
+}
+
sub primes {
my $optref = (ref $_[0] eq 'HASH') ? shift : {};
croak "no parameters to primes" unless scalar @_ > 0;
@@ -247,6 +274,15 @@ sub primes {
# Ignore method options in this code
+ $low = $low->numify if ref($low) =~ /^Math::Big/ && $low <= ~0;
+ $high = $high->numify if ref($high) =~ /^Math::Big/ && $high <= ~0;
+ # At some point even the pretty-fast pure perl sieve is going to be a
+ # dog, and we should move to trials. This is typical with a small range
+ # on a large base. More thought on the switchover should be done.
+ return trial_primes($low, $high) if ref($low) =~ /^Math::Big/
+ || ref($high) =~ /^Math::Big/
+ || ($low > 1_000_000_000_000 && ($high-$low) < int($low/1_000_000));
+
push @$sref, 2 if ($low <= 2) && ($high >= 2);
push @$sref, 3 if ($low <= 3) && ($high >= 3);
push @$sref, 5 if ($low <= 5) && ($high >= 5);
@@ -332,6 +368,21 @@ sub prime_count {
croak "Input must be a positive integer" unless _is_positive_int($low)
&& _is_positive_int($high);
+if (0) {
+ if ($low > ~0) {
+ croak "Input out of range" if !defined $Math::BigInt::VERSION;
+ $low = Math::BigInt->new($low) unless ref($low) =~ /^Math::Big/;
+ } else {
+ $low = $low->numify if $low < ~0 && ref($low) =~ /^Math::Big/;
+ }
+ if ($high > ~0) {
+ croak "Input out of range" if !defined $Math::BigInt::VERSION;
+ $high = Math::BigInt->new($high) unless ref($high) =~ /^Math::Big/;
+ } else {
+ $high = $high->numify if $high < ~0 && ref($high) =~ /^Math::Big/;
+ }
+}
+
my $count = 0;
$count++ if ($low <= 2) && ($high >= 2); # Count 2
@@ -341,8 +392,20 @@ sub prime_count {
$high-- if ($high % 2) == 0; # Make high go to odd number.
return $count if $low > $high;
- my $sieveref;
+ if ( ref($low) =~ /^Math::Big/ || ref($high) =~ /^Math::Big/
+ || $high > 16_000_000_000
+ || ($high-$low) < int($low/1_000_000) ) {
+ # Too big to sieve.
+ my $count = 0;
+ my $curprime = next_prime($low-1);
+ while ($curprime <= $high) {
+ $count++;
+ $curprime = next_prime($curprime);
+ }
+ return $count;
+ }
+ my $sieveref;
if ($low == 3) {
$sieveref = _sieve_erat($high);
} else {
@@ -464,6 +527,13 @@ sub miller_rabin {
croak "Input must be a positive integer" unless _is_positive_int($n);
croak "No bases given to miller_rabin" unless @bases;
+ if ($n > ~0) {
+ croak "Input out of range" if !defined $Math::BigInt::VERSION;
+ $n = Math::BigInt->new($n) unless ref($n) =~ /^Math::Big/;
+ } else {
+ $n = $n->numify if $n < ~0 && ref($n) =~ /^Math::Big/;
+ }
+
return 0 if ($n == 0) || ($n == 1);
return 2 if ($n == 2) || ($n == 3);
return 0 if ($n % 2) == 0;
@@ -522,6 +592,13 @@ sub miller_rabin {
sub is_prob_prime {
my($n) = @_;
+ if ($n > ~0) {
+ croak "Input out of range" if !defined $Math::BigInt::VERSION;
+ $n = Math::BigInt->new($n) unless ref($n) =~ /^Math::Big/;
+ } else {
+ $n = $n->numify if $n < ~0 && ref($n) =~ /^Math::Big/;
+ }
+
return 2 if ($n == 2) || ($n == 3) || ($n == 5) || ($n == 7);
return 0 if ($n < 2) || (($n % 2) == 0) || (($n % 3) == 0) || (($n % 5) == 0) || (($n % 7) == 0);
return 2 if $n < 121;
@@ -832,6 +909,8 @@ sub ExponentialIntegral {
croak "Invalid input to ExponentialIntegral: x must be != 0" if $x == 0;
+ $x = new Math::BigFloat "$x" if defined $Math::BigFloat::VERSION && ref($x) ne 'Math::BigFloat';
+
my $val; # The result from one of the four methods
if ($x < -1) {
@@ -979,6 +1058,8 @@ sub RiemannR {
croak "Invalid input to ReimannR: x must be > 0" if $x <= 0;
+ $x = new Math::BigFloat "$x" if defined $Math::BigFloat::VERSION && ref($x) ne 'Math::BigFloat';
+
$y = 1.0-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
my $flogx = log($x);
my $part_term = 1.0;
@@ -1011,7 +1092,7 @@ Math::Prime::Util::PP - Pure Perl version of Math::Prime::Util
=head1 VERSION
-Version 0.09
+Version 0.10
=head1 SYNOPSIS
diff --git a/lib/Math/Prime/Util/PrimeArray.pm b/lib/Math/Prime/Util/PrimeArray.pm
index e39e1d0..31eca43 100644
--- a/lib/Math/Prime/Util/PrimeArray.pm
+++ b/lib/Math/Prime/Util/PrimeArray.pm
@@ -4,7 +4,7 @@ use warnings;
BEGIN {
$Math::Prime::Util::PrimeArray::AUTHORITY = 'cpan:DANAJ';
- $Math::Prime::Util::PrimeArray::VERSION = '0.08';
+ $Math::Prime::Util::PrimeArray::VERSION = '0.10';
}
# parent is cleaner, and in the Perl 5.10.1 / 5.12.0 core, but not earlier.
@@ -109,7 +109,7 @@ Math::Prime::Util::PrimeArray - A tied array for primes
=head1 VERSION
-Version 0.08
+Version 0.10
=head1 SYNOPSIS
diff --git a/t/02-can.t b/t/02-can.t
index c2dd783..9d64513 100644
--- a/t/02-can.t
+++ b/t/02-can.t
@@ -6,6 +6,7 @@ use Math::Prime::Util;
use Test::More tests => 1;
my @functions = qw(
+ prime_get_config
prime_precalc prime_memfree
is_prime is_prob_prime miller_rabin
primes
@@ -13,7 +14,7 @@ my @functions = qw(
prime_count prime_count_lower prime_count_upper prime_count_approx
nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
random_prime random_ndigit_prime
- factor all_factors
+ factor all_factors moebius euler_phi
ExponentialIntegral LogarithmicIntegral RiemannR
);
can_ok( 'Math::Prime::Util', @functions);
diff --git a/t/81-bignum.t b/t/81-bignum.t
new file mode 100644
index 0000000..2ed5b02
--- /dev/null
+++ b/t/81-bignum.t
@@ -0,0 +1,163 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+use bigint; # <--------------- We're testing large numbers here: > 2^64
+
+my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
+use Test::More;
+
+my @primes = qw/
+777777777777777777777767 777777777777777777777787
+877777777777777777777753 877777777777777777777871
+87777777777777777777777795577 890745785790123461234805903467891234681243
+618970019642690137449562111
+/;
+push @primes, 531137992816767098689588206552468627329593117727031923199444138200403559860852242739162502265229285668889329486246501015346579337652707239409519978766587351943831270835393219031728127 if $extra;
+
+my @composites = qw/
+777777777777777777777777 877777777777777777777777
+87777777777777777777777795475 890745785790123461234805903467891234681234
+/;
+
+# pseudoprimes to various small prime bases
+# These must not be themselves prime, as we're going to primality test them.
+my %pseudoprimes = (
+ 75792980677 => [ qw/2/ ],
+ 21652684502221 => [ qw/2 7 37 61 9375/ ],
+ 3825123056546413051 => [ qw/2 3 5 7 11 13 17 19 23 29 31 325 9375/ ],
+ 318665857834031151167461 => [ qw/2 3 5 7 11 13 17 19 23 29 31 37 325 9375/ ],
+ 3317044064679887385961981 => [ qw/2 3 5 7 11 13 17 19 23 29 31 37 73 325 9375/ ],
+ 6003094289670105800312596501 => [ qw/2 3 5 7 11 13 17 19 23 29 31 37 61 325 9375/ ],
+ 59276361075595573263446330101 => [ qw/2 3 5 7 11 13 17 19 23 29 31 37 325 9375/ ],
+ 564132928021909221014087501701 => [ qw/2 3 5 7 11 13 17 19 23 29 31 37 325 9375/ ],
+ #1543267864443420616877677640751301 => [ qw/2 3 5 7 11 13 17 19 23 29 31 37 61 325 9375/ ],
+);
+my $num_pseudoprime_tests = 0;
+foreach my $psrp (keys %pseudoprimes) {
+ push @composites, $psrp;
+ $num_pseudoprime_tests += scalar @{$pseudoprimes{$psrp}};
+}
+
+my %factors = (
+ 1234567890 => [2, 3, 3, 5, 3607, 3803],
+ 190128090927491 => [61, 73, 196291, 217517],
+ 23489223467134234890234680 => [2, 2, 2, 5, 4073, 4283, 33662485846146713],
+ #7674353466844111807691499613711 => [11783, 12239, 18869, 22277, 37861, 55163, 60617],
+);
+
+my %allfactors = (
+ #7674353466844111807691499613711 => [11783,12239,18869,22277,37861,55163,60617,144212137,222333427,230937691,262489891,272648203,420344713,446116163,463380779,649985629,675139957,714250111,714399209,741891463,843429497,1040870647,1143782173,1228866151,1350364909,2088526343,2295020237,3343815571,2721138813053,3212613775949,4952921753279,5144598942407,5460015718957,7955174113331,8417765879647,8741707108529,8743531918951,9938129763151,10322733613783,12264578833601,12739215848633,134771853 [...]
+ 23489223467134234890234680 => [2,4,5,8,10,20,40,4073,4283,8146,8566,16292,17132,20365,21415,32584,34264,40730,42830,81460,85660,162920,171320,17444659,34889318,69778636,87223295,139557272,174446590,348893180,697786360,33662485846146713,67324971692293426,134649943384586852,168312429230733565,269299886769173704,336624858461467130,673249716922934260,1346499433845868520,137107304851355562049,144176426879046371779,274214609702711124098,288352853758092743558,548429219405422248196,57670570751 [...]
+);
+
+plan tests => 0 +
+ 2*(@primes + @composites) +
+ 1 +
+ 2 +
+ 1 +
+ $num_pseudoprime_tests +
+ 5 + # PC lower, upper, approx
+ scalar(keys %factors) +
+ scalar(keys %allfactors) +
+ 2 + # moebius, euler_phi
+ 0;
+
+use Math::Prime::Util qw/
+ 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
+/;
+
+ *primes = \&Math::Prime::Util::PP::primes;
+
+ *prime_count = \&Math::Prime::Util::PP::prime_count;
+ #*nth_prime = \&Math::Prime::Util::PP::nth_prime;
+
+ *is_prime = \&Math::Prime::Util::PP::is_prime;
+ *next_prime = \&Math::Prime::Util::PP::next_prime;
+ *prev_prime = \&Math::Prime::Util::PP::prev_prime;
+
+ *miller_rabin = \&Math::Prime::Util::PP::miller_rabin;
+ *is_prob_prime = \&Math::Prime::Util::PP::is_prob_prime;
+
+###############################################################################
+
+foreach my $n (@primes) {
+ ok( is_prime($n), "$n is prime" );
+ ok( is_prob_prime($n), "$n is probably prime");
+}
+foreach my $n (@composites) {
+ ok( !is_prime($n), "$n is not prime" );
+ ok( !is_prob_prime($n), "$n is not probably prime");
+}
+
+###############################################################################
+
+is_deeply( primes(2**66, 2**66+100), [73786976294838206473,73786976294838206549], "primes( 2^66, 2^66 + 100 )" );
+
+###############################################################################
+
+is( next_prime(777777777777777777777777), 777777777777777777777787, "next_prime(777777777777777777777777)");
+is( prev_prime(777777777777777777777777), 777777777777777777777767, "prev_prime(777777777777777777777777)");
+
+###############################################################################
+
+# Testing prime_count only on a small range -- it would take a very long time
+# otherwise.
+
+is( prime_count(877777777777777777777752, 877777777777777777777872), 2, "prime_count(87..7752, 87..7872)");
+
+###############################################################################
+
+# Testing nth_prime would be far too time consuming.
+
+###############################################################################
+
+while (my($psrp, $baseref) = each (%pseudoprimes)) {
+ foreach my $base (@$baseref) {
+ ok( miller_rabin($psrp, $base), "$psrp is a strong pseudoprime to base $base" );
+ }
+}
+
+###############################################################################
+
+{
+ # See: http://www.mersenneforum.org/showpost.php?p=206983&postcount=25
+ my $n = 31415926535897932384626433;
+ cmp_ok( prime_count_lower($n), '<=', 544551456594153032339707, "PC lower (high)" );
+ cmp_ok( prime_count_lower($n), '>=', 544503356940764609324440, "PC lower (low)" );
+ cmp_ok( prime_count_upper($n), '>=', 544551456620339227350566, "PC upper (low)" );
+ cmp_ok( prime_count_upper($n), '<=', 544613583498498996743730, "PC upper (high)" );
+ # TODO: Need to improve accuracy for this
+ ok( abs(prime_count_approx($n) - 544551456607147153724423) < 50_000_000, "PC approx" );
+}
+
+###############################################################################
+
+while (my($n, $factors) = each(%factors)) {
+ is_deeply( [factor($n)], $factors, "factor($n)" );
+}
+while (my($n, $allfactors) = each(%allfactors)) {
+ is_deeply( [all_factors($n)], $allfactors, "all_factors($n)" );
+}
+
+###############################################################################
+
+is( moebius(618970019642690137449562110), -1, "moebius(618970019642690137449562110)" );
+is( euler_phi(618970019642690137449562110), 145857122964987051805507584, "euler_phi(618970019642690137449562110)" );
+
+# ExponentialIntegral
+# LogarithmicIntegral
+# RiemannR
+
+###############################################################################
--
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