This is an automated email from the git hooks/post-receive script. ppm-guest pushed a commit to annotated tag v0.37 in repository libmath-prime-util-perl.
commit a17cb94897aeadca5faf79401af5752aa3e44581 Author: Dana Jacobsen <d...@acm.org> Date: Mon Jan 20 15:19:44 2014 -0800 Load PP and Math::BigInt only when used. More work needed --- Changes | 8 + MANIFEST | 1 + TODO | 9 + XS.xs | 2 + lib/Math/Prime/Util.pm | 1064 ++++--------------------------- lib/Math/Prime/Util/PP.pm | 73 ++- lib/Math/Prime/Util/PrimalityProving.pm | 1 + lib/Math/Prime/Util/RandomPrimes.pm | 1012 +++++++++++++++++++++++++++++ t/13-primecount.t | 4 +- t/23-primality-proofs.t | 1 + t/70-rt-bignum.t | 1 + t/81-bignum.t | 2 +- 12 files changed, 1243 insertions(+), 935 deletions(-) diff --git a/Changes b/Changes index f8f8deb..d81e0b2 100644 --- a/Changes +++ b/Changes @@ -7,6 +7,14 @@ Revision history for Perl module Math::Prime::Util - Simplified primes(). No longer takes an optional hashref as first arg, which was awkward and never documented. + - Dynamically loads the PP code and Math::BigInt only when needed. This + removes a lot of bloat for the usual cases: + 2.0 MB perl -E 'say 1' + 4.2 MB Math::Prime::XS + 4.6 MB MPU 0.37 + 9.6 MB MPU 0.36 + 10.0 MB MPU 0.35 + 0.36 2014-01-13 [API Changes] diff --git a/MANIFEST b/MANIFEST index b4757c8..d1a8f64 100644 --- a/MANIFEST +++ b/MANIFEST @@ -9,6 +9,7 @@ lib/Math/Prime/Util/ZetaBigFloat.pm lib/Math/Prime/Util/ECAffinePoint.pm lib/Math/Prime/Util/ECProjectivePoint.pm lib/Math/Prime/Util/PrimalityProving.pm +lib/Math/Prime/Util/RandomPrimes.pm LICENSE Makefile.PL MANIFEST diff --git a/TODO b/TODO index 71523b8..0c41e69 100644 --- a/TODO +++ b/TODO @@ -78,3 +78,12 @@ - Ensure a fast path for Math::GMP from MPU -> MPU:GMP -> GMP, and back. - znlog bignum tests, znlog better implementation + +- PP cleanup: + 1) Put front ends for all functions in a PPFE module. + PPFE uses PP. + PPFE does validation. + PP does no validation. + 2) No-XS segment will require PPFE and do aliases. + 3) XS code will call PP only. We've already validated input. + 4) Move some more of the functions to XS. diff --git a/XS.xs b/XS.xs index 9ae27f8..cb234f6 100644 --- a/XS.xs +++ b/XS.xs @@ -176,6 +176,8 @@ static int _vcallsubn(pTHX_ I32 flags, I32 stashflags, const char* name, int nar GV ** gvp = (GV**)hv_fetch(MY_CXT.MPUGMP,name,namelen,0); if (gvp) gv = *gvp; } + if (!gv && (stashflags & VCALL_PP)) + perl_require_pv("Math/Prime/Util/PP.pm"); if (!gv) { GV ** gvp = (GV**)hv_fetch(stashflags & VCALL_PP? MY_CXT.MPUPP : MY_CXT.MPUroot, name,namelen,0); if (gvp) gv = *gvp; diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm index 53ef540..d2f5ec6 100644 --- a/lib/Math/Prime/Util.pm +++ b/lib/Math/Prime/Util.pm @@ -5,17 +5,9 @@ use Carp qw/croak confess carp/; BEGIN { $Math::Prime::Util::AUTHORITY = 'cpan:DANAJ'; - $Math::Prime::Util::VERSION = '0.36'; + $Math::Prime::Util::VERSION = '0.37'; } -BEGIN { - # If they have used Math::BigInt already, make sure we don't change the - # back end. If they have not, try to get one of the fast ones. - do { require Math::BigInt; Math::BigInt->import(try=>"GMP,Pari"); } - unless defined $Math::BigInt::VERSION; -} - - # parent is cleaner, and in the Perl 5.10.1 / 5.12.0 core, but not earlier. # use parent qw( Exporter ); use base qw( Exporter ); @@ -81,11 +73,9 @@ BEGIN { use constant MPU_MAXDIGITS => MPU_32BIT ? 10 : 20; use constant MPU_MAXPRIME => MPU_32BIT ? 4294967291 : 18446744073709551557; use constant MPU_MAXPRIMEIDX => MPU_32BIT ? 203280221 : 425656284035217743; + use constant MPU_MAXNATIVE => OLD_PERL_VERSION ? 999999999999999 : ~0; use constant UVPACKLET => MPU_32BIT ? 'L' : 'Q'; - # Load PP code. Nothing exported. - require Math::Prime::Util::PP; Math::Prime::Util::PP->import(); - eval { return 0 if defined $ENV{MPU_NO_XS} && $ENV{MPU_NO_XS} == 1; require XSLoader; @@ -101,6 +91,9 @@ BEGIN { $_Config{'xs'} = 0; $_Config{'maxbits'} = MPU_MAXBITS; + # Load PP code now. Nothing is exported. + require Math::Prime::Util::PP; Math::Prime::Util::PP->import(); + *_validate_num = \&Math::Prime::Util::PP::_validate_num; *is_prime = \&Math::Prime::Util::PP::is_prime; *is_prob_prime = \&Math::Prime::Util::PP::is_prob_prime; @@ -218,7 +211,6 @@ sub prime_set_config { } elsif ($param eq 'irand') { croak "irand must supply a sub" unless (!defined $value) || (ref($value) eq 'CODE'); $_Config{'irand'} = $value; - _clear_randf(); # Force a new randf to be generated } elsif ($param =~ /^(assume[_ ]?)?[ge]?rh$/ || $param =~ /riemann\s*h/) { $_Config{'assume_rh'} = ($value) ? 1 : 0; } elsif ($param eq 'verbose') { @@ -236,13 +228,57 @@ sub prime_set_config { 1; } - sub _bigint_to_int { return (OLD_PERL_VERSION) ? unpack(UVPACKLET,pack(UVPACKLET,$_[0]->bstr)) : int($_[0]->bstr); } +sub _to_bigint { + do { require Math::BigInt; Math::BigInt->import(try=>"GMP,Pari"); } + unless defined $Math::BigInt::VERSION; + return Math::BigInt->new("$_[0]"); +} +sub _reftyped { + my $ref0 = ref($_[0]); + if ($ref0) { + return ($ref0 eq ref($_[1])) ? $_[1] : $ref0->new("$_[1]"); + } + my $strn = "$_[1]"; + return $_[1] if $strn <= ~0; + do { require Math::BigInt; Math::BigInt->import(try=>"GMP,Pari"); } + unless defined $Math::BigInt::VERSION; + return Math::BigInt->new($strn); +} + -*_validate_positive_integer = \&Math::Prime::Util::PP::_validate_positive_integer; +#*_validate_positive_integer = \&Math::Prime::Util::PP::_validate_positive_integer; +sub _validate_positive_integer { + my($n, $min, $max) = @_; + croak "Parameter must be defined" if !defined $n; + if (ref($n) eq 'CODE') { + $_[0] = $_[0]->(); + $n = $_[0]; + } + if (ref($n) eq 'Math::BigInt') { + croak "Parameter '$n' must be a positive integer" + if $n->sign() ne '+' || !$n->is_int(); + $_[0] = _bigint_to_int($_[0]) + if $n <= (OLD_PERL_VERSION ? 562949953421312 : ''.~0); + } else { + my $strn = "$n"; + croak "Parameter '$strn' must be a positive integer" + if $strn =~ tr/0123456789//c && $strn !~ /^\+?\d+$/; + if ($n <= (OLD_PERL_VERSION ? 562949953421312 : ''.~0)) { + $_[0] = $strn if ref($n); + } else { + #$_[0] = Math::BigInt->new($strn) + $_[0] = _to_bigint($strn); + } + } + $_[0]->upgrade(undef) if ref($_[0]) && $_[0]->upgrade(); + croak "Parameter '$_[0]' must be >= $min" if defined $min && $_[0] < $min; + croak "Parameter '$_[0]' must be <= $max" if defined $max && $_[0] > $max; + 1; +} sub _upgrade_to_float { do { require Math::BigFloat; Math::BigFloat->import(); } @@ -250,6 +286,7 @@ sub _upgrade_to_float { return Math::BigFloat->new($_[0]); } +# TODO: Remove these, push more funcs into XS 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, @@ -299,6 +336,7 @@ sub primes { } return $sref; } + require Math::Prime::Util::PP; return Math::Prime::Util::PP::primes($low,$high); } @@ -330,884 +368,89 @@ sub primes { return $sref; } - -# For random primes, there are two good papers that should be examined: -# -# "Fast Generation of Prime Numbers and Secure Public-Key -# Cryptographic Parameters" by Ueli M. Maurer, 1995 -# http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.26.2151 -# related discussions: -# http://www.daimi.au.dk/~ivan/provableprimesproject.pdf -# Handbook of Applied Cryptography by Menezes, et al. -# -# "Close to Uniform Prime Number Generation With Fewer Random Bits" -# by Pierre-Alain Fouque and Mehdi Tibouchi, 2011 -# http://eprint.iacr.org/2011/481 -# -# Some things to note: -# -# 1) Joye and Paillier have patents on their methods. Never use them. -# -# 2) The easy method of next_prime(random number), known as PRIMEINC, is -# fast but gives a terrible distribution. It has a positive bias and -# most importantly the probability for a prime is proportional to its -# gap, which makes a terrible distribution (some numbers in the range -# will be thousands of times more likely than others). -# -# We use: -# TRIVIAL range within native integer size (2^32 or 2^64) -# FTA1 random_nbit_prime with 65+ bits -# INVA1 other ranges with 65+ bit range -# where -# TRIVIAL = monte-carlo method or equivalent, perfect uniformity. -# FTA1 = Fouque/Tibouchi A1, very close to uniform -# INVA1 = inverted FTA1, less uniform but works with arbitrary ranges -# -# The random_maurer_prime function uses Maurer's FastPrime algorithm. -# -# If Math::Prime::Util::GMP is installed, these functions will be many times -# faster than other methods (e.g. Math::Pari monte-carlo or Crypt::Primes). -# -# Timings on x86_64, with Math::BigInt::GMP and Math::Random::ISAAC::XS. -# -# random_nbit_prime random_maurer_prime -# n-bits no GMP w/ MPU::GMP no GMP w/ MPU::GMP -# ---------- -------- ----------- -------- ----------- -# 24-bit 22uS same same same -# 64-bit 94uS same same same -# 128-bit 0.017s 0.0020s 0.098s 0.056s -# 256-bit 0.033s 0.0033s 0.25s 0.15s -# 512-bit 0.066s 0.0093s 0.65s 0.30s -# 1024-bit 0.16s 0.060s 1.3s 0.94s -# 2048-bit 0.83s 0.5s 3.2s 3.1s -# 4096-bit 6.6s 4.0s 23s 12.0s -# -# Writing these entirely in GMP has a problem, which is that we want to use -# a user-supplied rand function, which means a lot of callbacks. One -# possibility is to, if they do not supply a rand function, use the GMP MT -# function with an appropriate seed. -# -# Random timings for 10M calls: -# 1.92 system rand -# 2.62 Math::Random::MT::Auto -# 12.0 Math::Random::Secure w/ISAAC::XS -# 12.6 Bytes::Random::Secure OO w/ISAAC::XS <==== our -# 31.1 Bytes::Random::Secure OO <==== default -# 44.5 Bytes::Random::Secure function w/ISAAC::XS -# 44.8 Math::Random::Secure -# 71.5 Bytes::Random::Secure function -# 1840 Crypt::Random -# -# time perl -E 'sub irand {int(rand(4294967296));} irand() for 1..10000000;' -# time perl -E 'use Math::Random::MT::Auto qw/irand/; irand() for 1..10000000;' -# time perl -E 'use Math::Random::Secure qw/irand/; irand() for 1..10000000;' -# time perl -E 'use Bytes::Random::Secure qw/random_bytes/; sub irand {return unpack("L",random_bytes(4));} irand() for 1..10000000;' -# time perl -E 'use Bytes::Random::Secure; my $rng = Bytes::Random::Secure->new(); sub irand {return $rng->irand;} irand() for 1..10000000;' -# time perl -E 'use Crypt::Random qw/makerandom/; sub irand {makerandom(Size=>32, Uniform=>1, Strength=>0)} irand() for 1..100_000;' -# > haveged daemon running to stop /dev/random blocking -# > Both BRS and CR have more features that this isn't measuring. -# -# To verify distribution: -# perl -Iblib/lib -Iblib/arch -MMath::Prime::Util=:all -E 'my %freq; $n=1000000; $freq{random_nbit_prime(6)}++ for (1..$n); printf("%4d %6.3f%%\n", $_, 100.0*$freq{$_}/$n) for sort {$a<=>$b} keys %freq;' -# perl -Iblib/lib -Iblib/arch -MMath::Prime::Util=:all -E 'my %freq; $n=1000000; $freq{random_prime(1260437,1260733)}++ for (1..$n); printf("%4d %6.3f%%\n", $_, 100.0*$freq{$_}/$n) for sort {$a<=>$b} keys %freq;' - -{ - # These are much faster than straightforward trial division when n is big. - # You'll want to first do a test up to and including 23. - my @_big_gcd; - my $_big_gcd_top = 20046; - my $_big_gcd_use = -1; - sub _make_big_gcds { - return if $_big_gcd_use >= 0; - if ($_HAVE_GMP) { - $_big_gcd_use = 0; - return; - } - if (Math::BigInt->config()->{lib} !~ /^Math::BigInt::(GMP|Pari)/) { - $_big_gcd_use = 0; - return; - } - $_big_gcd_use = 1; - my $p0 = primorial(Math::BigInt->new( 520)); - my $p1 = primorial(Math::BigInt->new(2052)); - my $p2 = primorial(Math::BigInt->new(6028)); - my $p3 = primorial(Math::BigInt->new($_big_gcd_top)); - $_big_gcd[0] = $p0->bdiv(223092870)->bfloor->as_int; - $_big_gcd[1] = $p1->bdiv($p0)->bfloor->as_int; - $_big_gcd[2] = $p2->bdiv($p1)->bfloor->as_int; - $_big_gcd[3] = $p3->bdiv($p2)->bfloor->as_int; - } - - # Returns a function that will get a uniform random number - # between 0 and $max inclusive. $max can be a bigint. - my $_BRS; - my $_RANDF; - my $_RANDF_NBIT; - sub _set_randf { - # First define a function $irandf that returns a 32-bit integer. This - # corresponds to the irand function of many CPAN modules: - # Math::Random::MT - # Math::Random::ISAAC - # Math::Random::Xorshift - # Math::Random::Secure - # (but not Math::Random::MT::Auto which will return 64-bits) - my $irandf = $_Config{'irand'}; - if (!defined $irandf) { # Default irand: BRS nonblocking - require Bytes::Random::Secure; - $_BRS = Bytes::Random::Secure->new(NonBlocking=>1) unless defined $_BRS; - $_RANDF_NBIT = sub { - my($bits) = @_; - return 0 if $bits <= 0; - return ($_BRS->irand() >> (32-$bits)) - if $bits <= 32; - return ((($_BRS->irand() << 32) + $_BRS->irand()) >> (64-$bits)) - if $bits <= 64 && ~0 > 4294967295; - my $bytes = int(($bits+7)/8); - my $n = Math::BigInt->from_hex('0x' . $_BRS->bytes_hex($bytes)); - $n->brsft( 8*$bytes - $bits ) unless ($bits % 8) == 0; - return $n; - }; - $_RANDF = sub { - my($max) = @_; - my $range = $max+1; - my $U; - if (ref($range) eq 'Math::BigInt') { - my $bits = length($range->as_bin) - 2; # bits in range - my $bytes = 1 + int(($bits+7)/8); # extra byte to reduce ave. loops - my $rmax = Math::BigInt->bone->blsft($bytes*8)->bdec(); - my $overflow = $rmax - ($rmax % $range); - do { - $U = Math::BigInt->from_hex('0x' . $_BRS->bytes_hex($bytes)); - } while $U >= $overflow; - } elsif ($range <= 4294967295) { - my $overflow = (OLD_PERL_VERSION) ? 4294967295-(4294967295.0%$range) - : 4294967295-(4294967295 % $range); - do { - $U = $_BRS->irand(); - } while $U >= $overflow; - } else { - croak "randf given max out of bounds: $max" if $range > ~0; - my $overflow = 18446744073709551615 - (18446744073709551615 % $range); - do { - $U = ($_BRS->irand() << 32) + $_BRS->irand(); - } while $U >= $overflow; - } - return $U % $range; - }; - } else { # Custom irand - $_RANDF_NBIT = sub { - my($bits) = @_; - return 0 if $bits <= 0; - return ($irandf->() >> (32-$bits)) - if $bits <= 32; - return ((($irandf->() << 32) + $irandf->()) >> (64-$bits)) - if $bits <= 64 && MPU_64BIT; - my $words = int(($bits+31)/32); - my $n = Math::BigInt->from_hex - ("0x" . join '', map { sprintf("%08X", $irandf->()) } 1 .. $words ); - $n->brsft( 32*$words - $bits ) unless ($bits % 32) == 0; - return $n; - }; - $_RANDF = sub { - my($max) = @_; - return 0 if $max <= 0; - my $range = $max+1; - my $U; - if (ref($range) eq 'Math::BigInt') { - my $zero = $range->copy->bzero; - my $rbits = length($range->as_bin) - 2; # bits in range - my $rwords = int($rbits/32) + (($rbits % 32) ? 1 : 0); - my $rmax = Math::BigInt->bone->blsft($rwords*32)->bdec(); - my $overflow = $rmax - ($rmax % $range); - do { - $U = $range->copy->from_hex - ("0x" . join '', map { sprintf("%08X", $irandf->()) } 1 .. $rwords); - } while $U >= $overflow; - } elsif ($range <= 4294967295) { - my $overflow = 4294967295 - (4294967295 % $range); - do { - $U = $irandf->(); - } while $U >= $overflow; - } else { - croak "randf given max out of bounds: $max" if $range > ~0; - my $overflow = 18446744073709551615 - (18446744073709551615 % $range); - do { - $U = ($irandf->() << 32) + $irandf->(); - } while $U >= $overflow; - } - return $U % $range; - }; - } - } - sub _clear_randf { - undef $_RANDF; - undef $_RANDF_NBIT; - undef $_BRS; - } - sub _get_randf { - _set_randf() unless defined $_RANDF; - return $_RANDF; - } - sub _get_randf_nbit { - _set_randf() unless defined $_RANDF_NBIT; - return $_RANDF_NBIT; - } - - # Sub to call with low and high already primes and verified range. - my $_random_prime = sub { - my($low,$high) = @_; - my $prime; - - _set_randf() unless defined $_RANDF; - - #{ my $bsize = 100; my @bins; my $counts = 10000000; - # for my $c (1..$counts) { $bins[ $irandf->($bsize-1) ]++; } - # for my $b (0..$bsize) {printf("%4d %8.5f%%\n", $b, $bins[$b]/$counts);} } - - # low and high are both odds, and low < high. - - # This is fast for small values, low memory, perfectly uniform, and - # consumes the minimum amount of randomness needed. But it isn't feasible - # with large values. Also note that low must be a prime. - if ($high <= 262144 && $high <= $_XS_MAXVAL) { - my $li = prime_count(2, $low); - my $irange = prime_count($low, $high); - my $rand = $_RANDF->($irange-1); - return nth_prime($li + $rand); - } - - $low-- if $low == 2; # Low of 2 becomes 1 for our program. - # Math::BigInt::GMP's RT 71548 will wreak havoc if we don't do this. - $low = Math::BigInt->new("$low") if ref($high) eq 'Math::BigInt'; - confess "Invalid _random_prime parameters: $low, $high" if ($low % 2) == 0 || ($high % 2) == 0; - - # We're going to look at the odd numbers only. - my $oddrange = (($high - $low) >> 1) + 1; - - croak "Large random primes not supported on old Perl" - if OLD_PERL_VERSION && MPU_64BIT && $oddrange > 4294967295; - - # If $low is large (e.g. >10 digits) and $range is small (say ~10k), it - # 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 the range is reasonably small, generate using simple Monte Carlo - # method (aka the 'trivial' method). Completely uniform. - if ($oddrange < MPU_MAXPARAM) { - my $loop_limit = 2000 * 1000; # To protect against broken rand - if ($low > 11) { - while ($loop_limit-- > 0) { - $prime = $low + 2 * $_RANDF->($oddrange-1); - next if !($prime % 3) || !($prime % 5) || !($prime % 7) || !($prime % 11); - return $prime if is_prob_prime($prime); - } - } else { - while ($loop_limit-- > 0) { - $prime = $low + 2 * $_RANDF->($oddrange-1); - next if $prime > 11 && (!($prime % 3) || !($prime % 5) || !($prime % 7) || !($prime % 11)); - return 2 if $prime == 1; # Remember the special case for 2. - return $prime if is_prob_prime($prime); - } - } - croak "Random function broken?"; - } - - # We have an ocean of range, and a teaspoon to hold randomness. - - # 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). Similarly trying to simply generate - # upper bits is full of ways to trip up and get non-uniform results. - # - # What I'm doing here is: - # - # 1) divide the range into semi-evenly sized partitions, where each part - # is as close to $rand_max_val as we can. - # 2) randomly select one of the partitions. - # 3) iterate choosing random values within the partition. - # - # The downside is that we're skewing a _lot_ farther from uniformity than - # we'd like. Imagine we started at 0 with 1e18 partitions of size 100k - # each. - # Probability of '5' being returned = - # 1.04e-22 = 1e-18 (chose first partition) * 1/9592 (chose '5') - # Probability of '100003' being returned = - # 1.19e-22 = 1e-18 (chose second partition) * 1/8392 (chose '100003') - # Probability of '99999999999999999999977' being returned = - # 5.20e-22 = 1e-18 (chose last partition) * 1/1922 (chose '99...77') - # So the primes in the last partition will show up 5x more often. - # The partitions are selected uniformly, and the primes within are selected - # uniformly, but the number of primes in each bucket is _not_ uniform. - # Their individual probability of being selected is the probability of the - # partition (uniform) times the probability of being selected inside the - # partition (uniform with respect to all other primes in the same - # partition, but each partition is different and skewed). - # - # Partitions are typically much larger than 100k, but with a huge range - # we still see this (e.g. ~3x from 0-10^30, ~10x from 0-10^100). - # - # When selecting n-bit or n-digit primes, this effect is MUCH smaller, as - # the skew becomes approx lg(2^n) / lg(2^(n-1)) which is pretty close to 1. - # - # - # Another idea I'd like to try sometime is: - # pclo = prime_count_lower(low); - # pchi = prime_count_upper(high); - # do { - # $nth = random selection between pclo and pchi - # $prguess = nth_prime_approx($nth); - # } while ($prguess >= low) && ($prguess <= high); - # monte carlo select a prime in $prguess-2**24 to $prguess+2**24 - # which accounts for the prime distribution. - - my($binsize, $nparts); - my $rand_part_size = 1 << (MPU_64BIT ? 32 : 31); - if (ref($oddrange) eq 'Math::BigInt') { - # Go to some trouble here because some systems are wonky, such as - # giving us +a/+b = -r. Also note the quotes for the bigint argument. - # Without that, Math::BigInt::GMP can return garbage. - my($nbins, $rem); - ($nbins, $rem) = $oddrange->copy->bdiv( "$rand_part_size" ); - $nbins++ if $rem > 0; - $nbins = $nbins->as_int(); - ($binsize,$rem) = $oddrange->copy->bdiv($nbins); - $binsize++ if $rem > 0; - $binsize = $binsize->as_int(); - $nparts = $oddrange->copy->bdiv($binsize)->as_int(); - $low = $high->copy->bzero->badd($low) if ref($low) ne 'Math::BigInt'; - } else { - my $nbins = int($oddrange / $rand_part_size); - $nbins++ if $nbins * $rand_part_size != $oddrange; - $binsize = int($oddrange / $nbins); - $binsize++ if $binsize * $nbins != $oddrange; - $nparts = int($oddrange/$binsize); - } - $nparts-- if ($nparts * $binsize) == $oddrange; - - my $rpart = $_RANDF->($nparts); - - my $primelow = $low + 2 * $binsize * $rpart; - my $partsize = ($rpart < $nparts) ? $binsize - : $oddrange - ($nparts * $binsize); - $partsize = _bigint_to_int($partsize) if ref($partsize) eq 'Math::BigInt'; - #warn "range $oddrange = $nparts * $binsize + ", $oddrange - ($nparts * $binsize), "\n"; - #warn " chose part $rpart size $partsize\n"; - #warn " primelow is $low + 2 * $binsize * $rpart = $primelow\n"; - #die "Result could be too large" if ($primelow + 2*($partsize-1)) > $high; - - # Generate random numbers in the interval until one is prime. - my $loop_limit = 2000 * 1000; # To protect against broken rand - - # Simply things for non-bigints. - if (ref($low) ne 'Math::BigInt') { - while ($loop_limit-- > 0) { - my $rand = $_RANDF->($partsize-1); - $prime = $primelow + $rand + $rand; - croak "random prime failure, $prime > $high" if $prime > $high; - if ($prime <= 23) { - $prime = 2 if $prime == 1; # special case for low = 2 - next unless (0,0,1,1,0,1,0,1,0,0,0,1,0,1,0,0,0,1,0,1,0,0,0,1)[$prime]; - return $prime; - } - next if !($prime % 3) || !($prime % 5) || !($prime % 7) || !($prime % 11); - # It looks promising. Check it. - next unless is_prob_prime($prime); - return $prime; - } - croak "Random function broken?"; - } - - # By checking a wheel 30 mod, we can skip anything that would be a multiple - # of 2, 3, or 5, without even having to create the bigint prime. - my @w30 = (1,0,5,4,3,2,1,0,3,2,1,0,1,0,3,2,1,0,1,0,3,2,1,0,5,4,3,2,1,0); - my $primelow30 = $primelow % 30; - $primelow30 = _bigint_to_int($primelow30) if ref($primelow30) eq 'Math::BigInt'; - - # Big GCD's are hugely fast with GMP or Pari, but super slow with Calc. - _make_big_gcds() if $_big_gcd_use < 0; - - while ($loop_limit-- > 0) { - my $rand = $_RANDF->($partsize-1); - # Check wheel-30 mod - my $rand30 = $rand % 30; - next if $w30[($primelow30 + 2*$rand30) % 30] - && ($rand > 3 || $primelow > 5); - # Construct prime - $prime = $primelow + $rand + $rand; - croak "random prime failure, $prime > $high" if $prime > $high; - if ($prime <= 23) { - $prime = 2 if $prime == 1; # special case for low = 2 - next unless (0,0,1,1,0,1,0,1,0,0,0,1,0,1,0,0,0,1,0,1,0,0,0,1)[$prime]; - return $prime; - } - # With GMP, the fastest thing to do is check primality. - if ($_HAVE_GMP) { - next unless Math::Prime::Util::GMP::is_prime($prime); - return $prime; - } - # No MPU:GMP, so primality checking is slow. Skip some composites here. - next unless Math::BigInt::bgcd($prime, 7436429) == 1; - if ($_big_gcd_use && $prime > $_big_gcd_top) { - next unless Math::BigInt::bgcd($prime, $_big_gcd[0]) == 1; - next unless Math::BigInt::bgcd($prime, $_big_gcd[1]) == 1; - next unless Math::BigInt::bgcd($prime, $_big_gcd[2]) == 1; - next unless Math::BigInt::bgcd($prime, $_big_gcd[3]) == 1; - } - # It looks promising. Check it. - next unless is_prob_prime($prime); - return $prime; - } - croak "Random function broken?"; - }; - - # Cache of tight bounds for each digit. Helps performance a lot. - my @_random_ndigit_ranges = (undef, [2,7], [11,97] ); - my @_random_nbit_ranges = (undef, undef, [2,3],[5,7] ); - my %_random_cache_small; - - # For fixed small ranges with XS, e.g. 6-digit, 18-bit - sub _random_xscount_prime { - my($low,$high) = @_; - my($istart, $irange); - my $cachearef = $_random_cache_small{$low,$high}; - if (defined $cachearef) { - ($istart, $irange) = @$cachearef; - } else { - my $beg = ($low <= 2) ? 2 : next_prime($low-1); - my $end = ($high < ~0) ? prev_prime($high + 1) : prev_prime($high); - ($istart, $irange) = ( prime_count(2, $beg), prime_count($beg, $end) ); - $_random_cache_small{$low,$high} = [$istart, $irange]; - } - _set_randf() unless defined $_RANDF; - my $rand = $_RANDF->($irange-1); - return nth_prime($istart + $rand); - } - - sub random_prime { - my $low = (@_ == 2) ? shift : 2; - my $high = shift; +sub random_prime { + my($low,$high) = @_; + if (scalar @_ > 1) { _validate_num($low) || _validate_positive_integer($low); _validate_num($high) || _validate_positive_integer($high); - - # Tighten the range to the nearest prime. - $low = ($low <= 2) ? 2 : next_prime($low-1); - $high = ($high < ~0) ? prev_prime($high + 1) : prev_prime($high); - return $low if ($low == $high) && is_prob_prime($low); - return if $low >= $high; - - # At this point low and high are both primes, and low < high. - return $_random_prime->($low, $high); - } - - sub random_ndigit_prime { - my($digits) = @_; - _validate_num($digits, 1) || _validate_positive_integer($digits, 1); - - return _random_xscount_prime( int(10 ** ($digits-1)), int(10 ** $digits) ) - if $digits <= 6 && int(10**$digits) <= $_XS_MAXVAL; - - my $bigdigits = $digits >= MPU_MAXDIGITS; - if ($bigdigits && $_Config{'nobigint'}) { - _validate_positive_integer($digits, 1, MPU_MAXDIGITS); - # Special case for nobigint and threshold digits - if (!defined $_random_ndigit_ranges[$digits]) { - my $low = int(10 ** ($digits-1)); - my $high = ~0; - $_random_ndigit_ranges[$digits] = [next_prime($low),prev_prime($high)]; - } - } - - if (!defined $_random_ndigit_ranges[$digits]) { - if ($bigdigits) { - my $low = Math::BigInt->new('10')->bpow($digits-1); - my $high = Math::BigInt->new('10')->bpow($digits); - # Just pull the range in to the nearest odd. - $_random_ndigit_ranges[$digits] = [$low+1, $high-1]; - } else { - my $low = int(10 ** ($digits-1)); - my $high = int(10 ** $digits); - # Note: Perl 5.6.2 cannot represent 10**15 as an integer, so things - # will crash all over the place if you try. We can stringify it, but - # will just fail tests later. - $_random_ndigit_ranges[$digits] = [next_prime($low),prev_prime($high)]; - } - } - my ($low, $high) = @{$_random_ndigit_ranges[$digits]}; - return $_random_prime->($low, $high); - } - - my @_random_nbit_m; - my @_random_nbit_lambda; - my @_random_nbit_arange; - - sub random_nbit_prime { - my($bits) = @_; - _validate_num($bits, 2) || _validate_positive_integer($bits, 2); - - _set_randf() unless defined $_RANDF_NBIT; - - # Very small size, use the nth-prime method - if ($bits <= 18 && int(2**$bits) <= $_XS_MAXVAL) { - if ($bits <= 4) { - return (2,3)[$_RANDF_NBIT->(1)] if $bits == 2; - return (5,7)[$_RANDF_NBIT->(1)] if $bits == 3; - return (11,13)[$_RANDF_NBIT->(1)] if $bits == 4; - } - return _random_xscount_prime( 1 << ($bits-1), 1 << $bits ); - } - - croak "Mid-size random primes not supported on broken old Perl" - if OLD_PERL_VERSION && MPU_64BIT && $bits > 49 && $bits <= 64; - - # Fouque and Tibouchi (2011) Algorithm 1 (basic) - # Modified to make sure the nth bit is always set. - # - # Example for random_nbit_prime(512) on 64-bit Perl: - # p: 1aaaaaaaabbbbbbbbbbbbbbbbbbbb1 - # ^^ ^ ^--- Trailing 1 so p is odd - # || +--- 512-63-2 = 447 lower bits selected before loop - # |+--- 63 upper bits selected in loop, repeated until p is prime - # +--- upper bit is 1 so we generate an n-bit prime - # total: 1 + 63 + 447 + 1 = 512 bits - # - # Algorithm 2 is implemented in a previous commit on github. The problem - # is that it doesn't set the nth bit, and making that change requires a - # modification of the algorithm. It was not a lot faster than this A1 - # with the native int trial division. If the irandf function was very - # slow, then A2 would look more promising. - # - if (1 && $bits > 64) { - my $l = (MPU_64BIT && $bits > 79) ? 63 : 31; - $l = 49 if $l == 63 && OLD_PERL_VERSION; # Fix for broken Perl 5.6 - $l = $bits-2 if $bits-2 < $l; - - my $brand = $_RANDF_NBIT->($bits-$l-2); - $brand = Math::BigInt->new("$brand") unless ref($brand) eq 'Math::BigInt'; - my $b = $brand->blsft(1)->binc(); - - # Precalculate some modulii so we can do trial division on native int - # 9699690 = 2*3*5*7*11*13*17*19, so later operations can be native ints - my @premod; - my $bpremod = _bigint_to_int($b->copy->bmod(9699690)); - my $twopremod = _bigint_to_int(Math::BigInt->new(2)->bmodpow($bits-$l-1, 9699690)); - foreach my $zi (0 .. 19-1) { - foreach my $pm (3, 5, 7, 11, 13, 17, 19) { - next if $zi >= $pm || defined $premod[$pm]; - $premod[$pm] = $zi if ( ($twopremod*$zi+$bpremod) % $pm) == 0; - } - } - _make_big_gcds() if $_big_gcd_use < 0; - my $loop_limit = 1_000_000; - while ($loop_limit-- > 0) { - my $a = (1 << $l) + $_RANDF_NBIT->($l); - # $a % s == $premod[s] => $p % s == 0 => p will be composite - next if $a % 3 == $premod[ 3] || $a % 5 == $premod[ 5] - || $a % 7 == $premod[ 7] || $a % 11 == $premod[11] - || $a % 13 == $premod[13] || $a % 17 == $premod[17] - || $a % 19 == $premod[19]; - my $p = Math::BigInt->new("$a")->blsft($bits-$l-1)->badd($b); - #die " $a $b $p" if $a % 11 == $premod[11] && $p % 11 != 0; - #die "!$a $b $p" if $a % 11 != $premod[11] && $p % 11 == 0; - if ($_HAVE_GMP) { - next unless Math::Prime::Util::GMP::is_prime($p); - } else { - next unless Math::BigInt::bgcd($p, 1348781387) == 1; # 23-43 - if ($_big_gcd_use && $p > $_big_gcd_top) { - next unless Math::BigInt::bgcd($p, $_big_gcd[0]) == 1; - next unless Math::BigInt::bgcd($p, $_big_gcd[1]) == 1; - next unless Math::BigInt::bgcd($p, $_big_gcd[2]) == 1; - next unless Math::BigInt::bgcd($p, $_big_gcd[3]) == 1; - } - # We know we don't have GMP and are > 2^64, so skip all the middle. - #next unless is_prob_prime($p); - #next unless Math::Prime::Util::PP::is_strong_pseudoprime($p, 2); - #next unless Math::Prime::Util::PP::is_extra_strong_lucas_pseudoprime($p); - next unless Math::Prime::Util::PP::is_bpsw_prime($p); - } - return $p; - } - croak "Random function broken?"; - } - - # The Trivial method. Great uniformity, and fine for small sizes. It - # gets very slow as the bit size increases, but that is why we have the - # method above for bigints. - if (1) { - - my $loop_limit = 2_000_000; - if ($bits > MPU_MAXBITS) { - my $p = Math::BigInt->bone->blsft($bits-1)->binc(); - while ($loop_limit-- > 0) { - my $n = Math::BigInt->new(''.$_RANDF_NBIT->($bits-2))->blsft(1)->badd($p); - return $n if is_prob_prime($n); - } - } else { - my $p = (1 << ($bits-1)) + 1; - while ($loop_limit-- > 0) { - my $n = $p + ($_RANDF_NBIT->($bits-2) << 1); - return $n if is_prob_prime($n); - } - } - croak "Random function broken?"; - - } else { - - # Send through the generic random_prime function. Decently fast, but - # quite a bit slower than the F&T A1 method above. - if (!defined $_random_nbit_ranges[$bits]) { - if ($bits > MPU_MAXBITS) { - my $low = Math::BigInt->new('2')->bpow($bits-1); - my $high = Math::BigInt->new('2')->bpow($bits); - # Don't pull the range in to primes, just odds - $_random_nbit_ranges[$bits] = [$low+1, $high-1]; - } else { - my $low = 1 << ($bits-1); - my $high = ($bits == MPU_MAXBITS) - ? ~0-1 - : ~0 >> (MPU_MAXBITS - $bits); - $_random_nbit_ranges[$bits] = [next_prime($low-1),prev_prime($high+1)]; - # Example: bits = 7. - # low = 1<<6 = 64. next_prime(64-1) = 67 - # high = ~0 >> (64-7) = 127. prev_prime(127+1) = 127 - } - } - my ($low, $high) = @{$_random_nbit_ranges[$bits]}; - return $_random_prime->($low, $high); - - } - } - - sub random_maurer_prime { - my $k = shift; - _validate_num($k, 2) || _validate_positive_integer($k, 2); - if ($k <= MPU_MAXBITS && !OLD_PERL_VERSION) { - return random_nbit_prime($k); - } - my ($n, $cert) = random_maurer_prime_with_cert($k); - croak "maurer prime $n failed certificate verification!" - unless verify_prime($cert); - return $n; + } else { + ($low,$high) = (2, $low); + _validate_num($high) || _validate_positive_integer($high); } + require Math::Prime::Util::RandomPrimes; + return Math::Prime::Util::RandomPrimes::random_prime($low,$high); +} - sub random_maurer_prime_with_cert { - my($k) = @_; - _validate_num($k, 2) || _validate_positive_integer($k, 2); - # This should never happen. Trap now to prevent infinite loop. - croak "number of bits must not be a bigint" if ref($k) eq 'Math::BigInt'; - - # Results for random_nbit_prime are proven for all native bit sizes. - my $p0 = MPU_MAXBITS; - $p0 = 49 if OLD_PERL_VERSION && MPU_MAXBITS > 49; - - if ($k <= $p0) { - my $n = random_nbit_prime($k); - my ($isp, $cert) = is_provable_prime_with_cert($n); - croak "small nbit prime could not be proven" if $isp != 2; - return ($n, $cert); - } +sub random_ndigit_prime { + my($digits) = @_; + _validate_num($digits, 1) || _validate_positive_integer($digits, 1); + require Math::Prime::Util::RandomPrimes; + return Math::Prime::Util::RandomPrimes::random_ndigit_prime($digits); +} - # Set verbose to 3 to get pretty output like Crypt::Primes - my $verbose = $_Config{'verbose'}; - local $| = 1 if $verbose > 2; - - do { require Math::BigFloat; Math::BigFloat->import(); } - if !defined $Math::BigFloat::VERSION; - - # Ignore Maurer's g and c that controls how much trial division is done. - my $r = Math::BigFloat->new("0.5"); # relative size of the prime q - my $m = 20; # makes sure R is big enough - _set_randf() unless defined $_RANDF; - - # Generate a random prime q of size $r*$k, where $r >= 0.5. Try to - # cleverly select r to match the size of a typical random factor. - if ($k > 2*$m) { - do { - my $s = Math::BigFloat->new($_RANDF->(2147483647))->bdiv(2147483648); - $r = Math::BigFloat->new(2)->bpow($s-1); - } while ($k*$r >= $k-$m); - } +sub random_nbit_prime { + my($bits) = @_; + _validate_num($bits, 2) || _validate_positive_integer($bits, 2); + require Math::Prime::Util::RandomPrimes; + return Math::Prime::Util::RandomPrimes::random_nbit_prime($bits); +} - # I've seen +0, +1, and +2 here. Maurer uses +0. Menezes uses +1. - my ($q, $qcert) = random_maurer_prime_with_cert( ($r * $k)->bfloor->binc ); - $q = Math::BigInt->new("$q") unless ref($q) eq 'Math::BigInt'; - my $I = Math::BigInt->new(2)->bpow($k-2)->bdiv($q)->bfloor->as_int(); - print "r = $r k = $k q = $q I = $I\n" if $verbose && $verbose != 3; - $qcert = ($q < Math::BigInt->new("18446744073709551615")) - ? "" : _strip_proof_header($qcert); - - # Big GCD's are hugely fast with GMP or Pari, but super slow with Calc. - _make_big_gcds() if $_big_gcd_use < 0; - my $ONE = Math::BigInt->bone; - my $TWO = $ONE->copy->binc; - - my $loop_limit = 1_000_000 + $k * 1_000; - while ($loop_limit-- > 0) { - # R is a random number between $I+1 and 2*$I - #my $R = $I + 1 + $_RANDF->( $I - 1 ); - my $R = $I->copy->binc->badd( $_RANDF->($I->copy->bdec) ); - #my $n = 2 * $R * $q + 1; - my $nm1 = $TWO->copy->bmul($R)->bmul($q); - my $n = $nm1->copy->binc; - # We constructed a promising looking $n. Now test it. - print "." if $verbose > 2; - if ($_HAVE_GMP) { - # MPU::GMP::is_prob_prime has fast tests built in. - next unless Math::Prime::Util::GMP::is_prob_prime($n); - } else { - # No GMP, so first do trial divisions, then a SPSP test. - next unless Math::BigInt::bgcd($n, 111546435)->is_one; - if ($_big_gcd_use && $n > $_big_gcd_top) { - next unless Math::BigInt::bgcd($n, $_big_gcd[0])->is_one; - next unless Math::BigInt::bgcd($n, $_big_gcd[1])->is_one; - next unless Math::BigInt::bgcd($n, $_big_gcd[2])->is_one; - next unless Math::BigInt::bgcd($n, $_big_gcd[3])->is_one; - } - print "+" if $verbose > 2; - next unless is_strong_pseudoprime($n, 3); - } - print "*" if $verbose > 2; - - # We could pick a random generator by doing: - # Step 1: pick a random r - # Step 2: compute g = r^((n-1)/q) mod p - # Step 3: if g == 1, goto Step 1. - # Note that n = 2*R*q+1, hence the exponent is 2*R. - - # We could set r = 0.3333 earlier, then use BLS75 theorem 5 here. - # The chain would be shorter, requiring less overall work for - # large inputs. Maurer's paper discusses the idea. - - # Use BLS75 theorem 3. This is easier and possibly faster than - # BLS75 theorem 4 (Pocklington) used by Maurer's paper. - - # Check conditions -- these should be redundant. - my $m = $TWO * $R; - if (! ($q->is_odd && $q > 2 && $m > 0 && - $m * $q + $ONE == $n && $TWO*$q+$ONE > $n->copy->bsqrt()) ) { - carp "Maurer prime failed BLS75 theorem 3 conditions. Retry."; - next; - } - # Find a suitable a. Move on if one isn't found quickly. - foreach my $trya (2, 3, 5, 7, 11, 13) { - my $a = Math::BigInt->new($trya); - # m/2 = R (n-1)/2 = (2*R*q)/2 = R*q - next unless $a->copy->bmodpow($R, $n) != $nm1; - next unless $a->copy->bmodpow($R*$q, $n) == $nm1; - print "($k)" if $verbose > 2; - croak "Maurer prime $n=2*$R*$q+1 failed BPSW" unless is_prob_prime($n); - my $cert = "[MPU - Primality Certificate]\nVersion 1.0\n\n" . - "Proof for:\nN $n\n\n" . - "Type BLS3\nN $n\nQ $q\nA $a\n" . - $qcert; - return ($n, $cert); - } - # Didn't pass the selected a values. Try another R. - } - croak "Failure in random_maurer_prime, could not find a prime\n"; - } # End of random_maurer_prime - - # Gordon's algorithm for generating a strong prime. - sub random_strong_prime { - my($t) = @_; - _validate_num($t, 128) || _validate_positive_integer($t, 128); - croak "Random strong primes must be >= 173 bits on old Perl" - if OLD_PERL_VERSION && MPU_64BIT && $t < 173; - - _set_randf() unless defined $_RANDF; - - my $l = (($t+1) >> 1) - 2; - my $lp = int($t/2) - 20; - my $lpp = $l - 20; - while (1) { - my $qp = random_nbit_prime($lp); - my $qpp = random_nbit_prime($lpp); - $qp = Math::BigInt->new("$qp") unless ref($qp) eq 'Math::BigInt'; - $qpp = Math::BigInt->new("$qpp") unless ref($qpp) eq 'Math::BigInt'; - my ($il, $rem) = Math::BigInt->new(2)->bpow($l-1)->bdec()->bdiv(2*$qpp); - $il++ if $rem > 0; - $il = $il->as_int(); - my $iu = Math::BigInt->new(2)->bpow($l)->bsub(2)->bdiv(2*$qpp)->as_int(); - my $istart = $il + $_RANDF->($iu - $il); - for (my $i = $istart; $i <= $iu; $i++) { # Search for q - my $q = 2 * $i * $qpp + 1; - next unless is_prob_prime($q); - my $pp = $qp->copy->bmodpow($q-2, $q)->bmul(2)->bmul($qp)->bdec(); - my ($jl, $rem) = Math::BigInt->new(2)->bpow($t-1)->bsub($pp)->bdiv(2*$q*$qp); - $jl++ if $rem > 0; - $jl = $jl->as_int(); - my $ju = Math::BigInt->new(2)->bpow($t)->bdec()->bsub($pp)->bdiv(2*$q*$qp)->as_int(); - my $jstart = $jl + $_RANDF->($ju - $jl); - for (my $j = $jstart; $j <= $ju; $j++) { # Search for p - my $p = $pp + 2 * $j * $q * $qp; - return $p if is_prob_prime($p); - } - } - } - } +sub random_maurer_prime { + my($bits) = @_; + _validate_num($bits, 2) || _validate_positive_integer($bits, 2); + require Math::Prime::Util::RandomPrimes; + return Math::Prime::Util::RandomPrimes::random_maurer_prime($bits); +} - sub random_proven_prime { - my $k = shift; - my ($n, $cert) = random_proven_prime_with_cert($k); - croak "maurer prime $n failed certificate verification!" - unless verify_prime($cert); - return $n; - } +sub random_maurer_prime_with_cert { + my($bits) = @_; + _validate_num($bits, 2) || _validate_positive_integer($bits, 2); + require Math::Prime::Util::RandomPrimes; + return Math::Prime::Util::RandomPrimes::random_maurer_prime_with_cert($bits); +} - sub random_proven_prime_with_cert { - my $k = shift; - _validate_num($k, 2) || _validate_positive_integer($k, 2); +sub random_strong_prime { + my($bits) = @_; + _validate_num($bits, 128) || _validate_positive_integer($bits, 128); + require Math::Prime::Util::RandomPrimes; + return Math::Prime::Util::RandomPrimes::random_strong_prime($bits); +} - if ($_HAVE_GMP && $k <= 450) { - my $n = random_nbit_prime($k); - my ($isp, $cert) = is_provable_prime_with_cert($n); - croak "small nbit prime could not be proven" if $isp != 2; - return ($n, $cert); - } - return random_maurer_prime_with_cert($k); - } +sub random_proven_prime { + my($bits) = @_; + _validate_num($bits, 2) || _validate_positive_integer($bits, 2); + require Math::Prime::Util::RandomPrimes; + return Math::Prime::Util::RandomPrimes::random_proven_prime($bits); +} -} # end of the random prime section +sub random_proven_prime_with_cert { + my($bits) = @_; + _validate_num($bits, 2) || _validate_positive_integer($bits, 2); + require Math::Prime::Util::RandomPrimes; + return Math::Prime::Util::RandomPrimes::random_proven_prime_with_cert($bits); +} sub primorial { my($n) = @_; _validate_num($n) || _validate_positive_integer($n); - return Math::BigInt->new(''.Math::Prime::Util::GMP::primorial($n)) - if $_HAVE_GMP && defined &Math::Prime::Util::GMP::primorial; - - my $max = (MPU_32BIT) ? 29 : (OLD_PERL_VERSION) ? 43 : 53; - my $pn = (ref($_[0]) eq 'Math::BigInt') ? $_[0]->copy->bone() - : ($n >= $max) ? Math::BigInt->bone() - : 1; - if (ref($pn) eq 'Math::BigInt') { - my $start = 2; - if ($n >= 97) { - $start = 101; - $pn->bdec->badd(Math::BigInt->new("2305567963945518424753102147331756070")); - } - my @plist = @{primes($start,$n)}; - while (@plist > 2 && $plist[2] < 1625) { - $pn->bmul( Math::BigInt->new(shift(@plist)*shift(@plist)*shift(@plist)) ); - } - while (@plist > 1 && $plist[1] < 65536) { - $pn->bmul( Math::BigInt->new(shift(@plist)*shift(@plist)) ); - } - $pn->bmul($_) for @plist; - } else { - forprimes { $pn *= $_ } $n; + if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::primorial) { + return _reftyped($_[0], Math::Prime::Util::GMP::primorial($n)); } - return $pn; + require Math::Prime::Util::PP; + return Math::Prime::Util::PP::primorial($n); } sub pn_primorial { my($n) = @_; + _validate_num($n) || _validate_positive_integer($n); if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::pn_primorial) { - _validate_num($n) || _validate_positive_integer($n); - return Math::BigInt->new(''.Math::Prime::Util::GMP::pn_primorial($n)) + return _reftyped($_[0], Math::Prime::Util::GMP::pn_primorial($n)); } - return primorial(nth_prime($n)); + require Math::Prime::Util::PP; + return Math::Prime::Util::PP::primorial(nth_prime($n)); } sub consecutive_integer_lcm { @@ -1215,23 +458,11 @@ sub consecutive_integer_lcm { _validate_num($n) || _validate_positive_integer($n); return 0 if $n < 1; - my $max = (MPU_32BIT) ? 22 : (OLD_PERL_VERSION) ? 37 : 46; - if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::consecutive_integer_lcm) { - my $clcm = Math::Prime::Util::GMP::consecutive_integer_lcm($n); - return ($n < $max) ? int($clcm) : Math::BigInt->new("$clcm"); + return _reftyped($_[0],Math::Prime::Util::GMP::consecutive_integer_lcm($n)); } - - my $pn = (ref($_[0]) eq 'Math::BigInt') ? $_[0]->copy->bone() - : ($n >= $max) ? Math::BigInt->bone() - : 1; - forprimes { - my($p_power, $pmin) = ($_, int($n/$_)); - $p_power *= $_ while $p_power <= $pmin; - $pn *= $p_power; - } $n; - - return $pn; + require Math::Prime::Util::PP; + return Math::Prime::Util::PP::consecutive_integer_lcm($n); } # A008683 Moebius function mu(n) @@ -1239,6 +470,7 @@ sub consecutive_integer_lcm { sub _generic_moebius { my($n, $nend) = @_; return 0 if defined $n && $n < 0; + require Math::Prime::Util::PP; _validate_num($n) || _validate_positive_integer($n); return Math::Prime::Util::PP::moebius($n) if !defined $nend; _validate_num($nend) || _validate_positive_integer($nend); @@ -1282,6 +514,7 @@ sub _generic_mertens { sub _generic_euler_phi { my($n, $nend) = @_; return 0 if defined $n && $n < 0; + require Math::Prime::Util::PP; _validate_num($n) || _validate_positive_integer($n); return Math::Prime::Util::PP::euler_phi($n) if !defined $nend; _validate_num($nend) || _validate_positive_integer($nend); @@ -1291,6 +524,7 @@ sub _generic_euler_phi { sub _generic_divisor_sum { my($n) = @_; _validate_num($n) || _validate_positive_integer($n); + require Math::Prime::Util::PP; return Math::Prime::Util::PP::divisor_sum(@_); } @@ -1354,6 +588,7 @@ sub prime_iterator { } elsif ($_HAVE_GMP) { return sub { $p = $p-$p+Math::Prime::Util::GMP::next_prime($p); return $p;}; } else { + require Math::Prime::Util::PP; return sub { $p = Math::Prime::Util::PP::next_prime($p); return $p; } } } @@ -1399,22 +634,12 @@ sub partitions { return 1 if defined $n && $n <= 0; _validate_num($n) || _validate_positive_integer($n); - return Math::BigInt->new(''.Math::Prime::Util::GMP::partitions($n)) - if $_HAVE_GMP && defined &Math::Prime::Util::GMP::partitions; - - my $d = int(sqrt($n+1)); - my @pent = (1, map { (($_*(3*$_+1))>>1, (($_+1)*(3*$_+2))>>1) } 1 .. $d); - my @part = (Math::BigInt->bone); - foreach my $j (scalar @part .. $n) { - my ($psum1, $psum2, $k) = (Math::BigInt->bzero, Math::BigInt->bzero, 1); - foreach my $p (@pent) { - last if $p > $j; - if ((++$k) & 2) { $psum1->badd( $part[ $j - $p ] ); } - else { $psum2->badd( $part[ $j - $p ] ); } - } - $part[$j] = $psum1 - $psum2; + if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::partitions) { + return _reftyped($_[0],Math::Prime::Util::GMP::partitions($n)); } - return $part[$n]; + + require Math::Prime::Util::PP; + return Math::Prime::Util::PP::partitions($n); } sub chebyshev_theta { @@ -1510,11 +735,10 @@ sub _generic_next_prime { _validate_num($n) || _validate_positive_integer($n); if ($_HAVE_GMP) { - my $r = Math::Prime::Util::GMP::next_prime($n); - return (ref($n) eq 'Math::BigInt' || $n >= MPU_MAXPRIME) - ? Math::BigInt->new("$r") : int($r); + return _reftyped($_[0], Math::Prime::Util::GMP::next_prime($n)); } + require Math::Prime::Util::PP; return Math::Prime::Util::PP::next_prime($_[0]); } @@ -1523,11 +747,10 @@ sub _generic_prev_prime { _validate_num($n) || _validate_positive_integer($n); if ($_HAVE_GMP) { - my $r = Math::Prime::Util::GMP::prev_prime($n); - return (ref($n) eq 'Math::BigInt' && $r > MPU_MAXPRIME) - ? Math::BigInt->new("$r") : int($r); + return _reftyped($_[0], Math::Prime::Util::GMP::prev_prime($n)); } + require Math::Prime::Util::PP; return Math::Prime::Util::PP::prev_prime($_[0]); } @@ -1541,6 +764,7 @@ sub _generic_kronecker { return Math::BigInt->new(''.Math::Prime::Util::GMP::kronecker($a,$b)) if $_HAVE_GMP && defined &Math::Prime::Util::GMP::kronecker; + require Math::Prime::Util::PP; return Math::Prime::Util::PP::kronecker(@_); } @@ -1561,6 +785,7 @@ sub _generic_prime_count { && ( (ref($high) eq 'Math::BigInt') || (($high-$low) < int($low/1_000_000)) ); + require Math::Prime::Util::PP; return Math::Prime::Util::PP::prime_count($low,$high); } @@ -1579,6 +804,7 @@ sub _generic_factor { return @factors; } + require Math::Prime::Util::PP; return Math::Prime::Util::PP::factor($n); } @@ -1644,6 +870,7 @@ sub lucas_sequence { return map { ($_ > ''.~0) ? Math::BigInt->new(''.$_) : $_ } Math::Prime::Util::GMP::lucas_sequence($n, $P, $Q, $k); } + require Math::Prime::Util::PP; return map { ($_ <= ''.~0) ? _bigint_to_int($_) : $_ } Math::Prime::Util::PP::lucas_sequence($n, $P, $Q, $k); } @@ -1656,30 +883,13 @@ sub miller_rabin_random { return 1 if $k <= 0; return (is_prob_prime($n) > 0) if $n < 100; return 0 unless $n & 1; + return Math::Prime::Util::GMP::miller_rabin_random($n, $k) if $_HAVE_GMP && defined &Math::Prime::Util::GMP::miller_rabin_random; - # Testing this many bases is silly, but let's pretend they have some - # good reason. A composite n > 9 must have at least n/4 witnesses, - # hence we need to check only floor(3/4)+1 at most. We could improve - # this is $_Config{'assume_rh'} is true, to 1 .. 2(logn)^2. - if ($k >= int(3*$n/4)) { - return is_strong_pseudoprime($n, 2 .. int(3*$n/4)+1+2 ); - } - - my $brange = $n-3; - my $irandf = _get_randf(); - # Do one first before doing batches - return 0 unless is_strong_pseudoprime($n, $irandf->($brange)+2 ); - $k--; - while ($k > 0) { - my $nbases = ($k >= 20) ? 20 : $k; - my @bases = map { $irandf->($brange)+2 } 1..$nbases; - return 0 unless is_strong_pseudoprime($n, @bases); - $k -= $nbases; - } - 1; + require Math::Prime::Util::RandomPrimes; + return Math::Prime::Util::RandomPrimes::miller_rabin_random($n, $k, $seed); } @@ -2084,6 +1294,7 @@ sub RiemannZeta { return _XS_RiemannZeta($n) if !defined $bignum::VERSION && ref($n) ne 'Math::BigFloat' && $n <= $_XS_MAXVAL; + require Math::Prime::Util::PP; return Math::Prime::Util::PP::RiemannZeta($n); } @@ -2093,6 +1304,7 @@ sub RiemannR { return _XS_RiemannR($n) if !defined $bignum::VERSION && ref($n) ne 'Math::BigFloat' && $n <= $_XS_MAXVAL; + require Math::Prime::Util::PP; return Math::Prime::Util::PP::RiemannR($n); } @@ -2105,6 +1317,7 @@ sub ExponentialIntegral { return _XS_ExponentialIntegral($n) if !defined $bignum::VERSION && ref($n) ne 'Math::BigFloat' && $_Config{'xs'}; + require Math::Prime::Util::PP; return Math::Prime::Util::PP::ExponentialIntegral($n); } @@ -2121,6 +1334,7 @@ sub LogarithmicIntegral { return _XS_LogarithmicIntegral($n); } + require Math::Prime::Util::PP; return Math::Prime::Util::PP::LogarithmicIntegral($n); } diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm index db699cc..e72ae6d 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.36'; + $Math::Prime::Util::PP::VERSION = '0.37'; } BEGIN { @@ -72,8 +72,8 @@ sub _is_positive_int { } sub _bigint_to_int { - return (OLD_PERL_VERSION) ? unpack(UVPACKLET,pack(UVPACKLET,$_[0]->bstr)) - : int($_[0]->bstr); + return (OLD_PERL_VERSION) ? unpack(UVPACKLET,pack(UVPACKLET,"$_[0]")) + : int("$_[0]"); } sub _validate_num { @@ -491,6 +491,65 @@ sub prev_prime { #$d*30+$m; } +sub partitions { + my $n = shift; + + my $d = int(sqrt($n+1)); + my @pent = (1, map { (($_*(3*$_+1))>>1, (($_+1)*(3*$_+2))>>1) } 1 .. $d); + my @part = (Math::BigInt->bone); + foreach my $j (scalar @part .. $n) { + my ($psum1, $psum2, $k) = (Math::BigInt->bzero, Math::BigInt->bzero, 1); + foreach my $p (@pent) { + last if $p > $j; + if ((++$k) & 2) { $psum1->badd( $part[ $j - $p ] ); } + else { $psum2->badd( $part[ $j - $p ] ); } + } + $part[$j] = $psum1 - $psum2; + } + return $part[$n]; +} + +sub primorial { + my $n = shift; + + my $max = (MPU_32BIT) ? 29 : (OLD_PERL_VERSION) ? 43 : 53; + my $pn = (ref($_[0]) eq 'Math::BigInt') ? $_[0]->copy->bone() + : ($n >= $max) ? Math::BigInt->bone() + : 1; + if (ref($pn) eq 'Math::BigInt') { + my $start = 2; + if ($n >= 97) { + $start = 101; + $pn->bdec->badd(Math::BigInt->new("2305567963945518424753102147331756070")); + } + my @plist = @{primes($start,$n)}; + while (@plist > 2 && $plist[2] < 1625) { + $pn->bmul( Math::BigInt->new(shift(@plist)*shift(@plist)*shift(@plist)) ); + } + while (@plist > 1 && $plist[1] < 65536) { + $pn->bmul( Math::BigInt->new(shift(@plist)*shift(@plist)) ); + } + $pn->bmul($_) for @plist; + } else { + foreach my $p (@{primes($n)}) { $pn *= $p; } + } + return $pn; +} + +sub consecutive_integer_lcm { + my $n = shift; + + my $max = (MPU_32BIT) ? 22 : (OLD_PERL_VERSION) ? 37 : 46; + my $pn = ref($n) ? ref($n)->new(1) : ($n >= $max) ? Math::BigInt->bone() : 1; + foreach my $p (@{primes($n)}) { + my($p_power, $pmin) = ($p, int($n/$p)); + $p_power *= $p while $p_power <= $pmin; + $pn *= $p_power; + } + $pn = _bigint_to_int($pn) if $pn <= ''.~0; + return $pn; +} + sub jordan_totient { my($k, $n) = @_; _validate_num($k) || _validate_positive_integer($k); @@ -2370,10 +2429,8 @@ sub ecm_factor { # } #} - if (!defined $Math::Prime::Util::ECProjectivePoint::VERSION) { - eval { require Math::Prime::Util::ECProjectivePoint; 1; } - or do { croak "Cannot load Math::Prime::Util::ECProjectivePoint"; }; - } + require Math::Prime::Util::ECProjectivePoint; + require Math::Prime::Util::RandomPrimes; # With multiple curves, it's better to get all the primes at once. # The downside is this can kill memory with a very large B1. @@ -2385,7 +2442,7 @@ sub ecm_factor { $q = $k; } my @b2primes = ($B2 > $B1) ? @{primes($B1+1, $B2)} : (); - my $irandf = Math::Prime::Util::_get_randf(); + my $irandf = Math::Prime::Util::RandomPrimes::get_randf(); foreach my $curve (1 .. $ncurves) { my $sigma = $irandf->($n-1-6) + 6; diff --git a/lib/Math/Prime/Util/PrimalityProving.pm b/lib/Math/Prime/Util/PrimalityProving.pm index 1e2cfd5..8f17af5 100644 --- a/lib/Math/Prime/Util/PrimalityProving.pm +++ b/lib/Math/Prime/Util/PrimalityProving.pm @@ -117,6 +117,7 @@ sub primality_proof_bls75 { return @composite if ($n & 1) == 0; return @composite if is_strong_pseudoprime($n,2,15,325) == 0; + require Math::Prime::Util::PP; $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt'; my $nm1 = $n->copy->bdec; my $ONE = $nm1->copy->bone; diff --git a/lib/Math/Prime/Util/RandomPrimes.pm b/lib/Math/Prime/Util/RandomPrimes.pm new file mode 100644 index 0000000..057ba97 --- /dev/null +++ b/lib/Math/Prime/Util/RandomPrimes.pm @@ -0,0 +1,1012 @@ +package Math::Prime::Util::RandomPrimes; +use strict; +use warnings; +use Carp qw/carp croak confess/; +use Math::Prime::Util qw/ prime_get_config + verify_prime + is_provable_prime_with_cert + primorial prime_count nth_prime + is_prob_prime is_strong_pseudoprime + next_prime prev_prime + /; + +BEGIN { + $Math::Prime::Util::RandomPrimes::AUTHORITY = 'cpan:DANAJ'; + $Math::Prime::Util::RandomPrimes::VERSION = '0.37'; +} + +BEGIN { + do { require Math::BigInt; Math::BigInt->import(try=>"GMP,Pari"); } + unless defined $Math::BigInt::VERSION; + + use constant OLD_PERL_VERSION=> $] < 5.008; + use constant MPU_MAXBITS => (~0 == 4294967295) ? 32 : 64; + use constant MPU_64BIT => MPU_MAXBITS == 64; + use constant MPU_32BIT => MPU_MAXBITS == 32; + use constant MPU_MAXPARAM => MPU_32BIT ? 4294967295 : 18446744073709551615; + use constant MPU_MAXDIGITS => MPU_32BIT ? 10 : 20; + use constant MPU_USE_XS => prime_get_config->{'xs'}; + use constant MPU_USE_GMP => prime_get_config->{'gmp'}; + + *_bigint_to_int = \&Math::Prime::Util::_bigint_to_int; +} + +################################################################################ + +# These are much faster than straightforward trial division when n is big. +# You'll want to first do a test up to and including 23. +my @_big_gcd; +my $_big_gcd_top = 20046; +my $_big_gcd_use = -1; +sub _make_big_gcds { + return if $_big_gcd_use >= 0; + if (prime_get_config->{'gmp'}) { + $_big_gcd_use = 0; + return; + } + if (Math::BigInt->config()->{lib} !~ /^Math::BigInt::(GMP|Pari)/) { + $_big_gcd_use = 0; + return; + } + $_big_gcd_use = 1; + my $p0 = primorial(Math::BigInt->new( 520)); + my $p1 = primorial(Math::BigInt->new(2052)); + my $p2 = primorial(Math::BigInt->new(6028)); + my $p3 = primorial(Math::BigInt->new($_big_gcd_top)); + $_big_gcd[0] = $p0->bdiv(223092870)->bfloor->as_int; + $_big_gcd[1] = $p1->bdiv($p0)->bfloor->as_int; + $_big_gcd[2] = $p2->bdiv($p1)->bfloor->as_int; + $_big_gcd[3] = $p3->bdiv($p2)->bfloor->as_int; +} + +################################################################################ + +# Returns a function that will get a uniform random number +# between 0 and $max inclusive. $max can be a bigint. +my $_IRANDF; +my $_BRS; +my $_RANDF; +my $_RANDF_NBIT; +sub _set_randf { + # First define a function $irandf that returns a 32-bit integer. This + # corresponds to the irand function of many CPAN modules: + # Math::Random::MT + # Math::Random::ISAAC + # Math::Random::Xorshift + # Math::Random::Secure + # (but not Math::Random::MT::Auto which will return 64-bits) + my $irandf = prime_get_config->{'irand'}; + if ( ( defined $_IRANDF && !defined $irandf) || + (!defined $_IRANDF && defined $irandf) || + ( defined $_IRANDF && defined $irandf && $_IRANDF != $irandf) ) { + undef $_RANDF; + undef $_RANDF_NBIT; + $_IRANDF = $irandf; + } + return if defined $_RANDF; + + if (!defined $_IRANDF) { # Default irand: BRS nonblocking + require Bytes::Random::Secure; + $_BRS = Bytes::Random::Secure->new(NonBlocking=>1) unless defined $_BRS; + $_RANDF_NBIT = sub { + my($bits) = int("$_[0]"); + return 0 if $bits <= 0; + return ($_BRS->irand() >> (32-$bits)) + if $bits <= 32; + return ((($_BRS->irand() << 32) + $_BRS->irand()) >> (64-$bits)) + if $bits <= 64 && ~0 > 4294967295; + my $bytes = int(($bits+7)/8); + my $n = Math::BigInt->from_hex('0x' . $_BRS->bytes_hex($bytes)); + $n->brsft( 8*$bytes - $bits ) unless ($bits % 8) == 0; + return $n; + }; + $_RANDF = sub { + my($max) = @_; + my $range = $max+1; + my $U; + if (ref($range) eq 'Math::BigInt') { + my $bits = length($range->as_bin) - 2; # bits in range + my $bytes = 1 + int(($bits+7)/8); # extra byte to reduce ave. loops + my $rmax = Math::BigInt->bone->blsft($bytes*8)->bdec(); + my $overflow = $rmax - ($rmax % $range); + do { + $U = Math::BigInt->from_hex('0x' . $_BRS->bytes_hex($bytes)); + } while $U >= $overflow; + } elsif ($range <= 4294967295) { + my $overflow = (OLD_PERL_VERSION) ? 4294967295-(4294967295.0%$range) + : 4294967295-(4294967295 % $range); + do { + $U = $_BRS->irand(); + } while $U >= $overflow; + } else { + croak "randf given max out of bounds: $max" if $range > ~0; + my $overflow = 18446744073709551615 - (18446744073709551615 % $range); + do { + $U = ($_BRS->irand() << 32) + $_BRS->irand(); + } while $U >= $overflow; + } + return $U % $range; + }; + } else { # Custom irand + $_RANDF_NBIT = sub { + my($bits) = int("$_[0]"); + return 0 if $bits <= 0; + return ($_IRANDF->() >> (32-$bits)) + if $bits <= 32; + return ((($_IRANDF->() << 32) + $_IRANDF->()) >> (64-$bits)) + if $bits <= 64 && MPU_64BIT; + my $words = int(($bits+31)/32); + my $n = Math::BigInt->from_hex + ("0x" . join '', map { sprintf("%08X", $_IRANDF->()) } 1 .. $words ); + $n->brsft( 32*$words - $bits ) unless ($bits % 32) == 0; + return $n; + }; + $_RANDF = sub { + my($max) = @_; + return 0 if $max <= 0; + my $range = $max+1; + my $U; + if (ref($range) eq 'Math::BigInt') { + my $zero = $range->copy->bzero; + my $rbits = length($range->as_bin) - 2; # bits in range + my $rwords = int($rbits/32) + (($rbits % 32) ? 1 : 0); + my $rmax = Math::BigInt->bone->blsft($rwords*32)->bdec(); + my $overflow = $rmax - ($rmax % $range); + do { + $U = $range->copy->from_hex + ("0x" . join '', map { sprintf("%08X", $_IRANDF->()) } 1 .. $rwords); + } while $U >= $overflow; + } elsif ($range <= 4294967295) { + my $overflow = 4294967295 - (4294967295 % $range); + do { + $U = $_IRANDF->(); + } while $U >= $overflow; + } else { + croak "randf given max out of bounds: $max" if $range > ~0; + my $overflow = 18446744073709551615 - (18446744073709551615 % $range); + do { + $U = ($_IRANDF->() << 32) + $_IRANDF->(); + } while $U >= $overflow; + } + return $U % $range; + }; + } +} + +sub get_randf { + _set_randf(); + return $_RANDF; +} +sub get_randf_nbit { + _set_randf(); + return $_RANDF_NBIT; +} + +################################################################################ + + + +# For random primes, there are two good papers that should be examined: +# +# "Fast Generation of Prime Numbers and Secure Public-Key +# Cryptographic Parameters" by Ueli M. Maurer, 1995 +# http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.26.2151 +# related discussions: +# http://www.daimi.au.dk/~ivan/provableprimesproject.pdf +# Handbook of Applied Cryptography by Menezes, et al. +# +# "Close to Uniform Prime Number Generation With Fewer Random Bits" +# by Pierre-Alain Fouque and Mehdi Tibouchi, 2011 +# http://eprint.iacr.org/2011/481 +# +# Some things to note: +# +# 1) Joye and Paillier have patents on their methods. Never use them. +# +# 2) The easy method of next_prime(random number), known as PRIMEINC, is +# fast but gives a terrible distribution. It has a positive bias and +# most importantly the probability for a prime is proportional to its +# gap, which makes a terrible distribution (some numbers in the range +# will be thousands of times more likely than others). +# +# We use: +# TRIVIAL range within native integer size (2^32 or 2^64) +# FTA1 random_nbit_prime with 65+ bits +# INVA1 other ranges with 65+ bit range +# where +# TRIVIAL = monte-carlo method or equivalent, perfect uniformity. +# FTA1 = Fouque/Tibouchi A1, very close to uniform +# INVA1 = inverted FTA1, less uniform but works with arbitrary ranges +# +# The random_maurer_prime function uses Maurer's FastPrime algorithm. +# +# If Math::Prime::Util::GMP is installed, these functions will be many times +# faster than other methods (e.g. Math::Pari monte-carlo or Crypt::Primes). +# +# Timings on x86_64, with Math::BigInt::GMP and Math::Random::ISAAC::XS. +# +# random_nbit_prime random_maurer_prime +# n-bits no GMP w/ MPU::GMP no GMP w/ MPU::GMP +# ---------- -------- ----------- -------- ----------- +# 24-bit 22uS same same same +# 64-bit 94uS same same same +# 128-bit 0.017s 0.0020s 0.098s 0.056s +# 256-bit 0.033s 0.0033s 0.25s 0.15s +# 512-bit 0.066s 0.0093s 0.65s 0.30s +# 1024-bit 0.16s 0.060s 1.3s 0.94s +# 2048-bit 0.83s 0.5s 3.2s 3.1s +# 4096-bit 6.6s 4.0s 23s 12.0s +# +# Writing these entirely in GMP has a problem, which is that we want to use +# a user-supplied rand function, which means a lot of callbacks. One +# possibility is to, if they do not supply a rand function, use the GMP MT +# function with an appropriate seed. +# +# Random timings for 10M calls: +# 1.92 system rand +# 2.62 Math::Random::MT::Auto +# 12.0 Math::Random::Secure w/ISAAC::XS +# 12.6 Bytes::Random::Secure OO w/ISAAC::XS <==== our +# 31.1 Bytes::Random::Secure OO <==== default +# 44.5 Bytes::Random::Secure function w/ISAAC::XS +# 44.8 Math::Random::Secure +# 71.5 Bytes::Random::Secure function +# 1840 Crypt::Random +# +# time perl -E 'sub irand {int(rand(4294967296));} irand() for 1..10000000;' +# time perl -E 'use Math::Random::MT::Auto qw/irand/; irand() for 1..10000000;' +# time perl -E 'use Math::Random::Secure qw/irand/; irand() for 1..10000000;' +# time perl -E 'use Bytes::Random::Secure qw/random_bytes/; sub irand {return unpack("L",random_bytes(4));} irand() for 1..10000000;' +# time perl -E 'use Bytes::Random::Secure; my $rng = Bytes::Random::Secure->new(); sub irand {return $rng->irand;} irand() for 1..10000000;' +# time perl -E 'use Crypt::Random qw/makerandom/; sub irand {makerandom(Size=>32, Uniform=>1, Strength=>0)} irand() for 1..100_000;' +# > haveged daemon running to stop /dev/random blocking +# > Both BRS and CR have more features that this isn't measuring. +# +# To verify distribution: +# perl -Iblib/lib -Iblib/arch -MMath::Prime::Util=:all -E 'my %freq; $n=1000000; $freq{random_nbit_prime(6)}++ for (1..$n); printf("%4d %6.3f%%\n", $_, 100.0*$freq{$_}/$n) for sort {$a<=>$b} keys %freq;' +# perl -Iblib/lib -Iblib/arch -MMath::Prime::Util=:all -E 'my %freq; $n=1000000; $freq{random_prime(1260437,1260733)}++ for (1..$n); printf("%4d %6.3f%%\n", $_, 100.0*$freq{$_}/$n) for sort {$a<=>$b} keys %freq;' + +# Sub to call with low and high already primes and verified range. +my $_random_prime = sub { + my($low,$high) = @_; + my $prime; + + _set_randf(); + + #{ my $bsize = 100; my @bins; my $counts = 10000000; + # for my $c (1..$counts) { $bins[ $_IRANDF->($bsize-1) ]++; } + # for my $b (0..$bsize) {printf("%4d %8.5f%%\n", $b, $bins[$b]/$counts);} } + + # low and high are both odds, and low < high. + + # This is fast for small values, low memory, perfectly uniform, and + # consumes the minimum amount of randomness needed. But it isn't feasible + # with large values. Also note that low must be a prime. + if ($high <= 262144 && MPU_USE_XS) { + my $li = prime_count(2, $low); + my $irange = prime_count($low, $high); + my $rand = $_RANDF->($irange-1); + return nth_prime($li + $rand); + } + + $low-- if $low == 2; # Low of 2 becomes 1 for our program. + # Math::BigInt::GMP's RT 71548 will wreak havoc if we don't do this. + $low = Math::BigInt->new("$low") if ref($high) eq 'Math::BigInt'; + confess "Invalid _random_prime parameters: $low, $high" if ($low % 2) == 0 || ($high % 2) == 0; + + # We're going to look at the odd numbers only. + my $oddrange = (($high - $low) >> 1) + 1; + + croak "Large random primes not supported on old Perl" + if OLD_PERL_VERSION && MPU_64BIT && $oddrange > 4294967295; + + # If $low is large (e.g. >10 digits) and $range is small (say ~10k), it + # 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 the range is reasonably small, generate using simple Monte Carlo + # method (aka the 'trivial' method). Completely uniform. + if ($oddrange < MPU_MAXPARAM) { + my $loop_limit = 2000 * 1000; # To protect against broken rand + if ($low > 11) { + while ($loop_limit-- > 0) { + $prime = $low + 2 * $_RANDF->($oddrange-1); + next if !($prime % 3) || !($prime % 5) || !($prime % 7) || !($prime % 11); + return $prime if is_prob_prime($prime); + } + } else { + while ($loop_limit-- > 0) { + $prime = $low + 2 * $_RANDF->($oddrange-1); + next if $prime > 11 && (!($prime % 3) || !($prime % 5) || !($prime % 7) || !($prime % 11)); + return 2 if $prime == 1; # Remember the special case for 2. + return $prime if is_prob_prime($prime); + } + } + croak "Random function broken?"; + } + + # We have an ocean of range, and a teaspoon to hold randomness. + + # 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). Similarly trying to simply generate + # upper bits is full of ways to trip up and get non-uniform results. + # + # What I'm doing here is: + # + # 1) divide the range into semi-evenly sized partitions, where each part + # is as close to $rand_max_val as we can. + # 2) randomly select one of the partitions. + # 3) iterate choosing random values within the partition. + # + # The downside is that we're skewing a _lot_ farther from uniformity than + # we'd like. Imagine we started at 0 with 1e18 partitions of size 100k + # each. + # Probability of '5' being returned = + # 1.04e-22 = 1e-18 (chose first partition) * 1/9592 (chose '5') + # Probability of '100003' being returned = + # 1.19e-22 = 1e-18 (chose second partition) * 1/8392 (chose '100003') + # Probability of '99999999999999999999977' being returned = + # 5.20e-22 = 1e-18 (chose last partition) * 1/1922 (chose '99...77') + # So the primes in the last partition will show up 5x more often. + # The partitions are selected uniformly, and the primes within are selected + # uniformly, but the number of primes in each bucket is _not_ uniform. + # Their individual probability of being selected is the probability of the + # partition (uniform) times the probability of being selected inside the + # partition (uniform with respect to all other primes in the same + # partition, but each partition is different and skewed). + # + # Partitions are typically much larger than 100k, but with a huge range + # we still see this (e.g. ~3x from 0-10^30, ~10x from 0-10^100). + # + # When selecting n-bit or n-digit primes, this effect is MUCH smaller, as + # the skew becomes approx lg(2^n) / lg(2^(n-1)) which is pretty close to 1. + # + # + # Another idea I'd like to try sometime is: + # pclo = prime_count_lower(low); + # pchi = prime_count_upper(high); + # do { + # $nth = random selection between pclo and pchi + # $prguess = nth_prime_approx($nth); + # } while ($prguess >= low) && ($prguess <= high); + # monte carlo select a prime in $prguess-2**24 to $prguess+2**24 + # which accounts for the prime distribution. + + my($binsize, $nparts); + my $rand_part_size = 1 << (MPU_64BIT ? 32 : 31); + if (ref($oddrange) eq 'Math::BigInt') { + # Go to some trouble here because some systems are wonky, such as + # giving us +a/+b = -r. Also note the quotes for the bigint argument. + # Without that, Math::BigInt::GMP can return garbage. + my($nbins, $rem); + ($nbins, $rem) = $oddrange->copy->bdiv( "$rand_part_size" ); + $nbins++ if $rem > 0; + $nbins = $nbins->as_int(); + ($binsize,$rem) = $oddrange->copy->bdiv($nbins); + $binsize++ if $rem > 0; + $binsize = $binsize->as_int(); + $nparts = $oddrange->copy->bdiv($binsize)->as_int(); + $low = $high->copy->bzero->badd($low) if ref($low) ne 'Math::BigInt'; + } else { + my $nbins = int($oddrange / $rand_part_size); + $nbins++ if $nbins * $rand_part_size != $oddrange; + $binsize = int($oddrange / $nbins); + $binsize++ if $binsize * $nbins != $oddrange; + $nparts = int($oddrange/$binsize); + } + $nparts-- if ($nparts * $binsize) == $oddrange; + + my $rpart = $_RANDF->($nparts); + + my $primelow = $low + 2 * $binsize * $rpart; + my $partsize = ($rpart < $nparts) ? $binsize + : $oddrange - ($nparts * $binsize); + $partsize = _bigint_to_int($partsize) if ref($partsize) eq 'Math::BigInt'; + #warn "range $oddrange = $nparts * $binsize + ", $oddrange - ($nparts * $binsize), "\n"; + #warn " chose part $rpart size $partsize\n"; + #warn " primelow is $low + 2 * $binsize * $rpart = $primelow\n"; + #die "Result could be too large" if ($primelow + 2*($partsize-1)) > $high; + + # Generate random numbers in the interval until one is prime. + my $loop_limit = 2000 * 1000; # To protect against broken rand + + # Simply things for non-bigints. + if (ref($low) ne 'Math::BigInt') { + while ($loop_limit-- > 0) { + my $rand = $_RANDF->($partsize-1); + $prime = $primelow + $rand + $rand; + croak "random prime failure, $prime > $high" if $prime > $high; + if ($prime <= 23) { + $prime = 2 if $prime == 1; # special case for low = 2 + next unless (0,0,1,1,0,1,0,1,0,0,0,1,0,1,0,0,0,1,0,1,0,0,0,1)[$prime]; + return $prime; + } + next if !($prime % 3) || !($prime % 5) || !($prime % 7) || !($prime % 11); + # It looks promising. Check it. + next unless is_prob_prime($prime); + return $prime; + } + croak "Random function broken?"; + } + + # By checking a wheel 30 mod, we can skip anything that would be a multiple + # of 2, 3, or 5, without even having to create the bigint prime. + my @w30 = (1,0,5,4,3,2,1,0,3,2,1,0,1,0,3,2,1,0,1,0,3,2,1,0,5,4,3,2,1,0); + my $primelow30 = $primelow % 30; + $primelow30 = _bigint_to_int($primelow30) if ref($primelow30) eq 'Math::BigInt'; + + # Big GCD's are hugely fast with GMP or Pari, but super slow with Calc. + _make_big_gcds() if $_big_gcd_use < 0; + + while ($loop_limit-- > 0) { + my $rand = $_RANDF->($partsize-1); + # Check wheel-30 mod + my $rand30 = $rand % 30; + next if $w30[($primelow30 + 2*$rand30) % 30] + && ($rand > 3 || $primelow > 5); + # Construct prime + $prime = $primelow + $rand + $rand; + croak "random prime failure, $prime > $high" if $prime > $high; + if ($prime <= 23) { + $prime = 2 if $prime == 1; # special case for low = 2 + next unless (0,0,1,1,0,1,0,1,0,0,0,1,0,1,0,0,0,1,0,1,0,0,0,1)[$prime]; + return $prime; + } + # With GMP, the fastest thing to do is check primality. + if (MPU_USE_GMP) { + next unless Math::Prime::Util::GMP::is_prime($prime); + return $prime; + } + # No MPU:GMP, so primality checking is slow. Skip some composites here. + next unless Math::BigInt::bgcd($prime, 7436429) == 1; + if ($_big_gcd_use && $prime > $_big_gcd_top) { + next unless Math::BigInt::bgcd($prime, $_big_gcd[0]) == 1; + next unless Math::BigInt::bgcd($prime, $_big_gcd[1]) == 1; + next unless Math::BigInt::bgcd($prime, $_big_gcd[2]) == 1; + next unless Math::BigInt::bgcd($prime, $_big_gcd[3]) == 1; + } + # It looks promising. Check it. + next unless is_prob_prime($prime); + return $prime; + } + croak "Random function broken?"; +}; + +# Cache of tight bounds for each digit. Helps performance a lot. +my @_random_ndigit_ranges = (undef, [2,7], [11,97] ); +my @_random_nbit_ranges = (undef, undef, [2,3],[5,7] ); +my %_random_cache_small; + +# For fixed small ranges with XS, e.g. 6-digit, 18-bit +sub _random_xscount_prime { + my($low,$high) = @_; + my($istart, $irange); + my $cachearef = $_random_cache_small{$low,$high}; + if (defined $cachearef) { + ($istart, $irange) = @$cachearef; + } else { + my $beg = ($low <= 2) ? 2 : next_prime($low-1); + my $end = ($high < ~0) ? prev_prime($high + 1) : prev_prime($high); + ($istart, $irange) = ( prime_count(2, $beg), prime_count($beg, $end) ); + $_random_cache_small{$low,$high} = [$istart, $irange]; + } + _set_randf(); + my $rand = $_RANDF->($irange-1); + return nth_prime($istart + $rand); +} + +sub random_prime { + my($low,$high) = @_; + + # Tighten the range to the nearest prime. + $low = ($low <= 2) ? 2 : next_prime($low-1); + # TODO: if high is bigint, we should do high++? + $high = ($high < ~0) ? prev_prime($high + 1) : prev_prime($high); + return $low if ($low == $high) && is_prob_prime($low); + return if $low >= $high; + + # At this point low and high are both primes, and low < high. + return $_random_prime->($low, $high); +} + +sub random_ndigit_prime { + my($digits) = @_; + croak "random_ndigit_prime, digits must be >= 1" unless $digits >= 1; + + return _random_xscount_prime( int(10 ** ($digits-1)), int(10 ** $digits) ) + if $digits <= 6 && MPU_USE_XS; + + my $bigdigits = $digits >= MPU_MAXDIGITS; + if ($bigdigits && prime_get_config->{'nobigint'}) { + croak "random_ndigit_prime with -nobigint, digits out of range" + if $digits > MPU_MAXDIGITS; + # Special case for nobigint and threshold digits + if (!defined $_random_ndigit_ranges[$digits]) { + my $low = int(10 ** ($digits-1)); + my $high = ~0; + $_random_ndigit_ranges[$digits] = [next_prime($low),prev_prime($high)]; + } + } + + if (!defined $_random_ndigit_ranges[$digits]) { + if ($bigdigits) { + my $low = Math::BigInt->new('10')->bpow($digits-1); + my $high = Math::BigInt->new('10')->bpow($digits); + # Just pull the range in to the nearest odd. + $_random_ndigit_ranges[$digits] = [$low+1, $high-1]; + } else { + my $low = int(10 ** ($digits-1)); + my $high = int(10 ** $digits); + # Note: Perl 5.6.2 cannot represent 10**15 as an integer, so things + # will crash all over the place if you try. We can stringify it, but + # will just fail tests later. + $_random_ndigit_ranges[$digits] = [next_prime($low),prev_prime($high)]; + } + } + my ($low, $high) = @{$_random_ndigit_ranges[$digits]}; + return $_random_prime->($low, $high); +} + +my @_random_nbit_m; +my @_random_nbit_lambda; +my @_random_nbit_arange; + +sub random_nbit_prime { + my($bits) = @_; + croak "random_nbit_prime, bits must be >= 2" unless $bits >= 2; + + _set_randf(); + + # Very small size, use the nth-prime method + if ($bits <= 18 && MPU_USE_XS) { + if ($bits <= 4) { + return (2,3)[$_RANDF_NBIT->(1)] if $bits == 2; + return (5,7)[$_RANDF_NBIT->(1)] if $bits == 3; + return (11,13)[$_RANDF_NBIT->(1)] if $bits == 4; + } + return _random_xscount_prime( 1 << ($bits-1), 1 << $bits ); + } + + croak "Mid-size random primes not supported on broken old Perl" + if OLD_PERL_VERSION && MPU_64BIT && $bits > 49 && $bits <= 64; + + # Fouque and Tibouchi (2011) Algorithm 1 (basic) + # Modified to make sure the nth bit is always set. + # + # Example for random_nbit_prime(512) on 64-bit Perl: + # p: 1aaaaaaaabbbbbbbbbbbbbbbbbbbb1 + # ^^ ^ ^--- Trailing 1 so p is odd + # || +--- 512-63-2 = 447 lower bits selected before loop + # |+--- 63 upper bits selected in loop, repeated until p is prime + # +--- upper bit is 1 so we generate an n-bit prime + # total: 1 + 63 + 447 + 1 = 512 bits + # + # Algorithm 2 is implemented in a previous commit on github. The problem + # is that it doesn't set the nth bit, and making that change requires a + # modification of the algorithm. It was not a lot faster than this A1 + # with the native int trial division. If the irandf function was very + # slow, then A2 would look more promising. + # + if (1 && $bits > 64) { + my $l = (MPU_64BIT && $bits > 79) ? 63 : 31; + $l = 49 if $l == 63 && OLD_PERL_VERSION; # Fix for broken Perl 5.6 + $l = $bits-2 if $bits-2 < $l; + + my $brand = $_RANDF_NBIT->($bits-$l-2); + $brand = Math::BigInt->new("$brand") unless ref($brand) eq 'Math::BigInt'; + my $b = $brand->blsft(1)->binc(); + + # Precalculate some modulii so we can do trial division on native int + # 9699690 = 2*3*5*7*11*13*17*19, so later operations can be native ints + my @premod; + my $bpremod = _bigint_to_int($b->copy->bmod(9699690)); + my $twopremod = _bigint_to_int(Math::BigInt->new(2)->bmodpow($bits-$l-1, 9699690)); + foreach my $zi (0 .. 19-1) { + foreach my $pm (3, 5, 7, 11, 13, 17, 19) { + next if $zi >= $pm || defined $premod[$pm]; + $premod[$pm] = $zi if ( ($twopremod*$zi+$bpremod) % $pm) == 0; + } + } + _make_big_gcds() if $_big_gcd_use < 0; + if (!MPU_USE_GMP) { require Math::Prime::Util::PP; } + + my $loop_limit = 1_000_000; + while ($loop_limit-- > 0) { + my $a = (1 << $l) + $_RANDF_NBIT->($l); + # $a % s == $premod[s] => $p % s == 0 => p will be composite + next if $a % 3 == $premod[ 3] || $a % 5 == $premod[ 5] + || $a % 7 == $premod[ 7] || $a % 11 == $premod[11] + || $a % 13 == $premod[13] || $a % 17 == $premod[17] + || $a % 19 == $premod[19]; + my $p = Math::BigInt->new("$a")->blsft($bits-$l-1)->badd($b); + #die " $a $b $p" if $a % 11 == $premod[11] && $p % 11 != 0; + #die "!$a $b $p" if $a % 11 != $premod[11] && $p % 11 == 0; + if (MPU_USE_GMP) { + next unless Math::Prime::Util::GMP::is_prime($p); + } else { + next unless Math::BigInt::bgcd($p, 1348781387) == 1; # 23-43 + if ($_big_gcd_use && $p > $_big_gcd_top) { + next unless Math::BigInt::bgcd($p, $_big_gcd[0]) == 1; + next unless Math::BigInt::bgcd($p, $_big_gcd[1]) == 1; + next unless Math::BigInt::bgcd($p, $_big_gcd[2]) == 1; + next unless Math::BigInt::bgcd($p, $_big_gcd[3]) == 1; + } + # We know we don't have GMP and are > 2^64, so go directly to this. + next unless Math::Prime::Util::PP::is_bpsw_prime($p); + } + return $p; + } + croak "Random function broken?"; + } + + # The Trivial method. Great uniformity, and fine for small sizes. It + # gets very slow as the bit size increases, but that is why we have the + # method above for bigints. + if (1) { + + my $loop_limit = 2_000_000; + if ($bits > MPU_MAXBITS) { + my $p = Math::BigInt->bone->blsft($bits-1)->binc(); + while ($loop_limit-- > 0) { + my $n = Math::BigInt->new(''.$_RANDF_NBIT->($bits-2))->blsft(1)->badd($p); + return $n if is_prob_prime($n); + } + } else { + my $p = (1 << ($bits-1)) + 1; + while ($loop_limit-- > 0) { + my $n = $p + ($_RANDF_NBIT->($bits-2) << 1); + return $n if is_prob_prime($n); + } + } + croak "Random function broken?"; + + } else { + + # Send through the generic random_prime function. Decently fast, but + # quite a bit slower than the F&T A1 method above. + if (!defined $_random_nbit_ranges[$bits]) { + if ($bits > MPU_MAXBITS) { + my $low = Math::BigInt->new('2')->bpow($bits-1); + my $high = Math::BigInt->new('2')->bpow($bits); + # Don't pull the range in to primes, just odds + $_random_nbit_ranges[$bits] = [$low+1, $high-1]; + } else { + my $low = 1 << ($bits-1); + my $high = ($bits == MPU_MAXBITS) + ? ~0-1 + : ~0 >> (MPU_MAXBITS - $bits); + $_random_nbit_ranges[$bits] = [next_prime($low-1),prev_prime($high+1)]; + # Example: bits = 7. + # low = 1<<6 = 64. next_prime(64-1) = 67 + # high = ~0 >> (64-7) = 127. prev_prime(127+1) = 127 + } + } + my ($low, $high) = @{$_random_nbit_ranges[$bits]}; + return $_random_prime->($low, $high); + + } +} + +sub random_maurer_prime { + my $k = shift; + croak "random_maurer_prime, bits must be >= 2" unless $k >= 2; + + return random_nbit_prime($k) if $k <= MPU_MAXBITS && !OLD_PERL_VERSION; + + my ($n, $cert) = random_maurer_prime_with_cert($k); + croak "maurer prime $n failed certificate verification!" + unless verify_prime($cert); + return $n; +} + +sub random_maurer_prime_with_cert { + my $k = shift; + croak "random_maurer_prime, bits must be >= 2" unless $k >= 2; + + # This should never happen. Trap now to prevent infinite loop. + croak "number of bits must not be a bigint" if ref($k) eq 'Math::BigInt'; + + # Results for random_nbit_prime are proven for all native bit sizes. + my $p0 = MPU_MAXBITS; + $p0 = 49 if OLD_PERL_VERSION && MPU_MAXBITS > 49; + + if ($k <= $p0) { + my $n = random_nbit_prime($k); + my ($isp, $cert) = is_provable_prime_with_cert($n); + croak "small nbit prime could not be proven" if $isp != 2; + return ($n, $cert); + } + + # Set verbose to 3 to get pretty output like Crypt::Primes + my $verbose = prime_get_config->{'verbose'}; + local $| = 1 if $verbose > 2; + + do { require Math::BigFloat; Math::BigFloat->import(); } + if !defined $Math::BigFloat::VERSION; + + # Ignore Maurer's g and c that controls how much trial division is done. + my $r = Math::BigFloat->new("0.5"); # relative size of the prime q + my $m = 20; # makes sure R is big enough + _set_randf(); + + # Generate a random prime q of size $r*$k, where $r >= 0.5. Try to + # cleverly select r to match the size of a typical random factor. + if ($k > 2*$m) { + do { + my $s = Math::BigFloat->new($_RANDF->(2147483647))->bdiv(2147483648); + $r = Math::BigFloat->new(2)->bpow($s-1); + } while ($k*$r >= $k-$m); + } + + # I've seen +0, +1, and +2 here. Maurer uses +0. Menezes uses +1. + my ($q, $qcert) = random_maurer_prime_with_cert( ($r * $k)->bfloor->binc ); + $q = Math::BigInt->new("$q") unless ref($q) eq 'Math::BigInt'; + my $I = Math::BigInt->new(2)->bpow($k-2)->bdiv($q)->bfloor->as_int(); + print "r = $r k = $k q = $q I = $I\n" if $verbose && $verbose != 3; + $qcert = ($q < Math::BigInt->new("18446744073709551615")) + ? "" : _strip_proof_header($qcert); + + # Big GCD's are hugely fast with GMP or Pari, but super slow with Calc. + _make_big_gcds() if $_big_gcd_use < 0; + my $ONE = Math::BigInt->bone; + my $TWO = $ONE->copy->binc; + + my $loop_limit = 1_000_000 + $k * 1_000; + while ($loop_limit-- > 0) { + # R is a random number between $I+1 and 2*$I + #my $R = $I + 1 + $_RANDF->( $I - 1 ); + my $R = $I->copy->binc->badd( $_RANDF->($I->copy->bdec) ); + #my $n = 2 * $R * $q + 1; + my $nm1 = $TWO->copy->bmul($R)->bmul($q); + my $n = $nm1->copy->binc; + # We constructed a promising looking $n. Now test it. + print "." if $verbose > 2; + if (MPU_USE_GMP) { + # MPU::GMP::is_prob_prime has fast tests built in. + next unless Math::Prime::Util::GMP::is_prob_prime($n); + } else { + # No GMP, so first do trial divisions, then a SPSP test. + next unless Math::BigInt::bgcd($n, 111546435)->is_one; + if ($_big_gcd_use && $n > $_big_gcd_top) { + next unless Math::BigInt::bgcd($n, $_big_gcd[0])->is_one; + next unless Math::BigInt::bgcd($n, $_big_gcd[1])->is_one; + next unless Math::BigInt::bgcd($n, $_big_gcd[2])->is_one; + next unless Math::BigInt::bgcd($n, $_big_gcd[3])->is_one; + } + print "+" if $verbose > 2; + next unless is_strong_pseudoprime($n, 3); + } + print "*" if $verbose > 2; + + # We could pick a random generator by doing: + # Step 1: pick a random r + # Step 2: compute g = r^((n-1)/q) mod p + # Step 3: if g == 1, goto Step 1. + # Note that n = 2*R*q+1, hence the exponent is 2*R. + + # We could set r = 0.3333 earlier, then use BLS75 theorem 5 here. + # The chain would be shorter, requiring less overall work for + # large inputs. Maurer's paper discusses the idea. + + # Use BLS75 theorem 3. This is easier and possibly faster than + # BLS75 theorem 4 (Pocklington) used by Maurer's paper. + + # Check conditions -- these should be redundant. + my $m = $TWO * $R; + if (! ($q->is_odd && $q > 2 && $m > 0 && + $m * $q + $ONE == $n && $TWO*$q+$ONE > $n->copy->bsqrt()) ) { + carp "Maurer prime failed BLS75 theorem 3 conditions. Retry."; + next; + } + # Find a suitable a. Move on if one isn't found quickly. + foreach my $trya (2, 3, 5, 7, 11, 13) { + my $a = Math::BigInt->new($trya); + # m/2 = R (n-1)/2 = (2*R*q)/2 = R*q + next unless $a->copy->bmodpow($R, $n) != $nm1; + next unless $a->copy->bmodpow($R*$q, $n) == $nm1; + print "($k)" if $verbose > 2; + croak "Maurer prime $n=2*$R*$q+1 failed BPSW" unless is_prob_prime($n); + my $cert = "[MPU - Primality Certificate]\nVersion 1.0\n\n" . + "Proof for:\nN $n\n\n" . + "Type BLS3\nN $n\nQ $q\nA $a\n" . + $qcert; + return ($n, $cert); + } + # Didn't pass the selected a values. Try another R. + } + croak "Failure in random_maurer_prime, could not find a prime\n"; +} # End of random_maurer_prime + +# Gordon's algorithm for generating a strong prime. +sub random_strong_prime { + my $t = shift; + croak "random_strong_prime, bits must be >= 128" unless $t >= 128; + + croak "Random strong primes must be >= 173 bits on old Perl" + if OLD_PERL_VERSION && MPU_64BIT && $t < 173; + + _set_randf(); + + my $l = (($t+1) >> 1) - 2; + my $lp = int($t/2) - 20; + my $lpp = $l - 20; + while (1) { + my $qp = random_nbit_prime($lp); + my $qpp = random_nbit_prime($lpp); + $qp = Math::BigInt->new("$qp") unless ref($qp) eq 'Math::BigInt'; + $qpp = Math::BigInt->new("$qpp") unless ref($qpp) eq 'Math::BigInt'; + my ($il, $rem) = Math::BigInt->new(2)->bpow($l-1)->bdec()->bdiv(2*$qpp); + $il++ if $rem > 0; + $il = $il->as_int(); + my $iu = Math::BigInt->new(2)->bpow($l)->bsub(2)->bdiv(2*$qpp)->as_int(); + my $istart = $il + $_RANDF->($iu - $il); + for (my $i = $istart; $i <= $iu; $i++) { # Search for q + my $q = 2 * $i * $qpp + 1; + next unless is_prob_prime($q); + my $pp = $qp->copy->bmodpow($q-2, $q)->bmul(2)->bmul($qp)->bdec(); + my ($jl, $rem) = Math::BigInt->new(2)->bpow($t-1)->bsub($pp)->bdiv(2*$q*$qp); + $jl++ if $rem > 0; + $jl = $jl->as_int(); + my $ju = Math::BigInt->new(2)->bpow($t)->bdec()->bsub($pp)->bdiv(2*$q*$qp)->as_int(); + my $jstart = $jl + $_RANDF->($ju - $jl); + for (my $j = $jstart; $j <= $ju; $j++) { # Search for p + my $p = $pp + 2 * $j * $q * $qp; + return $p if is_prob_prime($p); + } + } + } +} + +sub random_proven_prime { + my $k = shift; + my ($n, $cert) = random_proven_prime_with_cert($k); + croak "random_proven_prime $n failed certificate verification!" + unless verify_prime($cert); + return $n; +} + +sub random_proven_prime_with_cert { + my $k = shift; + + if (prime_get_config->{'gmp'} && $k <= 450) { + my $n = random_nbit_prime($k); + my ($isp, $cert) = is_provable_prime_with_cert($n); + croak "small nbit prime could not be proven" if $isp != 2; + return ($n, $cert); + } + return random_maurer_prime_with_cert($k); +} + +sub miller_rabin_random { + my($n, $k, $seed) = @_; + + # Testing this many bases is silly, but let's pretend they have some + # good reason. A composite n > 9 must have at least n/4 witnesses, + # hence we need to check only floor(3/4)+1 at most. We could improve + # this is $_Config{'assume_rh'} is true, to 1 .. 2(logn)^2. + if ($k >= int(3*$n/4)) { + return is_strong_pseudoprime($n, 2 .. int(3*$n/4)+1+2 ); + } + + _set_randf(); + + my $brange = $n-3; + # Do one first before doing batches + return 0 unless is_strong_pseudoprime($n, $_RANDF->($brange)+2 ); + $k--; + while ($k > 0) { + my $nbases = ($k >= 20) ? 20 : $k; + my @bases = map { $_RANDF->($brange)+2 } 1..$nbases; + return 0 unless is_strong_pseudoprime($n, @bases); + $k -= $nbases; + } + 1; +} + +1; + +__END__ + + +# ABSTRACT: Generate random primes + +=pod + +=encoding utf8 + +=head1 NAME + +Math::Prime::Util::RandomPrimes - Generate random primes + + +=head1 VERSION + +Version 0.37 + + +=head1 SYNOPSIS + +=head1 DESCRIPTION + +Routines to generate random primes, including constructing proven primes. + + +=head1 RANDOM UTILITY FUNCTIONS + +=head2 get_randf + +Gets a subroutine that will produce random integers between 0 and C<n>, +inclusive. The argument C<n> can be a bigint. + +=head2 get_randf_nbit + +Gets a subroutine that will produce random integers between 0 and C<2^n-1>, +inclusive. + + +=head1 RANDOM PRIME FUNCTIONS + +=head2 random_prime + +Generate a random prime between C<low> and C<high>. If given one argument, +C<low> will be 2. + +=head2 random_ndigit_prime + +Generate a random prime with C<n> digits. C<n> must be at least 1. + +=head2 random_nbit_prime + +Generate a random prime with C<n> bits. C<n> must be at least 2. + +=head2 random_maurer_prime + +Construct a random provable prime of C<n> bits using Maurer's FastPrime +algorithm. C<n> must be at least 2. + +=head2 random_maurer_prime_with_cert + +Construct a random provable prime of C<n> bits using Maurer's FastPrime +algorithm. C<n> must be at least 2. Returns a list of two items: the +prime and the certificate. + +=head2 random_strong_prime + +Construct a random strong prime of C<n> bits. C<n> must be at least 128. + +=head2 random_proven_prime + +Generate or construct a random provable prime of C<n> bits. C<n> must +be at least 2. + +=head2 random_proven_prime_with_cert + +Generate or construct a random provable prime of C<n> bits. C<n> must +be at least 2. Returns a list of two items: the prime and the certificate. + + +=head1 RANDOM PRIMALITY FUNCTIONS + +=head2 miller_rabin_random + +Given a number C<n> and a number of tests to perform C<k>, this does C<k> +Miller-Rabin tests on C<n> using randomly selected bases. The return value +is 1 if all bases are a witness to to C<n>, or 0 if any of them fail. + +=head1 SEE ALSO + +L<Math::Prime::Util> + +=head1 AUTHORS + +Dana Jacobsen E<lt>d...@acm.orge<gt> + + +=head1 COPYRIGHT + +Copyright 2012-2013 by Dana Jacobsen E<lt>d...@acm.orge<gt> + +This program is free software; you can redistribute it and/or modify it under the same terms as Perl itself. + +=cut diff --git a/t/13-primecount.t b/t/13-primecount.t index 37fb1f4..2ba8afb 100644 --- a/t/13-primecount.t +++ b/t/13-primecount.t @@ -87,7 +87,7 @@ plan tests => 0 + 1 + $use64 * 3 * scalar(keys %pivals64) + scalar(keys %intervals) + 1 - + 4 + 2*$extra; # prime count specific methods + + 5 + 2*$extra; # prime count specific methods ok( eval { prime_count(13); 1; }, "prime_count in void context"); @@ -165,6 +165,8 @@ SKIP: { is(Math::Prime::Util::_XS_LMO_pi (66123456), 3903023,"XS LMO count"); is(Math::Prime::Util::_XS_segment_pi (66123456), 3903023,"XS segment count"); } + +require_ok 'Math::Prime::Util::PP'; is(Math::Prime::Util::PP::_lehmer_pi (1456789), 111119, "PP Lehmer count"); is(Math::Prime::Util::PP::_sieve_prime_count(145678), 13478, "PP sieve count"); if ($extra) { diff --git a/t/23-primality-proofs.t b/t/23-primality-proofs.t index 69866de..98a77f1 100644 --- a/t/23-primality-proofs.t +++ b/t/23-primality-proofs.t @@ -20,6 +20,7 @@ my $broken64 = (18446744073709550592 == ~0); # Do some tests only if: # EXTENDED_TESTING is on OR we have the GMP backend # Note that with Calc, these things are incredibly slow. +use Math::BigInt try=>"GMP,Pari"; my $doexpensive = 0 + ($extra || ( (!$use64 || !$broken64) && Math::BigInt->config()->{lib} eq 'Math::BigInt::GMP' )); my @plist = qw/20907001 809120722675364249/; diff --git a/t/70-rt-bignum.t b/t/70-rt-bignum.t index b0ec0ef..a3001ed 100644 --- a/t/70-rt-bignum.t +++ b/t/70-rt-bignum.t @@ -10,6 +10,7 @@ use warnings; # The second method in theory is all that is needed. use Math::Prime::Util qw/:all/; +use Math::Prime::Util::PP; use bignum; use Test::More tests => 2; diff --git a/t/81-bignum.t b/t/81-bignum.t index 6d90652..0e064cf 100644 --- a/t/81-bignum.t +++ b/t/81-bignum.t @@ -145,7 +145,7 @@ my $mpugmpver = $usegmp ? $Math::Prime::Util::GMP::VERSION : "<none>"; diag "BigInt $bignumver/$bigintver, lib: $bigintlib. MPU::GMP $mpugmpver\n"; # Turn off use of BRS - ECM tries to use this. -prime_set_config( irand => sub { int(rand(4294967296.0)) } ); +prime_set_config( irand => sub { int(rand(4294967296)) } ); ############################################################################### -- 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