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 | ||
41 | class ExpectedError (Exception): | |
42 | """ | |
43 | I represent an expected error, which should be reported in a friendly way. | |
44 | """ | |
45 | pass | |
46 | ||
47 | def prod(ff, one = 1): | |
48 | """ | |
49 | Return ONE times the product of the elements of FF. | |
50 | ||
51 | This is not done very efficiently. | |
52 | """ | |
53 | return reduce(lambda prod, f: prod*f, ff, one) | |
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 | |
242 | me._bits = n*[0] | |
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)) | |
606 | if VERBOSITY: print ';; %s = %d [%d]' % (label, p, nb) | |
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: | |
665 | raise ValueError('duplicate label `%s\'' % step.label) | |
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 | |
718 | except ExpectedError, e: | |
719 | raise ExpectedError('%s:%d: %s' % (file.name, lno, e.message)) | |
720 | return lastp | |
721 | ||
722 | ###-------------------------------------------------------------------------- | |
723 | ### Finding provable primes. | |
724 | ||
725 | class BasePrime (object): | |
726 | """ | |
727 | I represent a prime number which has been found and can be proven. | |
728 | ||
729 | This object can eventually be turned into a sequence of proof steps and | |
730 | added to a PrimeProof. This isn't done immediately, because some | |
731 | prime-search strategies want to build a pool of provable primes and will | |
732 | then select some subset of them to actually construct the number of final | |
733 | interest. This way, we avoid cluttering the output proof with proofs of | |
734 | uninteresting numbers. | |
735 | ||
736 | Protocol required. | |
737 | ||
738 | p The prime number in question. | |
739 | ||
740 | label(LABEL) Associate LABEL with this prime, and the corresponding proof | |
741 | step. A label can be set in the constructor, or later using | |
742 | this method. | |
743 | ||
744 | register(PP) Register the prime with a PrimeProof, adding any necessary | |
745 | proof steps. Returns the label of the proof step for this | |
746 | number. | |
747 | ||
748 | _mkstep(PP, **KW) | |
749 | Return a proof step for this prime. | |
750 | """ | |
751 | def __init__(me, label = None, *args, **kw): | |
752 | """Initialize a provable prime number object.""" | |
753 | super(BasePrime, me).__init__(*args, **kw) | |
754 | me._index = me._pp = None | |
755 | me._label = label | |
756 | def label(me, label): | |
757 | """Set this number's LABEL.""" | |
758 | me._label = label | |
759 | def register(me, pp): | |
760 | """ | |
761 | Register the prime's proof steps with PrimeProof PP. | |
762 | ||
763 | Return the final step's label. | |
764 | """ | |
765 | if me._pp is not None: | |
766 | assert me._pp == pp | |
767 | else: | |
768 | me._pp = pp | |
769 | me._index = pp.addstep(me._mkstep(pp, label = me._label)) | |
770 | ##try: me._index = pp.addstep(me._mkstep(pp, label = me._label)) | |
771 | ##except: raise RuntimeError('generated proof failed sanity check') | |
772 | return me._index | |
773 | ||
774 | class SmallPrime (BasePrime): | |
775 | """I represent a prime small enough to be checked in isolation.""" | |
776 | def __init__(me, p, *args, **kw): | |
777 | super(SmallPrime, me).__init__(*args, **kw) | |
778 | me.p = p | |
779 | def _mkstep(me, pp, **kw): | |
780 | return SmallStep(pp, me.p, **kw) | |
781 | ||
782 | class PockPrime (BasePrime): | |
783 | """I represent a prime proven using Pocklington's theorem.""" | |
784 | def __init__(me, p, a, qq, *args, **kw): | |
785 | super(PockPrime, me).__init__(*args, **kw) | |
786 | me.p = p | |
787 | me._a = a | |
788 | me._qq = qq | |
789 | def _mkstep(me, pp, **kw): | |
790 | return PockStep(pp, me._a, (me.p - 1)/prod((q.p for q in me._qq), 2), | |
791 | [q.register(pp) for q in me._qq], **kw) | |
792 | ||
793 | def gen_small(nbits, label = None, p = None): | |
794 | """ | |
795 | Return a new small prime. | |
796 | ||
797 | The prime will be exactly NBITS bits long. The proof step will have the | |
798 | given LABEL attached. Report progress to the ProgressReporter P. | |
799 | """ | |
800 | while True: | |
801 | ||
802 | ## Pick a random NBITS-bit number. | |
803 | n = C.rand.mp(nbits, 1) | |
804 | assert n.nbits == nbits | |
805 | ||
806 | ## If it's probably prime, then check it against the small primes we | |
807 | ## know. If it passes then we're done. Otherwise, try again. | |
808 | if n.primep(): | |
809 | for q in SIEVE.smallprimes(): | |
810 | if q*q > n: return SmallPrime(n, label = label) | |
811 | if n%q == 0: break | |
812 | ||
813 | def gen_pock(nbits, nsubbits = 0, label = None, p = ProgressReporter()): | |
814 | """ | |
815 | Return a new prime provable using Pocklington's theorem. | |
816 | ||
817 | The prime N will be exactly NBITS long, of the form N = 2 Q R + 1. If | |
818 | NSUBBITS is nonzero, then each prime factor of Q will be NSUBBITS bits | |
819 | long; otherwise a suitable default will be chosen. The proof step will | |
820 | have the given LABEL attached. Report progress to the ProgressReporter P. | |
821 | ||
822 | The prime numbers this function returns are a long way from being uniformly | |
823 | distributed. | |
824 | """ | |
825 | ||
826 | ## Pick a suitable value for NSUBBITS if we don't have one. | |
827 | if not nsubbits: | |
828 | ||
829 | ## This is remarkably tricky. Picking about 1/3 sqrt(NBITS) factors | |
830 | ## seems about right for large numbers, but there's serious trouble | |
831 | ## lurking for small sizes. | |
832 | nsubbits = int(3*M.sqrt(nbits)) | |
833 | if nbits < nsubbits + 3: nsubbits = nbits//2 + 1 | |
834 | if nbits == 2*nsubbits: nsubbits += 1 | |
835 | ||
836 | ## Figure out how many subgroups we'll need. | |
837 | npiece = ((nbits + 1)//2 + nsubbits - 1)//nsubbits | |
838 | p.push() | |
839 | ||
840 | ## Keep searching... | |
841 | while True: | |
842 | ||
843 | ## Come up with a collection of known prime factors. | |
844 | p.p('!'); qq = [gen(nsubbits, p = p) for i in xrange(npiece)] | |
845 | Q = prod(q.p for q in qq) | |
846 | ||
847 | ## Come up with bounds on the cofactor. If we're to have N = 2 Q R + 1, | |
848 | ## and 2^{B-1} <= N < 2^B, then we must have 2^{B-2}/Q <= R < 2^{B-1}/Q. | |
849 | Rbase = (C.MP(0).setbit(nbits - 2) + Q - 1)//Q | |
850 | Rwd = C.MP(0).setbit(nbits - 2)//Q | |
851 | ||
852 | ## Probe the available space of cofactors. If the space is kind of | |
853 | ## narrow, then we want to give up quickly if we're not finding anything | |
854 | ## suitable. | |
855 | step = 0 | |
856 | while step < Rwd: | |
857 | step += 1 | |
858 | ||
859 | ## Pick a random cofactor and examine the number we ended up with. | |
860 | ## Make sure it really does have the length we expect. | |
861 | R = C.rand.range(Rwd) + Rbase | |
862 | n = 2*Q*R + 1 | |
863 | assert n.nbits == nbits | |
864 | ||
865 | ## As a complication, if NPIECE is 1, it's just about possible that Q^2 | |
866 | ## <= n, in which case this isn't going to work. | |
867 | if Q*Q < n: continue | |
868 | ||
869 | ## If n has small factors, then pick another cofactor. | |
870 | if C.PrimeFilter.smallfactor(n) == C.PGEN_FAIL: continue | |
871 | ||
872 | ## Work through the small primes to find a suitable generator. The | |
873 | ## value 2 is almost always acceptable, so don't try too hard here. | |
874 | for a in I.islice(SIEVE.smallprimes(), 16): | |
875 | ||
876 | ## First, try the Fermat test. If that fails, then n is definitely | |
877 | ## composite. | |
878 | if pow(a, n - 1, n) != 1: p.p('.'); break | |
879 | p.p('*') | |
880 | ||
881 | ## Work through the subgroup orders, checking that suitable powers of | |
882 | ## a generate the necessary subgroups. | |
883 | for q in qq: | |
884 | if n.gcd(pow(a, (n - 1)/q.p, n) - 1) != 1: | |
885 | p.p('@'); ok = False; break | |
886 | else: | |
887 | ok = True | |
888 | ||
889 | ## we're all good. | |
890 | if ok: p.pop(); return PockPrime(n, a, qq, label = label) | |
891 | ||
892 | def gen(nbits, label = None, p = ProgressReporter()): | |
893 | """ | |
894 | Generate a prime number with NBITS bits. | |
895 | ||
896 | Give it the LABEL, and report progress to P. | |
897 | """ | |
898 | if SIEVE.limit >> (nbits + 1)/2: g = gen_small | |
899 | else: g = gen_pock | |
900 | return g(nbits, label = label, p = p) | |
901 | ||
902 | def gen_limlee(nbits, nsubbits, | |
903 | label = None, qlfmt = None, p = ProgressReporter()): | |
904 | """ | |
905 | Generate a Lim--Lee prime with NBITS bits. | |
906 | ||
907 | Let p be the prime. Then we'll have p = 2 q_0 q_1 ... q_k, with all q_i at | |
908 | least NSUBBITS bits long, and all but q_k exactly that long. | |
909 | ||
910 | The prime will be given the LABEL; progress is reported to P. The factors | |
911 | q_i will be labelled by filling in the `printf'-style format string QLFMT | |
912 | with the argument i. | |
913 | """ | |
914 | ||
915 | ## Figure out how many factors (p - 1)/2 will have. | |
916 | npiece = nbits//nsubbits | |
917 | if npiece < 2: raise ExpectedError('too few pieces') | |
918 | ||
919 | ## Decide how big to make the pool of factors. | |
920 | poolsz = max(3*npiece + 5, 25) # Heuristic from GnuPG | |
921 | ||
922 | ## Prepare for the main loop. | |
923 | disp = nstep = 0 | |
924 | qbig = None | |
925 | p.push() | |
926 | ||
927 | ## Try to make a prime. | |
928 | while True: | |
929 | p.p('!') | |
930 | ||
931 | ## Construct a pool of NSUBBITS-size primes. There's a problem with very | |
932 | ## small sizes: we might not be able to build a pool of distinct primes. | |
933 | pool = []; qmap = {} | |
934 | for i in xrange(poolsz): | |
935 | for j in xrange(64): | |
936 | q = gen(nsubbits, p = p) | |
937 | if q.p not in qmap: break | |
938 | else: | |
939 | raise ExpectedError('insufficient diversity') | |
940 | qmap[q.p] = q | |
941 | pool.append(q) | |
942 | ||
943 | ## Work through combinations of factors from the pool. | |
944 | for qq in combinations(npiece - 1, pool): | |
945 | ||
946 | ## Construct the product of the selected factors. | |
947 | qsmall = prod(q.p for q in qq) | |
948 | ||
949 | ## Maybe we'll need to replace the large factor. Try not to do this | |
950 | ## too often. DISP measures the large factor's performance at | |
951 | ## producing candidates with the right length. If it looks bad then | |
952 | ## we'll have to replace it. | |
953 | if 3*disp*disp > nstep*nstep: | |
954 | qbig = None | |
955 | if disp < 0: p.p('<') | |
956 | else: p.p('>') | |
957 | ||
958 | ## If we don't have a large factor, then make one. | |
959 | if qbig is None: | |
960 | qbig = gen(nbits - qsmall.nbits, p = p) | |
961 | disp = 0; nstep = 0 | |
962 | ||
963 | ## We have a candidate. Calculate it and make sure it has the right | |
964 | ## length. | |
965 | n = 2*qsmall*qbig.p + 1 | |
966 | nstep += 1 | |
967 | if n.nbits < nbits: disp -= 1 | |
968 | elif n.nbits > nbits: disp += 1 | |
969 | elif C.PrimeFilter.smallfactor(n) == C.PGEN_FAIL: pass | |
970 | else: | |
971 | ||
972 | ## The candidate has passed the small-primes test. Now check it | |
973 | ## against Pocklington. | |
974 | for a in I.islice(SIEVE.smallprimes(), 16): | |
975 | ||
976 | ## Fermat test. | |
977 | if pow(a, n - 1, n) != 1: p.p('.'); break | |
978 | p.p('*') | |
979 | ||
980 | ## Find a generator of a sufficiently large subgroup. | |
981 | if n.gcd(pow(a, (n - 1)/qbig.p, n) - 1) != 1: p.p('@'); continue | |
982 | ok = True | |
983 | for q in qq: | |
984 | if n.gcd(pow(a, (n - 1)/q.p, n) - 1) != 1: | |
985 | p.p('@'); ok = False; break | |
986 | ||
987 | ## We're done. | |
988 | if ok: | |
989 | ||
990 | ## Label the factors. | |
991 | qq.append(qbig) | |
992 | if qlfmt: | |
993 | for i, q in enumerate(qq): q.label(qlfmt % i) | |
994 | ||
995 | ## Return the number we found. | |
996 | p.pop(); return PockPrime(n, a, qq, label = label) | |
997 | ||
998 | ###-------------------------------------------------------------------------- | |
999 | ### Main program. | |
1000 | ||
1001 | def __main__(): | |
1002 | global VERBOSITY | |
1003 | ||
1004 | ## Prepare an option parser. | |
1005 | op = OP.OptionParser( | |
1006 | usage = '''\ | |
1007 | pock [-qv] CMD ARGS... | |
1008 | gen NBITS | |
1009 | ll NBITS NSUBBITS | |
1010 | [check] [FILE]''', | |
1011 | description = 'Generate or verify certified prime numbers.') | |
1012 | op.add_option('-v', '--verbose', dest = 'verbosity', | |
1013 | action = 'count', default = 1, | |
1014 | help = 'Print mysterious runes while looking for prime numbers.') | |
1015 | op.add_option('-q', '--quiet', dest = 'quietude', | |
1016 | action = 'count', default = 0, | |
1017 | help = 'be quiet while looking for prime numbers.') | |
1018 | op.add_option('-s', '--sievebits', dest = 'sievebits', | |
1019 | type = 'int', default = 32, | |
1020 | help = 'Size (in bits) of largest small prime.') | |
1021 | opts, argv = op.parse_args() | |
1022 | VERBOSITY = opts.verbosity - opts.quietude | |
1023 | p = ProgressReporter() | |
1024 | a = ArgFetcher(argv, op.error) | |
1025 | ||
1026 | ## Process arguments and do what the user asked. | |
1027 | w = a.arg() | |
1028 | ||
1029 | if w == 'gen': | |
1030 | ## Generate a prime with no special structure. | |
1031 | initsieve(opts.sievebits) | |
1032 | nbits = a.int(min = 4) | |
1033 | pp = PrimeProof() | |
1034 | p = gen(nbits, 'p', p = p) | |
1035 | p.register(pp) | |
1036 | pp.write(stdout) | |
1037 | ||
1038 | elif w == 'll': | |
1039 | ## Generate a Lim--Lee prime. | |
1040 | initsieve(opts.sievebits) | |
1041 | nbits = a.int(min = 4) | |
1042 | nsubbits = a.int(min = 4, max = nbits) | |
1043 | pp = PrimeProof() | |
1044 | p = gen_limlee(nbits, nsubbits, 'p', 'q_%d', p = p) | |
1045 | p.register(pp) | |
1046 | pp.write(stdout) | |
1047 | ||
1048 | elif w == 'check': | |
1049 | ## Check an existing certificate. | |
1050 | fn = a.arg(default = '-', must = False) | |
1051 | if fn == '-': f = stdin | |
1052 | else: f = open(fn, 'r') | |
1053 | pp = PrimeProof() | |
1054 | p = pp.read(f) | |
1055 | ||
1056 | else: | |
1057 | raise ExpectedError("unknown command `%s'" % w) | |
1058 | ||
1059 | if __name__ == '__main__': | |
1060 | prog = OS.path.basename(argv[0]) | |
1061 | try: __main__() | |
1062 | except ExpectedError, e: exit('%s: %s' % (prog, e.message)) | |
1063 | except IOError, e: exit('%s: %s' % (prog, e)) | |
1064 | ||
1065 | ###----- That's all, folks -------------------------------------------------- |