File Coverage

blib/lib/Math/Prime/Util/PrimalityProving.pm
Criterion Covered Total %
statement 497 567 87.6
branch 269 436 61.7
condition 59 144 40.9
subroutine 32 36 88.8
pod 4 4 100.0
total 861 1187 72.5


line stmt bran cond sub pod time code
1             package Math::Prime::Util::PrimalityProving;
2 5     5   1721 use strict;
  5         11  
  5         176  
3 5     5   30 use warnings;
  5         11  
  5         221  
4 5     5   28 use Carp qw/carp croak confess/;
  5         9  
  5         383  
5 5         31 use Math::Prime::Util qw/is_prob_prime is_strong_pseudoprime
6             is_provable_prime_with_cert
7             lucas_sequence
8             factor
9             prime_get_config
10 5     5   27 /;
  5         10  
11              
12             BEGIN {
13 5     5   18 $Math::Prime::Util::PrimalityProving::AUTHORITY = 'cpan:DANAJ';
14 5         211 $Math::Prime::Util::PrimalityProving::VERSION = '0.68';
15             }
16              
17             BEGIN {
18 5 50   5   29419 do { require Math::BigInt; Math::BigInt->import(try=>"GMP,Pari"); }
  0         0  
  0         0  
19             unless defined $Math::BigInt::VERSION;
20             }
21              
22             my $_smallval = Math::BigInt->new("18446744073709551615");
23             my $_maxint = Math::BigInt->new( (~0 > 4294967296 && $] < 5.008) ? "562949953421312" : ''.~0 );
24              
25              
26             ###############################################################################
27             # Pure Perl proofs
28             ###############################################################################
29              
30             my @_fsublist = (
31             sub { Math::Prime::Util::PP::pbrent_factor (shift, 32*1024, 1) },
32             sub { Math::Prime::Util::PP::pminus1_factor(shift, 1_000_000) },
33             sub { Math::Prime::Util::PP::ecm_factor (shift, 1_000, 5_000, 15) },
34             sub { Math::Prime::Util::PP::pbrent_factor (shift, 512*1024, 7) },
35             sub { Math::Prime::Util::PP::pminus1_factor(shift, 4_000_000) },
36             sub { Math::Prime::Util::PP::ecm_factor (shift, 10_000, 50_000, 10) },
37             sub { Math::Prime::Util::PP::pbrent_factor (shift, 512*1024, 11) },
38             sub { Math::Prime::Util::PP::pminus1_factor(shift,20_000_000) },
39             sub { Math::Prime::Util::PP::ecm_factor (shift, 100_000, 800_000, 10) },
40             sub { Math::Prime::Util::PP::pbrent_factor (shift, 2048*1024, 13) },
41             sub { Math::Prime::Util::PP::ecm_factor (shift, 1_000_000, 1_000_000, 20)},
42             sub { Math::Prime::Util::PP::pminus1_factor(shift, 100_000_000, 500_000_000)},
43             );
44              
45             sub _small_cert {
46 0     0   0 my $n = shift;
47 0 0       0 return '' unless is_prob_prime($n);
48 0         0 return join "\n", "[MPU - Primality Certificate]",
49             "Version 1.0",
50             "",
51             "Proof for:",
52             "N $n",
53             "",
54             "Type Small",
55             "N $n",
56             "";
57             }
58              
59             # For stripping off the header on certificates so they can be combined.
60             sub _strip_proof_header {
61 2     2   64 my $proof = shift;
62 2         17 $proof =~ s/^\[MPU - Primality Certificate\]\nVersion \S+\n+Proof for:\nN (\d+)\n+//ms;
63 2         10 return $proof;
64             }
65              
66              
67             sub primality_proof_lucas {
68 1     1 1 783 my ($n) = shift;
69 1         4 my @composite = (0, '');
70              
71             # Since this can take a very long time with a composite, try some easy cuts
72 1 50 33     8 return @composite if !defined $n || $n < 2;
73 1 50       4 return (2, _small_cert($n)) if $n < 4;
74 1 50       31 return @composite if is_strong_pseudoprime($n,2,15,325) == 0;
75              
76 1         3 my $nm1 = $n-1;
77 1         8 my @factors = factor($nm1);
78             { # remove duplicate factors and make a sorted array of bigints
79 1         2 my %uf;
  1         2  
80 1         6 undef @uf{@factors};
81 1         4 @factors = sort {$a<=>$b} map { Math::BigInt->new("$_") } keys %uf;
  5         128  
  4         130  
82             }
83 1         24 my $cert = "[MPU - Primality Certificate]\nVersion 1.0\n\nProof for:\nN $n\n\n";
84 1         4 $cert .= "Type Lucas\nN $n\n";
85 1         4 foreach my $i (1 .. scalar @factors) {
86 4         65 $cert .= "Q[$i] " . $factors[$i-1] . "\n";
87             }
88 1         19 for (my $a = 2; $a < $nm1; $a++) {
89 1         5 my $ap = Math::BigInt->new("$a");
90             # 1. a must be coprime to n
91 1 50       36 next unless Math::BigInt::bgcd($ap, $n) == 1;
92             # 2. a^(n-1) = 1 mod n.
93 1 50       295 next unless $ap->copy->bmodpow($nm1, $n) == 1;
94             # 3. a^((n-1)/f) != 1 mod n for all f.
95 4         704 next if (scalar grep { $_ == 1 }
96 1 50       775 map { $ap->copy->bmodpow(int($nm1/$_),$n); }
  4         2323  
97             @factors) > 0;
98             # Verify each factor and add to proof
99 1         84 my @fac_proofs;
100 1         3 foreach my $f (@factors) {
101 4         106 my ($isp, $fproof) = Math::Prime::Util::is_provable_prime_with_cert($f);
102 4 50       9 if ($isp != 2) {
103 0         0 carp "could not prove primality of $n.\n";
104 0         0 return (1, '');
105             }
106 4 50       10 push @fac_proofs, _strip_proof_header($fproof) if $f > $_smallval;
107             }
108 1         30 $cert .= "A $a\n";
109 1         4 foreach my $proof (@fac_proofs) {
110 0         0 $cert .= "\n$proof";
111             }
112 1         11 return (2, $cert);
113             }
114 0         0 return @composite;
115             }
116              
117             sub primality_proof_bls75 {
118 13     13 1 93 my ($n) = shift;
119 13         45 my @composite = (0, '');
120              
121             # Since this can take a very long time with a composite, try some easy tests
122 13 50 33     81 return @composite if !defined $n || $n < 2;
123 13 50       1174 return (2, _small_cert($n)) if $n < 4;
124 13 50       1003 return @composite if ($n & 1) == 0;
125 13 50       5008 return @composite if is_strong_pseudoprime($n,2,15,325) == 0;
126              
127 13         470 require Math::Prime::Util::PP;
128 13 100       66 $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
129 13         167 my $nm1 = $n->copy->bdec;
130 13         1130 my $ONE = $nm1->copy->bone;
131 13         1407 my $TWO = $ONE->copy->binc;
132 13         748 my $A = $ONE->copy; # factored part
133 13         220 my $B = $nm1->copy; # unfactored part
134 13         226 my @factors = ($TWO);
135 13 50       58 croak "BLS75 error: n-1 not even" unless $nm1->is_even();
136             {
137 13         199 while ($B->is_even) { $B->bdiv($TWO); $A->bmul($TWO); }
  13         41  
  19         482  
  19         3047  
138 13         776 my @tf;
139 13 100 66     57 if ($B <= $_maxint && prime_get_config->{'xs'}) {
140 2         8 @tf = Math::Prime::Util::trial_factor("$B", 20000);
141 2 100       141 pop @tf if $tf[-1] > 20000;
142             } else {
143 11         510 @tf = Math::Prime::Util::PP::trial_factor($B, 5000);
144 11 50       49 pop @tf if $tf[-1] > 5000;
145             }
146 13         451 foreach my $f (@tf) {
147 45 100       8643 next if $f == $factors[-1];
148 34         1578 push @factors, $f;
149 34         102 while (($B % $f) == 0) { $B /= $f; $A *= $f; }
  45         14770  
  45         9170  
150             }
151             }
152 13         4792 my @nstack;
153             # nstack should only hold composites
154 13 100       71 if ($B->is_one) {
    100          
155             # Completely factored. Nothing.
156             } elsif (is_prob_prime($B)) {
157 3         362 push @factors, $B;
158 3         8 $A *= $B; $B /= $B; # completely factored already
  3         338  
159             } else {
160 9         453 push @nstack, $B;
161             }
162 13         578 while (@nstack) {
163 14         72 my ($s,$r) = $B->copy->bdiv($A->copy->bmul($TWO));
164 14         4240 my $fpart = ($A+$ONE) * ($TWO*$A*$A + ($r-$ONE) * $A + $ONE);
165 14 100       8677 last if $n < $fpart;
166              
167 13         500 my $m = pop @nstack;
168             # Don't use bignum if it has gotten small enough.
169 13 100 100     80 $m = int($m->bstr) if ref($m) eq 'Math::BigInt' && $m <= $_maxint;
170             # Try to find factors of m, using the default set of factor subs.
171 13         517 my @ftry;
172 13         38 foreach my $sub (@_fsublist) {
173 13         49 @ftry = $sub->($m);
174 13 50       73 last if scalar @ftry >= 2;
175             }
176             # If we couldn't find a factor, skip it.
177 13 50       59 next unless scalar @ftry > 1;
178             # Process each factor
179 13         48 foreach my $f (@ftry) {
180 26 50 33     3405 croak "Invalid factoring: B=$B m=$m f=$f" if $f == 1 || $f == $m || !$B->copy->bmod($f)->is_zero;
      33        
181 26 100       7192 if (is_prob_prime($f)) {
182 21         785 push @factors, $f;
183 21         44 do { $B /= $f; $A *= $f; } while $B->copy->bmod($f)->is_zero;
  21         93  
  21         4610  
184             } else {
185 5         72 push @nstack, $f;
186             }
187             }
188             }
189             { # remove duplicate factors and make a sorted array of bigints
190 13         2976 my %uf = map { $_ => 1 } @factors;
  13         46  
  71         177  
191 13         679 @factors = sort {$a<=>$b} map { Math::BigInt->new("$_") } keys %uf;
  110         2425  
  71         2165  
192             }
193             # Just in case:
194 13         336 foreach my $f (@factors) {
195 71         5650 while ($B->copy->bmod($f)->is_zero) {
196 0         0 $B /= $f; $A *= $f;
  0         0  
197             }
198             }
199             # Did we factor enough?
200 13         1401 my ($s,$r) = $B->copy->bdiv($A->copy->bmul($TWO));
201 13         3383 my $fpart = ($A+$ONE) * ($TWO*$A*$A + ($r-$ONE) * $A + $ONE);
202 13 50       10140 return (1,'') if $n >= $fpart;
203             # Check we didn't mess up
204 13 50       539 croak "BLS75 error: $A * $B != $nm1" unless $A*$B == $nm1;
205 13 50       1631 croak "BLS75 error: $A not even" unless $A->is_even();
206 13 50       238 croak "BLS75 error: A and B not coprime" unless Math::BigInt::bgcd($A, $B)->is_one;
207              
208 13         3547 my $rtest = $r*$r - 8*$s;
209 13         3632 my $rtestroot = $rtest->copy->bsqrt;
210 13 50 66     1712 return @composite if $s != 0 && ($rtestroot*$rtestroot) == $rtest;
211              
212 13         2404 my $cert = "[MPU - Primality Certificate]\nVersion 1.0\n\nProof for:\nN $n\n\n";
213 13         493 $cert .= "Type BLS5\nN $n\n";
214 13         323 my $qnum = 0;
215 13         29 my $atext = '';
216 13         27 my @fac_proofs;
217 13         34 foreach my $f (@factors) {
218 71         2447 my $success = 0;
219 71 100       267 if ($qnum == 0) {
220 13 50       37 die "BLS5 Perl proof: Internal error, first factor not 2" unless $f == 2;
221             } else {
222 58         270 $cert .= "Q[$qnum] $f\n";
223             }
224 71         2612 my $nm1_div_f = $nm1 / $f;
225 71         15017 foreach my $a (2 .. 10000) {
226 81         202625 my $ap = Math::BigInt->new($a);
227 81 50       3253 next unless $ap->copy->bmodpow($nm1, $n)->is_one;
228 81 100       2245639 next unless Math::BigInt::bgcd($ap->copy->bmodpow($nm1_div_f, $n)->bdec, $n)->is_one;
229 71 100       1950781 $atext .= "A[$qnum] $a\n" unless $a == 2;
230 71         203 $success = 1;
231 71         311 last;
232             }
233 71         175 $qnum++;
234 71 50       185 return @composite unless $success;
235 71         506 my ($isp, $fproof) = is_provable_prime_with_cert($f);
236 71 50       246 if ($isp != 2) {
237 0         0 carp "could not prove primality of $n.\n";
238 0         0 return (1, '');
239             }
240 71 100       258 push @fac_proofs, _strip_proof_header($fproof) if $f > $_smallval;
241             }
242 13         520 $cert .= $atext;
243 13         45 $cert .= "----\n";
244 13         36 foreach my $proof (@fac_proofs) {
245 2         7 $cert .= "\n$proof";
246             }
247 13         315 return (2, $cert);
248             }
249              
250             ###############################################################################
251             # Convert certificates from old array format to new string format
252             ###############################################################################
253              
254             sub _convert_cert {
255 38     38   65 my $pdata = shift; # pdata is a ref
256              
257 38 50       98 return '' if scalar @$pdata == 0;
258 38         67 my $n = shift @$pdata;
259 38 100       100 if (length($n) == 1) {
260 3 100       15 return "Type Small\nN $n\n" if $n =~ /^[2357]$/;
261 2         5 return '';
262             }
263 35 50       164 $n = Math::BigInt->new("$n") if ref($n) ne 'Math::BigInt';
264 35 100       1742 return '' if $n->is_even;
265              
266 34 100       419 my $method = (scalar @$pdata > 0) ? shift @$pdata : 'BPSW';
267              
268 34 100       86 if ($method eq 'BPSW') {
269 3 100       8 return '' if $n > $_smallval;
270 2 50       62 return '' if is_prob_prime($n) != 2;
271 0         0 return "Type Small\nN $n\n";
272             }
273              
274 31 100 66     117 if ($method eq 'Pratt' || $method eq 'Lucas') {
275 12 50 33     80 if (scalar @$pdata != 2 || ref($$pdata[0]) ne 'ARRAY' || ref($$pdata[1]) eq 'ARRAY') {
      33        
276 0         0 carp "verify_prime: incorrect Pratt format, must have factors and a value\n";
277 0         0 return '';
278             }
279 12         20 my @factors = @{shift @$pdata};
  12         29  
280 12         21 my $a = shift @$pdata;
281 12         36 my $cert = "Type Lucas\nN $n\n";
282 12         248 foreach my $i (0 .. $#factors) {
283 49 100       90 my $f = (ref($factors[$i]) eq 'ARRAY') ? $factors[$i]->[0] : $factors[$i];
284 49         135 $cert .= sprintf("Q[%d] %s\n", $i+1, $f);
285             }
286 12         30 $cert .= "A $a\n\n";
287              
288 12         17 foreach my $farray (@factors) {
289 49 100       85 if (ref($farray) eq 'ARRAY') {
290 4         8 $cert .= _convert_cert($farray);
291             }
292             }
293 12         43 return $cert;
294             }
295 19 100       45 if ($method eq 'n-1') {
296 12 0 33     38 if (scalar @$pdata == 3 && ref($$pdata[0]) eq 'ARRAY' && $$pdata[0]->[0] =~ /^(B|T7|Theorem\s*7)$/i) {
      33        
297 0         0 croak "Unsupported BLS7 proof in conversion";
298             }
299 12 50 33     77 if (scalar @$pdata != 2 || ref($$pdata[0]) ne 'ARRAY' || ref($$pdata[1]) ne 'ARRAY') {
      33        
300 0         0 carp "verify_prime: incorrect n-1 format, must have factors and a values\n";
301 0         0 return '';
302             }
303 12         24 my @factors = @{shift @$pdata};
  12         28  
304 12         19 my @as = @{shift @$pdata};
  12         21  
305 12 50       34 if (scalar @factors != scalar @as) {
306 0         0 carp "verify_prime: incorrect n-1 format, must have a value for each factor\n";
307 0         0 return '';
308             }
309             # Make sure 2 is at the top
310 12         32 foreach my $i (1 .. $#factors) {
311 38 100       79 my $f = (ref($factors[$i]) eq 'ARRAY') ? $factors[$i]->[0] : $factors[$i];
312 38 50       79 if ($f == 2) {
313 0         0 my $tf = $factors[0]; $factors[0] = $factors[$i]; $factors[$i] = $tf;
  0         0  
  0         0  
314 0         0 my $ta = $as[0]; $as[0] = $as[$i]; $as[$i] = $ta;
  0         0  
  0         0  
315             }
316             }
317 12 100       58 return '' unless $factors[0] == 2;
318 11         37 my $cert = "Type BLS5\nN $n\n";
319 11         278 foreach my $i (1 .. $#factors) {
320 35 100       81 my $f = (ref($factors[$i]) eq 'ARRAY') ? $factors[$i]->[0] : $factors[$i];
321 35         97 $cert .= sprintf("Q[%d] %s\n", $i, $f);
322             }
323 11         27 foreach my $i (0 .. $#as) {
324 46 100       89 $cert .= sprintf("A[%d] %s\n", $i, $as[$i]) if $as[$i] != 2;
325             }
326 11         21 $cert .= "----\n\n";
327 11         21 foreach my $farray (@factors) {
328 46 100       80 if (ref($farray) eq 'ARRAY') {
329 3         38 $cert .= _convert_cert($farray);
330             }
331             }
332 11         79 return $cert;
333             }
334 7 50 33     28 if ($method eq 'ECPP' || $method eq 'AGKM') {
335 7 50       19 if (scalar @$pdata < 1) {
336 0         0 carp "verify_prime: incorrect AGKM format\n";
337 0         0 return '';
338             }
339 7         10 my $cert = '';
340 7         15 my $q = $n;
341 7         14 foreach my $block (@$pdata) {
342 14 50 33     70 if (ref($block) ne 'ARRAY' || scalar @$block != 6) {
343 0         0 carp "verify_prime: incorrect AGKM block format\n";
344 0         0 return '';
345             }
346 14         45 my($ni, $a, $b, $m, $qval, $P) = @$block;
347 14 50       43 if (Math::BigInt->new("$ni") != Math::BigInt->new("$q")) {
348 0         0 carp "verify_prime: incorrect AGKM block format: block n != q\n";
349 0         0 return '';
350             }
351 14 100       1546 $q = ref($qval) eq 'ARRAY' ? $qval->[0] : $qval;
352 14 50 33     63 if (ref($P) ne 'ARRAY' || scalar @$P != 2) {
353 0         0 carp "verify_prime: incorrect AGKM block point format\n";
354 0         0 return '';
355             }
356 14         24 my ($x, $y) = @{$P};
  14         28  
357 14         98 $cert .= "Type ECPP\nN $ni\nA $a\nB $b\nM $m\nQ $q\nX $x\nY $y\n\n";
358 14 100       45 if (ref($qval) eq 'ARRAY') {
359 1         3 $cert .= _convert_cert($qval);
360             }
361             }
362 7         25 return $cert;
363             }
364 0         0 carp "verify_prime: Unknown method: '$method'.\n";
365 0         0 return '';
366             }
367              
368              
369             sub convert_array_cert_to_string {
370 31     31 1 84 my @pdata = @_;
371              
372             # Convert reference input to array
373 31 100 66     174 @pdata = @{$pdata[0]} if scalar @pdata == 1 && ref($pdata[0]) eq 'ARRAY';
  29         83  
374              
375 31 100       92 return '' if scalar @pdata == 0;
376              
377 30         59 my $n = $pdata[0];
378 30         87 my $header = "[MPU - Primality Certificate]\nVersion 1.0\n\nProof for:\nN $n\n\n";
379              
380 30         82 my $cert = _convert_cert(\@pdata);
381 30 100       183 return '' if $cert eq '';
382 25         154 return $header . $cert;
383             }
384              
385              
386             ###############################################################################
387             # Verify certificate
388             ###############################################################################
389              
390             sub _primality_error ($) { ## no critic qw(ProhibitSubroutinePrototypes)
391 0 0   0   0 print "primality fail: $_[0]\n" if prime_get_config->{'verbose'};
392 0         0 return; # error in certificate
393             }
394              
395             sub _pfail ($) { ## no critic qw(ProhibitSubroutinePrototypes)
396 19 50   19   15833 print "primality fail: $_[0]\n" if prime_get_config->{'verbose'};
397 19         156 return; # Failed a condition
398             }
399              
400             sub _read_vars {
401 72     72   125 my $lines = shift;
402 72         131 my $type = shift;
403 72         168 my %vars = map { $_ => 1 } @_;
  178         499  
404 72         155 my %return;
405 72         218 while (scalar keys %vars) {
406 178         285 my $line = shift @$lines;
407 178 50       324 return _primality_error("end of file during type $type") unless defined $line;
408             # Skip comments and blank lines
409 178 50 33     688 next if $line =~ /^\s*#/ or $line =~ /^\s*$/;
410 178         297 chomp($line);
411 178 50       370 return _primality_error("Still missing values in type $type") if $line =~ /^Type /;
412 178 50       557 if ($line =~ /^(\S+)\s+(-?\d+)/) {
413 178         435 my ($var, $val) = ($1, $2);
414 178         292 $var =~ tr/a-z/A-Z/;
415 178 50       369 return _primality_error("Type $type: repeated or inappropriate var: $line") unless defined $vars{$var};
416 178         302 $return{$var} = $val;
417 178         448 delete $vars{$var};
418             } else {
419 0         0 return _primality_error("Unrecognized line: $line");
420             }
421             }
422             # Now return them in the order given, turned into bigints.
423 72         143 return map { Math::BigInt->new("$return{$_}") } @_;
  178         4411  
424             }
425              
426             sub _is_perfect_square {
427 4     4   2258 my($n) = @_;
428              
429 4 50       17 if (ref($n) eq 'Math::BigInt') {
430 4         14 my $mc = int(($n & 31)->bstr);
431 4 50 66     1230 if ($mc==0||$mc==1||$mc==4||$mc==9||$mc==16||$mc==17||$mc==25) {
      66        
      100        
      66        
      100        
      66        
432 4         12 my $sq = $n->copy->bsqrt->bfloor;
433 4         7773 $sq->bmul($sq);
434 4 50       632 return 1 if $sq == $n;
435             }
436             } else {
437 0         0 my $mc = $n & 31;
438 0 0 0     0 if ($mc==0||$mc==1||$mc==4||$mc==9||$mc==16||$mc==17||$mc==25) {
      0        
      0        
      0        
      0        
      0        
439 0         0 my $sq = int(sqrt($n));
440 0 0       0 return 1 if ($sq*$sq) == $n;
441             }
442             }
443 4         179 0;
444             }
445              
446             # Calculate Jacobi symbol (M|N)
447             sub _jacobi {
448 1     1   22 my($n, $m) = @_;
449 1 50 33     3 return 0 if $m <= 0 || ($m % 2) == 0;
450 1         434 my $j = 1;
451 1 50       4 if ($n < 0) {
452 1         154 $n = -$n;
453 1 50       30 $j = -$j if ($m % 4) == 3;
454             }
455             # Split loop so we can reduce n/m to non-bigints after first iteration.
456 1 50       227 if ($n != 0) {
457 1         152 while (($n % 2) == 0) {
458 3         1751 $n >>= 1;
459 3 50 33     424 $j = -$j if ($m % 8) == 3 || ($m % 8) == 5;
460             }
461 1         771 ($n, $m) = ($m, $n);
462 1 50 33     3 $j = -$j if ($n % 4) == 3 && ($m % 4) == 3;
463 1         227 $n = $n % $m;
464 1 50 33     87 $n = int($n->bstr) if ref($n) eq 'Math::BigInt' && $n <= $_maxint;
465 1 50 33     54 $m = int($m->bstr) if ref($m) eq 'Math::BigInt' && $m <= $_maxint;
466             }
467 1         42 while ($n != 0) {
468 0         0 while (($n % 2) == 0) {
469 0         0 $n >>= 1;
470 0 0 0     0 $j = -$j if ($m % 8) == 3 || ($m % 8) == 5;
471             }
472 0         0 ($n, $m) = ($m, $n);
473 0 0 0     0 $j = -$j if ($n % 4) == 3 && ($m % 4) == 3;
474 0         0 $n = $n % $m;
475             }
476 1 50       5 return ($m == 1) ? $j : 0;
477             }
478              
479             # Proof handlers (parse input and call verification)
480              
481             sub _prove_ecpp {
482 13     13   49 _verify_ecpp( _read_vars($_[0], 'ECPP', qw/N A B M Q X Y/) );
483             }
484             sub _prove_ecpp3 {
485 1     1   6 _verify_ecpp3( _read_vars($_[0], 'ECPP3', qw/N S R A B T/) );
486             }
487             sub _prove_ecpp4 {
488 0     0   0 _verify_ecpp4( _read_vars($_[0], 'ECPP4', qw/N S R J T/) );
489             }
490             sub _prove_bls15 {
491 1     1   5 _verify_bls15( _read_vars($_[0], 'BLS15', qw/N Q LP LQ/) );
492             }
493             sub _prove_bls3 {
494 5     5   28 _verify_bls3( _read_vars($_[0], 'BLS3', qw/N Q A/) );
495             }
496             sub _prove_pock {
497 5     5   26 _verify_pock( _read_vars($_[0], 'POCKLINGTON', qw/N Q A/) );
498             }
499             sub _prove_small {
500 3     3   13 _verify_small( _read_vars($_[0], 'Small', qw/N/) );
501             }
502             sub _prove_bls5 {
503 16     16   34 my $lines = shift;
504             # No good way to do this using read_vars
505 16         34 my ($n, @Q, @A);
506 16         28 my $index = 0;
507 16         60 $Q[0] = Math::BigInt->new(2); # 2 is implicit
508 16         592 while (1) {
509 124         3647 my $line = shift @$lines;
510 124 50       228 return _primality_error("end of file during type BLS5") unless defined $line;
511             # Skip comments and blank lines
512 124 50 33     475 next if $line =~ /^\s*#/ or $line =~ /^\s*$/;
513             # Stop when we see a line starting with -.
514 124 100       226 last if $line =~ /^-/;
515 108         140 chomp($line);
516 108 100       387 if ($line =~ /^N\s+(\d+)/) {
    100          
    50          
517 16 50       42 return _primality_error("BLS5: N redefined") if defined $n;
518 16         64 $n = Math::BigInt->new("$1");
519             } elsif ($line =~ /^Q\[(\d+)\]\s+(\d+)/) {
520 64         88 $index++;
521 64 50       145 return _primality_error("BLS5: Invalid index: $1") unless $1 == $index;
522 64         171 $Q[$1] = Math::BigInt->new("$2");
523             } elsif ($line =~ /^A\[(\d+)\]\s+(\d+)/) {
524 28 50 33     105 return _primality_error("BLS5: Invalid index: A[$1]") unless $1 >= 0 && $1 <= $index;
525 28         73 $A[$1] = Math::BigInt->new("$2");
526             } else {
527 0         0 return _primality_error("Unrecognized line: $line");
528             }
529             }
530 16         56 _verify_bls5($n, \@Q, \@A);
531             }
532              
533             sub _prove_lucas {
534 11     11   18 my $lines = shift;
535             # No good way to do this using read_vars
536 11         22 my ($n, @Q, $a);
537 11         13 my $index = 0;
538 11         19 while (1) {
539 67         1789 my $line = shift @$lines;
540 67 50       114 return _primality_error("end of file during type Lucas") unless defined $line;
541             # Skip comments and blank lines
542 67 50 33     241 next if $line =~ /^\s*#/ or $line =~ /^\s*$/;
543 67         90 chomp($line);
544 67 100       222 if ($line =~ /^N\s+(\d+)/) {
    100          
    50          
545 11 50       24 return _primality_error("Lucas: N redefined") if defined $n;
546 11         38 $n = Math::BigInt->new("$1");
547             } elsif ($line =~ /^Q\[(\d+)\]\s+(\d+)/) {
548 45         58 $index++;
549 45 50       96 return _primality_error("Lucas: Invalid index: $1") unless $1 == $index;
550 45         112 $Q[$1] = Math::BigInt->new("$2");
551             } elsif ($line =~ /^A\s+(\d+)/) {
552 11         38 $a = Math::BigInt->new("$1");
553 11         390 last;
554             } else {
555 0         0 return _primality_error("Unrecognized line: $line");
556             }
557             }
558 11         29 _verify_lucas($n, \@Q, $a);
559             }
560              
561             # Verification routines
562              
563             sub _verify_ecpp {
564 14     14   1258 my ($n, $a, $b, $m, $q, $x, $y) = @_;
565 14 50       43 return unless defined $n;
566 14 50       47 $a %= $n if $a < 0;
567 14 50       1958 $b %= $n if $b < 0;
568 14 50       1851 return _pfail "ECPP: $n failed N > 0" unless $n > 0;
569 14 100       1913 return _pfail "ECPP: $n failed gcd(N, 6) = 1" unless Math::BigInt::bgcd($n, 6) == 1;
570 13 100       3380 return _pfail "ECPP: $n failed gcd(4*a^3 + 27*b^2, N) = 1"
571             unless Math::BigInt::bgcd(4*$a*$a*$a+27*$b*$b,$n) == 1;
572 12 50       11821 return _pfail "ECPP: $n failed Y^2 = X^3 + A*X + B mod N"
573             unless ($y*$y) % $n == ($x*$x*$x + $a*$x + $b) % $n;
574 12 100       8439 return _pfail "ECPP: $n failed M >= N - 2*sqrt(N) + 1" unless $m >= $n + 1 - $n->copy->bmul(4)->bsqrt();
575 11 50       9171 return _pfail "ECPP: $n failed M <= N + 2*sqrt(N) + 1" unless $m <= $n + 1 + $n->copy->bmul(4)->bsqrt();
576 11 100       8644 return _pfail "ECPP: $n failed Q > (N^(1/4)+1)^2" unless $q > $n->copy->broot(4)->badd(1)->bpow(2);
577 10 50       13732 return _pfail "ECPP: $n failed Q < N" unless $q < $n;
578 10 50       274 return _pfail "ECPP: $n failed M != Q" unless $m != $q;
579 10         284 my ($mdivq, $rem) = $m->copy->bdiv($q);
580 10 100       1793 return _pfail "ECPP: $n failed Q divides M" unless $rem == 0;
581              
582             # Now verify the elliptic curve
583 9         1204 my $correct_point = 0;
584 9 50 33     36 if (prime_get_config->{'gmp'} && defined &Math::Prime::Util::GMP::_validate_ecpp_curve) {
585 0         0 $correct_point = Math::Prime::Util::GMP::_validate_ecpp_curve($a, $b, $n, $x, $y, $m, $q);
586             } else {
587 9 100       29 if (!defined $Math::Prime::Util::ECAffinePoint::VERSION) {
588 1         540 eval { require Math::Prime::Util::ECAffinePoint; 1; }
  1         5  
589 1 50       2 or do { die "Cannot load Math::Prime::Util::ECAffinePoint"; };
  0         0  
590             }
591 9         50 my $ECP = Math::Prime::Util::ECAffinePoint->new($a, $b, $n, $x, $y);
592             # Compute U = (m/q)P, check U != point at infinity
593 9         25 $ECP->mul( $m->copy->bdiv($q)->as_int );
594 9 50       48 if (!$ECP->is_infinity) {
595             # Compute V = qU, check V = point at infinity
596 9         134 $ECP->mul( $q );
597 9 50       44 $correct_point = 1 if $ECP->is_infinity;
598             }
599             }
600 9 50       280 return _pfail "ECPP: $n failed elliptic curve conditions" unless $correct_point;
601 9         66 ($n, $q);
602             }
603             sub _verify_ecpp3 {
604 1     1   55 my ($n, $s, $r, $a, $b, $t) = @_;
605 1 50       6 return unless defined $n;
606 1 50       6 return _pfail "ECPP3: $n failed |A| <= N/2" unless abs($a) <= $n/2;
607 1 50       303 return _pfail "ECPP3: $n failed |B| <= N/2" unless abs($b) <= $n/2;
608 1 50       243 return _pfail "ECPP3: $n failed T >= 0" unless $t >= 0;
609 1 50       169 return _pfail "ECPP3: $n failed T < N" unless $t < $n;
610 1         52 my $l = ($t*$t*$t + $a*$t + $b) % $n;
611 1         382 _verify_ecpp( $n,
612             ($a * $l*$l) % $n,
613             ($b * $l*$l*$l) % $n,
614             $r*$s,
615             $r,
616             ($t*$l) % $n,
617             ($l*$l) % $n );
618             }
619             sub _verify_ecpp4 {
620 0     0   0 my ($n, $s, $r, $j, $t) = @_;
621 0 0       0 return unless defined $n;
622 0 0       0 return _pfail "ECPP4: $n failed |J| <= N/2" unless abs($j) <= $n/2;
623 0 0       0 return _pfail "ECPP4: $n failed T >= 0" unless $t >= 0;
624 0 0       0 return _pfail "ECPP4: $n failed T < N" unless $t < $n;
625 0         0 my $a = 3 * $j * (1728 - $j);
626 0         0 my $b = 2 * $j * (1728 - $j) * (1728 - $j);
627 0         0 my $l = ($t*$t*$t + $a*$t + $b) % $n;
628 0         0 _verify_ecpp( $n,
629             ($a * $l*$l) % $n,
630             ($b * $l*$l*$l) % $n,
631             $r*$s,
632             $r,
633             ($t*$l) % $n,
634             ($l*$l) % $n );
635             }
636              
637             sub _verify_bls15 {
638 1     1   32 my ($n, $q, $lp, $lq) = @_;
639 1 50       4 return unless defined $n;
640 1 50       5 return _pfail "BLS15: $n failed Q odd" unless $q->is_odd();
641 1 50       16 return _pfail "BLS15: $n failed Q > 2" unless $q > 2;
642 1         128 my ($m, $rem) = ($n+1)->copy->bdiv($q);
643 1 50       421 return _pfail "BLS15: $n failed Q divides N+1" unless $rem == 0;
644 1 50       133 return _pfail "BLS15: $n failed MQ-1 = N" unless $m*$q-1 == $n;
645 1 50       297 return _pfail "BLS15: $n failed M > 0" unless $m > 0;
646 1 50       161 return _pfail "BLS15: $n failed 2Q-1 > sqrt(N)" unless 2*$q-1 > $n->copy->bsqrt();
647 1         1145 my $D = $lp*$lp - 4*$lq;
648 1 50       274 return _pfail "BLS15: $n failed D != 0" unless $D != 0;
649 1 50       127 return _pfail "BLS15: $n failed jacobi(D,N) = -1" unless _jacobi($D,$n) == -1;
650 1 50       4 return _pfail "BLS15: $n failed V_{m/2} mod N != 0"
651             unless (lucas_sequence($n, $lp, $lq, $m/2))[1] != 0;
652 1 50       227 return _pfail "BLS15: $n failed V_{(N+1)/2} mod N == 0"
653             unless (lucas_sequence($n, $lp, $lq, ($n+1)/2))[1] == 0;
654 1         207 ($n, $q);
655             }
656              
657             sub _verify_bls3 {
658 5     5   170 my ($n, $q, $a) = @_;
659 5 50       57 return unless defined $n;
660 5 50       26 return _pfail "BLS3: $n failed Q odd" unless $q->is_odd();
661 5 50       111 return _pfail "BLS3: $n failed Q > 2" unless $q > 2;
662 5         518 my ($m, $rem) = ($n-1)->copy->bdiv($q);
663 5 50       2215 return _pfail "BLS3: $n failed Q divides N-1" unless $rem == 0;
664 5 50       735 return _pfail "BLS3: $n failed MQ+1 = N" unless $m*$q+1 == $n;
665 5 50       1448 return _pfail "BLS3: $n failed M > 0" unless $m > 0;
666 5 50       755 return _pfail "BLS3: $n failed 2Q+1 > sqrt(n)" unless 2*$q+1 > $n->copy->bsqrt();
667 5 50       5310 return _pfail "BLS3: $n failed A^((N-1)/2) = N-1 mod N" unless $a->copy->bmodpow(($n-1)/2, $n) == $n-1;
668 5 50       122982 return _pfail "BLS3: $n failed A^(M/2) != N-1 mod N" unless $a->copy->bmodpow($m/2,$n) != $n-1;
669 5         51933 ($n, $q);
670             }
671             sub _verify_pock {
672 5     5   206 my ($n, $q, $a) = @_;
673 5 50       20 return unless defined $n;
674 5         21 my ($m, $rem) = ($n-1)->copy->bdiv($q);
675 5 50       2213 return _pfail "Pocklington: $n failed Q divides N-1" unless $rem == 0;
676 5 50       691 return _pfail "Pocklington: $n failed M is even" unless $m->is_even();
677 5 50       83 return _pfail "Pocklington: $n failed M > 0" unless $m > 0;
678 5 50       749 return _pfail "Pocklington: $n failed M < Q" unless $m < $q;
679 5 50       171 return _pfail "Pocklington: $n failed MQ+1 = N" unless $m*$q+1 == $n;
680 5 50       1405 return _pfail "Pocklington: $n failed A > 1" unless $a > 1;
681 5 50       448 return _pfail "Pocklington: $n failed A^(N-1) mod N = 1"
682             unless $a->copy->bmodpow($n-1, $n) == 1;
683 5 50       158144 return _pfail "Pocklington: $n failed gcd(A^M - 1, N) = 1"
684             unless Math::BigInt::bgcd($a->copy->bmodpow($m, $n)-1, $n) == 1;
685 5         106631 ($n, $q);
686             }
687             sub _verify_small {
688 3     3   117 my ($n) = @_;
689 3 50       7 return unless defined $n;
690 3 50       7 return _pfail "Small n $n is > 2^64\n" if $n > $_smallval;
691 3 50       117 return _pfail "Small n $n does not pass BPSW" unless is_prob_prime($n);
692 3         132 ($n);
693             }
694              
695             sub _verify_bls5 {
696 16     16   34 my ($n, $Qr, $Ar) = @_;
697 16 50       54 return unless defined $n;
698 16         25 my @Q = @{$Qr};
  16         46  
699 16         34 my @A = @{$Ar};
  16         32  
700 16         52 my $nm1 = $n - 1;
701 16         3119 my $F = Math::BigInt->bone;
702 16         416 my $R = $nm1->copy;
703 16         284 my $index = $#Q;
704 16         47 foreach my $i (0 .. $index) {
705 77 50       25676 return _primality_error "BLS5: $n failed Q[$i] doesn't exist" unless defined $Q[$i];
706 77 100       206 $A[$i] = Math::BigInt->new(2) unless defined $A[$i];
707 77 50       1606 return _pfail "BLS5: $n failed Q[$i] > 1" unless $Q[$i] > 1;
708 77 50       6634 return _pfail "BLS5: $n failed Q[$i] < N-1" unless $Q[$i] < $nm1;
709 77 50       2162 return _pfail "BLS5: $n failed A[$i] > 1" unless $A[$i] > 1;
710 77 50       6376 return _pfail "BLS5: $n failed A[$i] < N" unless $A[$i] < $n;
711 77 100       2154 return _pfail "BLS5: $n failed Q[$i] divides N-1" unless ($nm1 % $Q[$i]) == 0;
712 76         19959 while (($R % $Q[$i]) == 0) {
713 88         22886 $F *= $Q[$i];
714 88         5001 $R /= $Q[$i];
715             }
716             }
717 15 50       5069 die "BLS5: Internal error R != (N-1)/F\n" unless $R == $nm1/$F;
718 15 50       3769 return _pfail "BLS5: $n failed F is even" unless $F->is_even();
719 15 50       196 return _pfail "BLS5: $n failed gcd(F, R) = 1\n" unless Math::BigInt::bgcd($F,$R) == 1;
720 15         16277 my ($s, $r) = $R->copy->bdiv(2*$F);
721 15         4962 my $P = ($F+1) * (2 * $F * $F + ($r-1)*$F + 1);
722 15 100       14647 return _pfail "BLS5: $n failed n < P" unless $n < $P;
723 13 50 66     401 return _pfail "BLS5: $n failed s=0 OR r^2-8s not a perfect square"
724             unless $s == 0 or !_is_perfect_square($r*$r - 8*$s);
725 13         1176 foreach my $i (0 .. $index) {
726 68         34338384 my $a = $A[$i];
727 68         163 my $q = $Q[$i];
728 68 50       232 return _pfail "BLS5: $n failed A[i]^(N-1) mod N = 1"
729             unless $a->copy->bmodpow($nm1, $n)->is_one;
730 68 100       48507451 return _pfail "BLS5: $n failed gcd(A[i]^((N-1)/Q[i])-1, N) = 1"
731             unless Math::BigInt::bgcd($a->copy->bmodpow($nm1/$q, $n)->bdec, $n)->is_one;
732             }
733 12         2764348 ($n, @Q);
734             }
735              
736             sub _verify_lucas {
737 11     11   25 my ($n, $Qr, $a) = @_;
738 11 50       23 return unless defined $n;
739 11         19 my @Q = @{$Qr};
  11         26  
740 11         16 my $index = $#Q;
741 11 50       30 return _pfail "Lucas: $n failed A > 1" unless $a > 1;
742 11 100       998 return _pfail "Lucas: $n failed A < N" unless $a < $n;
743 10         332 my $nm1 = $n - 1;
744 10         1784 my $F = Math::BigInt->bone;
745 10         242 my $R = $nm1->copy;
746 10 50       161 return _pfail "Lucas: $n failed A^(N-1) mod N = 1"
747             unless $a->copy->bmodpow($nm1, $n) == 1;
748 10         37749 foreach my $i (1 .. $index) {
749 26 50       5002 return _primality_error "Lucas: $n failed Q[$i] doesn't exist" unless defined $Q[$i];
750 26 50       61 return _pfail "Lucas: $n failed Q[$i] > 1" unless $Q[$i] > 1;
751 26 50       2323 return _pfail "Lucas: $n failed Q[$i] < N-1" unless $Q[$i] < $nm1;
752 26 100       708 return _pfail "Lucas: $n failed Q[$i] divides N-1" unless ($nm1 % $Q[$i]) == 0;
753 23 100       5153 return _pfail "Lucas: $n failed A^((N-1)/Q[$i]) mod N != 1"
754             unless $a->copy->bmodpow($nm1/$Q[$i], $n) != 1;
755 22         60857 while (($R % $Q[$i]) == 0) {
756 29         6806 $F *= $Q[$i];
757 29         1398 $R /= $Q[$i];
758             }
759             }
760 6 100 66     1840 return _pfail("Lucas: $n failed N-1 has only factors Q") unless $R == 1 && $F == $nm1;
761 5         568 shift @Q; # Remove Q[0]
762 5         30 ($n, @Q);
763             }
764              
765             sub verify_cert {
766 50 100   50 1 170 my $cert = (@_ == 1) ? $_[0] : convert_array_cert_to_string(@_);
767 50 100       179 $cert = convert_array_cert_to_string($cert) if ref($cert) eq 'ARRAY';
768 50 100       164 return 0 if $cert eq '';
769              
770 44         74 my %parts; # Map of "N is prime if Q is prime"
771 44         370 my %proof_funcs = (
772             ECPP => \&_prove_ecpp, # Standard ECPP proof
773             ECPP3 => \&_prove_ecpp3, # Primo type 3
774             ECPP4 => \&_prove_ecpp4, # Primo type 4
775             BLS15 => \&_prove_bls15, # basic n+1, includes Primo type 2
776             BLS3 => \&_prove_bls3, # basic n-1
777             BLS5 => \&_prove_bls5, # much better n-1
778             SMALL => \&_prove_small, # n <= 2^64
779             POCKLINGTON => \&_prove_pock, # simple n-1, Primo type 1
780             LUCAS => \&_prove_lucas, # n-1 completely factored
781             );
782 44         81 my $base = 10;
783 44         74 my $cert_type = 'Unknown';
784 44         79 my $N;
785              
786 44         253 my @lines = split /^/, $cert;
787 44         92 my $lines = \@lines;
788 44         116 while (@$lines) {
789 287         1079 my $line = shift @$lines;
790 287 100 66     1393 next if $line =~ /^\s*#/ or $line =~ /^\s*$/; # Skip comments / blank lines
791 187         369 chomp($line);
792 187 100       502 if ($line =~ /^\[(\S+) - Primality Certificate\]/) {
793 44 50       157 if ($1 ne 'MPU') {
794 0         0 return _primality_error "Unknown certificate type: $1";
795             }
796 44         90 $cert_type = $1;
797 44         105 next;
798             }
799 143 100 33     861 if ( ($cert_type eq 'PRIMO' && $line =~ /^\[Candidate\]/) || ($cert_type eq 'MPU' && $line =~ /^Proof for:/) ) {
      66        
      66        
800 44 50       122 return _primality_error "Certificate with multiple N values" if defined $N;
801 44         126 ($N) = _read_vars($lines, 'Proof for', qw/N/);
802 44 100       2425 if (!is_prob_prime($N)) {
803 2         116 _pfail "N '$N' does not look prime.";
804 2         14 return 0;
805             }
806 42         4234 next;
807             }
808 99 50       253 if ($line =~ /^Base (\d+)/) {
809 0         0 $base = $1;
810 0 0       0 return _primality_error "Only base 10 supported, sorry" unless $base == 10;
811 0         0 next;
812             }
813 99 100       405 if ($line =~ /^Type (.*?)\s*$/) {
814 55 50       170 return _primality_error("Starting type without telling me the N value!") unless defined $N;
815 55         151 my $type = $1;
816 55         121 $type =~ tr/a-z/A-Z/;
817 55 50       172 error("Unknown type: $type") unless defined $proof_funcs{$type};
818 55         157 my ($n, @q) = $proof_funcs{$type}->($lines);
819 55 100       392 return 0 unless defined $n;
820 40         173 $parts{$n} = [@q];
821             }
822             }
823              
824 27 50       854 return _primality_error("No N") unless defined $N;
825 27         103 my @qs = ($N);
826 27         89 while (@qs) {
827 126         3348 my $q = shift @qs;
828             # Check that this q has a chain
829 126 100       257 if (!defined $parts{$q}) {
830 90 50       1584 if ($q > $_smallval) {
831 0         0 _primality_error "q value $q has no proof\n";
832 0         0 return 0;
833             }
834 90 100       2949 if (!is_prob_prime($q)) {
835 2         68 _pfail "Small n $q does not pass BPSW";
836 2         29 return 0;
837             }
838             } else {
839 36 50       879 die "Internal error: Invalid parts entry" if ref($parts{$q}) ne 'ARRAY';
840             # q is prime if all it's chains are prime.
841 36         829 push @qs, @{$parts{$q}};
  36         77  
842             }
843             }
844 25         1777 1;
845             }
846              
847             1;
848              
849             __END__