1package Math::Prime::Util::RandomPrimes; 2use strict; 3use warnings; 4use Carp qw/carp croak confess/; 5use Math::Prime::Util qw/ prime_get_config 6 verify_prime 7 is_provable_prime_with_cert 8 primorial prime_count nth_prime 9 is_prob_prime is_strong_pseudoprime 10 next_prime prev_prime 11 urandomb urandomm random_bytes 12 /; 13 14BEGIN { 15 $Math::Prime::Util::RandomPrimes::AUTHORITY = 'cpan:DANAJ'; 16 $Math::Prime::Util::RandomPrimes::VERSION = '0.73'; 17} 18 19BEGIN { 20 do { require Math::BigInt; Math::BigInt->import(try=>"GMP,Pari"); } 21 unless defined $Math::BigInt::VERSION; 22 23 use constant OLD_PERL_VERSION=> $] < 5.008; 24 use constant MPU_MAXBITS => (~0 == 4294967295) ? 32 : 64; 25 use constant MPU_64BIT => MPU_MAXBITS == 64; 26 use constant MPU_32BIT => MPU_MAXBITS == 32; 27 use constant MPU_MAXPARAM => MPU_32BIT ? 4294967295 : 18446744073709551615; 28 use constant MPU_MAXDIGITS => MPU_32BIT ? 10 : 20; 29 use constant MPU_USE_XS => prime_get_config->{'xs'}; 30 use constant MPU_USE_GMP => prime_get_config->{'gmp'}; 31 32 *_bigint_to_int = \&Math::Prime::Util::_bigint_to_int; 33} 34 35################################################################################ 36 37# These are much faster than straightforward trial division when n is big. 38# You'll want to first do a test up to and including 23. 39my @_big_gcd; 40my $_big_gcd_top = 20046; 41my $_big_gcd_use = -1; 42sub _make_big_gcds { 43 return if $_big_gcd_use >= 0; 44 if (prime_get_config->{'gmp'}) { 45 $_big_gcd_use = 0; 46 return; 47 } 48 if (Math::BigInt->config()->{lib} !~ /^Math::BigInt::(GMP|Pari)/) { 49 $_big_gcd_use = 0; 50 return; 51 } 52 $_big_gcd_use = 1; 53 my $p0 = primorial(Math::BigInt->new( 520)); 54 my $p1 = primorial(Math::BigInt->new(2052)); 55 my $p2 = primorial(Math::BigInt->new(6028)); 56 my $p3 = primorial(Math::BigInt->new($_big_gcd_top)); 57 $_big_gcd[0] = $p0->bdiv(223092870)->bfloor->as_int; 58 $_big_gcd[1] = $p1->bdiv($p0)->bfloor->as_int; 59 $_big_gcd[2] = $p2->bdiv($p1)->bfloor->as_int; 60 $_big_gcd[3] = $p3->bdiv($p2)->bfloor->as_int; 61} 62 63################################################################################ 64 65 66################################################################################ 67 68 69 70# For random primes, there are two good papers that should be examined: 71# 72# "Fast Generation of Prime Numbers and Secure Public-Key 73# Cryptographic Parameters" by Ueli M. Maurer, 1995 74# http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.26.2151 75# related discussions: 76# http://www.daimi.au.dk/~ivan/provableprimesproject.pdf 77# Handbook of Applied Cryptography by Menezes, et al. 78# 79# "Close to Uniform Prime Number Generation With Fewer Random Bits" 80# by Pierre-Alain Fouque and Mehdi Tibouchi, 2011 81# http://eprint.iacr.org/2011/481 82# 83# Some things to note: 84# 85# 1) Joye and Paillier have patents on their methods. Never use them. 86# 87# 2) The easy method of next_prime(random number), known as PRIMEINC, is 88# fast but gives a terrible distribution. It has a positive bias and 89# most importantly the probability for a prime is proportional to its 90# gap, meaning some numbers in the range will be thousands of times 91# more likely than others). On the contrary however, nobody has a way 92# to exploit this, and it's not-uncommon to see used. 93# 94# We use: 95# TRIVIAL range within native integer size (2^32 or 2^64) 96# FTA1 random_nbit_prime with 65+ bits 97# INVA1 other ranges with 65+ bit range 98# where 99# TRIVIAL = monte-carlo method or equivalent, perfect uniformity. 100# FTA1 = Fouque/Tibouchi A1, very close to uniform 101# INVA1 = inverted FTA1, less uniform but works with arbitrary ranges 102# 103# The random_maurer_prime function uses Maurer's FastPrime algorithm. 104# 105# If Math::Prime::Util::GMP is installed, these functions will be many times 106# faster than other methods (e.g. Math::Pari monte-carlo or Crypt::Primes). 107# 108# Timings on Macbook. 109# The "with GMP" numbers use Math::Prime::Util::GMP 0.44. 110# The "no GMP" numbers are with no Math::BigInt backend, so very slow in comparison. 111# If another backend was used (GMP, Pari, LTM) it would be more comparable. 112# 113# random_nbit_prime random_maurer_prime 114# n-bits no GMP w/ MPU::GMP no GMP w/ MPU::GMP 115# ---------- -------- ----------- -------- ----------- 116# 24-bit 1uS same same same 117# 64-bit 5uS same same same 118# 128-bit 0.12s 70uS 0.29s 166uS 119# 256-bit 0.66s 379uS 1.82s 800uS 120# 512-bit 7.8s 0.0022s 16.2s 0.0044s 121# 1024-bit ---- 0.019s ---- 0.037s 122# 2048-bit ---- 0.23s ---- 0.35s 123# 4096-bit ---- 2.4s ---- 5.2s 124# 125# Random timings for 10M calls on i4770K: 126# 0.39 Math::Random::MTwist 0.13 127# 0.41 ntheory <==== us 128# 0.89 system rand 129# 1.76 Math::Random::MT::Auto 130# 5.35 Bytes::Random::Secure OO w/ISAAC::XS 131# 7.43 Math::Random::Secure w/ISAAC::XS 132# 12.40 Math::Random::Secure 133# 12.78 Bytes::Random::Secure OO 134# 13.86 Bytes::Random::Secure function w/ISAAC::XS 135# 21.95 Bytes::Random::Secure function 136# 822.1 Crypt::Random 137# 138# time perl -E 'use Math::Random::MTwist "irand32"; irand32() for 1..10000000;' 139# time perl -E 'sub irand {int(rand(4294967296));} irand() for 1..10000000;' 140# time perl -E 'use Math::Random::MT::Auto; sub irand { Math::Random::MT::Auto::irand() & 0xFFFFFFFF } irand() for 1..10000000;' 141# time perl -E 'use Math::Random::Secure qw/irand/; irand() for 1..10000000;' 142# time perl -E 'use Bytes::Random::Secure qw/random_bytes/; sub irand {return unpack("L",random_bytes(4));} irand() for 1..10000000;' 143# time perl -E 'use Bytes::Random::Secure; my $rng = Bytes::Random::Secure->new(); sub irand {return $rng->irand;} irand() for 1..10000000;' 144# time perl -E 'use Crypt::Random qw/makerandom/; sub irand {makerandom(Size=>32, Uniform=>1, Strength=>0)} irand() for 1..100_000;' 145# > haveged daemon running to stop /dev/random blocking 146# > Both BRS and CR have more features that this isn't measuring. 147# 148# To verify distribution: 149# 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;' 150# 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;' 151 152# Sub to call with low and high already primes and verified range. 153my $_random_prime = sub { 154 my($low,$high) = @_; 155 my $prime; 156 157 #{ my $bsize = 100; my @bins; my $counts = 10000000; 158 # for my $c (1..$counts) { $bins[ $_IRANDF->($bsize-1) ]++; } 159 # for my $b (0..$bsize) {printf("%4d %8.5f%%\n", $b, $bins[$b]/$counts);} } 160 161 # low and high are both odds, and low < high. 162 163 # This is fast for small values, low memory, perfectly uniform, and 164 # consumes the minimum amount of randomness needed. But it isn't feasible 165 # with large values. Also note that low must be a prime. 166 if ($high <= 262144 && MPU_USE_XS) { 167 my $li = prime_count(2, $low); 168 my $irange = prime_count($low, $high); 169 my $rand = urandomm($irange); 170 return nth_prime($li + $rand); 171 } 172 173 $low-- if $low == 2; # Low of 2 becomes 1 for our program. 174 # Math::BigInt::GMP's RT 71548 will wreak havoc if we don't do this. 175 $low = Math::BigInt->new("$low") if ref($high) eq 'Math::BigInt'; 176 confess "Invalid _random_prime parameters: $low, $high" if ($low % 2) == 0 || ($high % 2) == 0; 177 178 # We're going to look at the odd numbers only. 179 my $oddrange = (($high - $low) >> 1) + 1; 180 181 croak "Large random primes not supported on old Perl" 182 if OLD_PERL_VERSION && MPU_64BIT && $oddrange > 4294967295; 183 184 # If $low is large (e.g. >10 digits) and $range is small (say ~10k), it 185 # would be fastest to call primes in the range and randomly pick one. I'm 186 # not implementing it now because it seems like a rare case. 187 188 # If the range is reasonably small, generate using simple Monte Carlo 189 # method (aka the 'trivial' method). Completely uniform. 190 if ($oddrange < MPU_MAXPARAM) { 191 my $loop_limit = 2000 * 1000; # To protect against broken rand 192 if ($low > 11) { 193 while ($loop_limit-- > 0) { 194 $prime = $low + 2 * urandomm($oddrange); 195 next if !($prime % 3) || !($prime % 5) || !($prime % 7) || !($prime % 11); 196 return $prime if is_prob_prime($prime); 197 } 198 } else { 199 while ($loop_limit-- > 0) { 200 $prime = $low + 2 * urandomm($oddrange); 201 next if $prime > 11 && (!($prime % 3) || !($prime % 5) || !($prime % 7) || !($prime % 11)); 202 return 2 if $prime == 1; # Remember the special case for 2. 203 return $prime if is_prob_prime($prime); 204 } 205 } 206 croak "Random function broken?"; 207 } 208 209 # We have an ocean of range, and a teaspoon to hold randomness. 210 211 # Since we have an arbitrary range and not a power of two, I don't see how 212 # Fouque's algorithm A1 could be used (where we generate lower bits and 213 # generate random sets of upper). Similarly trying to simply generate 214 # upper bits is full of ways to trip up and get non-uniform results. 215 # 216 # What I'm doing here is: 217 # 218 # 1) divide the range into semi-evenly sized partitions, where each part 219 # is as close to $rand_max_val as we can. 220 # 2) randomly select one of the partitions. 221 # 3) iterate choosing random values within the partition. 222 # 223 # The downside is that we're skewing a _lot_ farther from uniformity than 224 # we'd like. Imagine we started at 0 with 1e18 partitions of size 100k 225 # each. 226 # Probability of '5' being returned = 227 # 1.04e-22 = 1e-18 (chose first partition) * 1/9592 (chose '5') 228 # Probability of '100003' being returned = 229 # 1.19e-22 = 1e-18 (chose second partition) * 1/8392 (chose '100003') 230 # Probability of '99999999999999999999977' being returned = 231 # 5.20e-22 = 1e-18 (chose last partition) * 1/1922 (chose '99...77') 232 # So the primes in the last partition will show up 5x more often. 233 # The partitions are selected uniformly, and the primes within are selected 234 # uniformly, but the number of primes in each bucket is _not_ uniform. 235 # Their individual probability of being selected is the probability of the 236 # partition (uniform) times the probability of being selected inside the 237 # partition (uniform with respect to all other primes in the same 238 # partition, but each partition is different and skewed). 239 # 240 # Partitions are typically much larger than 100k, but with a huge range 241 # we still see this (e.g. ~3x from 0-10^30, ~10x from 0-10^100). 242 # 243 # When selecting n-bit or n-digit primes, this effect is MUCH smaller, as 244 # the skew becomes approx lg(2^n) / lg(2^(n-1)) which is pretty close to 1. 245 # 246 # 247 # Another idea I'd like to try sometime is: 248 # pclo = prime_count_lower(low); 249 # pchi = prime_count_upper(high); 250 # do { 251 # $nth = random selection between pclo and pchi 252 # $prguess = nth_prime_approx($nth); 253 # } while ($prguess >= low) && ($prguess <= high); 254 # monte carlo select a prime in $prguess-2**24 to $prguess+2**24 255 # which accounts for the prime distribution. 256 257 my($binsize, $nparts); 258 my $rand_part_size = 1 << (MPU_64BIT ? 32 : 31); 259 if (ref($oddrange) eq 'Math::BigInt') { 260 # Go to some trouble here because some systems are wonky, such as 261 # giving us +a/+b = -r. Also note the quotes for the bigint argument. 262 # Without that, Math::BigInt::GMP can return garbage. 263 my($nbins, $rem); 264 ($nbins, $rem) = $oddrange->copy->bdiv( "$rand_part_size" ); 265 $nbins++ if $rem > 0; 266 $nbins = $nbins->as_int(); 267 ($binsize,$rem) = $oddrange->copy->bdiv($nbins); 268 $binsize++ if $rem > 0; 269 $binsize = $binsize->as_int(); 270 $nparts = $oddrange->copy->bdiv($binsize)->as_int(); 271 $low = $high->copy->bzero->badd($low) if ref($low) ne 'Math::BigInt'; 272 } else { 273 my $nbins = int($oddrange / $rand_part_size); 274 $nbins++ if $nbins * $rand_part_size != $oddrange; 275 $binsize = int($oddrange / $nbins); 276 $binsize++ if $binsize * $nbins != $oddrange; 277 $nparts = int($oddrange/$binsize); 278 } 279 $nparts-- if ($nparts * $binsize) == $oddrange; 280 281 my $rpart = urandomm($nparts+1); 282 283 my $primelow = $low + 2 * $binsize * $rpart; 284 my $partsize = ($rpart < $nparts) ? $binsize 285 : $oddrange - ($nparts * $binsize); 286 $partsize = _bigint_to_int($partsize) if ref($partsize) eq 'Math::BigInt'; 287 #warn "range $oddrange = $nparts * $binsize + ", $oddrange - ($nparts * $binsize), "\n"; 288 #warn " chose part $rpart size $partsize\n"; 289 #warn " primelow is $low + 2 * $binsize * $rpart = $primelow\n"; 290 #die "Result could be too large" if ($primelow + 2*($partsize-1)) > $high; 291 292 # Generate random numbers in the interval until one is prime. 293 my $loop_limit = 2000 * 1000; # To protect against broken rand 294 295 # Simply things for non-bigints. 296 if (ref($low) ne 'Math::BigInt') { 297 while ($loop_limit-- > 0) { 298 my $rand = urandomm($partsize); 299 $prime = $primelow + $rand + $rand; 300 croak "random prime failure, $prime > $high" if $prime > $high; 301 if ($prime <= 23) { 302 $prime = 2 if $prime == 1; # special case for low = 2 303 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]; 304 return $prime; 305 } 306 next if !($prime % 3) || !($prime % 5) || !($prime % 7) || !($prime % 11); 307 # It looks promising. Check it. 308 next unless is_prob_prime($prime); 309 return $prime; 310 } 311 croak "Random function broken?"; 312 } 313 314 # By checking a wheel 30 mod, we can skip anything that would be a multiple 315 # of 2, 3, or 5, without even having to create the bigint prime. 316 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); 317 my $primelow30 = $primelow % 30; 318 $primelow30 = _bigint_to_int($primelow30) if ref($primelow30) eq 'Math::BigInt'; 319 320 # Big GCD's are hugely fast with GMP or Pari, but super slow with Calc. 321 _make_big_gcds() if $_big_gcd_use < 0; 322 323 while ($loop_limit-- > 0) { 324 my $rand = urandomm($partsize); 325 # Check wheel-30 mod 326 my $rand30 = $rand % 30; 327 next if $w30[($primelow30 + 2*$rand30) % 30] 328 && ($rand > 3 || $primelow > 5); 329 # Construct prime 330 $prime = $primelow + $rand + $rand; 331 croak "random prime failure, $prime > $high" if $prime > $high; 332 if ($prime <= 23) { 333 $prime = 2 if $prime == 1; # special case for low = 2 334 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]; 335 return $prime; 336 } 337 # With GMP, the fastest thing to do is check primality. 338 if (MPU_USE_GMP) { 339 next unless Math::Prime::Util::GMP::is_prime($prime); 340 return $prime; 341 } 342 # No MPU:GMP, so primality checking is slow. Skip some composites here. 343 next unless Math::BigInt::bgcd($prime, 7436429) == 1; 344 if ($_big_gcd_use && $prime > $_big_gcd_top) { 345 next unless Math::BigInt::bgcd($prime, $_big_gcd[0]) == 1; 346 next unless Math::BigInt::bgcd($prime, $_big_gcd[1]) == 1; 347 next unless Math::BigInt::bgcd($prime, $_big_gcd[2]) == 1; 348 next unless Math::BigInt::bgcd($prime, $_big_gcd[3]) == 1; 349 } 350 # It looks promising. Check it. 351 next unless is_prob_prime($prime); 352 return $prime; 353 } 354 croak "Random function broken?"; 355}; 356 357# Cache of tight bounds for each digit. Helps performance a lot. 358my @_random_ndigit_ranges = (undef, [2,7], [11,97] ); 359my @_random_nbit_ranges = (undef, undef, [2,3],[5,7] ); 360my %_random_cache_small; 361 362# For fixed small ranges with XS, e.g. 6-digit, 18-bit 363sub _random_xscount_prime { 364 my($low,$high) = @_; 365 my($istart, $irange); 366 my $cachearef = $_random_cache_small{$low,$high}; 367 if (defined $cachearef) { 368 ($istart, $irange) = @$cachearef; 369 } else { 370 my $beg = ($low <= 2) ? 2 : next_prime($low-1); 371 my $end = ($high < ~0) ? prev_prime($high + 1) : prev_prime($high); 372 ($istart, $irange) = ( prime_count(2, $beg), prime_count($beg, $end) ); 373 $_random_cache_small{$low,$high} = [$istart, $irange]; 374 } 375 my $rand = urandomm($irange); 376 return nth_prime($istart + $rand); 377} 378 379sub random_prime { 380 my($low,$high) = @_; 381 return if $high < 2 || $low > $high; 382 383 if ($high-$low > 1000000000) { 384 # Range is large, just make them odd if needed. 385 $low = 2 if $low < 2; 386 $low++ if $low > 2 && ($low % 2) == 0; 387 $high-- if ($high % 2) == 0; 388 } else { 389 # Tighten the range to the nearest prime. 390 $low = ($low <= 2) ? 2 : next_prime($low-1); 391 $high = ($high == ~0) ? prev_prime($high) : prev_prime($high + 1); 392 return $low if ($low == $high) && is_prob_prime($low); 393 return if $low >= $high; 394 # At this point low and high are both primes, and low < high. 395 } 396 397 # At this point low and high are both primes, and low < high. 398 return $_random_prime->($low, $high); 399} 400 401sub random_ndigit_prime { 402 my($digits) = @_; 403 croak "random_ndigit_prime, digits must be >= 1" unless $digits >= 1; 404 405 return _random_xscount_prime( int(10 ** ($digits-1)), int(10 ** $digits) ) 406 if $digits <= 6 && MPU_USE_XS; 407 408 my $bigdigits = $digits >= MPU_MAXDIGITS; 409 if ($bigdigits && prime_get_config->{'nobigint'}) { 410 croak "random_ndigit_prime with -nobigint, digits out of range" 411 if $digits > MPU_MAXDIGITS; 412 # Special case for nobigint and threshold digits 413 if (!defined $_random_ndigit_ranges[$digits]) { 414 my $low = int(10 ** ($digits-1)); 415 my $high = ~0; 416 $_random_ndigit_ranges[$digits] = [next_prime($low),prev_prime($high)]; 417 } 418 } 419 420 if (!defined $_random_ndigit_ranges[$digits]) { 421 if ($bigdigits) { 422 my $low = Math::BigInt->new('10')->bpow($digits-1); 423 my $high = Math::BigInt->new('10')->bpow($digits); 424 # Just pull the range in to the nearest odd. 425 $_random_ndigit_ranges[$digits] = [$low+1, $high-1]; 426 } else { 427 my $low = int(10 ** ($digits-1)); 428 my $high = int(10 ** $digits); 429 # Note: Perl 5.6.2 cannot represent 10**15 as an integer, so things 430 # will crash all over the place if you try. We can stringify it, but 431 # will just fail tests later. 432 $_random_ndigit_ranges[$digits] = [next_prime($low),prev_prime($high)]; 433 } 434 } 435 my ($low, $high) = @{$_random_ndigit_ranges[$digits]}; 436 return $_random_prime->($low, $high); 437} 438 439my @_random_nbit_m; 440my @_random_nbit_lambda; 441my @_random_nbit_arange; 442 443sub random_nbit_prime { 444 my($bits) = @_; 445 croak "random_nbit_prime, bits must be >= 2" unless $bits >= 2; 446 $bits = int("$bits"); 447 448 # Very small size, use the nth-prime method 449 if ($bits <= 20 && MPU_USE_XS) { 450 if ($bits <= 4) { 451 return (2,3)[urandomb(1)] if $bits == 2; 452 return (5,7)[urandomb(1)] if $bits == 3; 453 return (11,13)[urandomb(1)] if $bits == 4; 454 } 455 return _random_xscount_prime( 1 << ($bits-1), 1 << $bits ); 456 } 457 458 croak "Mid-size random primes not supported on broken old Perl" 459 if OLD_PERL_VERSION && MPU_64BIT && $bits > 49 && $bits <= 64; 460 461 # Fouque and Tibouchi (2011) Algorithm 1 (basic) 462 # Modified to make sure the nth bit is always set. 463 # 464 # Example for random_nbit_prime(512) on 64-bit Perl: 465 # p: 1aaaaaaaabbbbbbbbbbbbbbbbbbbb1 466 # ^^ ^ ^--- Trailing 1 so p is odd 467 # || +--- 512-63-2 = 447 lower bits selected before loop 468 # |+--- 63 upper bits selected in loop, repeated until p is prime 469 # +--- upper bit is 1 so we generate an n-bit prime 470 # total: 1 + 63 + 447 + 1 = 512 bits 471 # 472 # Algorithm 2 is implemented in a previous commit on github. The problem 473 # is that it doesn't set the nth bit, and making that change requires a 474 # modification of the algorithm. It was not a lot faster than this A1 475 # with the native int trial division. If the irandf function was very 476 # slow, then A2 would look more promising. 477 # 478 if (1 && $bits > 64) { 479 my $l = (MPU_64BIT && $bits > 79) ? 63 : 31; 480 $l = 49 if $l == 63 && OLD_PERL_VERSION; # Fix for broken Perl 5.6 481 $l = $bits-2 if $bits-2 < $l; 482 483 my $brand = urandomb($bits-$l-2); 484 $brand = Math::BigInt->new("$brand") unless ref($brand) eq 'Math::BigInt'; 485 my $b = $brand->blsft(1)->binc(); 486 487 # Precalculate some modulii so we can do trial division on native int 488 # 9699690 = 2*3*5*7*11*13*17*19, so later operations can be native ints 489 my @premod; 490 my $bpremod = _bigint_to_int($b->copy->bmod(9699690)); 491 my $twopremod = _bigint_to_int(Math::BigInt->new(2)->bmodpow($bits-$l-1, 9699690)); 492 foreach my $zi (0 .. 19-1) { 493 foreach my $pm (3, 5, 7, 11, 13, 17, 19) { 494 next if $zi >= $pm || defined $premod[$pm]; 495 $premod[$pm] = $zi if ( ($twopremod*$zi+$bpremod) % $pm) == 0; 496 } 497 } 498 _make_big_gcds() if $_big_gcd_use < 0; 499 if (!MPU_USE_GMP) { require Math::Prime::Util::PP; } 500 501 my $loop_limit = 1_000_000; 502 while ($loop_limit-- > 0) { 503 my $a = (1 << $l) + urandomb($l); 504 # $a % s == $premod[s] => $p % s == 0 => p will be composite 505 next if $a % 3 == $premod[ 3] || $a % 5 == $premod[ 5] 506 || $a % 7 == $premod[ 7] || $a % 11 == $premod[11] 507 || $a % 13 == $premod[13] || $a % 17 == $premod[17] 508 || $a % 19 == $premod[19]; 509 my $p = Math::BigInt->new("$a")->blsft($bits-$l-1)->badd($b); 510 #die " $a $b $p" if $a % 11 == $premod[11] && $p % 11 != 0; 511 #die "!$a $b $p" if $a % 11 != $premod[11] && $p % 11 == 0; 512 if (MPU_USE_GMP) { 513 next unless Math::Prime::Util::GMP::is_prime($p); 514 } else { 515 next unless Math::BigInt::bgcd($p, 1348781387) == 1; # 23-43 516 if ($_big_gcd_use && $p > $_big_gcd_top) { 517 next unless Math::BigInt::bgcd($p, $_big_gcd[0]) == 1; 518 next unless Math::BigInt::bgcd($p, $_big_gcd[1]) == 1; 519 next unless Math::BigInt::bgcd($p, $_big_gcd[2]) == 1; 520 next unless Math::BigInt::bgcd($p, $_big_gcd[3]) == 1; 521 } 522 # We know we don't have GMP and are > 2^64, so go directly to this. 523 next unless Math::Prime::Util::PP::is_bpsw_prime($p); 524 } 525 return $p; 526 } 527 croak "Random function broken?"; 528 } 529 530 # The Trivial method. Great uniformity, and fine for small sizes. It 531 # gets very slow as the bit size increases, but that is why we have the 532 # method above for bigints. 533 if (1) { 534 535 my $loop_limit = 2_000_000; 536 if ($bits > MPU_MAXBITS) { 537 my $p = Math::BigInt->bone->blsft($bits-1)->binc(); 538 while ($loop_limit-- > 0) { 539 my $n = Math::BigInt->new(''.urandomb($bits-2))->blsft(1)->badd($p); 540 return $n if is_prob_prime($n); 541 } 542 } else { 543 my $p = (1 << ($bits-1)) + 1; 544 while ($loop_limit-- > 0) { 545 my $n = $p + (urandomb($bits-2) << 1); 546 return $n if is_prob_prime($n); 547 } 548 } 549 croak "Random function broken?"; 550 551 } else { 552 553 # Send through the generic random_prime function. Decently fast, but 554 # quite a bit slower than the F&T A1 method above. 555 if (!defined $_random_nbit_ranges[$bits]) { 556 if ($bits > MPU_MAXBITS) { 557 my $low = Math::BigInt->new('2')->bpow($bits-1); 558 my $high = Math::BigInt->new('2')->bpow($bits); 559 # Don't pull the range in to primes, just odds 560 $_random_nbit_ranges[$bits] = [$low+1, $high-1]; 561 } else { 562 my $low = 1 << ($bits-1); 563 my $high = ($bits == MPU_MAXBITS) 564 ? ~0-1 565 : ~0 >> (MPU_MAXBITS - $bits); 566 $_random_nbit_ranges[$bits] = [next_prime($low-1),prev_prime($high+1)]; 567 # Example: bits = 7. 568 # low = 1<<6 = 64. next_prime(64-1) = 67 569 # high = ~0 >> (64-7) = 127. prev_prime(127+1) = 127 570 } 571 } 572 my ($low, $high) = @{$_random_nbit_ranges[$bits]}; 573 return $_random_prime->($low, $high); 574 575 } 576} 577 578 579# For stripping off the header on certificates so they can be combined. 580sub _strip_proof_header { 581 my $proof = shift; 582 $proof =~ s/^\[MPU - Primality Certificate\]\nVersion \S+\n+Proof for:\nN (\d+)\n+//ms; 583 return $proof; 584} 585 586 587sub random_maurer_prime { 588 my $k = shift; 589 croak "random_maurer_prime, bits must be >= 2" unless $k >= 2; 590 $k = int("$k"); 591 592 return random_nbit_prime($k) if $k <= MPU_MAXBITS && !OLD_PERL_VERSION; 593 594 my ($n, $cert) = random_maurer_prime_with_cert($k); 595 croak "maurer prime $n failed certificate verification!" 596 unless verify_prime($cert); 597 return $n; 598} 599 600sub random_maurer_prime_with_cert { 601 my $k = shift; 602 croak "random_maurer_prime, bits must be >= 2" unless $k >= 2; 603 $k = int("$k"); 604 605 # This should never happen. Trap now to prevent infinite loop. 606 croak "number of bits must not be a bigint" if ref($k) eq 'Math::BigInt'; 607 608 # Results for random_nbit_prime are proven for all native bit sizes. 609 my $p0 = MPU_MAXBITS; 610 $p0 = 49 if OLD_PERL_VERSION && MPU_MAXBITS > 49; 611 612 if ($k <= $p0) { 613 my $n = random_nbit_prime($k); 614 my ($isp, $cert) = is_provable_prime_with_cert($n); 615 croak "small nbit prime could not be proven" if $isp != 2; 616 return ($n, $cert); 617 } 618 619 # Set verbose to 3 to get pretty output like Crypt::Primes 620 my $verbose = prime_get_config->{'verbose'}; 621 local $| = 1 if $verbose > 2; 622 623 do { require Math::BigFloat; Math::BigFloat->import(); } 624 if !defined $Math::BigFloat::VERSION; 625 626 # Ignore Maurer's g and c that controls how much trial division is done. 627 my $r = Math::BigFloat->new("0.5"); # relative size of the prime q 628 my $m = 20; # makes sure R is big enough 629 630 # Generate a random prime q of size $r*$k, where $r >= 0.5. Try to 631 # cleverly select r to match the size of a typical random factor. 632 if ($k > 2*$m) { 633 do { 634 my $s = Math::Prime::Util::drand(); 635 $r = Math::BigFloat->new(2)->bpow($s-1); 636 } while ($k*$r >= $k-$m); 637 } 638 639 # I've seen +0, +1, and +2 here. Maurer uses +0. Menezes uses +1. 640 # We can use +1 because we're using BLS75 theorem 3 later. 641 my $smallk = int(($r * $k)->bfloor->bstr) + 1; 642 my ($q, $qcert) = random_maurer_prime_with_cert($smallk); 643 $q = Math::BigInt->new("$q") unless ref($q) eq 'Math::BigInt'; 644 my $I = Math::BigInt->new(2)->bpow($k-2)->bdiv($q)->bfloor->as_int(); 645 print "r = $r k = $k q = $q I = $I\n" if $verbose && $verbose != 3; 646 $qcert = ($q < Math::BigInt->new("18446744073709551615")) 647 ? "" : _strip_proof_header($qcert); 648 649 # Big GCD's are hugely fast with GMP or Pari, but super slow with Calc. 650 _make_big_gcds() if $_big_gcd_use < 0; 651 my $ONE = Math::BigInt->bone; 652 my $TWO = $ONE->copy->binc; 653 654 my $loop_limit = 1_000_000 + $k * 1_000; 655 while ($loop_limit-- > 0) { 656 # R is a random number between $I+1 and 2*$I 657 #my $R = $I + 1 + urandomm( $I ); 658 my $R = $I->copy->binc->badd( urandomm($I) ); 659 #my $n = 2 * $R * $q + 1; 660 my $nm1 = $TWO->copy->bmul($R)->bmul($q); 661 my $n = $nm1->copy->binc; 662 # We constructed a promising looking $n. Now test it. 663 print "." if $verbose > 2; 664 if (MPU_USE_GMP) { 665 # MPU::GMP::is_prob_prime has fast tests built in. 666 next unless Math::Prime::Util::GMP::is_prob_prime($n); 667 } else { 668 # No GMP, so first do trial divisions, then a SPSP test. 669 next unless Math::BigInt::bgcd($n, 111546435)->is_one; 670 if ($_big_gcd_use && $n > $_big_gcd_top) { 671 next unless Math::BigInt::bgcd($n, $_big_gcd[0])->is_one; 672 next unless Math::BigInt::bgcd($n, $_big_gcd[1])->is_one; 673 next unless Math::BigInt::bgcd($n, $_big_gcd[2])->is_one; 674 next unless Math::BigInt::bgcd($n, $_big_gcd[3])->is_one; 675 } 676 print "+" if $verbose > 2; 677 next unless is_strong_pseudoprime($n, 3); 678 } 679 print "*" if $verbose > 2; 680 681 # We could pick a random generator by doing: 682 # Step 1: pick a random r 683 # Step 2: compute g = r^((n-1)/q) mod p 684 # Step 3: if g == 1, goto Step 1. 685 # Note that n = 2*R*q+1, hence the exponent is 2*R. 686 687 # We could set r = 0.3333 earlier, then use BLS75 theorem 5 here. 688 # The chain would be shorter, requiring less overall work for 689 # large inputs. Maurer's paper discusses the idea. 690 691 # Use BLS75 theorem 3. This is easier and possibly faster than 692 # BLS75 theorem 4 (Pocklington) used by Maurer's paper. 693 694 # Check conditions -- these should be redundant. 695 my $m = $TWO * $R; 696 if (! ($q->is_odd && $q > 2 && $m > 0 && 697 $m * $q + $ONE == $n && $TWO*$q+$ONE > $n->copy->bsqrt()) ) { 698 carp "Maurer prime failed BLS75 theorem 3 conditions. Retry."; 699 next; 700 } 701 # Find a suitable a. Move on if one isn't found quickly. 702 foreach my $trya (2, 3, 5, 7, 11, 13) { 703 my $a = Math::BigInt->new($trya); 704 # m/2 = R (n-1)/2 = (2*R*q)/2 = R*q 705 next unless $a->copy->bmodpow($R, $n) != $nm1; 706 next unless $a->copy->bmodpow($R*$q, $n) == $nm1; 707 print "($k)" if $verbose > 2; 708 croak "Maurer prime $n=2*$R*$q+1 failed BPSW" unless is_prob_prime($n); 709 my $cert = "[MPU - Primality Certificate]\nVersion 1.0\n\n" . 710 "Proof for:\nN $n\n\n" . 711 "Type BLS3\nN $n\nQ $q\nA $a\n" . 712 $qcert; 713 return ($n, $cert); 714 } 715 # Didn't pass the selected a values. Try another R. 716 } 717 croak "Failure in random_maurer_prime, could not find a prime\n"; 718} # End of random_maurer_prime 719 720 721sub random_shawe_taylor_prime_with_cert { 722 my $k = shift; 723 724 my $seed = random_bytes(512/8); 725 726 my($status,$prime,$prime_seed,$prime_gen_counter,$cert) 727 = _ST_Random_prime($k, $seed); 728 croak "Shawe-Taylor random prime failure" unless $status; 729 croak "Shawe-Taylor random prime failure: prime $prime failed certificate verification!" unless verify_prime($cert); 730 731 return ($prime,$cert); 732} 733 734sub _seed_plus_one { 735 my($s) = @_; 736 for (my $i = length($s)-1; $i >= 0; $i--) { 737 vec($s, $i, 8)++; 738 last unless vec($s, $i, 8) == 0; 739 } 740 return $s; 741} 742 743sub _ST_Random_prime { # From FIPS 186-4 744 my($k, $input_seed) = @_; 745 croak "Shawe-Taylor random prime must have length >= 2" if $k < 2; 746 $k = int("$k"); 747 748 croak "Shawe-Taylor random prime, invalid input seed" 749 unless defined $input_seed && length($input_seed) >= 32; 750 751 if (!defined $Digest::SHA::VERSION) { 752 eval { require Digest::SHA; 753 my $version = $Digest::SHA::VERSION; 754 $version =~ s/[^\d.]//g; 755 $version >= 4.00; } 756 or do { croak "Must have Digest::SHA 4.00 or later"; }; 757 } 758 759 my $k2 = Math::BigInt->new(2)->bpow($k-1); 760 761 if ($k < 33) { 762 my $seed = $input_seed; 763 my $prime_gen_counter = 0; 764 my $kmask = 0xFFFFFFFF >> (32-$k); # Does the mod operation 765 my $kstencil = (1 << ($k-1)) | 1; # Sets high and low bits 766 while (1) { 767 my $seedp1 = _seed_plus_one($seed); 768 my $cvec = Digest::SHA::sha256($seed) ^ Digest::SHA::sha256($seedp1); 769 # my $c = Math::BigInt->from_hex('0x' . unpack("H*", $cvec)); 770 # $c = $k2 + ($c % $k2); 771 # $c = (2 * ($c >> 1)) + 1; 772 my($c) = unpack("N*", substr($cvec,-4,4)); 773 $c = ($c & $kmask) | $kstencil; 774 $prime_gen_counter++; 775 $seed = _seed_plus_one($seedp1); 776 my ($isp, $cert) = is_provable_prime_with_cert($c); 777 return (1,$c,$seed,$prime_gen_counter,$cert) if $isp; 778 return (0,0,0,0) if $prime_gen_counter > 10000 + 16*$k; 779 } 780 } 781 my($status,$c0,$seed,$prime_gen_counter,$cert) 782 = _ST_Random_prime( (($k+1)>>1)+1, $input_seed); 783 return (0,0,0,0) unless $status; 784 $cert = ($c0 < Math::BigInt->new("18446744073709551615")) 785 ? "" : _strip_proof_header($cert); 786 my $iterations = int(($k + 255) / 256) - 1; # SHA256 generates 256 bits 787 my $old_counter = $prime_gen_counter; 788 my $xstr = ''; 789 for my $i (0 .. $iterations) { 790 $xstr = Digest::SHA::sha256_hex($seed) . $xstr; 791 $seed = _seed_plus_one($seed); 792 } 793 my $x = Math::BigInt->from_hex('0x'.$xstr); 794 $x = $k2 + ($x % $k2); 795 my $t = ($x + 2*$c0 - 1) / (2*$c0); 796 _make_big_gcds() if $_big_gcd_use < 0; 797 while (1) { 798 if (2*$t*$c0 + 1 > 2*$k2) { $t = ($k2 + 2*$c0 - 1) / (2*$c0); } 799 my $c = 2*$t*$c0 + 1; 800 $prime_gen_counter++; 801 802 # Don't do the Pocklington check unless the candidate looks prime 803 my $looks_prime = 0; 804 if (MPU_USE_GMP) { 805 # MPU::GMP::is_prob_prime has fast tests built in. 806 $looks_prime = Math::Prime::Util::GMP::is_prob_prime($c); 807 } else { 808 # No GMP, so first do trial divisions, then a SPSP test. 809 $looks_prime = Math::BigInt::bgcd($c, 111546435)->is_one; 810 if ($looks_prime && $_big_gcd_use && $c > $_big_gcd_top) { 811 $looks_prime = Math::BigInt::bgcd($c, $_big_gcd[0])->is_one && 812 Math::BigInt::bgcd($c, $_big_gcd[1])->is_one && 813 Math::BigInt::bgcd($c, $_big_gcd[2])->is_one && 814 Math::BigInt::bgcd($c, $_big_gcd[3])->is_one; 815 } 816 $looks_prime = 0 if $looks_prime && !is_strong_pseudoprime($c, 3); 817 } 818 819 if ($looks_prime) { 820 # We could use a in (2,3,5,7,11,13), but pedantically use FIPS 186-4. 821 my $astr = ''; 822 for my $i (0 .. $iterations) { 823 $astr = Digest::SHA::sha256_hex($seed) . $astr; 824 $seed = _seed_plus_one($seed); 825 } 826 my $a = Math::BigInt->from_hex('0x'.$astr); 827 $a = ($a % ($c-3)) + 2; 828 my $z = $a->copy->bmodpow(2*$t,$c); 829 if (Math::BigInt::bgcd($z-1,$c)->is_one && $z->copy->bmodpow($c0,$c)->is_one) { 830 croak "Shawe-Taylor random prime failure at ($k): $c not prime" 831 unless is_prob_prime($c); 832 $cert = "[MPU - Primality Certificate]\nVersion 1.0\n\n" . 833 "Proof for:\nN $c\n\n" . 834 "Type Pocklington\nN $c\nQ $c0\nA $a\n" . 835 $cert; 836 return (1, $c, $seed, $prime_gen_counter, $cert); 837 } 838 } else { 839 # Update seed "as if" we performed the Pocklington check from FIPS 186-4 840 for my $i (0 .. $iterations) { 841 $seed = _seed_plus_one($seed); 842 } 843 } 844 return (0,0,0,0) if $prime_gen_counter > 10000 + 16*$k + $old_counter; 845 $t++; 846 } 847} 848 849 850# Gordon's algorithm for generating a strong prime. 851sub random_strong_prime { 852 my $t = shift; 853 croak "random_strong_prime, bits must be >= 128" unless $t >= 128; 854 $t = int("$t"); 855 856 croak "Random strong primes must be >= 173 bits on old Perl" 857 if OLD_PERL_VERSION && MPU_64BIT && $t < 173; 858 859 my $l = (($t+1) >> 1) - 2; 860 my $lp = int($t/2) - 20; 861 my $lpp = $l - 20; 862 while (1) { 863 my $qp = random_nbit_prime($lp); 864 my $qpp = random_nbit_prime($lpp); 865 $qp = Math::BigInt->new("$qp") unless ref($qp) eq 'Math::BigInt'; 866 $qpp = Math::BigInt->new("$qpp") unless ref($qpp) eq 'Math::BigInt'; 867 my ($il, $rem) = Math::BigInt->new(2)->bpow($l-1)->bdec()->bdiv(2*$qpp); 868 $il++ if $rem > 0; 869 $il = $il->as_int(); 870 my $iu = Math::BigInt->new(2)->bpow($l)->bsub(2)->bdiv(2*$qpp)->as_int(); 871 my $istart = $il + urandomm($iu - $il + 1); 872 for (my $i = $istart; $i <= $iu; $i++) { # Search for q 873 my $q = 2 * $i * $qpp + 1; 874 next unless is_prob_prime($q); 875 my $pp = $qp->copy->bmodpow($q-2, $q)->bmul(2)->bmul($qp)->bdec(); 876 my ($jl, $rem) = Math::BigInt->new(2)->bpow($t-1)->bsub($pp)->bdiv(2*$q*$qp); 877 $jl++ if $rem > 0; 878 $jl = $jl->as_int(); 879 my $ju = Math::BigInt->new(2)->bpow($t)->bdec()->bsub($pp)->bdiv(2*$q*$qp)->as_int(); 880 my $jstart = $jl + urandomm($ju - $jl + 1); 881 for (my $j = $jstart; $j <= $ju; $j++) { # Search for p 882 my $p = $pp + 2 * $j * $q * $qp; 883 return $p if is_prob_prime($p); 884 } 885 } 886 } 887} 888 889sub random_proven_prime { 890 my $k = shift; 891 my ($n, $cert) = random_proven_prime_with_cert($k); 892 croak "random_proven_prime $n failed certificate verification!" 893 unless verify_prime($cert); 894 return $n; 895} 896 897sub random_proven_prime_with_cert { 898 my $k = shift; 899 900 if (prime_get_config->{'gmp'} && $k <= 450) { 901 my $n = random_nbit_prime($k); 902 my ($isp, $cert) = is_provable_prime_with_cert($n); 903 croak "small nbit prime could not be proven" if $isp != 2; 904 return ($n, $cert); 905 } 906 return random_maurer_prime_with_cert($k); 907} 908 9091; 910 911__END__ 912 913 914# ABSTRACT: Generate random primes 915 916=pod 917 918=encoding utf8 919 920=head1 NAME 921 922Math::Prime::Util::RandomPrimes - Generate random primes 923 924 925=head1 VERSION 926 927Version 0.73 928 929 930=head1 SYNOPSIS 931 932=head1 DESCRIPTION 933 934Routines to generate random primes, including constructing proven primes. 935 936 937=head1 RANDOM PRIME FUNCTIONS 938 939=head2 random_prime 940 941Generate a random prime between C<low> and C<high>. If given one argument, 942C<low> will be 2. 943 944=head2 random_ndigit_prime 945 946Generate a random prime with C<n> digits. C<n> must be at least 1. 947 948=head2 random_nbit_prime 949 950Generate a random prime with C<n> bits. C<n> must be at least 2. 951 952=head2 random_maurer_prime 953 954Construct a random provable prime of C<n> bits using Maurer's FastPrime 955algorithm. C<n> must be at least 2. 956 957=head2 random_maurer_prime_with_cert 958 959Construct a random provable prime of C<n> bits using Maurer's FastPrime 960algorithm. C<n> must be at least 2. Returns a list of two items: the 961prime and the certificate. 962 963=head2 random_shawe_taylor_prime 964 965Construct a random provable prime of C<n> bits using Shawe-Taylor's 966algorithm. C<n> must be at least 2. The implementation is from 967FIPS 186-4 and uses SHA-256 with 512 bits of randomness. 968 969=head2 random_shawe_taylor_prime_with_cert 970 971Construct a random provable prime of C<n> bits using Shawe-Taylor's 972algorithm. C<n> must be at least 2. Returns a list of two items: the 973prime and the certificate. 974 975=head2 random_strong_prime 976 977Construct a random strong prime of C<n> bits. C<n> must be at least 128. 978 979=head2 random_proven_prime 980 981Generate or construct a random provable prime of C<n> bits. C<n> must 982be at least 2. 983 984=head2 random_proven_prime_with_cert 985 986Generate or construct a random provable prime of C<n> bits. C<n> must 987be at least 2. Returns a list of two items: the prime and the certificate. 988 989 990=head1 SEE ALSO 991 992L<Math::Prime::Util> 993 994=head1 AUTHORS 995 996Dana Jacobsen E<lt>dana@acm.orgE<gt> 997 998 999=head1 COPYRIGHT 1000 1001Copyright 2012-2013 by Dana Jacobsen E<lt>dana@acm.orgE<gt> 1002 1003This program is free software; you can redistribute it and/or modify it under the same terms as Perl itself. 1004 1005=cut 1006