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 <d...@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 _______________________________________________ Pkg-perl-cvs-commits mailing list Pkg-perl-cvs-commits@lists.alioth.debian.org http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits