[libmath-prime-util-perl] 35/59: Updates for portable bigints and random primes
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:58 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 b72adcd241d0214d2053c0217881b632d8312fcf
Author: Dana Jacobsen <dana at acm.org>
Date: Sat Jul 7 22:58:39 2012 -0600
Updates for portable bigints and random primes
---
Changes | 44 +++++++++++++++-------
Makefile.PL | 6 +--
lib/Math/Prime/Util.pm | 96 +++++++++++++++++++++++++++++++++--------------
lib/Math/Prime/Util/PP.pm | 16 +++++---
t/81-bignum.t | 69 +++++++++++++++++++++-------------
5 files changed, 153 insertions(+), 78 deletions(-)
diff --git a/Changes b/Changes
index 55ccac3..fc27e1c 100644
--- a/Changes
+++ b/Changes
@@ -1,20 +1,36 @@
Revision history for Perl extension Math::Prime::Util.
0.10 7 July 2012
- - 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.
- - many more bigint/bignum changes
- - factor always returns sorted results
- - added BPSW primality test for large is_prob_prime
- - miller_rabin() deprecated. Use is_strong_pseudoprime instead.
- - Add '-bigint' tag to turn off bigint support for some functions. In
- conjunction, add wrappers for all remaining functions so we support
- bigints everywhere.
+ - Add:
+ prime_get_config to get configuration options
+ is_strong_pseudoprime better name for miller_rabin
+ is_strong_lucas_pseudoprime strong lucas-selfridge psp test
+ random_nbit_prime for n-bit primes
+ random_maurer_prime provable n-bit primes
+ moebius Mo:bius function
+ euler_phi Euler's phi aka totient
+
+ - full bigint support for everything. Use '-nobigint' as an import to
+ shortcut straight to XS for better speed on some of the very fast functions.
+ This involved moving a lot of functions into Util.pm.
+
+ - added BPSW primality test for large (>2^64) is_prob_prime and is_prime.
+
+ - Add tests for pp and bignum, cleanup of many tests.
+
+ - Minor changes:
+ - Make miller_rabin return 0 if given even number.
+ - The XS miller_rabin code now works with large base > n.
+ - factor always returns sorted results
+ - miller_rabin() deprecated. Use is_strong_pseudoprime instead.
+
+ - We now should support most of the functionality of:
+ Math::Prime::XS (more functions)
+ Math::Prime::FastSieve
+ Math::Prime::TiedArray (much faster)
+ Math::Factor::XS (faster, missing multiplicity)
+ Math::Primality (more portable, fast for native, slow for bigint)
+ Crypt::Primes (more portable, much slower, no fancy options)
0.09 25 June 2012
- Pure Perl code added. Passes all tests. Used only if the XSLoader
diff --git a/Makefile.PL b/Makefile.PL
index b2f5988..fe742a3 100644
--- a/Makefile.PL
+++ b/Makefile.PL
@@ -23,9 +23,9 @@ WriteMakefile1(
'Carp' => 0,
'Tie::Array' => 0,
'base' => 0,
- 'bignum' => 0,
- 'Math::BigInt' => 0,
- 'Math::BigFloat' => 0,
+ 'bignum' => '0.23', # Used for bigint tests
+ 'Math::BigInt' => '1.88',
+ 'Math::BigFloat' => '1.59',
},
MIN_PERL_VERSION => 5.006002,
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 49803d4..5d8408a 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -309,12 +309,14 @@ sub primes {
my $irandf = (defined &::rand) ? sub { return int(::rand()*shift); }
: sub { return int(rand()*shift); };
# TODO: Look at RANDBITS if using system rand
+ my $rand_max_bits = 31;
+ my $rand_max_val = 1 << $rand_max_bits;
my $_rdata = 0;
my $_rbits = 0;
my $get_rand_bit = sub {
if ($_rbits == 0) {
- $_rdata = $irandf->(2147483648);
- $_rbits = 31;
+ $_rdata = $irandf->($rand_max_val);
+ $_rbits = $rand_max_bits;
}
my $r = $_rdata & 1;
$_rdata >>= 1;
@@ -325,19 +327,29 @@ sub primes {
# Returns a uniform number between [0,$range] inclusive.
my $get_rand_range = sub {
my($range) = @_;
- $range = int($range);
+ my $max = int($range) + 1;
my $offset = 0;
- my $nbits = 0;
- # TODO: consider range == 1. We'll get 0 three out of four times.
- while ($range > 0) {
- my $part = $range >> 1;
- $part++ if ($range & 1) && $get_rand_bit->();
+ while ($max > 1) {
+ if ($max <= $rand_max_val) { $offset += $irandf->($max); last; }
+ my $part = $max >> 1;
+ $part++ if ($max & 1) && $get_rand_bit->();
$offset += $part if $get_rand_bit->();
- $range >>= 1;
- $nbits++;
+ $max -= $part;
}
- wantarray ? ($offset, $nbits) : $offset;
+ $offset;
};
+ # The above routine isn't perfect, but it works pretty well. It's repeatedly
+ # partitioning the space into two pieces selected at random. For odd
+ # ranges the two edges are selected with slightly higher priority because
+ # we're approximating 1/r using powers of 2. The error rapidly reduces
+ # as r increases. By calling out to irandf when max is small enough we can
+ # make it basically go away.
+ #
+ # The other implementation choice I can think of is to call irandf a bunch
+ # of times to get a random number R >= r. Let m = int(R/r). If R < m*r
+ # then return R % m. Repeat otherwise. This description isn't quite right
+ # in that we want to generate R with at least as many random bits as r, not
+ # necessarily greater, and m is related to the bits in each.
# Sub to call with low and high already primes and verified range.
@@ -345,6 +357,10 @@ sub primes {
my($low,$high) = @_;
my $prime;
+ # { my $bsize = 100; my @bins; my $counts = 10000000;
+ # for my $c (1..$counts) { $bins[ $get_rand_range->($bsize) ]++; }
+ # for my $b (0..$bsize) {printf("%4d %8.5f%%\n", $b, $bins[$b]/$counts);}
+
# low and high are both primes, and low < high.
if ($high < 30000) {
@@ -367,7 +383,7 @@ sub primes {
# would be fastest to call primes in the range and randomly pick one. I'm
# not implementing it now because it seems like a rare case.
- if ($oddrange <= 2147483648) {
+ if ($oddrange <= $rand_max_val) {
# Our range is small enough we can just call rand once and be happy.
# Generate random numbers in the interval until one is prime.
my $loop_limit = 2000 * 1000; # To protect against broken rand
@@ -383,19 +399,20 @@ sub primes {
# We have an ocean of range, and a teaspoon to hold randomness.
- # The strategy I'm going to take is to randomly select the bottom
- # portion, leaving the top 2^31 bits for us to iterate over.
+ # Since we have an arbitrary range and not a power of two, I don't see how
+ # Fouque's algorithm A1 could be used (where we generate lower bits and
+ # generate random sets of upper). What I'm doing is pulling out 2^31 lower
+ # bits, then randomly select all the uppers. We iterate adding in the lower
+ # bits.
- my $srange = $oddrange >> 31;
- my ($offset, $nbits) = $get_rand_range->($srange);
+ my $srange = $oddrange - $rand_max_val - 1;
+ my $offset = $get_rand_range->($srange);
my $primelow = $low + 2 * $offset;
- my $uppermult = int($oddrange / $oddrange); # 1 in possible bigint
- $uppermult *= 2 for (1 .. $nbits);
# Generate random numbers in the interval until one is prime.
my $loop_limit = 2000 * 1000; # To protect against broken rand
while (1) {
- $prime = $primelow + ($irandf->(2147483648) * $uppermult);
+ $prime = $primelow + ( 2 * $irandf->($rand_max_val) );
die "$prime > $high" if $prime > $high;
croak "Random function broken?" if $loop_limit-- < 0;
if (ref($prime) eq 'Math::BigInt') {
@@ -492,12 +509,12 @@ sub primes {
if ($k > 2*$m) {
my $rbits = 0;
while ($rbits <= $m) {
- my $s = Math::BigFloat->new( $irandf->(2147483648) )->bdiv(2147483648);
+ my $s = Math::BigFloat->new( $irandf->($rand_max_val) )->bdiv($rand_max_val);
my $r = Math::BigFloat->new(2)->bpow($s-1);
$rbits = $k - ($r*$k);
}
}
- # I've seen +0, +1, and +2 here
+ # I've seen +0, +1, and +2 here. Menezes uses +1.
my $q = random_maurer_prime( ($r * $k)->bfloor + 1 );
#warn "B = $B r = $r k = $k q = $q\n";
my $I = Math::BigInt->new(2)->bpow($k-1)->bdiv(2 * $q)->bfloor;
@@ -525,17 +542,23 @@ sub primes {
next unless $b == 1;
#warn "$n passes a^n-1 == 1\n";
- # Crypt::Primes includes this:
- # next if ($q <= $n->copy->bpow(1/3));
- # but nothing else.
+ # We now get to choose between Maurer's original proposal:
+ # check gcd(a^((n-1)/q)-1,n)==1 for each factor q of n-1
+ # thusly:
- # Maurer's paper indicates we need to check gcd(a^((n-1)/q)-1,n)==1
- # for each factor q of n-1.
$b = $a->copy->bmodpow(2*$R, $n);
next unless Math::BigInt::bgcd($b-1, $n) == 1;
#warn "$n passes final gcd\n";
- # Verify with a BPSW test on the result.
+ # Or via a different method, where we check q >= n**1/3 and also do
+ # some tests on x & y from 2R = xq+y. Crypt::Primes does the q test
+ # but doesn't seem to do the x/y and perfect square portions.
+ # next if ($q <= $n->copy->bpow(1/3));
+ # next if ....
+
+ # Finally, verify with a BPSW test on the result. This will either,
+ # 1) save us from accidently outputing a non-prime due to some mistake
+ # 2) make history by finding the first known BPSW pseudo-prime
die "Maurer prime $n failed BPSW" unless is_prob_prime($n);
#warn " and passed BPSW.\n";
@@ -625,6 +648,16 @@ sub euler_phi {
$totient;
}
+# Omega function A001221. Don't export.
+sub omega {
+ my($n) = @_;
+ return 0 if defined $n && $n <= 1;
+ _validate_positive_integer($n);
+ my %factor_mult;
+ my @factors = grep { !$factor_mult{$_}++ } factor($n);
+ return scalar @factors;
+}
+
#############################################################################
# Front ends to functions.
@@ -803,7 +836,7 @@ sub prime_count_approx {
my $result = RiemannR($x) + 0.5;
- $result = Math::BigInt->new($result->bfloor->bstr()) if ref($result) eq 'Math::BigFloat';
+ return Math::BigInt->new($result->bfloor->bstr()) if ref($result) eq 'Math::BigFloat';
return int($result);
}
@@ -894,7 +927,7 @@ sub prime_count_upper {
# Old versions of Math::BigFloat will do the Wrong Thing with this.
#return int( ($x/$flogx) * (1.0 + 1.0/$flogx + $a/($flogx*$flogx)) + 1.0 );
my $result = ($x/$flogx) * (1.0 + 1.0/$flogx + $a/($flogx*$flogx)) + 1.0;
- $result = Math::BigInt->new($result->bfloor->bstr()) if ref($result) eq 'Math::BigFloat';
+ return Math::BigInt->new($result->bfloor->bstr()) if ref($result) eq 'Math::BigFloat';
return int($result);
}
@@ -1211,6 +1244,11 @@ Using the flag:
will turn off bigint support for those functions. Those functions will then
go directly to the XS versions, which will speed up very small inputs a B<lot>.
+Having run these functions on many versions of Perl, if you're using anything
+older than Perl 5.14, I would recommend you upgrade if you want bigint support.
+There are a lot of brittle behaviors on 5.12.4 and earlier.
+
+
=head1 FUNCTIONS
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index adcb29a..525bcb3 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -333,7 +333,7 @@ sub next_prime {
# my $m = $n - $d*30;
# See: int(9999999999999999403 / 30) => 333333333333333312 (off by 1)
my $m = $n % 30;
- my $d = int( ($n - $m) / 30 );
+ my $d = ($n - $m) / 30;
if ($m == 29) { $d++; $m = 1;} else { $m = $_nextwheel30[$m]; }
while (!_is_prime7($d*30+$m)) {
$m = $_nextwheel30[$m];
@@ -845,7 +845,7 @@ sub factor {
while ( ($n >= (31*31)) && !_is_prime7($n) ) {
my @ftry;
my $holf_rounds = 0;
- if ($n <= ~0) {
+ if ($n < $_half_word) {
$holf_rounds = 64*1024;
#warn "trying holf 64k on $n\n";
@ftry = holf_factor($n, $holf_rounds);
@@ -904,9 +904,13 @@ sub prho_factor {
$U = $n->copy->bzero->badd($U);
$V = $n->copy->bzero->badd($V);
for my $i (1 .. $rounds) {
- $U->bmuladd($U, $a); $U->bmod($n);
- $V->bmuladd($V, $a); $V->bmod($n);
- $V->bmuladd($V, $a); $V->bmod($n);
+ # Would use bmuladd here, but old Math::BigInt's barf with scalar $a.
+ #$U->bmuladd($U, $a); $U->bmod($n);
+ #$V->bmuladd($V, $a); $V->bmod($n);
+ #$V->bmuladd($V, $a); $V->bmod($n);
+ $U->bmul($U); $U->badd($a); $U->bmod($n);
+ $V->bmul($V); $V->badd($a); $V->bmod($n);
+ $V->bmul($V); $V->badd($a); $V->bmod($n);
my $f = Math::BigInt::bgcd( ($U > $V) ? $U-$V : $V-$U, $n);
if ($f == $n) {
last if $inloop++; # We've been here before
@@ -981,7 +985,7 @@ sub pbrent_factor {
$Xi = $n->copy->bzero->badd($Xi);
$Xm = $n->copy->bzero->badd($Xm);
for my $i (1 .. $rounds) {
- $Xi->bmuladd($Xi, $a); $Xi->bmod($n);
+ $Xi->bmul($Xi); $Xi->badd($a); $Xi->bmod($n);
my $f = Math::BigInt::bgcd( ($Xi > $Xm) ? $Xi-$Xm : $Xm-$Xi, $n);
if ( ($f != 1) && ($f != $n) ) {
my $f2 = $n->copy->bdiv($f);
diff --git a/t/81-bignum.t b/t/81-bignum.t
index 802d79f..1c318f6 100644
--- a/t/81-bignum.t
+++ b/t/81-bignum.t
@@ -1,9 +1,15 @@
#!/usr/bin/env perl
use strict;
use warnings;
-use bigint; # <--------------- We're testing large numbers here: > 2^64
+
+# If you're not using ancient perl 5.6.2 with super early releases of bigint,
+# then you can define bigint up here and lose all the stupid quotes
+# around every number.
+
my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
+my $broken64 = (18446744073709550592 == ~0);
+
use Test::More;
my @primes = qw/
@@ -12,7 +18,7 @@ my @primes = qw/
87777777777777777777777795577 890745785790123461234805903467891234681243
618970019642690137449562111
/;
-push @primes, 531137992816767098689588206552468627329593117727031923199444138200403559860852242739162502265229285668889329486246501015346579337652707239409519978766587351943831270835393219031728127 if $extra;
+push @primes, "531137992816767098689588206552468627329593117727031923199444138200403559860852242739162502265229285668889329486246501015346579337652707239409519978766587351943831270835393219031728127" if $extra;
my @composites = qw/
777777777777777777777777 877777777777777777777777
@@ -22,15 +28,15 @@ my @composites = qw/
# 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/ ],
+ '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) {
@@ -39,15 +45,15 @@ foreach my $psrp (keys %pseudoprimes) {
}
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],
+ '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 [...]
+ #'7674353466844111807691499613711' => [qw/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 1347 [...]
+ '23489223467134234890234680' => [qw/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 576705 [...]
);
plan tests => 0 +
@@ -62,6 +68,8 @@ plan tests => 0 +
2 + # moebius, euler_phi
0;
+use bigint; # <--------------- We're testing large numbers here: > 2^64
+
use Math::Prime::Util qw/
is_prob_prime
prime_count_lower
@@ -120,8 +128,11 @@ is( prime_count(877777777777777777777752, 877777777777777777777872), 2, "prime_c
###############################################################################
while (my($psrp, $baseref) = each (%pseudoprimes)) {
- foreach my $base (@$baseref) {
- ok( is_strong_pseudoprime($psrp, $base), "$psrp is a strong pseudoprime to base $base" );
+ SKIP: {
+ skip "Your 64-bit Perl is broken, skipping pseudoprime tests for $psrp", scalar @$baseref if $broken64 && $psrp == 3825123056546413051;
+ foreach my $base (@$baseref) {
+ ok( is_strong_pseudoprime($psrp, $base), "$psrp is a strong pseudoprime to base $base" );
+ }
}
}
@@ -140,17 +151,23 @@ while (my($psrp, $baseref) = each (%pseudoprimes)) {
###############################################################################
-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)" );
+SKIP: {
+ skip "Your 64-bit Perl is broken, skipping bignum factoring tests", scalar(keys %factors) + scalar(keys %allfactors) if $broken64;
+ 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)" );
+SKIP: {
+ skip "Your 64-bit Perl is broken, skipping moebius and euler_phi tests", 2 if $broken64;
+ is( moebius(618970019642690137449562110), -1, "moebius(618970019642690137449562110)" );
+ is( euler_phi(618970019642690137449562110), 145857122964987051805507584, "euler_phi(618970019642690137449562110)" );
+}
# ExponentialIntegral
# LogarithmicIntegral
--
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