| 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 = sievesz*[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 ExpectedError('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] [-s SIEVEBITS] 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 -------------------------------------------------- |