File Coverage

blib/lib/Math/Prime/Util/RandomPrimes.pm
Criterion Covered Total %
statement 263 343 76.6
branch 105 214 49.0
condition 46 83 55.4
subroutine 24 28 85.7
pod 9 9 100.0
total 447 677 66.0


line stmt bran cond sub pod time code
1             package Math::Prime::Util::RandomPrimes;
2 5     5   31 use strict;
  5         10  
  5         138  
3 5     5   21 use warnings;
  5         9  
  5         133  
4 5     5   22 use Carp qw/carp croak confess/;
  5         9  
  5         311  
5 5         38 use 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 5     5   28 /;
  5         9  
13              
14             BEGIN {
15 5     5   16 $Math::Prime::Util::RandomPrimes::AUTHORITY = 'cpan:DANAJ';
16 5         248 $Math::Prime::Util::RandomPrimes::VERSION = '0.69';
17             }
18              
19             BEGIN {
20 5 50   5   27 do { require Math::BigInt; Math::BigInt->import(try=>"GMP,Pari"); }
  0         0  
  0         0  
21             unless defined $Math::BigInt::VERSION;
22              
23 5     5   26 use constant OLD_PERL_VERSION=> $] < 5.008;
  5         12  
  5         398  
24 5     5   27 use constant MPU_MAXBITS => (~0 == 4294967295) ? 32 : 64;
  5         10  
  5         219  
25 5     5   25 use constant MPU_64BIT => MPU_MAXBITS == 64;
  5         8  
  5         204  
26 5     5   23 use constant MPU_32BIT => MPU_MAXBITS == 32;
  5         10  
  5         205  
27 5     5   26 use constant MPU_MAXPARAM => MPU_32BIT ? 4294967295 : 18446744073709551615;
  5         11  
  5         254  
28 5     5   28 use constant MPU_MAXDIGITS => MPU_32BIT ? 10 : 20;
  5         11  
  5         246  
29 5     5   25 use constant MPU_USE_XS => prime_get_config->{'xs'};
  5         12  
  5         23  
30 5     5   24 use constant MPU_USE_GMP => prime_get_config->{'gmp'};
  5         8  
  5         16  
31              
32 5         19880 *_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.
39             my @_big_gcd;
40             my $_big_gcd_top = 20046;
41             my $_big_gcd_use = -1;
42             sub _make_big_gcds {
43 3 50   3   10 return if $_big_gcd_use >= 0;
44 3 50       13 if (prime_get_config->{'gmp'}) {
45 0         0 $_big_gcd_use = 0;
46 0         0 return;
47             }
48 3 50       23 if (Math::BigInt->config()->{lib} !~ /^Math::BigInt::(GMP|Pari)/) {
49 3         181 $_big_gcd_use = 0;
50 3         8 return;
51             }
52 0         0 $_big_gcd_use = 1;
53 0         0 my $p0 = primorial(Math::BigInt->new( 520));
54 0         0 my $p1 = primorial(Math::BigInt->new(2052));
55 0         0 my $p2 = primorial(Math::BigInt->new(6028));
56 0         0 my $p3 = primorial(Math::BigInt->new($_big_gcd_top));
57 0         0 $_big_gcd[0] = $p0->bdiv(223092870)->bfloor->as_int;
58 0         0 $_big_gcd[1] = $p1->bdiv($p0)->bfloor->as_int;
59 0         0 $_big_gcd[2] = $p2->bdiv($p1)->bfloor->as_int;
60 0         0 $_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.
153             my $_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.
358             my @_random_ndigit_ranges = (undef, [2,7], [11,97] );
359             my @_random_nbit_ranges = (undef, undef, [2,3],[5,7] );
360             my %_random_cache_small;
361              
362             # For fixed small ranges with XS, e.g. 6-digit, 18-bit
363             sub _random_xscount_prime {
364 0     0   0 my($low,$high) = @_;
365 0         0 my($istart, $irange);
366 0         0 my $cachearef = $_random_cache_small{$low,$high};
367 0 0       0 if (defined $cachearef) {
368 0         0 ($istart, $irange) = @$cachearef;
369             } else {
370 0 0       0 my $beg = ($low <= 2) ? 2 : next_prime($low-1);
371 0 0       0 my $end = ($high < ~0) ? prev_prime($high + 1) : prev_prime($high);
372 0         0 ($istart, $irange) = ( prime_count(2, $beg), prime_count($beg, $end) );
373 0         0 $_random_cache_small{$low,$high} = [$istart, $irange];
374             }
375 0         0 my $rand = urandomm($irange);
376 0         0 return nth_prime($istart + $rand);
377             }
378              
379             sub random_prime {
380 2     2 1 7 my($low,$high) = @_;
381 2 50 33     7 return if $high < 2 || $low > $high;
382              
383 2 50       298 if ($high-$low > 1000000000) {
384             # Range is large, just make them odd if needed.
385 2 50       494 $low = 2 if $low < 2;
386 2 100 66     180 $low++ if $low > 2 && ($low % 2) == 0;
387 2 100       902 $high-- if ($high % 2) == 0;
388             } else {
389             # Tighten the range to the nearest prime.
390 0 0       0 $low = ($low <= 2) ? 2 : next_prime($low-1);
391 0 0       0 $high = ($high == ~0) ? prev_prime($high) : prev_prime($high + 1);
392 0 0 0     0 return $low if ($low == $high) && is_prob_prime($low);
393 0 0       0 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 2         968 return $_random_prime->($low, $high);
399             }
400              
401             sub random_ndigit_prime {
402 3     3 1 10 my($digits) = @_;
403 3 50       12 croak "random_ndigit_prime, digits must be >= 1" unless $digits >= 1;
404              
405 3 50 50     11 return _random_xscount_prime( int(10 ** ($digits-1)), int(10 ** $digits) )
406             if $digits <= 6 && MPU_USE_XS;
407              
408 3         6 my $bigdigits = $digits >= MPU_MAXDIGITS;
409 3 100 66     22 if ($bigdigits && prime_get_config->{'nobigint'}) {
410 1 50       3 croak "random_ndigit_prime with -nobigint, digits out of range"
411             if $digits > MPU_MAXDIGITS;
412             # Special case for nobigint and threshold digits
413 1 50       5 if (!defined $_random_ndigit_ranges[$digits]) {
414 1         4 my $low = int(10 ** ($digits-1));
415 1         2 my $high = ~0;
416 1         33 $_random_ndigit_ranges[$digits] = [next_prime($low),prev_prime($high)];
417             }
418             }
419              
420 3 100       15 if (!defined $_random_ndigit_ranges[$digits]) {
421 2 50       9 if ($bigdigits) {
422 2         11 my $low = Math::BigInt->new('10')->bpow($digits-1);
423 2         576 my $high = Math::BigInt->new('10')->bpow($digits);
424             # Just pull the range in to the nearest odd.
425 2         503 $_random_ndigit_ranges[$digits] = [$low+1, $high-1];
426             } else {
427 0         0 my $low = int(10 ** ($digits-1));
428 0         0 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 0         0 $_random_ndigit_ranges[$digits] = [next_prime($low),prev_prime($high)];
433             }
434             }
435 3         710 my ($low, $high) = @{$_random_ndigit_ranges[$digits]};
  3         9  
436 3         12 return $_random_prime->($low, $high);
437             }
438              
439             my @_random_nbit_m;
440             my @_random_nbit_lambda;
441             my @_random_nbit_arange;
442              
443             sub random_nbit_prime {
444 13     13 1 40 my($bits) = @_;
445 13 50       48 croak "random_nbit_prime, bits must be >= 2" unless $bits >= 2;
446 13         41 $bits = int("$bits");
447              
448             # Very small size, use the nth-prime method
449 13 50 50     52 if ($bits <= 20 && MPU_USE_XS) {
450 0 0       0 if ($bits <= 4) {
451 0 0       0 return (2,3)[urandomb(1)] if $bits == 2;
452 0 0       0 return (5,7)[urandomb(1)] if $bits == 3;
453 0 0       0 return (11,13)[urandomb(1)] if $bits == 4;
454             }
455 0         0 return _random_xscount_prime( 1 << ($bits-1), 1 << $bits );
456             }
457              
458 13         28 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 13 100       43 if (1 && $bits > 64) {
479 6 100       20 my $l = (MPU_64BIT && $bits > 79) ? 63 : 31;
480 6 50 100     33 $l = 49 if $l == 63 && OLD_PERL_VERSION; # Fix for broken Perl 5.6
481 6 50       22 $l = $bits-2 if $bits-2 < $l;
482              
483 6         45 my $brand = urandomb($bits-$l-2);
484 6 50       46 $brand = Math::BigInt->new("$brand") unless ref($brand) eq 'Math::BigInt';
485 6         322 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 6         1282 my @premod;
490 6         23 my $bpremod = _bigint_to_int($b->copy->bmod(9699690));
491 6         141 my $twopremod = _bigint_to_int(Math::BigInt->new(2)->bmodpow($bits-$l-1, 9699690));
492 6         106 foreach my $zi (0 .. 19-1) {
493 114         148 foreach my $pm (3, 5, 7, 11, 13, 17, 19) {
494 798 100 100     1553 next if $zi >= $pm || defined $premod[$pm];
495 282 100       540 $premod[$pm] = $zi if ( ($twopremod*$zi+$bpremod) % $pm) == 0;
496             }
497             }
498 6 100       23 _make_big_gcds() if $_big_gcd_use < 0;
499 6         12 if (!MPU_USE_GMP) { require Math::Prime::Util::PP; }
  6         37  
500              
501 6         16 my $loop_limit = 1_000_000;
502 6         25 while ($loop_limit-- > 0) {
503 203         11400 my $a = (1 << $l) + urandomb($l);
504             # $a % s == $premod[s] => $p % s == 0 => p will be composite
505 203 100 100     1607 next if $a % 3 == $premod[ 3] || $a % 5 == $premod[ 5]
      100        
      100        
      100        
      100        
      100        
506             || $a % 7 == $premod[ 7] || $a % 11 == $premod[11]
507             || $a % 13 == $premod[13] || $a % 17 == $premod[17]
508             || $a % 19 == $premod[19];
509 75         464 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 75         32370 if (MPU_USE_GMP) {
513             next unless Math::Prime::Util::GMP::is_prime($p);
514             } else {
515 75 100       301 next unless Math::BigInt::bgcd($p, 1348781387) == 1; # 23-43
516 59 50 33     38167 if ($_big_gcd_use && $p > $_big_gcd_top) {
517 0 0       0 next unless Math::BigInt::bgcd($p, $_big_gcd[0]) == 1;
518 0 0       0 next unless Math::BigInt::bgcd($p, $_big_gcd[1]) == 1;
519 0 0       0 next unless Math::BigInt::bgcd($p, $_big_gcd[2]) == 1;
520 0 0       0 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 59 100       223 next unless Math::Prime::Util::PP::is_bpsw_prime($p);
524             }
525 6         812 return $p;
526             }
527 0         0 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 7         16 if (1) {
534              
535 7         16 my $loop_limit = 2_000_000;
536 7 50       25 if ($bits > MPU_MAXBITS) {
537 0         0 my $p = Math::BigInt->bone->blsft($bits-1)->binc();
538 0         0 while ($loop_limit-- > 0) {
539 0         0 my $n = Math::BigInt->new(''.urandomb($bits-2))->blsft(1)->badd($p);
540 0 0       0 return $n if is_prob_prime($n);
541             }
542             } else {
543 7         21 my $p = (1 << ($bits-1)) + 1;
544 7         28 while ($loop_limit-- > 0) {
545 109         281 my $n = $p + (urandomb($bits-2) << 1);
546 109 100       413 return $n if is_prob_prime($n);
547             }
548             }
549 0         0 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.
580             sub _strip_proof_header {
581 0     0   0 my $proof = shift;
582 0         0 $proof =~ s/^\[MPU - Primality Certificate\]\nVersion \S+\n+Proof for:\nN (\d+)\n+//ms;
583 0         0 return $proof;
584             }
585              
586              
587             sub random_maurer_prime {
588 0     0 1 0 my $k = shift;
589 0 0       0 croak "random_maurer_prime, bits must be >= 2" unless $k >= 2;
590 0         0 $k = int("$k");
591              
592 0 0 0     0 return random_nbit_prime($k) if $k <= MPU_MAXBITS && !OLD_PERL_VERSION;
593              
594 0         0 my ($n, $cert) = random_maurer_prime_with_cert($k);
595 0 0       0 croak "maurer prime $n failed certificate verification!"
596             unless verify_prime($cert);
597 0         0 return $n;
598             }
599              
600             sub random_maurer_prime_with_cert {
601 10     10 1 44 my $k = shift;
602 10 50       46 croak "random_maurer_prime, bits must be >= 2" unless $k >= 2;
603 10         40 $k = int("$k");
604              
605             # This should never happen. Trap now to prevent infinite loop.
606 10 50       46 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 10         25 my $p0 = MPU_MAXBITS;
610 10         19 $p0 = 49 if OLD_PERL_VERSION && MPU_MAXBITS > 49;
611              
612 10 100       39 if ($k <= $p0) {
613 5         25 my $n = random_nbit_prime($k);
614 5         36 my ($isp, $cert) = is_provable_prime_with_cert($n);
615 5 50       21 croak "small nbit prime could not be proven" if $isp != 2;
616 5         16 return ($n, $cert);
617             }
618              
619             # Set verbose to 3 to get pretty output like Crypt::Primes
620 5         28 my $verbose = prime_get_config->{'verbose'};
621 5 50       33 local $| = 1 if $verbose > 2;
622              
623 5 100       22 do { require Math::BigFloat; Math::BigFloat->import(); }
  2         1826  
  2         46989  
624             if !defined $Math::BigFloat::VERSION;
625              
626             # Ignore Maurer's g and c that controls how much trial division is done.
627 5         1971 my $r = Math::BigFloat->new("0.5"); # relative size of the prime q
628 5         1572 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 5 50       26 if ($k > 2*$m) {
633 5         13 do {
634 6         50340 my $s = Math::Prime::Util::drand();
635 6         26 $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 5         356266 my $smallk = int(($r * $k)->bfloor->bstr) + 1;
642 5         2491 my ($q, $qcert) = random_maurer_prime_with_cert($smallk);
643 5 50       48 $q = Math::BigInt->new("$q") unless ref($q) eq 'Math::BigInt';
644 5         348 my $I = Math::BigInt->new(2)->bpow($k-2)->bdiv($q)->bfloor->as_int();
645 5 50 33     3464 print "r = $r k = $k q = $q I = $I\n" if $verbose && $verbose != 3;
646 5 50       23 $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 5 100       543 _make_big_gcds() if $_big_gcd_use < 0;
651 5         23 my $ONE = Math::BigInt->bone;
652 5         164 my $TWO = $ONE->copy->binc;
653              
654 5         299 my $loop_limit = 1_000_000 + $k * 1_000;
655 5         27 while ($loop_limit-- > 0) {
656             # R is a random number between $I+1 and 2*$I
657             #my $R = $I + 1 + urandomm( $I );
658 91         20724 my $R = $I->copy->binc->badd( urandomm($I) );
659             #my $n = 2 * $R * $q + 1;
660 91         22036 my $nm1 = $TWO->copy->bmul($R)->bmul($q);
661 91         13449 my $n = $nm1->copy->binc;
662             # We constructed a promising looking $n. Now test it.
663 91 50       4145 print "." if $verbose > 2;
664 91         134 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 91 100       273 next unless Math::BigInt::bgcd($n, 111546435)->is_one;
670 32 50 33     12332 if ($_big_gcd_use && $n > $_big_gcd_top) {
671 0 0       0 next unless Math::BigInt::bgcd($n, $_big_gcd[0])->is_one;
672 0 0       0 next unless Math::BigInt::bgcd($n, $_big_gcd[1])->is_one;
673 0 0       0 next unless Math::BigInt::bgcd($n, $_big_gcd[2])->is_one;
674 0 0       0 next unless Math::BigInt::bgcd($n, $_big_gcd[3])->is_one;
675             }
676 32 50       112 print "+" if $verbose > 2;
677 32 100       213 next unless is_strong_pseudoprime($n, 3);
678             }
679 5 50       27 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 5         30 my $m = $TWO * $R;
696 5 50 33     538 if (! ($q->is_odd && $q > 2 && $m > 0 &&
      33        
      33        
      33        
697             $m * $q + $ONE == $n && $TWO*$q+$ONE > $n->copy->bsqrt()) ) {
698 0         0 carp "Maurer prime failed BLS75 theorem 3 conditions. Retry.";
699 0         0 next;
700             }
701             # Find a suitable a. Move on if one isn't found quickly.
702 5         7997 foreach my $trya (2, 3, 5, 7, 11, 13) {
703 11         162596 my $a = Math::BigInt->new($trya);
704             # m/2 = R (n-1)/2 = (2*R*q)/2 = R*q
705 11 50       743 next unless $a->copy->bmodpow($R, $n) != $nm1;
706 11 100       87630 next unless $a->copy->bmodpow($R*$q, $n) == $nm1;
707 5 50       123282 print "($k)" if $verbose > 2;
708 5 50       57 croak "Maurer prime $n=2*$R*$q+1 failed BPSW" unless is_prob_prime($n);
709 5         767 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 5         555 return ($n, $cert);
714             }
715             # Didn't pass the selected a values. Try another R.
716             }
717 0         0 croak "Failure in random_maurer_prime, could not find a prime\n";
718             } # End of random_maurer_prime
719              
720              
721             sub random_shawe_taylor_prime_with_cert {
722 2     2 1 7 my $k = shift;
723              
724 2         65 my $seed = random_bytes(512/8);
725              
726 2         12 my($status,$prime,$prime_seed,$prime_gen_counter,$cert)
727             = _ST_Random_prime($k, $seed);
728 2 50       18 croak "Shawe-Taylor random prime failure" unless $status;
729 2 50       19 croak "Shawe-Taylor random prime failure: prime $prime failed certificate verification!" unless verify_prime($cert);
730              
731 2         21 return ($prime,$cert);
732             }
733              
734             sub _seed_plus_one {
735 57     57   115 my($s) = @_;
736 57         159 for (my $i = length($s)-1; $i >= 0; $i--) {
737 58         213 vec($s, $i, 8)++;
738 58 100       159 last unless vec($s, $i, 8) == 0;
739             }
740 57         145 return $s;
741             }
742              
743             sub _ST_Random_prime { # From FIPS 186-4
744 6     6   19 my($k, $input_seed) = @_;
745 6 50       18 croak "Shawe-Taylor random prime must have length >= 2" if $k < 2;
746 6         20 $k = int("$k");
747              
748 6 50 33     48 croak "Shawe-Taylor random prime, invalid input seed"
749             unless defined $input_seed && length($input_seed) >= 32;
750              
751 6 50       19 if (!defined $Digest::SHA::VERSION) {
752 0         0 eval { require Digest::SHA;
753 0         0 my $version = $Digest::SHA::VERSION;
754 0         0 $version =~ s/[^\d.]//g;
755 0         0 $version >= 4.00; }
756 0 0       0 or do { croak "Must have Digest::SHA 4.00 or later"; };
  0         0  
757             }
758              
759 6         23 my $k2 = Math::BigInt->new(2)->bpow($k-1);
760              
761 6 100       1622 if ($k < 33) {
762 2         5 my $seed = $input_seed;
763 2         6 my $prime_gen_counter = 0;
764 2         6 my $kmask = 0xFFFFFFFF >> (32-$k); # Does the mod operation
765 2         7 my $kstencil = (1 << ($k-1)) | 1; # Sets high and low bits
766 2         6 while (1) {
767 2         8 my $seedp1 = _seed_plus_one($seed);
768 2         42 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 2         23 my($c) = unpack("N*", substr($cvec,-4,4));
773 2         10 $c = ($c & $kmask) | $kstencil;
774 2         18 $prime_gen_counter++;
775 2         7 $seed = _seed_plus_one($seedp1);
776 2         10 my ($isp, $cert) = is_provable_prime_with_cert($c);
777 2 50       16 return (1,$c,$seed,$prime_gen_counter,$cert) if $isp;
778 0 0       0 return (0,0,0,0) if $prime_gen_counter > 10000 + 16*$k;
779             }
780             }
781 4         45 my($status,$c0,$seed,$prime_gen_counter,$cert)
782             = _ST_Random_prime( (($k+1)>>1)+1, $input_seed);
783 4 50       17 return (0,0,0,0) unless $status;
784 4 50       18 $cert = ($c0 < Math::BigInt->new("18446744073709551615"))
785             ? "" : _strip_proof_header($cert);
786 4         543 my $iterations = int(($k + 255) / 256) - 1; # SHA256 generates 256 bits
787 4         9 my $old_counter = $prime_gen_counter;
788 4         8 my $xstr = '';
789 4         13 for my $i (0 .. $iterations) {
790 4         67 $xstr = Digest::SHA::sha256_hex($seed) . $xstr;
791 4         14 $seed = _seed_plus_one($seed);
792             }
793 4         23 my $x = Math::BigInt->from_hex('0x'.$xstr);
794 4         3753 $x = $k2 + ($x % $k2);
795 4         1656 my $t = ($x + 2*$c0 - 1) / (2*$c0);
796 4 50       2667 _make_big_gcds() if $_big_gcd_use < 0;
797 4         13 while (1) {
798 49 50       1739 if (2*$t*$c0 + 1 > 2*$k2) { $t = ($k2 + 2*$c0 - 1) / (2*$c0); }
  0         0  
799 49         27176 my $c = 2*$t*$c0 + 1;
800 49         19005 $prime_gen_counter++;
801              
802             # Don't do the Pocklington check unless the candidate looks prime
803 49         81 my $looks_prime = 0;
804 49         65 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 49         124 $looks_prime = Math::BigInt::bgcd($c, 111546435)->is_one;
810 49 50 66     17266 if ($looks_prime && $_big_gcd_use && $c > $_big_gcd_top) {
      33        
811 0   0     0 $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 49 100 100     182 $looks_prime = 0 if $looks_prime && !is_strong_pseudoprime($c, 3);
817             }
818              
819 49 100       1051 if ($looks_prime) {
820             # We could use a in (2,3,5,7,11,13), but pedantically use FIPS 186-4.
821 4         14 my $astr = '';
822 4         14 for my $i (0 .. $iterations) {
823 4         57 $astr = Digest::SHA::sha256_hex($seed) . $astr;
824 4         18 $seed = _seed_plus_one($seed);
825             }
826 4         42 my $a = Math::BigInt->from_hex('0x'.$astr);
827 4         3742 $a = ($a % ($c-3)) + 2;
828 4         2951 my $z = $a->copy->bmodpow(2*$t,$c);
829 4 50 33     35034 if (Math::BigInt::bgcd($z-1,$c)->is_one && $z->copy->bmodpow($c0,$c)->is_one) {
830 4 50       48549 croak "Shawe-Taylor random prime failure at ($k): $c not prime"
831             unless is_prob_prime($c);
832 4         556 $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 4         435 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 45         103 for my $i (0 .. $iterations) {
841 45         86 $seed = _seed_plus_one($seed);
842             }
843             }
844 45 50       108 return (0,0,0,0) if $prime_gen_counter > 10000 + 16*$k + $old_counter;
845 45         119 $t++;
846             }
847             }
848              
849              
850             # Gordon's algorithm for generating a strong prime.
851             sub random_strong_prime {
852 1     1 1 3 my $t = shift;
853 1 50       3 croak "random_strong_prime, bits must be >= 128" unless $t >= 128;
854 1         3 $t = int("$t");
855              
856 1         2 croak "Random strong primes must be >= 173 bits on old Perl"
857             if OLD_PERL_VERSION && MPU_64BIT && $t < 173;
858              
859 1         3 my $l = (($t+1) >> 1) - 2;
860 1         4 my $lp = int($t/2) - 20;
861 1         4 my $lpp = $l - 20;
862 1         2 while (1) {
863 1         4 my $qp = random_nbit_prime($lp);
864 1         6 my $qpp = random_nbit_prime($lpp);
865 1 50       6 $qp = Math::BigInt->new("$qp") unless ref($qp) eq 'Math::BigInt';
866 1 50       5 $qpp = Math::BigInt->new("$qpp") unless ref($qpp) eq 'Math::BigInt';
867 1         5 my ($il, $rem) = Math::BigInt->new(2)->bpow($l-1)->bdec()->bdiv(2*$qpp);
868 1 50       764 $il++ if $rem > 0;
869 1         188 $il = $il->as_int();
870 1         23 my $iu = Math::BigInt->new(2)->bpow($l)->bsub(2)->bdiv(2*$qpp)->as_int();
871 1         802 my $istart = $il + urandomm($iu - $il + 1);
872 1         413 for (my $i = $istart; $i <= $iu; $i++) { # Search for q
873 4         3415 my $q = 2 * $i * $qpp + 1;
874 4 100       1440 next unless is_prob_prime($q);
875 1         139 my $pp = $qp->copy->bmodpow($q-2, $q)->bmul(2)->bmul($qp)->bdec();
876 1         31550 my ($jl, $rem) = Math::BigInt->new(2)->bpow($t-1)->bsub($pp)->bdiv(2*$q*$qp);
877 1 50       1171 $jl++ if $rem > 0;
878 1         224 $jl = $jl->as_int();
879 1         22 my $ju = Math::BigInt->new(2)->bpow($t)->bdec()->bsub($pp)->bdiv(2*$q*$qp)->as_int();
880 1         1075 my $jstart = $jl + urandomm($ju - $jl + 1);
881 1         419 for (my $j = $jstart; $j <= $ju; $j++) { # Search for p
882 130         144098 my $p = $pp + 2 * $j * $q * $qp;
883 130 100       65314 return $p if is_prob_prime($p);
884             }
885             }
886             }
887             }
888              
889             sub random_proven_prime {
890 0     0 1 0 my $k = shift;
891 0         0 my ($n, $cert) = random_proven_prime_with_cert($k);
892 0 0       0 croak "random_proven_prime $n failed certificate verification!"
893             unless verify_prime($cert);
894 0         0 return $n;
895             }
896              
897             sub random_proven_prime_with_cert {
898 1     1 1 2 my $k = shift;
899              
900 1 50 33     5 if (prime_get_config->{'gmp'} && $k <= 450) {
901 0         0 my $n = random_nbit_prime($k);
902 0         0 my ($isp, $cert) = is_provable_prime_with_cert($n);
903 0 0       0 croak "small nbit prime could not be proven" if $isp != 2;
904 0         0 return ($n, $cert);
905             }
906 1         7 return random_maurer_prime_with_cert($k);
907             }
908              
909             1;
910              
911             __END__