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 import sys
as SYS
; from sys
import argv
, stdin
, stdout
, stderr
38 if SYS
.version_info
>= (3,):
40 range = lambda *args
: list(xrange(*args
))
42 ###--------------------------------------------------------------------------
46 """Return the most recent exception object."""
47 return SYS
.exc_info()[1]
49 class ExpectedError (Exception):
51 I represent an expected error, which should be reported in a friendly way.
55 def prod(ff
, one
= 1):
56 """Return ONE times the product of the elements of FF."""
57 return C
.MPMul().factor(one
).factor(ff
).done()
59 def parse_label(line
):
61 Split LINE at an `=' character and return the left and right halves.
63 The returned pieces have leading and trailing whitespace trimmed.
66 if eq
< 0: raise ExpectedError('expected `LABEL = ...\'')
67 return line
[:eq
].strip(), line
[eq
+ 1:].strip()
70 l
= s
.split(',', n
- 1)
71 if n
is not None and len(l
) != n
:
72 raise ExpectedError('expected `,\'-separated list of %d items' % n
)
76 """Convert S to a integer."""
77 try: return C
.MP(s
, 0)
78 except TypeError: raise ExpectedError('invalid integer `%s\'' % s
)
82 class ProgressReporter (object):
84 I keep users amused while the program looks for large prime numbers.
86 My main strategy is the printing of incomprehensible runes. I can be
87 muffled by lowering by verbosity level.
89 Prime searches are recursive in nature. When a new recursive level is
90 started, call `push'; and call `pop' when the level is finished. This must
91 be done around the top level too.
94 """Initialize the ProgressReporter."""
99 Update our idea of whether we're active.
101 We don't write inscrutable runes when inactive. The current policy is to
102 write nothing if verbosity is zero, to write runes for the top level only
103 if verbosity is 1, and to write runes always if verbosity is higher than
106 me
._active
= VERBOSITY
>= 2 or (VERBOSITY
== 1 and me
._level
== 0)
108 """Push a new search level."""
111 if me
._level
> 0: me
.p('[')
114 """Pop a search level."""
115 if me
._level
> 0: me
.p(']')
120 """Print CH as a progress rune."""
121 if me
._active
: stderr
.write(ch
); stderr
.flush()
123 def combinations(r
, v
):
125 Return an iterator which yields all combinations of R elements from V.
127 V must be an indexable sequence. The each combination is returned as a
128 list, containing elements from V in their original order.
131 ## Set up the selection vector. C will contain the indices of the items of
132 ## V we've selected for the current combination. At all times, C contains
133 ## a strictly increasing sequence of integers in the interval [0, N).
139 ## Yield up the current combination.
140 vv
= [v
[i
] for i
in c
]
143 ## Now advance to the next one. Find the last index in C which we can
144 ## increment subject to the rules. As we iterate downwards, i will
145 ## contain the index into C, and j will be the maximum acceptable value
146 ## for the corresponding item. We'll step the last index until it
147 ## reaches the limit, and then step the next one down, resetting the last
152 ## If i is zero here, then we've advanced everything as far as it will
156 ## Move down to the next index.
159 ## If this index isn't at its maximum value, then we've found the place
163 ## Step this index on by one, and set the following indices to the
164 ## immediately following values.
166 while i
< r
: c
[i
] = j
; i
+= 1; j
+= 1
168 class ArgFetcher (object):
170 I return arguments from a list, reporting problems when they occur.
172 def __init__(me
, argv
, errfn
):
174 Initialize, returning successive arguments from ARGV.
176 Errors are reported to ERRFN.
182 def arg(me
, default
= None, must
= True):
184 Return the next argument.
186 If MUST is true, then report an error (to the ERRFN) if there are no more
187 arguments; otherwise, return the DEFAULT.
189 if me
._i
>= me
._argc
:
190 if must
: me
._errfn('missing argument')
192 arg
= me
._argv
[me
._i
]; me
._i
+= 1
194 def int(me
, default
= None, must
= True, min = None, max = None):
196 Return the next argument converted to an integer.
198 If MUST is true, then report an error (to the ERRFN) if there are no more
199 arguments; otherwise return the DEFAULT. Report an error if the next
200 argument is not a valid integer, or if the integer is beyond the MIN and
203 arg
= me
.arg(default
= None, must
= must
)
204 if arg
is None: return default
206 except ValueError: me
._errfn('bad integer')
207 if (min is not None and arg
< min) or (max is not None and arg
> max):
208 me
._errfn('out of range')
211 ###--------------------------------------------------------------------------
212 ### Sieving for small primes.
214 class Sieve (object):
216 I represent a collection of small primes, up to some chosen limit.
218 The limit is available as the `limit' attribute. Let L be this limit;
219 then, if N < L^2 is some composite, then N has at least one prime factor
223 ## Figure out the number of bits in a (nonnegative) primitive `int'. We'll
224 ## use a list of these as our sieve.
225 try: _MAX
= SYS
.maxint
226 except AttributeError: _MAX
= SYS
.maxsize
228 while 1 << (_NBIT
+ 1) < _MAX
: _NBIT
+= 1
230 def __init__(me
, limit
):
232 Initialize a sieve holding all primes below LIMIT.
235 ## The sieve is maintained in the `_bits' attribute. This is a list of
236 ## integers, used as a bitmask: let 2 < n < L be an odd integer; then bit
237 ## (n - 3)/2 will be clear iff n is prime. Let W be the value of
238 ## `_NBIT', above; then bit W i + j in the sieve is stored in bit j of
241 ## Store the limit for later inspection.
244 ## Calculate the size of sieve we'll need and initialize the bit list.
246 sievesz
= (n
+ me
._NBIT
- 1)//me
._NBIT
247 me
._sievemax
= sievesz
*me
._NBIT
248 me
._bits
= sievesz
*[0]
250 ## This is standard Sieve of Eratosthenes. For each index i: if
251 ## bit i is clear, then p = 2 i + 3 is prime, so set the bits
252 ## corresponding to each multiple of p, i.e., bits (k p - 3)/2 =
253 ## (2 k i + 3 - 3)/2 = k i for k > 1.
254 for i
in xrange(me
._sievemax
):
255 if me
._bitp(i
): i
+= 1; continue
258 for j
in xrange(i
+ p
, me
._sievemax
, p
): me
._setbit(j
)
261 def _bitp(me
, i
): i
, j
= divmod(i
, me
._NBIT
); return (me
._bits
[i
] >> j
)&1
262 def _setbit(me
, i
): i
, j
= divmod(i
, me
._NBIT
); me
._bits
[i
] |
= 1 << j
266 Return an iterator over the known small primes.
271 for j
in xrange(me
._NBIT
):
272 if not (b
&1): yield n
275 ## We generate the sieve on demand.
278 def initsieve(sievebits
):
282 Ensure that it can be used to check the primality of numbers up to (but not
283 including) 2^SIEVEBITS.
286 if SIEVE
is not None: raise ValueError('sieve already defined')
287 if sievebits
< 6: sievebits
= 6
288 SIEVE
= Sieve(1 << (sievebits
+ 1)//2)
290 ###--------------------------------------------------------------------------
291 ### Primality checking.
295 Check that P is a small prime.
297 If not, raise an `ExpectedError'. The `SIEVE' variable must have been
300 if p
< 2: raise ExpectedError('%d too small' % p
)
301 if SIEVE
.limit
*SIEVE
.limit
< p
:
302 raise ExpectedError('%d too large for small prime' % p
)
303 for q
in SIEVE
.smallprimes():
305 if p
%q
== 0: raise ExpectedError('%d divides %d' %
(q
, p
))
307 def pock_test(p
, a
, qq
):
309 Check that P is prime using Pocklington's criterion.
311 If not, raise an `ExpectedError'.
313 Let Q be the product of the elements of the sequence QQ. The test works as
314 follows. Suppose p is the smallest prime factor of P. If A^{P-1} /== 1
315 (mod P) then P is certainly composite (Fermat's test); otherwise, we have
316 established that the order of A in (Z/pZ)^* divides P - 1. Next, let t =
317 A^{(P-1)/q} for some prime factor q of Q, and let g = gcd(t - 1, P). If g
318 = P then the proof is inconclusive; if 1 < g < P then g is a nontrivial
319 factor of P, so P is composite; otherwise, t has order q in (Z/pZ)^*, so
320 (Z/pZ)^* contains a subgroup of size q, and therefore q divides p - 1. If
321 QQ is a sequence of distinct primes, and the preceding criterion holds for
322 all q in QQ, then Q divides p - 1. If Q^2 < P then the proof is
323 inconclusive; otherwise, let p' be any prime dividing P/p. Then p' >= p >
324 Q, so p p' > Q^2 > P; but p p' divides P, so this is a contradiction.
325 Therefore P/p has no prime factors, and P is prime.
328 ## We don't actually need the distinctness criterion. Suppose that q^e
329 ## divides Q. Then gcd(t - 1, P) = 1 implies that A^{(P-1)/q^{e-1}} has
330 ## order q^e in (Z/pZ)^*, which accounts for the multiplicity.
333 if p
< 2: raise ExpectedError('%d too small' % p
)
335 raise ExpectedError('too few Pocklington factors for %d' % p
)
336 if pow(a
, p
- 1, p
) != 1:
337 raise ExpectedError('%d is Fermat witness for %d' %
(a
, p
))
340 raise ExpectedError('duplicate Pocklington factor %d for %d' %
(q
, p
))
341 g
= p
.gcd(pow(a
, (p
- 1)/q
, p
) - 1)
343 raise ExpectedError('%d order not multiple of %d mod %d' %
(a
, q
, p
))
345 raise ExpectedError('%d divides %d' %
(g
, p
))
347 def ecpp_test(p
, a
, b
, x
, y
, qq
):
349 Check that P is prime using Goldwasser and Kilian's ECPP method.
351 If not, raise an `ExpectedError'.
353 Let Q be the product of the elements of the sequence QQ. Suppose p is the
354 smallest prime factor of P. Let g = gcd(4 A^3 + 27 B^2, P). If g = P then
355 the test is inconclusive; otherwise, if g /= 1 then g is a nontrivial
356 factor of P. Define E(GF(p)) = { (x, y) | y^2 = x^3 + A x + B } U { inf }
357 to be the elliptic curve over p with short-Weierstraß coefficients A and B;
358 we have just checked that this curve is not singular. If R = (X, Y) is not
359 a point on this curve, then the test is inconclusive. If Q R is not the
360 point at infinity, then the test fails; otherwise we deduce that P has
361 Q-torsion in E. Let S = (Q/q) R for some prime factor q of Q. If S is the
362 point at infinity then the test is inconclusive; otherwise, q divides the
363 order of S in E. If QQ is a sequence of distinct primes, and the preceding
364 criterion holds for all q in QQ, then Q divides the order of S. Therefore
365 #E(p) >= Q. If Q <= (qrrt(P) + 1)^2 then the test is inconclusive.
366 Otherwise, Hasse's theorem tells us that |p + 1 - #E(p)| <= 2 sqrt(p);
367 hence we must have p + 1 + 2 sqrt(p) = (sqrt(p) + 1)^2 >= #E(p) >= Q >
368 (qrrt(P) + 1)^2; so sqrt(p) + 1 > qrrt(P) + 1, i.e., p^2 > P. As for
369 Pocklington above, if p' is any prime factor of P/p, then p p' >= p^2 > P,
370 which is a contradiction, and we conclude that P is prime.
373 ## This isn't going to work if gcd(P, 6) /= 1: we're going to use the
374 ## large-characteristic addition formulae.
376 if g
!= 1: raise ExpectedError('%d divides %d' %
(g
, p
))
378 ## We want to check that Q > (qrrt(P) + 1)^2 iff sqrt(Q) > qrrt(P) + 1; but
379 ## calculating square roots is not enjoyable (partly because we have to
380 ## deal with the imprecision). Fortunately, some algebra will help: the
381 ## condition holds iff qrrt(P) < sqrt(Q) - 1 iff P < Q^2 - 4 Q sqrt(Q) +
382 ## 6 Q - 4 sqrt(Q) + 1 = Q (Q + 6) + 1 - 4 sqrt(Q) (Q + 1) iff Q (Q + 6) -
383 ## P + 1 > 4 sqrt(Q) (Q + 1) iff (Q (Q + 6) - P + 1)^2 > 16 Q (Q + 1)^2
385 t
, u
= Q
*(Q
+ 6) - p
+ 1, 4*(Q
+ 1)
386 if t
*t
<= Q
*u
*u
: raise ExpectedError('too few subgroups for ECPP')
388 ## Construct the curve.
389 E
= C
.PrimeField(p
).ec(a
, b
) # careful: may not be a prime!
391 ## Find the base point.
394 raise ExpectedError('(%d, %d) is not on the curve' %
(x
, y
))
396 ## Check that it has Q-torsion.
397 if Q
*R
: raise ExpectedError('(%d, %d) not a %d-torsion point' %
(x
, y
, Q
))
399 ## Now check the individual factors.
402 raise ExpectedError('duplicate ECPP factor %d for %d' %
(q
, p
))
405 raise ExpectedError('(%d, %d) order not a multiple of %d' %
(x
, y
, q
))
408 raise ExpectedError('%d divides %d' %
(g
, p
))
410 ###--------------------------------------------------------------------------
411 ### Proof steps and proofs.
413 class BaseStep (object):
415 I'm a step in a primality proof.
417 I assert that a particular number is prime, and can check this.
419 This class provides basic protocol for proof steps, mostly to do with
422 The step's label is kept in its `label' attribute. It can be set by the
423 constructor, and is `None' by default. Users can modify this attribute if
424 they like. Labels beginning `$' are assumed to be internal and
425 uninteresting; other labels cause `check' lines to be written to the output
426 listing the actual number of interest.
428 Protocol that proof steps should provide:
430 label A string labelling the proof step and the associated prime
433 p The prime number which this step proves to be prime.
435 check() Check that the proof step is actually correct, assuming that
436 any previous steps have already been verified.
438 out(FILE) Write an appropriate encoding of the proof step to the output
441 def __init__(me
, label
= None, *arg
, **kw
):
442 """Initialize a proof step, setting a default label if necessary."""
443 super(BaseStep
, me
).__init__(*arg
, **kw
)
447 Write the proof step to an output FILE.
449 Subclasses must implement a method `_out' which actually does the work.
450 Here, we write a `check' line to verify that the proof actually applies
451 to the number we wanted, if the label says that this is an interesting
455 if me
.label
is not None and not me
.label
.startswith('$'):
456 file.write('check %s, %d, %d\n' %
(me
.label
, me
.p
.nbits
, me
.p
))
458 class SmallStep (BaseStep
):
460 I represent a claim that a number is a small prime.
462 Such claims act as the base cases in a complicated primality proof. When
463 verifying, the claim is checked by trial division using a collection of
466 def __init__(me
, pp
, p
, *arg
, **kw
):
468 Initialize a small-prime step.
470 PP is the overall PrimeProof object of which this is a step; P is the
471 small number whose primality is asserted.
473 super(SmallStep
, me
).__init__(*arg
, **kw
)
476 """Check that the number is indeed a small prime."""
477 return small_test(me
.p
)
479 """Write a small-prime step to the FILE."""
480 file.write('small %s = %d\n' %
(me
.label
, me
.p
))
481 def __repr__(me
): return 'SmallStep(%d)' %
(me
.p
)
483 def parse(cls
, pp
, line
):
485 Parse a small-prime step from a LINE in a proof file.
487 SMALL-STEP ::= `small' LABEL `=' P
489 PP is a PrimeProof object holding the results from the previous steps.
491 if SIEVE
is None: raise ExpectedError('missing `sievebits\' line')
492 label
, p
= parse_label(line
)
493 return cls(pp
, conv_int(p
), label
= label
)
495 class PockStep (BaseStep
):
497 I represent a Pocklington certificate for a number.
499 The number is not explicitly represented in a proof file. See `pock_test'
500 for the underlying mathematics.
502 def __init__(me
, pp
, a
, R
, qqi
, *arg
, **kw
):
504 Inititialize a Pocklington step.
506 PP is the overall PrimeProof object of which this is a step; A is the
507 generator of a substantial subgroup of units; R is a cofactor; and QQI is
508 a sequence of labels for previous proof steps. If Q is the product of
509 the primes listed in QQI, then the number whose primality is asserted is
512 super(PockStep
, me
).__init__(*arg
, **kw
)
516 me
._qq
= [pp
.get_step(qi
).p
for qi
in qqi
]
517 me
.p
= prod(me
._qq
, 2*R
) + 1
519 """Verify a proof step based on Pocklington's theorem."""
520 return pock_test(me
.p
, me
._a
, me
._qq
)
522 """Write a Pocklington step to the FILE."""
523 file.write('pock %s = %d, %d, [%s]\n' % \
525 me
._R
, ', '.join('%s' % qi
for qi
in me
._qqi
)))
526 def __repr__(me
): return 'PockStep(%d, %d, %s)' %
(me
._a
, me
._R
, me
._qqi
)
528 def parse(cls
, pp
, line
):
530 Parse a Pocklington step from a LINE in a proof file.
532 POCK-STEP ::= `pock' LABEL `=' A `,' R `,' `[' Q-LIST `]'
533 Q-LIST ::= Q [`,' Q-LIST]
535 PP is a PrimeProof object holding the results from the previous steps.
537 label
, rest
= parse_label(line
)
538 a
, R
, qq
= parse_list(rest
, 3)
540 if not qq
.startswith('[') or not qq
.endswith(']'):
541 raise ExpectedError('missing `[...]\' around Pocklington factors')
542 return cls(pp
, conv_int(a
), conv_int(R
),
543 [q
.strip() for q
in qq
[1:-1].split(',')], label
= label
)
545 class ECPPStep (BaseStep
):
547 I represent a Goldwasser--Kilian ECPP certificate for a number.
549 def __init__(me
, pp
, p
, a
, b
, x
, y
, qqi
, *arg
, **kw
):
551 Inititialize an ECPP step.
553 PP is the overall PrimeProof object of which this is a step; P is the
554 number whose primality is asserted; A and B are the short Weierstraß
555 curve coefficients; X and Y are the base point coordinates; and QQI is a
556 sequence of labels for previous proof steps.
558 super(ECPPStep
, me
).__init__(*arg
, **kw
)
562 me
._qq
= [pp
.get_step(qi
).p
for qi
in qqi
]
565 """Verify a proof step based on Goldwasser and Kilian's theorem."""
566 return ecpp_test(me
.p
, me
._a
, me
._b
, me
._x
, me
._y
, me
._qq
)
568 """Write an ECPP step to the FILE."""
569 file.write('ecpp %s = %d, %d, %d, %d, %d, [%s]\n' % \
570 (me
.label
, me
.p
, me
._a
, me
._b
, me
._x
, me
._y
,
571 ', '.join('%s' % qi
for qi
in me
._qqi
)))
573 return 'ECPPstep(%d, %d, %d, %d, %d, %s)' % \
574 (me
.p
, me
._a
, me
._b
, me
._x
, me
._y
, me
._qqi
)
576 def parse(cls
, pp
, line
):
578 Parse an ECPP step from a LINE in a proof file.
580 ECPP-STEP ::= `ecpp' LABEL `=' P `,' A `,' B `,' X `,' Y `,'
582 Q-LIST ::= Q [`,' Q-LIST]
584 PP is a PrimeProof object holding the results from the previous steps.
586 label
, rest
= parse_label(line
)
587 p
, a
, b
, x
, y
, qq
= parse_list(rest
, 6)
589 if not qq
.startswith('[') or not qq
.endswith(']'):
590 raise ExpectedError('missing `[...]\' around ECPP factors')
591 return cls(pp
, conv_int(p
), conv_int(a
), conv_int(b
),
592 conv_int(x
), conv_int(y
),
593 [q
.strip() for q
in qq
[1:-1].split(',')], label
= label
)
597 Handle a `check' line in a proof file.
599 CHECK ::= `check' LABEL, B, N
601 Verify that the proof step with the given LABEL asserts the primality of
602 the integer N, and that 2^{B-1} <= N < 2^B.
604 label
, nb
, p
= parse_list(line
, 3)
605 label
, nb
, p
= label
.strip(), conv_int(nb
), conv_int(p
)
606 pi
= pp
.get_step(label
).p
608 raise ExpectedError('check failed: %s = %d /= %d' %
(label
, pi
, p
))
610 raise ExpectedError('check failed: nbits(%s) = %d /= %d' % \
611 (label
, p
.nbits
, nb
))
612 if VERBOSITY
: print(';; %s = %d [%d]' %
(label
, p
, nb
))
614 def setsievebits(pp
, line
):
616 Handle a `sievebits' line in a proof file.
618 SIEVEBITS ::= `sievebits' N
620 Ensure that the verifier is willing to accept small primes up to 2^N.
624 class PrimeProof (object):
626 I represent a proof of primality for one or more numbers.
628 I can encode my proof as a line-oriented text file, in a simple format, and
629 read such a proof back to check it.
632 ## A table to dispatch on keywords read from a file.
633 STEPMAP
= { 'small': SmallStep
.parse
,
634 'pock': PockStep
.parse
,
635 'ecpp': ECPPStep
.parse
,
636 'sievebits': setsievebits
,
641 Initialize a proof object.
643 me
._steps
= {} # Maps labels to steps.
644 me
._stepseq
= [] # Sequence of labels, in order.
645 me
._pmap
= {} # Maps primes to steps.
648 def addstep(me
, step
):
650 Add a new STEP to the proof.
652 The STEP may have a label already. If not, a new internal label is
653 chosen. The proof step is checked before being added to the proof. The
657 ## If there's already a step for this prime, and the new step doesn't
658 ## have a label, then return the old one instead.
659 if step
.label
is None:
660 try: return me
._pmap
[step
.p
]
661 except KeyError: pass
663 ## Make sure the step is actually correct.
666 ## Generate a label if the step doesn't have one already.
667 if step
.label
is None: step
.label
= '$t%d' % me
._i
; me
._i
+= 1
669 ## If the label is already taken then we have a problem.
670 if step
.label
in me
._steps
:
671 raise ExpectedError('duplicate label `%s\'' % step
.label
)
673 ## Store the proof step.
674 me
._pmap
[step
.p
] = step
.label
675 me
._steps
[step
.label
] = step
676 me
._stepseq
.append(step
.label
)
679 def get_step(me
, label
):
681 Check that LABEL labels a known step, and return that step.
683 try: return me
._steps
[label
]
684 except KeyError: raise ExpectedError('unknown label `%s\'' % label
)
688 Write the proof to the given FILE.
691 ## Prefix the main steps with a `sievebits' line.
692 file.write('sievebits %d\n' %
(2*(SIEVE
.limit
.bit_length() - 1)))
694 ## Write the steps out one by one.
695 for label
in me
._stepseq
: me
._steps
[label
].out(file)
699 Read a proof from a given FILE.
701 FILE ::= {STEP | CHECK | SIEVEBITS} [FILE]
702 STEP ::= SMALL-STEP | POCK-STEP
704 Comments (beginning `;') and blank lines are ignored. Other lines begin
708 for lno
, line
in enumerate(file, 1):
710 if line
.startswith(';'): continue
711 ww
= line
.split(None, 1)
714 if len(ww
) > 1: tail
= ww
[1]
717 try: op
= me
.STEPMAP
[w
]
719 raise ExpectedError('unrecognized keyword `%s\'' % w
)
724 except ExpectedError
:
725 raise ExpectedError('%s:%d: %s' %
726 (file.name
, lno
, _excval().message
))
729 ###--------------------------------------------------------------------------
730 ### Finding provable primes.
732 class BasePrime (object):
734 I represent a prime number which has been found and can be proven.
736 This object can eventually be turned into a sequence of proof steps and
737 added to a PrimeProof. This isn't done immediately, because some
738 prime-search strategies want to build a pool of provable primes and will
739 then select some subset of them to actually construct the number of final
740 interest. This way, we avoid cluttering the output proof with proofs of
741 uninteresting numbers.
745 p The prime number in question.
747 label(LABEL) Associate LABEL with this prime, and the corresponding proof
748 step. A label can be set in the constructor, or later using
751 register(PP) Register the prime with a PrimeProof, adding any necessary
752 proof steps. Returns the label of the proof step for this
756 Return a proof step for this prime.
758 def __init__(me
, label
= None, *args
, **kw
):
759 """Initialize a provable prime number object."""
760 super(BasePrime
, me
).__init__(*args
, **kw
)
761 me
._index
= me
._pp
= None
763 def label(me
, label
):
764 """Set this number's LABEL."""
766 def register(me
, pp
):
768 Register the prime's proof steps with PrimeProof PP.
770 Return the final step's label.
772 if me
._pp
is not None:
776 me
._index
= pp
.addstep(me
._mkstep(pp
, label
= me
._label
))
777 ##try: me._index = pp.addstep(me._mkstep(pp, label = me._label))
778 ##except: raise RuntimeError('generated proof failed sanity check')
781 class SmallPrime (BasePrime
):
782 """I represent a prime small enough to be checked in isolation."""
783 def __init__(me
, p
, *args
, **kw
):
784 super(SmallPrime
, me
).__init__(*args
, **kw
)
786 def _mkstep(me
, pp
, **kw
):
787 return SmallStep(pp
, me
.p
, **kw
)
789 class PockPrime (BasePrime
):
790 """I represent a prime proven using Pocklington's theorem."""
791 def __init__(me
, p
, a
, qq
, *args
, **kw
):
792 super(PockPrime
, me
).__init__(*args
, **kw
)
796 def _mkstep(me
, pp
, **kw
):
797 return PockStep(pp
, me
._a
, (me
.p
- 1)/prod((q
.p
for q
in me
._qq
), 2),
798 [q
.register(pp
) for q
in me
._qq
], **kw
)
800 def gen_small(nbits
, label
= None, p
= None):
802 Return a new small prime.
804 The prime will be exactly NBITS bits long. The proof step will have the
805 given LABEL attached. Report progress to the ProgressReporter P.
809 ## Pick a random NBITS-bit number.
810 n
= C
.rand
.mp(nbits
, 1)
811 assert n
.nbits
== nbits
813 ## If it's probably prime, then check it against the small primes we
814 ## know. If it passes then we're done. Otherwise, try again.
816 for q
in SIEVE
.smallprimes():
817 if q
*q
> n
: return SmallPrime(n
, label
= label
)
820 def gen_pock(nbits
, nsubbits
= 0, label
= None, p
= ProgressReporter()):
822 Return a new prime provable using Pocklington's theorem.
824 The prime N will be exactly NBITS long, of the form N = 2 Q R + 1. If
825 NSUBBITS is nonzero, then each prime factor of Q will be NSUBBITS bits
826 long; otherwise a suitable default will be chosen. The proof step will
827 have the given LABEL attached. Report progress to the ProgressReporter P.
829 The prime numbers this function returns are a long way from being uniformly
833 ## Pick a suitable value for NSUBBITS if we don't have one.
836 ## This is remarkably tricky. Picking about 1/3 sqrt(NBITS) factors
837 ## seems about right for large numbers, but there's serious trouble
838 ## lurking for small sizes.
839 nsubbits
= int(3*M
.sqrt(nbits
))
840 if nbits
< nsubbits
+ 3: nsubbits
= nbits
//2 + 1
841 if nbits
== 2*nsubbits
: nsubbits
+= 1
843 ## Figure out how many subgroups we'll need.
844 npiece
= ((nbits
+ 1)//2 + nsubbits
- 1)//nsubbits
850 ## Come up with a collection of known prime factors.
851 p
.p('!'); qq
= [gen(nsubbits
, p
= p
) for i
in xrange(npiece
)]
852 Q
= prod(q
.p
for q
in qq
)
854 ## Come up with bounds on the cofactor. If we're to have N = 2 Q R + 1,
855 ## and 2^{B-1} <= N < 2^B, then we must have 2^{B-2}/Q <= R < 2^{B-1}/Q.
856 Rbase
= (C
.MP(0).setbit(nbits
- 2) + Q
- 1)//Q
857 Rwd
= C
.MP(0).setbit(nbits
- 2)//Q
859 ## Probe the available space of cofactors. If the space is kind of
860 ## narrow, then we want to give up quickly if we're not finding anything
866 ## Pick a random cofactor and examine the number we ended up with.
867 ## Make sure it really does have the length we expect.
868 R
= C
.rand
.range(Rwd
) + Rbase
870 assert n
.nbits
== nbits
872 ## As a complication, if NPIECE is 1, it's just about possible that Q^2
873 ## <= n, in which case this isn't going to work.
876 ## If n has small factors, then pick another cofactor.
877 if C
.PrimeFilter
.smallfactor(n
) == C
.PGEN_FAIL
: continue
879 ## Work through the small primes to find a suitable generator. The
880 ## value 2 is almost always acceptable, so don't try too hard here.
881 for a
in I
.islice(SIEVE
.smallprimes(), 16):
883 ## First, try the Fermat test. If that fails, then n is definitely
885 if pow(a
, n
- 1, n
) != 1: p
.p('.'); break
888 ## Work through the subgroup orders, checking that suitable powers of
889 ## a generate the necessary subgroups.
891 if n
.gcd(pow(a
, (n
- 1)/q
.p
, n
) - 1) != 1:
892 p
.p('@'); ok
= False; break
897 if ok
: p
.pop(); return PockPrime(n
, a
, qq
, label
= label
)
899 def gen(nbits
, label
= None, p
= ProgressReporter()):
901 Generate a prime number with NBITS bits.
903 Give it the LABEL, and report progress to P.
905 if SIEVE
.limit
>> (nbits
+ 1)//2: g
= gen_small
907 return g(nbits
, label
= label
, p
= p
)
909 def gen_limlee(nbits
, nsubbits
,
910 label
= None, qlfmt
= None, p
= ProgressReporter()):
912 Generate a Lim--Lee prime with NBITS bits.
914 Let p be the prime. Then we'll have p = 2 q_0 q_1 ... q_k, with all q_i at
915 least NSUBBITS bits long, and all but q_k exactly that long.
917 The prime will be given the LABEL; progress is reported to P. The factors
918 q_i will be labelled by filling in the `printf'-style format string QLFMT
922 ## Figure out how many factors (p - 1)/2 will have.
923 npiece
= nbits
//nsubbits
924 if npiece
< 2: raise ExpectedError('too few pieces')
926 ## Decide how big to make the pool of factors.
927 poolsz
= max(3*npiece
+ 5, 25) # Heuristic from GnuPG
929 ## Prepare for the main loop.
934 ## Try to make a prime.
938 ## Construct a pool of NSUBBITS-size primes. There's a problem with very
939 ## small sizes: we might not be able to build a pool of distinct primes.
941 for i
in xrange(poolsz
):
943 q
= gen(nsubbits
, p
= p
)
944 if q
.p
not in qmap
: break
946 raise ExpectedError('insufficient diversity')
950 ## Work through combinations of factors from the pool.
951 for qq
in combinations(npiece
- 1, pool
):
953 ## Construct the product of the selected factors.
954 qsmall
= prod(q
.p
for q
in qq
)
956 ## Maybe we'll need to replace the large factor. Try not to do this
957 ## too often. DISP measures the large factor's performance at
958 ## producing candidates with the right length. If it looks bad then
959 ## we'll have to replace it.
960 if 3*disp
*disp
> nstep
*nstep
:
962 if disp
< 0: p
.p('<')
965 ## If we don't have a large factor, then make one.
967 qbig
= gen(nbits
- qsmall
.nbits
, p
= p
)
970 ## We have a candidate. Calculate it and make sure it has the right
972 n
= 2*qsmall
*qbig
.p
+ 1
974 if n
.nbits
< nbits
: disp
-= 1
975 elif n
.nbits
> nbits
: disp
+= 1
976 elif C
.PrimeFilter
.smallfactor(n
) == C
.PGEN_FAIL
: pass
979 ## The candidate has passed the small-primes test. Now check it
980 ## against Pocklington.
981 for a
in I
.islice(SIEVE
.smallprimes(), 16):
984 if pow(a
, n
- 1, n
) != 1: p
.p('.'); break
987 ## Find a generator of a sufficiently large subgroup.
988 if n
.gcd(pow(a
, (n
- 1)/qbig
.p
, n
) - 1) != 1: p
.p('@'); continue
991 if n
.gcd(pow(a
, (n
- 1)/q
.p
, n
) - 1) != 1:
992 p
.p('@'); ok
= False; break
997 ## Label the factors.
1000 for i
, q
in enumerate(qq
): q
.label(qlfmt % i
)
1002 ## Return the number we found.
1003 p
.pop(); return PockPrime(n
, a
, qq
, label
= label
)
1005 ###--------------------------------------------------------------------------
1011 ## Prepare an option parser.
1012 op
= OP
.OptionParser(
1014 pock [-qv] [-s SIEVEBITS] CMD ARGS...
1018 description
= 'Generate or verify certified prime numbers.')
1019 op
.add_option('-v', '--verbose', dest
= 'verbosity',
1020 action
= 'count', default
= 1,
1021 help = 'print mysterious runes while looking for prime numbers')
1022 op
.add_option('-q', '--quiet', dest
= 'quietude',
1023 action
= 'count', default
= 0,
1024 help = 'be quiet while looking for prime numbers')
1025 op
.add_option('-s', '--sievebits', dest
= 'sievebits',
1026 type = 'int', default
= 32,
1027 help = 'size (in bits) of largest small prime')
1028 opts
, argv
= op
.parse_args()
1029 VERBOSITY
= opts
.verbosity
- opts
.quietude
1030 p
= ProgressReporter()
1031 a
= ArgFetcher(argv
, op
.error
)
1033 ## Process arguments and do what the user asked.
1037 ## Generate a prime with no special structure.
1038 initsieve(opts
.sievebits
)
1039 nbits
= a
.int(min = 4)
1041 p
= gen(nbits
, 'p', p
= p
)
1046 ## Generate a Lim--Lee prime.
1047 initsieve(opts
.sievebits
)
1048 nbits
= a
.int(min = 4)
1049 nsubbits
= a
.int(min = 4, max = nbits
)
1051 p
= gen_limlee(nbits
, nsubbits
, 'p', 'q_%d', p
= p
)
1056 ## Check an existing certificate.
1057 fn
= a
.arg(default
= '-', must
= False)
1058 if fn
== '-': f
= stdin
1059 else: f
= open(fn
, 'r')
1064 raise ExpectedError("unknown command `%s'" % w
)
1066 if __name__
== '__main__':
1067 prog
= OS
.path
.basename(argv
[0])
1069 except ExpectedError
: exit('%s: %s' %
(prog
, _excval().message
))
1070 except IOError: exit('%s: %s' %
(prog
, _excval()))
1072 ###----- That's all, folks --------------------------------------------------