[libmath-prime-util-perl] 10/33: Load PP and Math::BigInt only when used. More work needed
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:51:41 UTC 2015
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 <dana at 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>dana at acm.orgE<gt>
+
+
+=head1 COPYRIGHT
+
+Copyright 2012-2013 by Dana Jacobsen E<lt>dana at 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
More information about the Pkg-perl-cvs-commits
mailing list