Commit | Line | Data |
---|---|---|
7de4d3fe MW |
1 | #! /usr/bin/python |
2 | ### -*- mode: python, coding: utf-8 -*- | |
3 | ### | |
4 | ### Tool for generating and verifying primality certificates | |
5 | ### | |
6 | ### (c) 2017 Straylight/Edgeware | |
7 | ### | |
8 | ||
9 | ###----- Licensing notice --------------------------------------------------- | |
10 | ### | |
11 | ### This file is part of the Python interface to Catacomb. | |
12 | ### | |
13 | ### Catacomb/Python is free software; you can redistribute it and/or modify | |
14 | ### it under the terms of the GNU General Public License as published by | |
15 | ### the Free Software Foundation; either version 2 of the License, or | |
16 | ### (at your option) any later version. | |
17 | ### | |
18 | ### Catacomb/Python is distributed in the hope that it will be useful, | |
19 | ### but WITHOUT ANY WARRANTY; without even the implied warranty of | |
20 | ### MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
21 | ### GNU General Public License for more details. | |
22 | ### | |
23 | ### You should have received a copy of the GNU General Public License | |
24 | ### along with Catacomb/Python; if not, write to the Free Software Foundation, | |
25 | ### Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. | |
26 | ||
27 | ###-------------------------------------------------------------------------- | |
28 | ### Imported modules. | |
29 | ||
30 | from sys import argv, stdin, stdout, stderr | |
31 | import os as OS | |
32 | import itertools as I | |
33 | import math as M | |
34 | import optparse as OP | |
35 | ||
36 | import catacomb as C | |
37 | ||
38 | ###-------------------------------------------------------------------------- | |
39 | ### Utilities. | |
40 | ||
f6d012db MW |
41 | def _excval(): |
42 | """Return the most recent exception object.""" | |
43 | return SYS.exc_info()[1] | |
44 | ||
7de4d3fe MW |
45 | class ExpectedError (Exception): |
46 | """ | |
47 | I represent an expected error, which should be reported in a friendly way. | |
48 | """ | |
49 | pass | |
50 | ||
51 | def prod(ff, one = 1): | |
41bfac1d MW |
52 | """Return ONE times the product of the elements of FF.""" |
53 | return C.MPMul().factor(one).factor(ff).done() | |
7de4d3fe MW |
54 | |
55 | def parse_label(line): | |
56 | """ | |
57 | Split LINE at an `=' character and return the left and right halves. | |
58 | ||
59 | The returned pieces have leading and trailing whitespace trimmed. | |
60 | """ | |
61 | eq = line.find('=') | |
62 | if eq < 0: raise ExpectedError('expected `LABEL = ...\'') | |
63 | return line[:eq].strip(), line[eq + 1:].strip() | |
64 | ||
65 | def parse_list(s, n): | |
66 | l = s.split(',', n - 1) | |
67 | if n is not None and len(l) != n: | |
68 | raise ExpectedError('expected `,\'-separated list of %d items' % n) | |
69 | return l | |
70 | ||
71 | def conv_int(s): | |
72 | """Convert S to a integer.""" | |
73 | try: return C.MP(s, 0) | |
74 | except TypeError: raise ExpectedError('invalid integer `%s\'' % s) | |
75 | ||
76 | VERBOSITY = 1 | |
77 | ||
78 | class ProgressReporter (object): | |
79 | """ | |
80 | I keep users amused while the program looks for large prime numbers. | |
81 | ||
82 | My main strategy is the printing of incomprehensible runes. I can be | |
83 | muffled by lowering by verbosity level. | |
84 | ||
85 | Prime searches are recursive in nature. When a new recursive level is | |
86 | started, call `push'; and call `pop' when the level is finished. This must | |
87 | be done around the top level too. | |
88 | """ | |
89 | def __init__(me): | |
90 | """Initialize the ProgressReporter.""" | |
91 | me._level = -1 | |
92 | me._update() | |
93 | def _update(me): | |
94 | """ | |
95 | Update our idea of whether we're active. | |
96 | ||
97 | We don't write inscrutable runes when inactive. The current policy is to | |
98 | write nothing if verbosity is zero, to write runes for the top level only | |
99 | if verbosity is 1, and to write runes always if verbosity is higher than | |
100 | that. | |
101 | """ | |
102 | me._active = VERBOSITY >= 2 or (VERBOSITY == 1 and me._level == 0) | |
103 | def push(me): | |
104 | """Push a new search level.""" | |
105 | me._level += 1 | |
106 | me._update() | |
107 | if me._level > 0: me.p('[') | |
108 | else: me.p(';; ') | |
109 | def pop(me): | |
110 | """Pop a search level.""" | |
111 | if me._level > 0: me.p(']') | |
112 | else: me.p('\n') | |
113 | me._level -= 1 | |
114 | me._update() | |
115 | def p(me, ch): | |
116 | """Print CH as a progress rune.""" | |
117 | if me._active: stderr.write(ch); stderr.flush() | |
118 | ||
119 | def combinations(r, v): | |
120 | """ | |
121 | Return an iterator which yields all combinations of R elements from V. | |
122 | ||
123 | V must be an indexable sequence. The each combination is returned as a | |
124 | list, containing elements from V in their original order. | |
125 | """ | |
126 | ||
127 | ## Set up the selection vector. C will contain the indices of the items of | |
128 | ## V we've selected for the current combination. At all times, C contains | |
129 | ## a strictly increasing sequence of integers in the interval [0, N). | |
130 | n = len(v) | |
131 | c = range(r) | |
132 | ||
133 | while True: | |
134 | ||
135 | ## Yield up the current combination. | |
136 | vv = [v[i] for i in c] | |
137 | yield vv | |
138 | ||
139 | ## Now advance to the next one. Find the last index in C which we can | |
140 | ## increment subject to the rules. As we iterate downwards, i will | |
141 | ## contain the index into C, and j will be the maximum acceptable value | |
142 | ## for the corresponding item. We'll step the last index until it | |
143 | ## reaches the limit, and then step the next one down, resetting the last | |
144 | ## index, and so on. | |
145 | i, j = r, n | |
146 | while True: | |
147 | ||
148 | ## If i is zero here, then we've advanced everything as far as it will | |
149 | ## go. We're done. | |
150 | if i == 0: return | |
151 | ||
152 | ## Move down to the next index. | |
153 | i -= 1; j -= 1 | |
154 | ||
155 | ## If this index isn't at its maximum value, then we've found the place | |
156 | ## to step. | |
157 | if c[i] != j: break | |
158 | ||
159 | ## Step this index on by one, and set the following indices to the | |
160 | ## immediately following values. | |
161 | j = c[i] + 1 | |
162 | while i < r: c[i] = j; i += 1; j += 1 | |
163 | ||
164 | class ArgFetcher (object): | |
165 | """ | |
166 | I return arguments from a list, reporting problems when they occur. | |
167 | """ | |
168 | def __init__(me, argv, errfn): | |
169 | """ | |
170 | Initialize, returning successive arguments from ARGV. | |
171 | ||
172 | Errors are reported to ERRFN. | |
173 | """ | |
174 | me._argv = argv | |
175 | me._argc = len(argv) | |
176 | me._errfn = errfn | |
177 | me._i = 0 | |
178 | def arg(me, default = None, must = True): | |
179 | """ | |
180 | Return the next argument. | |
181 | ||
182 | If MUST is true, then report an error (to the ERRFN) if there are no more | |
183 | arguments; otherwise, return the DEFAULT. | |
184 | """ | |
185 | if me._i >= me._argc: | |
186 | if must: me._errfn('missing argument') | |
187 | return default | |
188 | arg = me._argv[me._i]; me._i += 1 | |
189 | return arg | |
190 | def int(me, default = None, must = True, min = None, max = None): | |
191 | """ | |
192 | Return the next argument converted to an integer. | |
193 | ||
194 | If MUST is true, then report an error (to the ERRFN) if there are no more | |
195 | arguments; otherwise return the DEFAULT. Report an error if the next | |
196 | argument is not a valid integer, or if the integer is beyond the MIN and | |
197 | MAX bounds. | |
198 | """ | |
199 | arg = me.arg(default = None, must = must) | |
200 | if arg is None: return default | |
201 | try: arg = int(arg) | |
202 | except ValueError: me._errfn('bad integer') | |
203 | if (min is not None and arg < min) or (max is not None and arg > max): | |
204 | me._errfn('out of range') | |
205 | return arg | |
206 | ||
207 | ###-------------------------------------------------------------------------- | |
208 | ### Sieving for small primes. | |
209 | ||
210 | class Sieve (object): | |
211 | """ | |
212 | I represent a collection of small primes, up to some chosen limit. | |
213 | ||
214 | The limit is available as the `limit' attribute. Let L be this limit; | |
215 | then, if N < L^2 is some composite, then N has at least one prime factor | |
216 | less than L. | |
217 | """ | |
218 | ||
219 | ## Figure out the number of bits in a (nonnegative) primitive `int'. We'll | |
220 | ## use a list of these as our sieve. | |
221 | _NBIT = 15 | |
222 | while type(1 << (_NBIT + 1)) == int: _NBIT += 1 | |
223 | ||
224 | def __init__(me, limit): | |
225 | """ | |
226 | Initialize a sieve holding all primes below LIMIT. | |
227 | """ | |
228 | ||
229 | ## The sieve is maintained in the `_bits' attribute. This is a list of | |
230 | ## integers, used as a bitmask: let 2 < n < L be an odd integer; then bit | |
231 | ## (n - 3)/2 will be clear iff n is prime. Let W be the value of | |
232 | ## `_NBIT', above; then bit W i + j in the sieve is stored in bit j of | |
233 | ## `_bits[i]'. | |
234 | ||
235 | ## Store the limit for later inspection. | |
236 | me.limit = limit | |
237 | ||
238 | ## Calculate the size of sieve we'll need and initialize the bit list. | |
239 | n = (limit - 2)/2 | |
240 | sievesz = (n + me._NBIT - 1)/me._NBIT | |
241 | me._sievemax = sievesz*me._NBIT | |
537eab73 | 242 | me._bits = sievesz*[0] |
7de4d3fe MW |
243 | |
244 | ## This is standard Sieve of Eratosthenes. For each index i: if | |
245 | ## bit i is clear, then p = 2 i + 3 is prime, so set the bits | |
246 | ## corresponding to each multiple of p, i.e., bits (k p - 3)/2 = | |
247 | ## (2 k i + 3 - 3)/2 = k i for k > 1. | |
248 | for i in xrange(me._sievemax): | |
249 | if me._bitp(i): i += 1; continue | |
250 | p = 2*i + 3 | |
251 | if p >= limit: break | |
252 | for j in xrange(i + p, me._sievemax, p): me._setbit(j) | |
253 | i += 1 | |
254 | ||
255 | def _bitp(me, i): i, j = divmod(i, me._NBIT); return (me._bits[i] >> j)&1 | |
256 | def _setbit(me, i): i, j = divmod(i, me._NBIT); me._bits[i] |= 1 << j | |
257 | ||
258 | def smallprimes(me): | |
259 | """ | |
260 | Return an iterator over the known small primes. | |
261 | """ | |
262 | yield 2 | |
263 | n = 3 | |
264 | for b in me._bits: | |
265 | for j in xrange(me._NBIT): | |
266 | if not (b&1): yield n | |
267 | b >>= 1; n += 2 | |
268 | ||
269 | ## We generate the sieve on demand. | |
270 | SIEVE = None | |
271 | ||
272 | def initsieve(sievebits): | |
273 | """ | |
274 | Generate the sieve. | |
275 | ||
276 | Ensure that it can be used to check the primality of numbers up to (but not | |
277 | including) 2^SIEVEBITS. | |
278 | """ | |
279 | global SIEVE | |
280 | if SIEVE is not None: raise ValueError('sieve already defined') | |
281 | if sievebits < 6: sievebits = 6 | |
282 | SIEVE = Sieve(1 << (sievebits + 1)/2) | |
283 | ||
284 | ###-------------------------------------------------------------------------- | |
285 | ### Primality checking. | |
286 | ||
287 | def small_test(p): | |
288 | """ | |
289 | Check that P is a small prime. | |
290 | ||
291 | If not, raise an `ExpectedError'. The `SIEVE' variable must have been | |
292 | initialized. | |
293 | """ | |
294 | if p < 2: raise ExpectedError('%d too small' % p) | |
295 | if SIEVE.limit*SIEVE.limit < p: | |
296 | raise ExpectedError('%d too large for small prime' % p) | |
297 | for q in SIEVE.smallprimes(): | |
298 | if q*q > p: return | |
299 | if p%q == 0: raise ExpectedError('%d divides %d' % (q, p)) | |
300 | ||
301 | def pock_test(p, a, qq): | |
302 | """ | |
303 | Check that P is prime using Pocklington's criterion. | |
304 | ||
305 | If not, raise an `ExpectedError'. | |
306 | ||
307 | Let Q be the product of the elements of the sequence QQ. The test works as | |
308 | follows. Suppose p is the smallest prime factor of P. If A^{P-1} /== 1 | |
309 | (mod P) then P is certainly composite (Fermat's test); otherwise, we have | |
310 | established that the order of A in (Z/pZ)^* divides P - 1. Next, let t = | |
311 | A^{(P-1)/q} for some prime factor q of Q, and let g = gcd(t - 1, P). If g | |
312 | = P then the proof is inconclusive; if 1 < g < P then g is a nontrivial | |
313 | factor of P, so P is composite; otherwise, t has order q in (Z/pZ)^*, so | |
314 | (Z/pZ)^* contains a subgroup of size q, and therefore q divides p - 1. If | |
315 | QQ is a sequence of distinct primes, and the preceding criterion holds for | |
316 | all q in QQ, then Q divides p - 1. If Q^2 < P then the proof is | |
317 | inconclusive; otherwise, let p' be any prime dividing P/p. Then p' >= p > | |
318 | Q, so p p' > Q^2 > P; but p p' divides P, so this is a contradiction. | |
319 | Therefore P/p has no prime factors, and P is prime. | |
320 | """ | |
321 | ||
322 | ## We don't actually need the distinctness criterion. Suppose that q^e | |
323 | ## divides Q. Then gcd(t - 1, P) = 1 implies that A^{(P-1)/q^{e-1}} has | |
324 | ## order q^e in (Z/pZ)^*, which accounts for the multiplicity. | |
325 | ||
326 | Q = prod(qq) | |
327 | if p < 2: raise ExpectedError('%d too small' % p) | |
328 | if Q*Q <= p: | |
329 | raise ExpectedError('too few Pocklington factors for %d' % p) | |
330 | if pow(a, p - 1, p) != 1: | |
331 | raise ExpectedError('%d is Fermat witness for %d' % (a, p)) | |
332 | for q in qq: | |
333 | if Q%(q*q) == 0: | |
334 | raise ExpectedError('duplicate Pocklington factor %d for %d' % (q, p)) | |
335 | g = p.gcd(pow(a, (p - 1)/q, p) - 1) | |
336 | if g == p: | |
337 | raise ExpectedError('%d order not multiple of %d mod %d' % (a, q, p)) | |
338 | elif g != 1: | |
339 | raise ExpectedError('%d divides %d' % (g, p)) | |
340 | ||
341 | def ecpp_test(p, a, b, x, y, qq): | |
342 | """ | |
343 | Check that P is prime using Goldwasser and Kilian's ECPP method. | |
344 | ||
345 | If not, raise an `ExpectedError'. | |
346 | ||
347 | Let Q be the product of the elements of the sequence QQ. Suppose p is the | |
348 | smallest prime factor of P. Let g = gcd(4 A^3 + 27 B^2, P). If g = P then | |
349 | the test is inconclusive; otherwise, if g /= 1 then g is a nontrivial | |
350 | factor of P. Define E(GF(p)) = { (x, y) | y^2 = x^3 + A x + B } U { inf } | |
351 | to be the elliptic curve over p with short-Weierstraß coefficients A and B; | |
352 | we have just checked that this curve is not singular. If R = (X, Y) is not | |
353 | a point on this curve, then the test is inconclusive. If Q R is not the | |
354 | point at infinity, then the test fails; otherwise we deduce that P has | |
355 | Q-torsion in E. Let S = (Q/q) R for some prime factor q of Q. If S is the | |
356 | point at infinity then the test is inconclusive; otherwise, q divides the | |
357 | order of S in E. If QQ is a sequence of distinct primes, and the preceding | |
358 | criterion holds for all q in QQ, then Q divides the order of S. Therefore | |
359 | #E(p) >= Q. If Q <= (qrrt(P) + 1)^2 then the test is inconclusive. | |
360 | Otherwise, Hasse's theorem tells us that |p + 1 - #E(p)| <= 2 sqrt(p); | |
361 | hence we must have p + 1 + 2 sqrt(p) = (sqrt(p) + 1)^2 >= #E(p) >= Q > | |
362 | (qrrt(P) + 1)^2; so sqrt(p) + 1 > qrrt(P) + 1, i.e., p^2 > P. As for | |
363 | Pocklington above, if p' is any prime factor of P/p, then p p' >= p^2 > P, | |
364 | which is a contradiction, and we conclude that P is prime. | |
365 | """ | |
366 | ||
367 | ## This isn't going to work if gcd(P, 6) /= 1: we're going to use the | |
368 | ## large-characteristic addition formulae. | |
369 | g = p.gcd(6) | |
370 | if g != 1: raise ExpectedError('%d divides %d' % (g, p)) | |
371 | ||
372 | ## We want to check that Q > (qrrt(P) + 1)^2 iff sqrt(Q) > qrrt(P) + 1; but | |
373 | ## calculating square roots is not enjoyable (partly because we have to | |
374 | ## deal with the imprecision). Fortunately, some algebra will help: the | |
375 | ## condition holds iff qrrt(P) < sqrt(Q) - 1 iff P < Q^2 - 4 Q sqrt(Q) + | |
376 | ## 6 Q - 4 sqrt(Q) + 1 = Q (Q + 6) + 1 - 4 sqrt(Q) (Q + 1) iff Q (Q + 6) - | |
377 | ## P + 1 > 4 sqrt(Q) (Q + 1) iff (Q (Q + 6) - P + 1)^2 > 16 Q (Q + 1)^2 | |
378 | Q = prod(qq) | |
379 | t, u = Q*(Q + 6) - p + 1, 4*(Q + 1) | |
380 | if t*t <= Q*u*u: raise ExpectedError('too few subgroups for ECPP') | |
381 | ||
382 | ## Construct the curve. | |
383 | E = C.PrimeField(p).ec(a, b) # careful: may not be a prime! | |
384 | ||
385 | ## Find the base point. | |
386 | R = E(x, y) | |
387 | if not R.oncurvep(): | |
388 | raise ExpectedError('(%d, %d) is not on the curve' % (x, y)) | |
389 | ||
390 | ## Check that it has Q-torsion. | |
391 | if Q*R: raise ExpectedError('(%d, %d) not a %d-torsion point' % (x, y, Q)) | |
392 | ||
393 | ## Now check the individual factors. | |
394 | for q in qq: | |
395 | if Q%(q*q) == 0: | |
396 | raise ExpectedError('duplicate ECPP factor %d for %d' % (q, p)) | |
397 | S = (Q/q)*R | |
398 | if not S: | |
399 | raise ExpectedError('(%d, %d) order not a multiple of %d' % (x, y, q)) | |
400 | g = p.gcd(S._z) | |
401 | if g != 1: | |
402 | raise ExpectedError('%d divides %d' % (g, p)) | |
403 | ||
404 | ###-------------------------------------------------------------------------- | |
405 | ### Proof steps and proofs. | |
406 | ||
407 | class BaseStep (object): | |
408 | """ | |
409 | I'm a step in a primality proof. | |
410 | ||
411 | I assert that a particular number is prime, and can check this. | |
412 | ||
413 | This class provides basic protocol for proof steps, mostly to do with | |
414 | handling labels. | |
415 | ||
416 | The step's label is kept in its `label' attribute. It can be set by the | |
417 | constructor, and is `None' by default. Users can modify this attribute if | |
418 | they like. Labels beginning `$' are assumed to be internal and | |
419 | uninteresting; other labels cause `check' lines to be written to the output | |
420 | listing the actual number of interest. | |
421 | ||
422 | Protocol that proof steps should provide: | |
423 | ||
424 | label A string labelling the proof step and the associated prime | |
425 | number. | |
426 | ||
427 | p The prime number which this step proves to be prime. | |
428 | ||
429 | check() Check that the proof step is actually correct, assuming that | |
430 | any previous steps have already been verified. | |
431 | ||
432 | out(FILE) Write an appropriate encoding of the proof step to the output | |
433 | FILE. | |
434 | """ | |
435 | def __init__(me, label = None, *arg, **kw): | |
436 | """Initialize a proof step, setting a default label if necessary.""" | |
437 | super(BaseStep, me).__init__(*arg, **kw) | |
438 | me.label = label | |
439 | def out(me, file): | |
440 | """ | |
441 | Write the proof step to an output FILE. | |
442 | ||
443 | Subclasses must implement a method `_out' which actually does the work. | |
444 | Here, we write a `check' line to verify that the proof actually applies | |
445 | to the number we wanted, if the label says that this is an interesting | |
446 | step. | |
447 | """ | |
448 | me._out(file) | |
449 | if me.label is not None and not me.label.startswith('$'): | |
450 | file.write('check %s, %d, %d\n' % (me.label, me.p.nbits, me.p)) | |
451 | ||
452 | class SmallStep (BaseStep): | |
453 | """ | |
454 | I represent a claim that a number is a small prime. | |
455 | ||
456 | Such claims act as the base cases in a complicated primality proof. When | |
457 | verifying, the claim is checked by trial division using a collection of | |
458 | known small primes. | |
459 | """ | |
460 | def __init__(me, pp, p, *arg, **kw): | |
461 | """ | |
462 | Initialize a small-prime step. | |
463 | ||
464 | PP is the overall PrimeProof object of which this is a step; P is the | |
465 | small number whose primality is asserted. | |
466 | """ | |
467 | super(SmallStep, me).__init__(*arg, **kw) | |
468 | me.p = p | |
469 | def check(me): | |
470 | """Check that the number is indeed a small prime.""" | |
471 | return small_test(me.p) | |
472 | def _out(me, file): | |
473 | """Write a small-prime step to the FILE.""" | |
474 | file.write('small %s = %d\n' % (me.label, me.p)) | |
475 | def __repr__(me): return 'SmallStep(%d)' % (me.p) | |
476 | @classmethod | |
477 | def parse(cls, pp, line): | |
478 | """ | |
479 | Parse a small-prime step from a LINE in a proof file. | |
480 | ||
481 | SMALL-STEP ::= `small' LABEL `=' P | |
482 | ||
483 | PP is a PrimeProof object holding the results from the previous steps. | |
484 | """ | |
485 | if SIEVE is None: raise ExpectedError('missing `sievebits\' line') | |
486 | label, p = parse_label(line) | |
487 | return cls(pp, conv_int(p), label = label) | |
488 | ||
489 | class PockStep (BaseStep): | |
490 | """ | |
491 | I represent a Pocklington certificate for a number. | |
492 | ||
493 | The number is not explicitly represented in a proof file. See `pock_test' | |
494 | for the underlying mathematics. | |
495 | """ | |
496 | def __init__(me, pp, a, R, qqi, *arg, **kw): | |
497 | """ | |
498 | Inititialize a Pocklington step. | |
499 | ||
500 | PP is the overall PrimeProof object of which this is a step; A is the | |
501 | generator of a substantial subgroup of units; R is a cofactor; and QQI is | |
502 | a sequence of labels for previous proof steps. If Q is the product of | |
503 | the primes listed in QQI, then the number whose primality is asserted is | |
504 | 2 Q R + 1. | |
505 | """ | |
506 | super(PockStep, me).__init__(*arg, **kw) | |
507 | me._a = a | |
508 | me._R = R | |
509 | me._qqi = qqi | |
510 | me._qq = [pp.get_step(qi).p for qi in qqi] | |
511 | me.p = prod(me._qq, 2*R) + 1 | |
512 | def check(me): | |
513 | """Verify a proof step based on Pocklington's theorem.""" | |
514 | return pock_test(me.p, me._a, me._qq) | |
515 | def _out(me, file): | |
516 | """Write a Pocklington step to the FILE.""" | |
517 | file.write('pock %s = %d, %d, [%s]\n' % \ | |
518 | (me.label, me._a, | |
519 | me._R, ', '.join('%s' % qi for qi in me._qqi))) | |
520 | def __repr__(me): return 'PockStep(%d, %d, %s)' % (me._a, me._R, me._qqi) | |
521 | @classmethod | |
522 | def parse(cls, pp, line): | |
523 | """ | |
524 | Parse a Pocklington step from a LINE in a proof file. | |
525 | ||
526 | POCK-STEP ::= `pock' LABEL `=' A `,' R `,' `[' Q-LIST `]' | |
527 | Q-LIST ::= Q [`,' Q-LIST] | |
528 | ||
529 | PP is a PrimeProof object holding the results from the previous steps. | |
530 | """ | |
531 | label, rest = parse_label(line) | |
532 | a, R, qq = parse_list(rest, 3) | |
533 | qq = qq.strip() | |
534 | if not qq.startswith('[') or not qq.endswith(']'): | |
535 | raise ExpectedError('missing `[...]\' around Pocklington factors') | |
536 | return cls(pp, conv_int(a), conv_int(R), | |
537 | [q.strip() for q in qq[1:-1].split(',')], label = label) | |
538 | ||
539 | class ECPPStep (BaseStep): | |
540 | """ | |
541 | I represent a Goldwasser--Kilian ECPP certificate for a number. | |
542 | """ | |
543 | def __init__(me, pp, p, a, b, x, y, qqi, *arg, **kw): | |
544 | """ | |
545 | Inititialize an ECPP step. | |
546 | ||
547 | PP is the overall PrimeProof object of which this is a step; P is the | |
548 | number whose primality is asserted; A and B are the short Weierstraß | |
549 | curve coefficients; X and Y are the base point coordinates; and QQI is a | |
550 | sequence of labels for previous proof steps. | |
551 | """ | |
552 | super(ECPPStep, me).__init__(*arg, **kw) | |
553 | me._a, me._b = a, b | |
554 | me._x, me._y = x, y | |
555 | me._qqi = qqi | |
556 | me._qq = [pp.get_step(qi).p for qi in qqi] | |
557 | me.p = p | |
558 | def check(me): | |
559 | """Verify a proof step based on Goldwasser and Kilian's theorem.""" | |
560 | return ecpp_test(me.p, me._a, me._b, me._x, me._y, me._qq) | |
561 | def _out(me, file): | |
562 | """Write an ECPP step to the FILE.""" | |
563 | file.write('ecpp %s = %d, %d, %d, %d, %d, [%s]\n' % \ | |
564 | (me.label, me.p, me._a, me._b, me._x, me._y, | |
565 | ', '.join('%s' % qi for qi in me._qqi))) | |
566 | def __repr__(me): | |
567 | return 'ECPPstep(%d, %d, %d, %d, %d, %s)' % \ | |
568 | (me.p, me._a, me._b, me._x, me._y, me._qqi) | |
569 | @classmethod | |
570 | def parse(cls, pp, line): | |
571 | """ | |
572 | Parse an ECPP step from a LINE in a proof file. | |
573 | ||
574 | ECPP-STEP ::= `ecpp' LABEL `=' P `,' A `,' B `,' X `,' Y `,' | |
575 | `[' Q-LIST `]' | |
576 | Q-LIST ::= Q [`,' Q-LIST] | |
577 | ||
578 | PP is a PrimeProof object holding the results from the previous steps. | |
579 | """ | |
580 | label, rest = parse_label(line) | |
581 | p, a, b, x, y, qq = parse_list(rest, 6) | |
582 | qq = qq.strip() | |
583 | if not qq.startswith('[') or not qq.endswith(']'): | |
584 | raise ExpectedError('missing `[...]\' around ECPP factors') | |
585 | return cls(pp, conv_int(p), conv_int(a), conv_int(b), | |
586 | conv_int(x), conv_int(y), | |
587 | [q.strip() for q in qq[1:-1].split(',')], label = label) | |
588 | ||
589 | def check(pp, line): | |
590 | """ | |
591 | Handle a `check' line in a proof file. | |
592 | ||
593 | CHECK ::= `check' LABEL, B, N | |
594 | ||
595 | Verify that the proof step with the given LABEL asserts the primality of | |
596 | the integer N, and that 2^{B-1} <= N < 2^B. | |
597 | """ | |
598 | label, nb, p = parse_list(line, 3) | |
599 | label, nb, p = label.strip(), conv_int(nb), conv_int(p) | |
600 | pi = pp.get_step(label).p | |
601 | if pi != p: | |
602 | raise ExpectedError('check failed: %s = %d /= %d' % (label, pi, p)) | |
603 | if p.nbits != nb: | |
604 | raise ExpectedError('check failed: nbits(%s) = %d /= %d' % \ | |
605 | (label, p.nbits, nb)) | |
bbb113f6 | 606 | if VERBOSITY: print(';; %s = %d [%d]' % (label, p, nb)) |
7de4d3fe MW |
607 | |
608 | def setsievebits(pp, line): | |
609 | """ | |
610 | Handle a `sievebits' line in a proof file. | |
611 | ||
612 | SIEVEBITS ::= `sievebits' N | |
613 | ||
614 | Ensure that the verifier is willing to accept small primes up to 2^N. | |
615 | """ | |
616 | initsieve(int(line)) | |
617 | ||
618 | class PrimeProof (object): | |
619 | """ | |
620 | I represent a proof of primality for one or more numbers. | |
621 | ||
622 | I can encode my proof as a line-oriented text file, in a simple format, and | |
623 | read such a proof back to check it. | |
624 | """ | |
625 | ||
626 | ## A table to dispatch on keywords read from a file. | |
627 | STEPMAP = { 'small': SmallStep.parse, | |
628 | 'pock': PockStep.parse, | |
629 | 'ecpp': ECPPStep.parse, | |
630 | 'sievebits': setsievebits, | |
631 | 'check': check } | |
632 | ||
633 | def __init__(me): | |
634 | """ | |
635 | Initialize a proof object. | |
636 | """ | |
637 | me._steps = {} # Maps labels to steps. | |
638 | me._stepseq = [] # Sequence of labels, in order. | |
639 | me._pmap = {} # Maps primes to steps. | |
640 | me._i = 0 | |
641 | ||
642 | def addstep(me, step): | |
643 | """ | |
644 | Add a new STEP to the proof. | |
645 | ||
646 | The STEP may have a label already. If not, a new internal label is | |
647 | chosen. The proof step is checked before being added to the proof. The | |
648 | label is returned. | |
649 | """ | |
650 | ||
651 | ## If there's already a step for this prime, and the new step doesn't | |
652 | ## have a label, then return the old one instead. | |
653 | if step.label is None: | |
654 | try: return me._pmap[step.p] | |
655 | except KeyError: pass | |
656 | ||
657 | ## Make sure the step is actually correct. | |
658 | step.check() | |
659 | ||
660 | ## Generate a label if the step doesn't have one already. | |
661 | if step.label is None: step.label = '$t%d' % me._i; me._i += 1 | |
662 | ||
663 | ## If the label is already taken then we have a problem. | |
664 | if step.label in me._steps: | |
31b5a1fe | 665 | raise ExpectedError('duplicate label `%s\'' % step.label) |
7de4d3fe MW |
666 | |
667 | ## Store the proof step. | |
668 | me._pmap[step.p] = step.label | |
669 | me._steps[step.label] = step | |
670 | me._stepseq.append(step.label) | |
671 | return step.label | |
672 | ||
673 | def get_step(me, label): | |
674 | """ | |
675 | Check that LABEL labels a known step, and return that step. | |
676 | """ | |
677 | try: return me._steps[label] | |
678 | except KeyError: raise ExpectedError('unknown label `%s\'' % label) | |
679 | ||
680 | def write(me, file): | |
681 | """ | |
682 | Write the proof to the given FILE. | |
683 | """ | |
684 | ||
685 | ## Prefix the main steps with a `sievebits' line. | |
686 | file.write('sievebits %d\n' % (2*(SIEVE.limit.bit_length() - 1))) | |
687 | ||
688 | ## Write the steps out one by one. | |
689 | for label in me._stepseq: me._steps[label].out(file) | |
690 | ||
691 | def read(me, file): | |
692 | """ | |
693 | Read a proof from a given FILE. | |
694 | ||
695 | FILE ::= {STEP | CHECK | SIEVEBITS} [FILE] | |
696 | STEP ::= SMALL-STEP | POCK-STEP | |
697 | ||
698 | Comments (beginning `;') and blank lines are ignored. Other lines begin | |
699 | with a keyword. | |
700 | """ | |
701 | lastp = None | |
702 | for lno, line in enumerate(file, 1): | |
703 | line = line.strip() | |
704 | if line.startswith(';'): continue | |
705 | ww = line.split(None, 1) | |
706 | if not ww: continue | |
707 | w = ww[0] | |
708 | if len(ww) > 1: tail = ww[1] | |
709 | else: tail = '' | |
710 | try: | |
711 | try: op = me.STEPMAP[w] | |
712 | except KeyError: | |
713 | raise ExpectedError('unrecognized keyword `%s\'' % w) | |
714 | step = op(me, tail) | |
715 | if step is not None: | |
716 | me.addstep(step) | |
717 | lastp = step.p | |
f6d012db MW |
718 | except ExpectedError: |
719 | raise ExpectedError('%s:%d: %s' % | |
720 | (file.name, lno, _excval().message)) | |
7de4d3fe MW |
721 | return lastp |
722 | ||
723 | ###-------------------------------------------------------------------------- | |
724 | ### Finding provable primes. | |
725 | ||
726 | class BasePrime (object): | |
727 | """ | |
728 | I represent a prime number which has been found and can be proven. | |
729 | ||
730 | This object can eventually be turned into a sequence of proof steps and | |
731 | added to a PrimeProof. This isn't done immediately, because some | |
732 | prime-search strategies want to build a pool of provable primes and will | |
733 | then select some subset of them to actually construct the number of final | |
734 | interest. This way, we avoid cluttering the output proof with proofs of | |
735 | uninteresting numbers. | |
736 | ||
737 | Protocol required. | |
738 | ||
739 | p The prime number in question. | |
740 | ||
741 | label(LABEL) Associate LABEL with this prime, and the corresponding proof | |
742 | step. A label can be set in the constructor, or later using | |
743 | this method. | |
744 | ||
745 | register(PP) Register the prime with a PrimeProof, adding any necessary | |
746 | proof steps. Returns the label of the proof step for this | |
747 | number. | |
748 | ||
749 | _mkstep(PP, **KW) | |
750 | Return a proof step for this prime. | |
751 | """ | |
752 | def __init__(me, label = None, *args, **kw): | |
753 | """Initialize a provable prime number object.""" | |
754 | super(BasePrime, me).__init__(*args, **kw) | |
755 | me._index = me._pp = None | |
756 | me._label = label | |
757 | def label(me, label): | |
758 | """Set this number's LABEL.""" | |
759 | me._label = label | |
760 | def register(me, pp): | |
761 | """ | |
762 | Register the prime's proof steps with PrimeProof PP. | |
763 | ||
764 | Return the final step's label. | |
765 | """ | |
766 | if me._pp is not None: | |
767 | assert me._pp == pp | |
768 | else: | |
769 | me._pp = pp | |
770 | me._index = pp.addstep(me._mkstep(pp, label = me._label)) | |
771 | ##try: me._index = pp.addstep(me._mkstep(pp, label = me._label)) | |
772 | ##except: raise RuntimeError('generated proof failed sanity check') | |
773 | return me._index | |
774 | ||
775 | class SmallPrime (BasePrime): | |
776 | """I represent a prime small enough to be checked in isolation.""" | |
777 | def __init__(me, p, *args, **kw): | |
778 | super(SmallPrime, me).__init__(*args, **kw) | |
779 | me.p = p | |
780 | def _mkstep(me, pp, **kw): | |
781 | return SmallStep(pp, me.p, **kw) | |
782 | ||
783 | class PockPrime (BasePrime): | |
784 | """I represent a prime proven using Pocklington's theorem.""" | |
785 | def __init__(me, p, a, qq, *args, **kw): | |
786 | super(PockPrime, me).__init__(*args, **kw) | |
787 | me.p = p | |
788 | me._a = a | |
789 | me._qq = qq | |
790 | def _mkstep(me, pp, **kw): | |
791 | return PockStep(pp, me._a, (me.p - 1)/prod((q.p for q in me._qq), 2), | |
792 | [q.register(pp) for q in me._qq], **kw) | |
793 | ||
794 | def gen_small(nbits, label = None, p = None): | |
795 | """ | |
796 | Return a new small prime. | |
797 | ||
798 | The prime will be exactly NBITS bits long. The proof step will have the | |
799 | given LABEL attached. Report progress to the ProgressReporter P. | |
800 | """ | |
801 | while True: | |
802 | ||
803 | ## Pick a random NBITS-bit number. | |
804 | n = C.rand.mp(nbits, 1) | |
805 | assert n.nbits == nbits | |
806 | ||
807 | ## If it's probably prime, then check it against the small primes we | |
808 | ## know. If it passes then we're done. Otherwise, try again. | |
809 | if n.primep(): | |
810 | for q in SIEVE.smallprimes(): | |
811 | if q*q > n: return SmallPrime(n, label = label) | |
812 | if n%q == 0: break | |
813 | ||
814 | def gen_pock(nbits, nsubbits = 0, label = None, p = ProgressReporter()): | |
815 | """ | |
816 | Return a new prime provable using Pocklington's theorem. | |
817 | ||
818 | The prime N will be exactly NBITS long, of the form N = 2 Q R + 1. If | |
819 | NSUBBITS is nonzero, then each prime factor of Q will be NSUBBITS bits | |
820 | long; otherwise a suitable default will be chosen. The proof step will | |
821 | have the given LABEL attached. Report progress to the ProgressReporter P. | |
822 | ||
823 | The prime numbers this function returns are a long way from being uniformly | |
824 | distributed. | |
825 | """ | |
826 | ||
827 | ## Pick a suitable value for NSUBBITS if we don't have one. | |
828 | if not nsubbits: | |
829 | ||
830 | ## This is remarkably tricky. Picking about 1/3 sqrt(NBITS) factors | |
831 | ## seems about right for large numbers, but there's serious trouble | |
832 | ## lurking for small sizes. | |
833 | nsubbits = int(3*M.sqrt(nbits)) | |
834 | if nbits < nsubbits + 3: nsubbits = nbits//2 + 1 | |
835 | if nbits == 2*nsubbits: nsubbits += 1 | |
836 | ||
837 | ## Figure out how many subgroups we'll need. | |
838 | npiece = ((nbits + 1)//2 + nsubbits - 1)//nsubbits | |
839 | p.push() | |
840 | ||
841 | ## Keep searching... | |
842 | while True: | |
843 | ||
844 | ## Come up with a collection of known prime factors. | |
845 | p.p('!'); qq = [gen(nsubbits, p = p) for i in xrange(npiece)] | |
846 | Q = prod(q.p for q in qq) | |
847 | ||
848 | ## Come up with bounds on the cofactor. If we're to have N = 2 Q R + 1, | |
849 | ## and 2^{B-1} <= N < 2^B, then we must have 2^{B-2}/Q <= R < 2^{B-1}/Q. | |
850 | Rbase = (C.MP(0).setbit(nbits - 2) + Q - 1)//Q | |
851 | Rwd = C.MP(0).setbit(nbits - 2)//Q | |
852 | ||
853 | ## Probe the available space of cofactors. If the space is kind of | |
854 | ## narrow, then we want to give up quickly if we're not finding anything | |
855 | ## suitable. | |
856 | step = 0 | |
857 | while step < Rwd: | |
858 | step += 1 | |
859 | ||
860 | ## Pick a random cofactor and examine the number we ended up with. | |
861 | ## Make sure it really does have the length we expect. | |
862 | R = C.rand.range(Rwd) + Rbase | |
863 | n = 2*Q*R + 1 | |
864 | assert n.nbits == nbits | |
865 | ||
866 | ## As a complication, if NPIECE is 1, it's just about possible that Q^2 | |
867 | ## <= n, in which case this isn't going to work. | |
868 | if Q*Q < n: continue | |
869 | ||
870 | ## If n has small factors, then pick another cofactor. | |
871 | if C.PrimeFilter.smallfactor(n) == C.PGEN_FAIL: continue | |
872 | ||
873 | ## Work through the small primes to find a suitable generator. The | |
874 | ## value 2 is almost always acceptable, so don't try too hard here. | |
875 | for a in I.islice(SIEVE.smallprimes(), 16): | |
876 | ||
877 | ## First, try the Fermat test. If that fails, then n is definitely | |
878 | ## composite. | |
879 | if pow(a, n - 1, n) != 1: p.p('.'); break | |
880 | p.p('*') | |
881 | ||
882 | ## Work through the subgroup orders, checking that suitable powers of | |
883 | ## a generate the necessary subgroups. | |
884 | for q in qq: | |
885 | if n.gcd(pow(a, (n - 1)/q.p, n) - 1) != 1: | |
886 | p.p('@'); ok = False; break | |
887 | else: | |
888 | ok = True | |
889 | ||
890 | ## we're all good. | |
891 | if ok: p.pop(); return PockPrime(n, a, qq, label = label) | |
892 | ||
893 | def gen(nbits, label = None, p = ProgressReporter()): | |
894 | """ | |
895 | Generate a prime number with NBITS bits. | |
896 | ||
897 | Give it the LABEL, and report progress to P. | |
898 | """ | |
899 | if SIEVE.limit >> (nbits + 1)/2: g = gen_small | |
900 | else: g = gen_pock | |
901 | return g(nbits, label = label, p = p) | |
902 | ||
903 | def gen_limlee(nbits, nsubbits, | |
904 | label = None, qlfmt = None, p = ProgressReporter()): | |
905 | """ | |
906 | Generate a Lim--Lee prime with NBITS bits. | |
907 | ||
908 | Let p be the prime. Then we'll have p = 2 q_0 q_1 ... q_k, with all q_i at | |
909 | least NSUBBITS bits long, and all but q_k exactly that long. | |
910 | ||
911 | The prime will be given the LABEL; progress is reported to P. The factors | |
912 | q_i will be labelled by filling in the `printf'-style format string QLFMT | |
913 | with the argument i. | |
914 | """ | |
915 | ||
916 | ## Figure out how many factors (p - 1)/2 will have. | |
917 | npiece = nbits//nsubbits | |
918 | if npiece < 2: raise ExpectedError('too few pieces') | |
919 | ||
920 | ## Decide how big to make the pool of factors. | |
921 | poolsz = max(3*npiece + 5, 25) # Heuristic from GnuPG | |
922 | ||
923 | ## Prepare for the main loop. | |
924 | disp = nstep = 0 | |
925 | qbig = None | |
926 | p.push() | |
927 | ||
928 | ## Try to make a prime. | |
929 | while True: | |
930 | p.p('!') | |
931 | ||
932 | ## Construct a pool of NSUBBITS-size primes. There's a problem with very | |
933 | ## small sizes: we might not be able to build a pool of distinct primes. | |
934 | pool = []; qmap = {} | |
935 | for i in xrange(poolsz): | |
936 | for j in xrange(64): | |
937 | q = gen(nsubbits, p = p) | |
938 | if q.p not in qmap: break | |
939 | else: | |
940 | raise ExpectedError('insufficient diversity') | |
941 | qmap[q.p] = q | |
942 | pool.append(q) | |
943 | ||
944 | ## Work through combinations of factors from the pool. | |
945 | for qq in combinations(npiece - 1, pool): | |
946 | ||
947 | ## Construct the product of the selected factors. | |
948 | qsmall = prod(q.p for q in qq) | |
949 | ||
950 | ## Maybe we'll need to replace the large factor. Try not to do this | |
951 | ## too often. DISP measures the large factor's performance at | |
952 | ## producing candidates with the right length. If it looks bad then | |
953 | ## we'll have to replace it. | |
954 | if 3*disp*disp > nstep*nstep: | |
955 | qbig = None | |
956 | if disp < 0: p.p('<') | |
957 | else: p.p('>') | |
958 | ||
959 | ## If we don't have a large factor, then make one. | |
960 | if qbig is None: | |
961 | qbig = gen(nbits - qsmall.nbits, p = p) | |
962 | disp = 0; nstep = 0 | |
963 | ||
964 | ## We have a candidate. Calculate it and make sure it has the right | |
965 | ## length. | |
966 | n = 2*qsmall*qbig.p + 1 | |
967 | nstep += 1 | |
968 | if n.nbits < nbits: disp -= 1 | |
969 | elif n.nbits > nbits: disp += 1 | |
970 | elif C.PrimeFilter.smallfactor(n) == C.PGEN_FAIL: pass | |
971 | else: | |
972 | ||
973 | ## The candidate has passed the small-primes test. Now check it | |
974 | ## against Pocklington. | |
975 | for a in I.islice(SIEVE.smallprimes(), 16): | |
976 | ||
977 | ## Fermat test. | |
978 | if pow(a, n - 1, n) != 1: p.p('.'); break | |
979 | p.p('*') | |
980 | ||
981 | ## Find a generator of a sufficiently large subgroup. | |
982 | if n.gcd(pow(a, (n - 1)/qbig.p, n) - 1) != 1: p.p('@'); continue | |
983 | ok = True | |
984 | for q in qq: | |
985 | if n.gcd(pow(a, (n - 1)/q.p, n) - 1) != 1: | |
986 | p.p('@'); ok = False; break | |
987 | ||
988 | ## We're done. | |
989 | if ok: | |
990 | ||
991 | ## Label the factors. | |
992 | qq.append(qbig) | |
993 | if qlfmt: | |
994 | for i, q in enumerate(qq): q.label(qlfmt % i) | |
995 | ||
996 | ## Return the number we found. | |
997 | p.pop(); return PockPrime(n, a, qq, label = label) | |
998 | ||
999 | ###-------------------------------------------------------------------------- | |
1000 | ### Main program. | |
1001 | ||
1002 | def __main__(): | |
1003 | global VERBOSITY | |
1004 | ||
1005 | ## Prepare an option parser. | |
1006 | op = OP.OptionParser( | |
1007 | usage = '''\ | |
ba2dd9e2 | 1008 | pock [-qv] [-s SIEVEBITS] CMD ARGS... |
7de4d3fe MW |
1009 | gen NBITS |
1010 | ll NBITS NSUBBITS | |
596f3d96 | 1011 | check [FILE]''', |
7de4d3fe MW |
1012 | description = 'Generate or verify certified prime numbers.') |
1013 | op.add_option('-v', '--verbose', dest = 'verbosity', | |
1014 | action = 'count', default = 1, | |
fa965026 | 1015 | help = 'print mysterious runes while looking for prime numbers') |
7de4d3fe MW |
1016 | op.add_option('-q', '--quiet', dest = 'quietude', |
1017 | action = 'count', default = 0, | |
fa965026 | 1018 | help = 'be quiet while looking for prime numbers') |
7de4d3fe MW |
1019 | op.add_option('-s', '--sievebits', dest = 'sievebits', |
1020 | type = 'int', default = 32, | |
fa965026 | 1021 | help = 'size (in bits) of largest small prime') |
7de4d3fe MW |
1022 | opts, argv = op.parse_args() |
1023 | VERBOSITY = opts.verbosity - opts.quietude | |
1024 | p = ProgressReporter() | |
1025 | a = ArgFetcher(argv, op.error) | |
1026 | ||
1027 | ## Process arguments and do what the user asked. | |
1028 | w = a.arg() | |
1029 | ||
1030 | if w == 'gen': | |
1031 | ## Generate a prime with no special structure. | |
1032 | initsieve(opts.sievebits) | |
1033 | nbits = a.int(min = 4) | |
1034 | pp = PrimeProof() | |
1035 | p = gen(nbits, 'p', p = p) | |
1036 | p.register(pp) | |
1037 | pp.write(stdout) | |
1038 | ||
1039 | elif w == 'll': | |
1040 | ## Generate a Lim--Lee prime. | |
1041 | initsieve(opts.sievebits) | |
1042 | nbits = a.int(min = 4) | |
1043 | nsubbits = a.int(min = 4, max = nbits) | |
1044 | pp = PrimeProof() | |
1045 | p = gen_limlee(nbits, nsubbits, 'p', 'q_%d', p = p) | |
1046 | p.register(pp) | |
1047 | pp.write(stdout) | |
1048 | ||
1049 | elif w == 'check': | |
1050 | ## Check an existing certificate. | |
1051 | fn = a.arg(default = '-', must = False) | |
1052 | if fn == '-': f = stdin | |
1053 | else: f = open(fn, 'r') | |
1054 | pp = PrimeProof() | |
1055 | p = pp.read(f) | |
1056 | ||
1057 | else: | |
1058 | raise ExpectedError("unknown command `%s'" % w) | |
1059 | ||
1060 | if __name__ == '__main__': | |
1061 | prog = OS.path.basename(argv[0]) | |
1062 | try: __main__() | |
f6d012db MW |
1063 | except ExpectedError: exit('%s: %s' % (prog, _excval().message)) |
1064 | except IOError: exit('%s: %s' % (prog, _excval())) | |
7de4d3fe MW |
1065 | |
1066 | ###----- That's all, folks -------------------------------------------------- |