2 ### -*- mode: python, coding: utf-8 -*-
4 ### Tool for generating and verifying primality certificates
6 ### (c) 2017 Straylight/Edgeware
9 ###----- Licensing notice ---------------------------------------------------
11 ### This file is part of the Python interface to Catacomb.
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.
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.
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.
27 ###--------------------------------------------------------------------------
30 from sys
import argv
, stdin
, stdout
, stderr
38 ###--------------------------------------------------------------------------
42 """Return the most recent exception object."""
43 return SYS
.exc_info()[1]
45 class ExpectedError (Exception):
47 I represent an expected error, which should be reported in a friendly way.
51 def prod(ff
, one
= 1):
52 """Return ONE times the product of the elements of FF."""
53 return C
.MPMul().factor(one
).factor(ff
).done()
55 def parse_label(line
):
57 Split LINE at an `=' character and return the left and right halves.
59 The returned pieces have leading and trailing whitespace trimmed.
62 if eq
< 0: raise ExpectedError('expected `LABEL = ...\'')
63 return line
[:eq
].strip(), line
[eq
+ 1:].strip()
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
)
72 """Convert S to a integer."""
73 try: return C
.MP(s
, 0)
74 except TypeError: raise ExpectedError('invalid integer `%s\'' % s
)
78 class ProgressReporter (object):
80 I keep users amused while the program looks for large prime numbers.
82 My main strategy is the printing of incomprehensible runes. I can be
83 muffled by lowering by verbosity level.
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.
90 """Initialize the ProgressReporter."""
95 Update our idea of whether we're active.
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
102 me
._active
= VERBOSITY
>= 2 or (VERBOSITY
== 1 and me
._level
== 0)
104 """Push a new search level."""
107 if me
._level
> 0: me
.p('[')
110 """Pop a search level."""
111 if me
._level
> 0: me
.p(']')
116 """Print CH as a progress rune."""
117 if me
._active
: stderr
.write(ch
); stderr
.flush()
119 def combinations(r
, v
):
121 Return an iterator which yields all combinations of R elements from V.
123 V must be an indexable sequence. The each combination is returned as a
124 list, containing elements from V in their original order.
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).
135 ## Yield up the current combination.
136 vv
= [v
[i
] for i
in c
]
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
148 ## If i is zero here, then we've advanced everything as far as it will
152 ## Move down to the next index.
155 ## If this index isn't at its maximum value, then we've found the place
159 ## Step this index on by one, and set the following indices to the
160 ## immediately following values.
162 while i
< r
: c
[i
] = j
; i
+= 1; j
+= 1
164 class ArgFetcher (object):
166 I return arguments from a list, reporting problems when they occur.
168 def __init__(me
, argv
, errfn
):
170 Initialize, returning successive arguments from ARGV.
172 Errors are reported to ERRFN.
178 def arg(me
, default
= None, must
= True):
180 Return the next argument.
182 If MUST is true, then report an error (to the ERRFN) if there are no more
183 arguments; otherwise, return the DEFAULT.
185 if me
._i
>= me
._argc
:
186 if must
: me
._errfn('missing argument')
188 arg
= me
._argv
[me
._i
]; me
._i
+= 1
190 def int(me
, default
= None, must
= True, min = None, max = None):
192 Return the next argument converted to an integer.
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
199 arg
= me
.arg(default
= None, must
= must
)
200 if arg
is None: return default
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')
207 ###--------------------------------------------------------------------------
208 ### Sieving for small primes.
210 class Sieve (object):
212 I represent a collection of small primes, up to some chosen limit.
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
219 ## Figure out the number of bits in a (nonnegative) primitive `int'. We'll
220 ## use a list of these as our sieve.
221 try: _MAX
= SYS
.maxint
222 except AttributeError: _MAX
= SYS
.maxsize
224 while 1 << (_NBIT
+ 1) < _MAX
: _NBIT
+= 1
226 def __init__(me
, limit
):
228 Initialize a sieve holding all primes below LIMIT.
231 ## The sieve is maintained in the `_bits' attribute. This is a list of
232 ## integers, used as a bitmask: let 2 < n < L be an odd integer; then bit
233 ## (n - 3)/2 will be clear iff n is prime. Let W be the value of
234 ## `_NBIT', above; then bit W i + j in the sieve is stored in bit j of
237 ## Store the limit for later inspection.
240 ## Calculate the size of sieve we'll need and initialize the bit list.
242 sievesz
= (n
+ me
._NBIT
- 1)//me
._NBIT
243 me
._sievemax
= sievesz
*me
._NBIT
244 me
._bits
= sievesz
*[0]
246 ## This is standard Sieve of Eratosthenes. For each index i: if
247 ## bit i is clear, then p = 2 i + 3 is prime, so set the bits
248 ## corresponding to each multiple of p, i.e., bits (k p - 3)/2 =
249 ## (2 k i + 3 - 3)/2 = k i for k > 1.
250 for i
in xrange(me
._sievemax
):
251 if me
._bitp(i
): i
+= 1; continue
254 for j
in xrange(i
+ p
, me
._sievemax
, p
): me
._setbit(j
)
257 def _bitp(me
, i
): i
, j
= divmod(i
, me
._NBIT
); return (me
._bits
[i
] >> j
)&1
258 def _setbit(me
, i
): i
, j
= divmod(i
, me
._NBIT
); me
._bits
[i
] |
= 1 << j
262 Return an iterator over the known small primes.
267 for j
in xrange(me
._NBIT
):
268 if not (b
&1): yield n
271 ## We generate the sieve on demand.
274 def initsieve(sievebits
):
278 Ensure that it can be used to check the primality of numbers up to (but not
279 including) 2^SIEVEBITS.
282 if SIEVE
is not None: raise ValueError('sieve already defined')
283 if sievebits
< 6: sievebits
= 6
284 SIEVE
= Sieve(1 << (sievebits
+ 1)//2)
286 ###--------------------------------------------------------------------------
287 ### Primality checking.
291 Check that P is a small prime.
293 If not, raise an `ExpectedError'. The `SIEVE' variable must have been
296 if p
< 2: raise ExpectedError('%d too small' % p
)
297 if SIEVE
.limit
*SIEVE
.limit
< p
:
298 raise ExpectedError('%d too large for small prime' % p
)
299 for q
in SIEVE
.smallprimes():
301 if p
%q
== 0: raise ExpectedError('%d divides %d' %
(q
, p
))
303 def pock_test(p
, a
, qq
):
305 Check that P is prime using Pocklington's criterion.
307 If not, raise an `ExpectedError'.
309 Let Q be the product of the elements of the sequence QQ. The test works as
310 follows. Suppose p is the smallest prime factor of P. If A^{P-1} /== 1
311 (mod P) then P is certainly composite (Fermat's test); otherwise, we have
312 established that the order of A in (Z/pZ)^* divides P - 1. Next, let t =
313 A^{(P-1)/q} for some prime factor q of Q, and let g = gcd(t - 1, P). If g
314 = P then the proof is inconclusive; if 1 < g < P then g is a nontrivial
315 factor of P, so P is composite; otherwise, t has order q in (Z/pZ)^*, so
316 (Z/pZ)^* contains a subgroup of size q, and therefore q divides p - 1. If
317 QQ is a sequence of distinct primes, and the preceding criterion holds for
318 all q in QQ, then Q divides p - 1. If Q^2 < P then the proof is
319 inconclusive; otherwise, let p' be any prime dividing P/p. Then p' >= p >
320 Q, so p p' > Q^2 > P; but p p' divides P, so this is a contradiction.
321 Therefore P/p has no prime factors, and P is prime.
324 ## We don't actually need the distinctness criterion. Suppose that q^e
325 ## divides Q. Then gcd(t - 1, P) = 1 implies that A^{(P-1)/q^{e-1}} has
326 ## order q^e in (Z/pZ)^*, which accounts for the multiplicity.
329 if p
< 2: raise ExpectedError('%d too small' % p
)
331 raise ExpectedError('too few Pocklington factors for %d' % p
)
332 if pow(a
, p
- 1, p
) != 1:
333 raise ExpectedError('%d is Fermat witness for %d' %
(a
, p
))
336 raise ExpectedError('duplicate Pocklington factor %d for %d' %
(q
, p
))
337 g
= p
.gcd(pow(a
, (p
- 1)/q
, p
) - 1)
339 raise ExpectedError('%d order not multiple of %d mod %d' %
(a
, q
, p
))
341 raise ExpectedError('%d divides %d' %
(g
, p
))
343 def ecpp_test(p
, a
, b
, x
, y
, qq
):
345 Check that P is prime using Goldwasser and Kilian's ECPP method.
347 If not, raise an `ExpectedError'.
349 Let Q be the product of the elements of the sequence QQ. Suppose p is the
350 smallest prime factor of P. Let g = gcd(4 A^3 + 27 B^2, P). If g = P then
351 the test is inconclusive; otherwise, if g /= 1 then g is a nontrivial
352 factor of P. Define E(GF(p)) = { (x, y) | y^2 = x^3 + A x + B } U { inf }
353 to be the elliptic curve over p with short-Weierstraß coefficients A and B;
354 we have just checked that this curve is not singular. If R = (X, Y) is not
355 a point on this curve, then the test is inconclusive. If Q R is not the
356 point at infinity, then the test fails; otherwise we deduce that P has
357 Q-torsion in E. Let S = (Q/q) R for some prime factor q of Q. If S is the
358 point at infinity then the test is inconclusive; otherwise, q divides the
359 order of S in E. If QQ is a sequence of distinct primes, and the preceding
360 criterion holds for all q in QQ, then Q divides the order of S. Therefore
361 #E(p) >= Q. If Q <= (qrrt(P) + 1)^2 then the test is inconclusive.
362 Otherwise, Hasse's theorem tells us that |p + 1 - #E(p)| <= 2 sqrt(p);
363 hence we must have p + 1 + 2 sqrt(p) = (sqrt(p) + 1)^2 >= #E(p) >= Q >
364 (qrrt(P) + 1)^2; so sqrt(p) + 1 > qrrt(P) + 1, i.e., p^2 > P. As for
365 Pocklington above, if p' is any prime factor of P/p, then p p' >= p^2 > P,
366 which is a contradiction, and we conclude that P is prime.
369 ## This isn't going to work if gcd(P, 6) /= 1: we're going to use the
370 ## large-characteristic addition formulae.
372 if g
!= 1: raise ExpectedError('%d divides %d' %
(g
, p
))
374 ## We want to check that Q > (qrrt(P) + 1)^2 iff sqrt(Q) > qrrt(P) + 1; but
375 ## calculating square roots is not enjoyable (partly because we have to
376 ## deal with the imprecision). Fortunately, some algebra will help: the
377 ## condition holds iff qrrt(P) < sqrt(Q) - 1 iff P < Q^2 - 4 Q sqrt(Q) +
378 ## 6 Q - 4 sqrt(Q) + 1 = Q (Q + 6) + 1 - 4 sqrt(Q) (Q + 1) iff Q (Q + 6) -
379 ## P + 1 > 4 sqrt(Q) (Q + 1) iff (Q (Q + 6) - P + 1)^2 > 16 Q (Q + 1)^2
381 t
, u
= Q
*(Q
+ 6) - p
+ 1, 4*(Q
+ 1)
382 if t
*t
<= Q
*u
*u
: raise ExpectedError('too few subgroups for ECPP')
384 ## Construct the curve.
385 E
= C
.PrimeField(p
).ec(a
, b
) # careful: may not be a prime!
387 ## Find the base point.
390 raise ExpectedError('(%d, %d) is not on the curve' %
(x
, y
))
392 ## Check that it has Q-torsion.
393 if Q
*R
: raise ExpectedError('(%d, %d) not a %d-torsion point' %
(x
, y
, Q
))
395 ## Now check the individual factors.
398 raise ExpectedError('duplicate ECPP factor %d for %d' %
(q
, p
))
401 raise ExpectedError('(%d, %d) order not a multiple of %d' %
(x
, y
, q
))
404 raise ExpectedError('%d divides %d' %
(g
, p
))
406 ###--------------------------------------------------------------------------
407 ### Proof steps and proofs.
409 class BaseStep (object):
411 I'm a step in a primality proof.
413 I assert that a particular number is prime, and can check this.
415 This class provides basic protocol for proof steps, mostly to do with
418 The step's label is kept in its `label' attribute. It can be set by the
419 constructor, and is `None' by default. Users can modify this attribute if
420 they like. Labels beginning `$' are assumed to be internal and
421 uninteresting; other labels cause `check' lines to be written to the output
422 listing the actual number of interest.
424 Protocol that proof steps should provide:
426 label A string labelling the proof step and the associated prime
429 p The prime number which this step proves to be prime.
431 check() Check that the proof step is actually correct, assuming that
432 any previous steps have already been verified.
434 out(FILE) Write an appropriate encoding of the proof step to the output
437 def __init__(me
, label
= None, *arg
, **kw
):
438 """Initialize a proof step, setting a default label if necessary."""
439 super(BaseStep
, me
).__init__(*arg
, **kw
)
443 Write the proof step to an output FILE.
445 Subclasses must implement a method `_out' which actually does the work.
446 Here, we write a `check' line to verify that the proof actually applies
447 to the number we wanted, if the label says that this is an interesting
451 if me
.label
is not None and not me
.label
.startswith('$'):
452 file.write('check %s, %d, %d\n' %
(me
.label
, me
.p
.nbits
, me
.p
))
454 class SmallStep (BaseStep
):
456 I represent a claim that a number is a small prime.
458 Such claims act as the base cases in a complicated primality proof. When
459 verifying, the claim is checked by trial division using a collection of
462 def __init__(me
, pp
, p
, *arg
, **kw
):
464 Initialize a small-prime step.
466 PP is the overall PrimeProof object of which this is a step; P is the
467 small number whose primality is asserted.
469 super(SmallStep
, me
).__init__(*arg
, **kw
)
472 """Check that the number is indeed a small prime."""
473 return small_test(me
.p
)
475 """Write a small-prime step to the FILE."""
476 file.write('small %s = %d\n' %
(me
.label
, me
.p
))
477 def __repr__(me
): return 'SmallStep(%d)' %
(me
.p
)
479 def parse(cls
, pp
, line
):
481 Parse a small-prime step from a LINE in a proof file.
483 SMALL-STEP ::= `small' LABEL `=' P
485 PP is a PrimeProof object holding the results from the previous steps.
487 if SIEVE
is None: raise ExpectedError('missing `sievebits\' line')
488 label
, p
= parse_label(line
)
489 return cls(pp
, conv_int(p
), label
= label
)
491 class PockStep (BaseStep
):
493 I represent a Pocklington certificate for a number.
495 The number is not explicitly represented in a proof file. See `pock_test'
496 for the underlying mathematics.
498 def __init__(me
, pp
, a
, R
, qqi
, *arg
, **kw
):
500 Inititialize a Pocklington step.
502 PP is the overall PrimeProof object of which this is a step; A is the
503 generator of a substantial subgroup of units; R is a cofactor; and QQI is
504 a sequence of labels for previous proof steps. If Q is the product of
505 the primes listed in QQI, then the number whose primality is asserted is
508 super(PockStep
, me
).__init__(*arg
, **kw
)
512 me
._qq
= [pp
.get_step(qi
).p
for qi
in qqi
]
513 me
.p
= prod(me
._qq
, 2*R
) + 1
515 """Verify a proof step based on Pocklington's theorem."""
516 return pock_test(me
.p
, me
._a
, me
._qq
)
518 """Write a Pocklington step to the FILE."""
519 file.write('pock %s = %d, %d, [%s]\n' % \
521 me
._R
, ', '.join('%s' % qi
for qi
in me
._qqi
)))
522 def __repr__(me
): return 'PockStep(%d, %d, %s)' %
(me
._a
, me
._R
, me
._qqi
)
524 def parse(cls
, pp
, line
):
526 Parse a Pocklington step from a LINE in a proof file.
528 POCK-STEP ::= `pock' LABEL `=' A `,' R `,' `[' Q-LIST `]'
529 Q-LIST ::= Q [`,' Q-LIST]
531 PP is a PrimeProof object holding the results from the previous steps.
533 label
, rest
= parse_label(line
)
534 a
, R
, qq
= parse_list(rest
, 3)
536 if not qq
.startswith('[') or not qq
.endswith(']'):
537 raise ExpectedError('missing `[...]\' around Pocklington factors')
538 return cls(pp
, conv_int(a
), conv_int(R
),
539 [q
.strip() for q
in qq
[1:-1].split(',')], label
= label
)
541 class ECPPStep (BaseStep
):
543 I represent a Goldwasser--Kilian ECPP certificate for a number.
545 def __init__(me
, pp
, p
, a
, b
, x
, y
, qqi
, *arg
, **kw
):
547 Inititialize an ECPP step.
549 PP is the overall PrimeProof object of which this is a step; P is the
550 number whose primality is asserted; A and B are the short Weierstraß
551 curve coefficients; X and Y are the base point coordinates; and QQI is a
552 sequence of labels for previous proof steps.
554 super(ECPPStep
, me
).__init__(*arg
, **kw
)
558 me
._qq
= [pp
.get_step(qi
).p
for qi
in qqi
]
561 """Verify a proof step based on Goldwasser and Kilian's theorem."""
562 return ecpp_test(me
.p
, me
._a
, me
._b
, me
._x
, me
._y
, me
._qq
)
564 """Write an ECPP step to the FILE."""
565 file.write('ecpp %s = %d, %d, %d, %d, %d, [%s]\n' % \
566 (me
.label
, me
.p
, me
._a
, me
._b
, me
._x
, me
._y
,
567 ', '.join('%s' % qi
for qi
in me
._qqi
)))
569 return 'ECPPstep(%d, %d, %d, %d, %d, %s)' % \
570 (me
.p
, me
._a
, me
._b
, me
._x
, me
._y
, me
._qqi
)
572 def parse(cls
, pp
, line
):
574 Parse an ECPP step from a LINE in a proof file.
576 ECPP-STEP ::= `ecpp' LABEL `=' P `,' A `,' B `,' X `,' Y `,'
578 Q-LIST ::= Q [`,' Q-LIST]
580 PP is a PrimeProof object holding the results from the previous steps.
582 label
, rest
= parse_label(line
)
583 p
, a
, b
, x
, y
, qq
= parse_list(rest
, 6)
585 if not qq
.startswith('[') or not qq
.endswith(']'):
586 raise ExpectedError('missing `[...]\' around ECPP factors')
587 return cls(pp
, conv_int(p
), conv_int(a
), conv_int(b
),
588 conv_int(x
), conv_int(y
),
589 [q
.strip() for q
in qq
[1:-1].split(',')], label
= label
)
593 Handle a `check' line in a proof file.
595 CHECK ::= `check' LABEL, B, N
597 Verify that the proof step with the given LABEL asserts the primality of
598 the integer N, and that 2^{B-1} <= N < 2^B.
600 label
, nb
, p
= parse_list(line
, 3)
601 label
, nb
, p
= label
.strip(), conv_int(nb
), conv_int(p
)
602 pi
= pp
.get_step(label
).p
604 raise ExpectedError('check failed: %s = %d /= %d' %
(label
, pi
, p
))
606 raise ExpectedError('check failed: nbits(%s) = %d /= %d' % \
607 (label
, p
.nbits
, nb
))
608 if VERBOSITY
: print(';; %s = %d [%d]' %
(label
, p
, nb
))
610 def setsievebits(pp
, line
):
612 Handle a `sievebits' line in a proof file.
614 SIEVEBITS ::= `sievebits' N
616 Ensure that the verifier is willing to accept small primes up to 2^N.
620 class PrimeProof (object):
622 I represent a proof of primality for one or more numbers.
624 I can encode my proof as a line-oriented text file, in a simple format, and
625 read such a proof back to check it.
628 ## A table to dispatch on keywords read from a file.
629 STEPMAP
= { 'small': SmallStep
.parse
,
630 'pock': PockStep
.parse
,
631 'ecpp': ECPPStep
.parse
,
632 'sievebits': setsievebits
,
637 Initialize a proof object.
639 me
._steps
= {} # Maps labels to steps.
640 me
._stepseq
= [] # Sequence of labels, in order.
641 me
._pmap
= {} # Maps primes to steps.
644 def addstep(me
, step
):
646 Add a new STEP to the proof.
648 The STEP may have a label already. If not, a new internal label is
649 chosen. The proof step is checked before being added to the proof. The
653 ## If there's already a step for this prime, and the new step doesn't
654 ## have a label, then return the old one instead.
655 if step
.label
is None:
656 try: return me
._pmap
[step
.p
]
657 except KeyError: pass
659 ## Make sure the step is actually correct.
662 ## Generate a label if the step doesn't have one already.
663 if step
.label
is None: step
.label
= '$t%d' % me
._i
; me
._i
+= 1
665 ## If the label is already taken then we have a problem.
666 if step
.label
in me
._steps
:
667 raise ExpectedError('duplicate label `%s\'' % step
.label
)
669 ## Store the proof step.
670 me
._pmap
[step
.p
] = step
.label
671 me
._steps
[step
.label
] = step
672 me
._stepseq
.append(step
.label
)
675 def get_step(me
, label
):
677 Check that LABEL labels a known step, and return that step.
679 try: return me
._steps
[label
]
680 except KeyError: raise ExpectedError('unknown label `%s\'' % label
)
684 Write the proof to the given FILE.
687 ## Prefix the main steps with a `sievebits' line.
688 file.write('sievebits %d\n' %
(2*(SIEVE
.limit
.bit_length() - 1)))
690 ## Write the steps out one by one.
691 for label
in me
._stepseq
: me
._steps
[label
].out(file)
695 Read a proof from a given FILE.
697 FILE ::= {STEP | CHECK | SIEVEBITS} [FILE]
698 STEP ::= SMALL-STEP | POCK-STEP
700 Comments (beginning `;') and blank lines are ignored. Other lines begin
704 for lno
, line
in enumerate(file, 1):
706 if line
.startswith(';'): continue
707 ww
= line
.split(None, 1)
710 if len(ww
) > 1: tail
= ww
[1]
713 try: op
= me
.STEPMAP
[w
]
715 raise ExpectedError('unrecognized keyword `%s\'' % w
)
720 except ExpectedError
:
721 raise ExpectedError('%s:%d: %s' %
722 (file.name
, lno
, _excval().message
))
725 ###--------------------------------------------------------------------------
726 ### Finding provable primes.
728 class BasePrime (object):
730 I represent a prime number which has been found and can be proven.
732 This object can eventually be turned into a sequence of proof steps and
733 added to a PrimeProof. This isn't done immediately, because some
734 prime-search strategies want to build a pool of provable primes and will
735 then select some subset of them to actually construct the number of final
736 interest. This way, we avoid cluttering the output proof with proofs of
737 uninteresting numbers.
741 p The prime number in question.
743 label(LABEL) Associate LABEL with this prime, and the corresponding proof
744 step. A label can be set in the constructor, or later using
747 register(PP) Register the prime with a PrimeProof, adding any necessary
748 proof steps. Returns the label of the proof step for this
752 Return a proof step for this prime.
754 def __init__(me
, label
= None, *args
, **kw
):
755 """Initialize a provable prime number object."""
756 super(BasePrime
, me
).__init__(*args
, **kw
)
757 me
._index
= me
._pp
= None
759 def label(me
, label
):
760 """Set this number's LABEL."""
762 def register(me
, pp
):
764 Register the prime's proof steps with PrimeProof PP.
766 Return the final step's label.
768 if me
._pp
is not None:
772 me
._index
= pp
.addstep(me
._mkstep(pp
, label
= me
._label
))
773 ##try: me._index = pp.addstep(me._mkstep(pp, label = me._label))
774 ##except: raise RuntimeError('generated proof failed sanity check')
777 class SmallPrime (BasePrime
):
778 """I represent a prime small enough to be checked in isolation."""
779 def __init__(me
, p
, *args
, **kw
):
780 super(SmallPrime
, me
).__init__(*args
, **kw
)
782 def _mkstep(me
, pp
, **kw
):
783 return SmallStep(pp
, me
.p
, **kw
)
785 class PockPrime (BasePrime
):
786 """I represent a prime proven using Pocklington's theorem."""
787 def __init__(me
, p
, a
, qq
, *args
, **kw
):
788 super(PockPrime
, me
).__init__(*args
, **kw
)
792 def _mkstep(me
, pp
, **kw
):
793 return PockStep(pp
, me
._a
, (me
.p
- 1)/prod((q
.p
for q
in me
._qq
), 2),
794 [q
.register(pp
) for q
in me
._qq
], **kw
)
796 def gen_small(nbits
, label
= None, p
= None):
798 Return a new small prime.
800 The prime will be exactly NBITS bits long. The proof step will have the
801 given LABEL attached. Report progress to the ProgressReporter P.
805 ## Pick a random NBITS-bit number.
806 n
= C
.rand
.mp(nbits
, 1)
807 assert n
.nbits
== nbits
809 ## If it's probably prime, then check it against the small primes we
810 ## know. If it passes then we're done. Otherwise, try again.
812 for q
in SIEVE
.smallprimes():
813 if q
*q
> n
: return SmallPrime(n
, label
= label
)
816 def gen_pock(nbits
, nsubbits
= 0, label
= None, p
= ProgressReporter()):
818 Return a new prime provable using Pocklington's theorem.
820 The prime N will be exactly NBITS long, of the form N = 2 Q R + 1. If
821 NSUBBITS is nonzero, then each prime factor of Q will be NSUBBITS bits
822 long; otherwise a suitable default will be chosen. The proof step will
823 have the given LABEL attached. Report progress to the ProgressReporter P.
825 The prime numbers this function returns are a long way from being uniformly
829 ## Pick a suitable value for NSUBBITS if we don't have one.
832 ## This is remarkably tricky. Picking about 1/3 sqrt(NBITS) factors
833 ## seems about right for large numbers, but there's serious trouble
834 ## lurking for small sizes.
835 nsubbits
= int(3*M
.sqrt(nbits
))
836 if nbits
< nsubbits
+ 3: nsubbits
= nbits
//2 + 1
837 if nbits
== 2*nsubbits
: nsubbits
+= 1
839 ## Figure out how many subgroups we'll need.
840 npiece
= ((nbits
+ 1)//2 + nsubbits
- 1)//nsubbits
846 ## Come up with a collection of known prime factors.
847 p
.p('!'); qq
= [gen(nsubbits
, p
= p
) for i
in xrange(npiece
)]
848 Q
= prod(q
.p
for q
in qq
)
850 ## Come up with bounds on the cofactor. If we're to have N = 2 Q R + 1,
851 ## and 2^{B-1} <= N < 2^B, then we must have 2^{B-2}/Q <= R < 2^{B-1}/Q.
852 Rbase
= (C
.MP(0).setbit(nbits
- 2) + Q
- 1)//Q
853 Rwd
= C
.MP(0).setbit(nbits
- 2)//Q
855 ## Probe the available space of cofactors. If the space is kind of
856 ## narrow, then we want to give up quickly if we're not finding anything
862 ## Pick a random cofactor and examine the number we ended up with.
863 ## Make sure it really does have the length we expect.
864 R
= C
.rand
.range(Rwd
) + Rbase
866 assert n
.nbits
== nbits
868 ## As a complication, if NPIECE is 1, it's just about possible that Q^2
869 ## <= n, in which case this isn't going to work.
872 ## If n has small factors, then pick another cofactor.
873 if C
.PrimeFilter
.smallfactor(n
) == C
.PGEN_FAIL
: continue
875 ## Work through the small primes to find a suitable generator. The
876 ## value 2 is almost always acceptable, so don't try too hard here.
877 for a
in I
.islice(SIEVE
.smallprimes(), 16):
879 ## First, try the Fermat test. If that fails, then n is definitely
881 if pow(a
, n
- 1, n
) != 1: p
.p('.'); break
884 ## Work through the subgroup orders, checking that suitable powers of
885 ## a generate the necessary subgroups.
887 if n
.gcd(pow(a
, (n
- 1)/q
.p
, n
) - 1) != 1:
888 p
.p('@'); ok
= False; break
893 if ok
: p
.pop(); return PockPrime(n
, a
, qq
, label
= label
)
895 def gen(nbits
, label
= None, p
= ProgressReporter()):
897 Generate a prime number with NBITS bits.
899 Give it the LABEL, and report progress to P.
901 if SIEVE
.limit
>> (nbits
+ 1)//2: g
= gen_small
903 return g(nbits
, label
= label
, p
= p
)
905 def gen_limlee(nbits
, nsubbits
,
906 label
= None, qlfmt
= None, p
= ProgressReporter()):
908 Generate a Lim--Lee prime with NBITS bits.
910 Let p be the prime. Then we'll have p = 2 q_0 q_1 ... q_k, with all q_i at
911 least NSUBBITS bits long, and all but q_k exactly that long.
913 The prime will be given the LABEL; progress is reported to P. The factors
914 q_i will be labelled by filling in the `printf'-style format string QLFMT
918 ## Figure out how many factors (p - 1)/2 will have.
919 npiece
= nbits
//nsubbits
920 if npiece
< 2: raise ExpectedError('too few pieces')
922 ## Decide how big to make the pool of factors.
923 poolsz
= max(3*npiece
+ 5, 25) # Heuristic from GnuPG
925 ## Prepare for the main loop.
930 ## Try to make a prime.
934 ## Construct a pool of NSUBBITS-size primes. There's a problem with very
935 ## small sizes: we might not be able to build a pool of distinct primes.
937 for i
in xrange(poolsz
):
939 q
= gen(nsubbits
, p
= p
)
940 if q
.p
not in qmap
: break
942 raise ExpectedError('insufficient diversity')
946 ## Work through combinations of factors from the pool.
947 for qq
in combinations(npiece
- 1, pool
):
949 ## Construct the product of the selected factors.
950 qsmall
= prod(q
.p
for q
in qq
)
952 ## Maybe we'll need to replace the large factor. Try not to do this
953 ## too often. DISP measures the large factor's performance at
954 ## producing candidates with the right length. If it looks bad then
955 ## we'll have to replace it.
956 if 3*disp
*disp
> nstep
*nstep
:
958 if disp
< 0: p
.p('<')
961 ## If we don't have a large factor, then make one.
963 qbig
= gen(nbits
- qsmall
.nbits
, p
= p
)
966 ## We have a candidate. Calculate it and make sure it has the right
968 n
= 2*qsmall
*qbig
.p
+ 1
970 if n
.nbits
< nbits
: disp
-= 1
971 elif n
.nbits
> nbits
: disp
+= 1
972 elif C
.PrimeFilter
.smallfactor(n
) == C
.PGEN_FAIL
: pass
975 ## The candidate has passed the small-primes test. Now check it
976 ## against Pocklington.
977 for a
in I
.islice(SIEVE
.smallprimes(), 16):
980 if pow(a
, n
- 1, n
) != 1: p
.p('.'); break
983 ## Find a generator of a sufficiently large subgroup.
984 if n
.gcd(pow(a
, (n
- 1)/qbig
.p
, n
) - 1) != 1: p
.p('@'); continue
987 if n
.gcd(pow(a
, (n
- 1)/q
.p
, n
) - 1) != 1:
988 p
.p('@'); ok
= False; break
993 ## Label the factors.
996 for i
, q
in enumerate(qq
): q
.label(qlfmt % i
)
998 ## Return the number we found.
999 p
.pop(); return PockPrime(n
, a
, qq
, label
= label
)
1001 ###--------------------------------------------------------------------------
1007 ## Prepare an option parser.
1008 op
= OP
.OptionParser(
1010 pock [-qv] [-s SIEVEBITS] CMD ARGS...
1014 description
= 'Generate or verify certified prime numbers.')
1015 op
.add_option('-v', '--verbose', dest
= 'verbosity',
1016 action
= 'count', default
= 1,
1017 help = 'print mysterious runes while looking for prime numbers')
1018 op
.add_option('-q', '--quiet', dest
= 'quietude',
1019 action
= 'count', default
= 0,
1020 help = 'be quiet while looking for prime numbers')
1021 op
.add_option('-s', '--sievebits', dest
= 'sievebits',
1022 type = 'int', default
= 32,
1023 help = 'size (in bits) of largest small prime')
1024 opts
, argv
= op
.parse_args()
1025 VERBOSITY
= opts
.verbosity
- opts
.quietude
1026 p
= ProgressReporter()
1027 a
= ArgFetcher(argv
, op
.error
)
1029 ## Process arguments and do what the user asked.
1033 ## Generate a prime with no special structure.
1034 initsieve(opts
.sievebits
)
1035 nbits
= a
.int(min = 4)
1037 p
= gen(nbits
, 'p', p
= p
)
1042 ## Generate a Lim--Lee prime.
1043 initsieve(opts
.sievebits
)
1044 nbits
= a
.int(min = 4)
1045 nsubbits
= a
.int(min = 4, max = nbits
)
1047 p
= gen_limlee(nbits
, nsubbits
, 'p', 'q_%d', p
= p
)
1052 ## Check an existing certificate.
1053 fn
= a
.arg(default
= '-', must
= False)
1054 if fn
== '-': f
= stdin
1055 else: f
= open(fn
, 'r')
1060 raise ExpectedError("unknown command `%s'" % w
)
1062 if __name__
== '__main__':
1063 prog
= OS
.path
.basename(argv
[0])
1065 except ExpectedError
: exit('%s: %s' %
(prog
, _excval().message
))
1066 except IOError: exit('%s: %s' %
(prog
, _excval()))
1068 ###----- That's all, folks --------------------------------------------------