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):
53 Return ONE times the product of the elements of FF.
55 This is not done very efficiently.
57 return reduce(lambda prod
, f
: prod
*f
, ff
, one
)
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.
226 while type(1 << (_NBIT
+ 1)) == int: _NBIT
+= 1
228 def __init__(me
, limit
):
230 Initialize a sieve holding all primes below LIMIT.
233 ## The sieve is maintained in the `_bits' attribute. This is a list of
234 ## integers, used as a bitmask: let 2 < n < L be an odd integer; then bit
235 ## (n - 3)/2 will be clear iff n is prime. Let W be the value of
236 ## `_NBIT', above; then bit W i + j in the sieve is stored in bit j of
239 ## Store the limit for later inspection.
242 ## Calculate the size of sieve we'll need and initialize the bit list.
244 sievesz
= (n
+ me
._NBIT
- 1)/me
._NBIT
245 me
._sievemax
= sievesz
*me
._NBIT
246 me
._bits
= sievesz
*[0]
248 ## This is standard Sieve of Eratosthenes. For each index i: if
249 ## bit i is clear, then p = 2 i + 3 is prime, so set the bits
250 ## corresponding to each multiple of p, i.e., bits (k p - 3)/2 =
251 ## (2 k i + 3 - 3)/2 = k i for k > 1.
252 for i
in xrange(me
._sievemax
):
253 if me
._bitp(i
): i
+= 1; continue
256 for j
in xrange(i
+ p
, me
._sievemax
, p
): me
._setbit(j
)
259 def _bitp(me
, i
): i
, j
= divmod(i
, me
._NBIT
); return (me
._bits
[i
] >> j
)&1
260 def _setbit(me
, i
): i
, j
= divmod(i
, me
._NBIT
); me
._bits
[i
] |
= 1 << j
264 Return an iterator over the known small primes.
269 for j
in xrange(me
._NBIT
):
270 if not (b
&1): yield n
273 ## We generate the sieve on demand.
276 def initsieve(sievebits
):
280 Ensure that it can be used to check the primality of numbers up to (but not
281 including) 2^SIEVEBITS.
284 if SIEVE
is not None: raise ValueError('sieve already defined')
285 if sievebits
< 6: sievebits
= 6
286 SIEVE
= Sieve(1 << (sievebits
+ 1)/2)
288 ###--------------------------------------------------------------------------
289 ### Primality checking.
293 Check that P is a small prime.
295 If not, raise an `ExpectedError'. The `SIEVE' variable must have been
298 if p
< 2: raise ExpectedError('%d too small' % p
)
299 if SIEVE
.limit
*SIEVE
.limit
< p
:
300 raise ExpectedError('%d too large for small prime' % p
)
301 for q
in SIEVE
.smallprimes():
303 if p
%q
== 0: raise ExpectedError('%d divides %d' %
(q
, p
))
305 def pock_test(p
, a
, qq
):
307 Check that P is prime using Pocklington's criterion.
309 If not, raise an `ExpectedError'.
311 Let Q be the product of the elements of the sequence QQ. The test works as
312 follows. Suppose p is the smallest prime factor of P. If A^{P-1} /== 1
313 (mod P) then P is certainly composite (Fermat's test); otherwise, we have
314 established that the order of A in (Z/pZ)^* divides P - 1. Next, let t =
315 A^{(P-1)/q} for some prime factor q of Q, and let g = gcd(t - 1, P). If g
316 = P then the proof is inconclusive; if 1 < g < P then g is a nontrivial
317 factor of P, so P is composite; otherwise, t has order q in (Z/pZ)^*, so
318 (Z/pZ)^* contains a subgroup of size q, and therefore q divides p - 1. If
319 QQ is a sequence of distinct primes, and the preceding criterion holds for
320 all q in QQ, then Q divides p - 1. If Q^2 < P then the proof is
321 inconclusive; otherwise, let p' be any prime dividing P/p. Then p' >= p >
322 Q, so p p' > Q^2 > P; but p p' divides P, so this is a contradiction.
323 Therefore P/p has no prime factors, and P is prime.
326 ## We don't actually need the distinctness criterion. Suppose that q^e
327 ## divides Q. Then gcd(t - 1, P) = 1 implies that A^{(P-1)/q^{e-1}} has
328 ## order q^e in (Z/pZ)^*, which accounts for the multiplicity.
331 if p
< 2: raise ExpectedError('%d too small' % p
)
333 raise ExpectedError('too few Pocklington factors for %d' % p
)
334 if pow(a
, p
- 1, p
) != 1:
335 raise ExpectedError('%d is Fermat witness for %d' %
(a
, p
))
338 raise ExpectedError('duplicate Pocklington factor %d for %d' %
(q
, p
))
339 g
= p
.gcd(pow(a
, (p
- 1)/q
, p
) - 1)
341 raise ExpectedError('%d order not multiple of %d mod %d' %
(a
, q
, p
))
343 raise ExpectedError('%d divides %d' %
(g
, p
))
345 def ecpp_test(p
, a
, b
, x
, y
, qq
):
347 Check that P is prime using Goldwasser and Kilian's ECPP method.
349 If not, raise an `ExpectedError'.
351 Let Q be the product of the elements of the sequence QQ. Suppose p is the
352 smallest prime factor of P. Let g = gcd(4 A^3 + 27 B^2, P). If g = P then
353 the test is inconclusive; otherwise, if g /= 1 then g is a nontrivial
354 factor of P. Define E(GF(p)) = { (x, y) | y^2 = x^3 + A x + B } U { inf }
355 to be the elliptic curve over p with short-Weierstraß coefficients A and B;
356 we have just checked that this curve is not singular. If R = (X, Y) is not
357 a point on this curve, then the test is inconclusive. If Q R is not the
358 point at infinity, then the test fails; otherwise we deduce that P has
359 Q-torsion in E. Let S = (Q/q) R for some prime factor q of Q. If S is the
360 point at infinity then the test is inconclusive; otherwise, q divides the
361 order of S in E. If QQ is a sequence of distinct primes, and the preceding
362 criterion holds for all q in QQ, then Q divides the order of S. Therefore
363 #E(p) >= Q. If Q <= (qrrt(P) + 1)^2 then the test is inconclusive.
364 Otherwise, Hasse's theorem tells us that |p + 1 - #E(p)| <= 2 sqrt(p);
365 hence we must have p + 1 + 2 sqrt(p) = (sqrt(p) + 1)^2 >= #E(p) >= Q >
366 (qrrt(P) + 1)^2; so sqrt(p) + 1 > qrrt(P) + 1, i.e., p^2 > P. As for
367 Pocklington above, if p' is any prime factor of P/p, then p p' >= p^2 > P,
368 which is a contradiction, and we conclude that P is prime.
371 ## This isn't going to work if gcd(P, 6) /= 1: we're going to use the
372 ## large-characteristic addition formulae.
374 if g
!= 1: raise ExpectedError('%d divides %d' %
(g
, p
))
376 ## We want to check that Q > (qrrt(P) + 1)^2 iff sqrt(Q) > qrrt(P) + 1; but
377 ## calculating square roots is not enjoyable (partly because we have to
378 ## deal with the imprecision). Fortunately, some algebra will help: the
379 ## condition holds iff qrrt(P) < sqrt(Q) - 1 iff P < Q^2 - 4 Q sqrt(Q) +
380 ## 6 Q - 4 sqrt(Q) + 1 = Q (Q + 6) + 1 - 4 sqrt(Q) (Q + 1) iff Q (Q + 6) -
381 ## P + 1 > 4 sqrt(Q) (Q + 1) iff (Q (Q + 6) - P + 1)^2 > 16 Q (Q + 1)^2
383 t
, u
= Q
*(Q
+ 6) - p
+ 1, 4*(Q
+ 1)
384 if t
*t
<= Q
*u
*u
: raise ExpectedError('too few subgroups for ECPP')
386 ## Construct the curve.
387 E
= C
.PrimeField(p
).ec(a
, b
) # careful: may not be a prime!
389 ## Find the base point.
392 raise ExpectedError('(%d, %d) is not on the curve' %
(x
, y
))
394 ## Check that it has Q-torsion.
395 if Q
*R
: raise ExpectedError('(%d, %d) not a %d-torsion point' %
(x
, y
, Q
))
397 ## Now check the individual factors.
400 raise ExpectedError('duplicate ECPP factor %d for %d' %
(q
, p
))
403 raise ExpectedError('(%d, %d) order not a multiple of %d' %
(x
, y
, q
))
406 raise ExpectedError('%d divides %d' %
(g
, p
))
408 ###--------------------------------------------------------------------------
409 ### Proof steps and proofs.
411 class BaseStep (object):
413 I'm a step in a primality proof.
415 I assert that a particular number is prime, and can check this.
417 This class provides basic protocol for proof steps, mostly to do with
420 The step's label is kept in its `label' attribute. It can be set by the
421 constructor, and is `None' by default. Users can modify this attribute if
422 they like. Labels beginning `$' are assumed to be internal and
423 uninteresting; other labels cause `check' lines to be written to the output
424 listing the actual number of interest.
426 Protocol that proof steps should provide:
428 label A string labelling the proof step and the associated prime
431 p The prime number which this step proves to be prime.
433 check() Check that the proof step is actually correct, assuming that
434 any previous steps have already been verified.
436 out(FILE) Write an appropriate encoding of the proof step to the output
439 def __init__(me
, label
= None, *arg
, **kw
):
440 """Initialize a proof step, setting a default label if necessary."""
441 super(BaseStep
, me
).__init__(*arg
, **kw
)
445 Write the proof step to an output FILE.
447 Subclasses must implement a method `_out' which actually does the work.
448 Here, we write a `check' line to verify that the proof actually applies
449 to the number we wanted, if the label says that this is an interesting
453 if me
.label
is not None and not me
.label
.startswith('$'):
454 file.write('check %s, %d, %d\n' %
(me
.label
, me
.p
.nbits
, me
.p
))
456 class SmallStep (BaseStep
):
458 I represent a claim that a number is a small prime.
460 Such claims act as the base cases in a complicated primality proof. When
461 verifying, the claim is checked by trial division using a collection of
464 def __init__(me
, pp
, p
, *arg
, **kw
):
466 Initialize a small-prime step.
468 PP is the overall PrimeProof object of which this is a step; P is the
469 small number whose primality is asserted.
471 super(SmallStep
, me
).__init__(*arg
, **kw
)
474 """Check that the number is indeed a small prime."""
475 return small_test(me
.p
)
477 """Write a small-prime step to the FILE."""
478 file.write('small %s = %d\n' %
(me
.label
, me
.p
))
479 def __repr__(me
): return 'SmallStep(%d)' %
(me
.p
)
481 def parse(cls
, pp
, line
):
483 Parse a small-prime step from a LINE in a proof file.
485 SMALL-STEP ::= `small' LABEL `=' P
487 PP is a PrimeProof object holding the results from the previous steps.
489 if SIEVE
is None: raise ExpectedError('missing `sievebits\' line')
490 label
, p
= parse_label(line
)
491 return cls(pp
, conv_int(p
), label
= label
)
493 class PockStep (BaseStep
):
495 I represent a Pocklington certificate for a number.
497 The number is not explicitly represented in a proof file. See `pock_test'
498 for the underlying mathematics.
500 def __init__(me
, pp
, a
, R
, qqi
, *arg
, **kw
):
502 Inititialize a Pocklington step.
504 PP is the overall PrimeProof object of which this is a step; A is the
505 generator of a substantial subgroup of units; R is a cofactor; and QQI is
506 a sequence of labels for previous proof steps. If Q is the product of
507 the primes listed in QQI, then the number whose primality is asserted is
510 super(PockStep
, me
).__init__(*arg
, **kw
)
514 me
._qq
= [pp
.get_step(qi
).p
for qi
in qqi
]
515 me
.p
= prod(me
._qq
, 2*R
) + 1
517 """Verify a proof step based on Pocklington's theorem."""
518 return pock_test(me
.p
, me
._a
, me
._qq
)
520 """Write a Pocklington step to the FILE."""
521 file.write('pock %s = %d, %d, [%s]\n' % \
523 me
._R
, ', '.join('%s' % qi
for qi
in me
._qqi
)))
524 def __repr__(me
): return 'PockStep(%d, %d, %s)' %
(me
._a
, me
._R
, me
._qqi
)
526 def parse(cls
, pp
, line
):
528 Parse a Pocklington step from a LINE in a proof file.
530 POCK-STEP ::= `pock' LABEL `=' A `,' R `,' `[' Q-LIST `]'
531 Q-LIST ::= Q [`,' Q-LIST]
533 PP is a PrimeProof object holding the results from the previous steps.
535 label
, rest
= parse_label(line
)
536 a
, R
, qq
= parse_list(rest
, 3)
538 if not qq
.startswith('[') or not qq
.endswith(']'):
539 raise ExpectedError('missing `[...]\' around Pocklington factors')
540 return cls(pp
, conv_int(a
), conv_int(R
),
541 [q
.strip() for q
in qq
[1:-1].split(',')], label
= label
)
543 class ECPPStep (BaseStep
):
545 I represent a Goldwasser--Kilian ECPP certificate for a number.
547 def __init__(me
, pp
, p
, a
, b
, x
, y
, qqi
, *arg
, **kw
):
549 Inititialize an ECPP step.
551 PP is the overall PrimeProof object of which this is a step; P is the
552 number whose primality is asserted; A and B are the short Weierstraß
553 curve coefficients; X and Y are the base point coordinates; and QQI is a
554 sequence of labels for previous proof steps.
556 super(ECPPStep
, me
).__init__(*arg
, **kw
)
560 me
._qq
= [pp
.get_step(qi
).p
for qi
in qqi
]
563 """Verify a proof step based on Goldwasser and Kilian's theorem."""
564 return ecpp_test(me
.p
, me
._a
, me
._b
, me
._x
, me
._y
, me
._qq
)
566 """Write an ECPP step to the FILE."""
567 file.write('ecpp %s = %d, %d, %d, %d, %d, [%s]\n' % \
568 (me
.label
, me
.p
, me
._a
, me
._b
, me
._x
, me
._y
,
569 ', '.join('%s' % qi
for qi
in me
._qqi
)))
571 return 'ECPPstep(%d, %d, %d, %d, %d, %s)' % \
572 (me
.p
, me
._a
, me
._b
, me
._x
, me
._y
, me
._qqi
)
574 def parse(cls
, pp
, line
):
576 Parse an ECPP step from a LINE in a proof file.
578 ECPP-STEP ::= `ecpp' LABEL `=' P `,' A `,' B `,' X `,' Y `,'
580 Q-LIST ::= Q [`,' Q-LIST]
582 PP is a PrimeProof object holding the results from the previous steps.
584 label
, rest
= parse_label(line
)
585 p
, a
, b
, x
, y
, qq
= parse_list(rest
, 6)
587 if not qq
.startswith('[') or not qq
.endswith(']'):
588 raise ExpectedError('missing `[...]\' around ECPP factors')
589 return cls(pp
, conv_int(p
), conv_int(a
), conv_int(b
),
590 conv_int(x
), conv_int(y
),
591 [q
.strip() for q
in qq
[1:-1].split(',')], label
= label
)
595 Handle a `check' line in a proof file.
597 CHECK ::= `check' LABEL, B, N
599 Verify that the proof step with the given LABEL asserts the primality of
600 the integer N, and that 2^{B-1} <= N < 2^B.
602 label
, nb
, p
= parse_list(line
, 3)
603 label
, nb
, p
= label
.strip(), conv_int(nb
), conv_int(p
)
604 pi
= pp
.get_step(label
).p
606 raise ExpectedError('check failed: %s = %d /= %d' %
(label
, pi
, p
))
608 raise ExpectedError('check failed: nbits(%s) = %d /= %d' % \
609 (label
, p
.nbits
, nb
))
610 if VERBOSITY
: print ';; %s = %d [%d]' %
(label
, p
, nb
)
612 def setsievebits(pp
, line
):
614 Handle a `sievebits' line in a proof file.
616 SIEVEBITS ::= `sievebits' N
618 Ensure that the verifier is willing to accept small primes up to 2^N.
622 class PrimeProof (object):
624 I represent a proof of primality for one or more numbers.
626 I can encode my proof as a line-oriented text file, in a simple format, and
627 read such a proof back to check it.
630 ## A table to dispatch on keywords read from a file.
631 STEPMAP
= { 'small': SmallStep
.parse
,
632 'pock': PockStep
.parse
,
633 'ecpp': ECPPStep
.parse
,
634 'sievebits': setsievebits
,
639 Initialize a proof object.
641 me
._steps
= {} # Maps labels to steps.
642 me
._stepseq
= [] # Sequence of labels, in order.
643 me
._pmap
= {} # Maps primes to steps.
646 def addstep(me
, step
):
648 Add a new STEP to the proof.
650 The STEP may have a label already. If not, a new internal label is
651 chosen. The proof step is checked before being added to the proof. The
655 ## If there's already a step for this prime, and the new step doesn't
656 ## have a label, then return the old one instead.
657 if step
.label
is None:
658 try: return me
._pmap
[step
.p
]
659 except KeyError: pass
661 ## Make sure the step is actually correct.
664 ## Generate a label if the step doesn't have one already.
665 if step
.label
is None: step
.label
= '$t%d' % me
._i
; me
._i
+= 1
667 ## If the label is already taken then we have a problem.
668 if step
.label
in me
._steps
:
669 raise ExpectedError('duplicate label `%s\'' % step
.label
)
671 ## Store the proof step.
672 me
._pmap
[step
.p
] = step
.label
673 me
._steps
[step
.label
] = step
674 me
._stepseq
.append(step
.label
)
677 def get_step(me
, label
):
679 Check that LABEL labels a known step, and return that step.
681 try: return me
._steps
[label
]
682 except KeyError: raise ExpectedError('unknown label `%s\'' % label
)
686 Write the proof to the given FILE.
689 ## Prefix the main steps with a `sievebits' line.
690 file.write('sievebits %d\n' %
(2*(SIEVE
.limit
.bit_length() - 1)))
692 ## Write the steps out one by one.
693 for label
in me
._stepseq
: me
._steps
[label
].out(file)
697 Read a proof from a given FILE.
699 FILE ::= {STEP | CHECK | SIEVEBITS} [FILE]
700 STEP ::= SMALL-STEP | POCK-STEP
702 Comments (beginning `;') and blank lines are ignored. Other lines begin
706 for lno
, line
in enumerate(file, 1):
708 if line
.startswith(';'): continue
709 ww
= line
.split(None, 1)
712 if len(ww
) > 1: tail
= ww
[1]
715 try: op
= me
.STEPMAP
[w
]
717 raise ExpectedError('unrecognized keyword `%s\'' % w
)
722 except ExpectedError
:
723 raise ExpectedError('%s:%d: %s' %
724 (file.name
, lno
, _excval().message
))
727 ###--------------------------------------------------------------------------
728 ### Finding provable primes.
730 class BasePrime (object):
732 I represent a prime number which has been found and can be proven.
734 This object can eventually be turned into a sequence of proof steps and
735 added to a PrimeProof. This isn't done immediately, because some
736 prime-search strategies want to build a pool of provable primes and will
737 then select some subset of them to actually construct the number of final
738 interest. This way, we avoid cluttering the output proof with proofs of
739 uninteresting numbers.
743 p The prime number in question.
745 label(LABEL) Associate LABEL with this prime, and the corresponding proof
746 step. A label can be set in the constructor, or later using
749 register(PP) Register the prime with a PrimeProof, adding any necessary
750 proof steps. Returns the label of the proof step for this
754 Return a proof step for this prime.
756 def __init__(me
, label
= None, *args
, **kw
):
757 """Initialize a provable prime number object."""
758 super(BasePrime
, me
).__init__(*args
, **kw
)
759 me
._index
= me
._pp
= None
761 def label(me
, label
):
762 """Set this number's LABEL."""
764 def register(me
, pp
):
766 Register the prime's proof steps with PrimeProof PP.
768 Return the final step's label.
770 if me
._pp
is not None:
774 me
._index
= pp
.addstep(me
._mkstep(pp
, label
= me
._label
))
775 ##try: me._index = pp.addstep(me._mkstep(pp, label = me._label))
776 ##except: raise RuntimeError('generated proof failed sanity check')
779 class SmallPrime (BasePrime
):
780 """I represent a prime small enough to be checked in isolation."""
781 def __init__(me
, p
, *args
, **kw
):
782 super(SmallPrime
, me
).__init__(*args
, **kw
)
784 def _mkstep(me
, pp
, **kw
):
785 return SmallStep(pp
, me
.p
, **kw
)
787 class PockPrime (BasePrime
):
788 """I represent a prime proven using Pocklington's theorem."""
789 def __init__(me
, p
, a
, qq
, *args
, **kw
):
790 super(PockPrime
, me
).__init__(*args
, **kw
)
794 def _mkstep(me
, pp
, **kw
):
795 return PockStep(pp
, me
._a
, (me
.p
- 1)/prod((q
.p
for q
in me
._qq
), 2),
796 [q
.register(pp
) for q
in me
._qq
], **kw
)
798 def gen_small(nbits
, label
= None, p
= None):
800 Return a new small prime.
802 The prime will be exactly NBITS bits long. The proof step will have the
803 given LABEL attached. Report progress to the ProgressReporter P.
807 ## Pick a random NBITS-bit number.
808 n
= C
.rand
.mp(nbits
, 1)
809 assert n
.nbits
== nbits
811 ## If it's probably prime, then check it against the small primes we
812 ## know. If it passes then we're done. Otherwise, try again.
814 for q
in SIEVE
.smallprimes():
815 if q
*q
> n
: return SmallPrime(n
, label
= label
)
818 def gen_pock(nbits
, nsubbits
= 0, label
= None, p
= ProgressReporter()):
820 Return a new prime provable using Pocklington's theorem.
822 The prime N will be exactly NBITS long, of the form N = 2 Q R + 1. If
823 NSUBBITS is nonzero, then each prime factor of Q will be NSUBBITS bits
824 long; otherwise a suitable default will be chosen. The proof step will
825 have the given LABEL attached. Report progress to the ProgressReporter P.
827 The prime numbers this function returns are a long way from being uniformly
831 ## Pick a suitable value for NSUBBITS if we don't have one.
834 ## This is remarkably tricky. Picking about 1/3 sqrt(NBITS) factors
835 ## seems about right for large numbers, but there's serious trouble
836 ## lurking for small sizes.
837 nsubbits
= int(3*M
.sqrt(nbits
))
838 if nbits
< nsubbits
+ 3: nsubbits
= nbits
//2 + 1
839 if nbits
== 2*nsubbits
: nsubbits
+= 1
841 ## Figure out how many subgroups we'll need.
842 npiece
= ((nbits
+ 1)//2 + nsubbits
- 1)//nsubbits
848 ## Come up with a collection of known prime factors.
849 p
.p('!'); qq
= [gen(nsubbits
, p
= p
) for i
in xrange(npiece
)]
850 Q
= prod(q
.p
for q
in qq
)
852 ## Come up with bounds on the cofactor. If we're to have N = 2 Q R + 1,
853 ## and 2^{B-1} <= N < 2^B, then we must have 2^{B-2}/Q <= R < 2^{B-1}/Q.
854 Rbase
= (C
.MP(0).setbit(nbits
- 2) + Q
- 1)//Q
855 Rwd
= C
.MP(0).setbit(nbits
- 2)//Q
857 ## Probe the available space of cofactors. If the space is kind of
858 ## narrow, then we want to give up quickly if we're not finding anything
864 ## Pick a random cofactor and examine the number we ended up with.
865 ## Make sure it really does have the length we expect.
866 R
= C
.rand
.range(Rwd
) + Rbase
868 assert n
.nbits
== nbits
870 ## As a complication, if NPIECE is 1, it's just about possible that Q^2
871 ## <= n, in which case this isn't going to work.
874 ## If n has small factors, then pick another cofactor.
875 if C
.PrimeFilter
.smallfactor(n
) == C
.PGEN_FAIL
: continue
877 ## Work through the small primes to find a suitable generator. The
878 ## value 2 is almost always acceptable, so don't try too hard here.
879 for a
in I
.islice(SIEVE
.smallprimes(), 16):
881 ## First, try the Fermat test. If that fails, then n is definitely
883 if pow(a
, n
- 1, n
) != 1: p
.p('.'); break
886 ## Work through the subgroup orders, checking that suitable powers of
887 ## a generate the necessary subgroups.
889 if n
.gcd(pow(a
, (n
- 1)/q
.p
, n
) - 1) != 1:
890 p
.p('@'); ok
= False; break
895 if ok
: p
.pop(); return PockPrime(n
, a
, qq
, label
= label
)
897 def gen(nbits
, label
= None, p
= ProgressReporter()):
899 Generate a prime number with NBITS bits.
901 Give it the LABEL, and report progress to P.
903 if SIEVE
.limit
>> (nbits
+ 1)/2: g
= gen_small
905 return g(nbits
, label
= label
, p
= p
)
907 def gen_limlee(nbits
, nsubbits
,
908 label
= None, qlfmt
= None, p
= ProgressReporter()):
910 Generate a Lim--Lee prime with NBITS bits.
912 Let p be the prime. Then we'll have p = 2 q_0 q_1 ... q_k, with all q_i at
913 least NSUBBITS bits long, and all but q_k exactly that long.
915 The prime will be given the LABEL; progress is reported to P. The factors
916 q_i will be labelled by filling in the `printf'-style format string QLFMT
920 ## Figure out how many factors (p - 1)/2 will have.
921 npiece
= nbits
//nsubbits
922 if npiece
< 2: raise ExpectedError('too few pieces')
924 ## Decide how big to make the pool of factors.
925 poolsz
= max(3*npiece
+ 5, 25) # Heuristic from GnuPG
927 ## Prepare for the main loop.
932 ## Try to make a prime.
936 ## Construct a pool of NSUBBITS-size primes. There's a problem with very
937 ## small sizes: we might not be able to build a pool of distinct primes.
939 for i
in xrange(poolsz
):
941 q
= gen(nsubbits
, p
= p
)
942 if q
.p
not in qmap
: break
944 raise ExpectedError('insufficient diversity')
948 ## Work through combinations of factors from the pool.
949 for qq
in combinations(npiece
- 1, pool
):
951 ## Construct the product of the selected factors.
952 qsmall
= prod(q
.p
for q
in qq
)
954 ## Maybe we'll need to replace the large factor. Try not to do this
955 ## too often. DISP measures the large factor's performance at
956 ## producing candidates with the right length. If it looks bad then
957 ## we'll have to replace it.
958 if 3*disp
*disp
> nstep
*nstep
:
960 if disp
< 0: p
.p('<')
963 ## If we don't have a large factor, then make one.
965 qbig
= gen(nbits
- qsmall
.nbits
, p
= p
)
968 ## We have a candidate. Calculate it and make sure it has the right
970 n
= 2*qsmall
*qbig
.p
+ 1
972 if n
.nbits
< nbits
: disp
-= 1
973 elif n
.nbits
> nbits
: disp
+= 1
974 elif C
.PrimeFilter
.smallfactor(n
) == C
.PGEN_FAIL
: pass
977 ## The candidate has passed the small-primes test. Now check it
978 ## against Pocklington.
979 for a
in I
.islice(SIEVE
.smallprimes(), 16):
982 if pow(a
, n
- 1, n
) != 1: p
.p('.'); break
985 ## Find a generator of a sufficiently large subgroup.
986 if n
.gcd(pow(a
, (n
- 1)/qbig
.p
, n
) - 1) != 1: p
.p('@'); continue
989 if n
.gcd(pow(a
, (n
- 1)/q
.p
, n
) - 1) != 1:
990 p
.p('@'); ok
= False; break
995 ## Label the factors.
998 for i
, q
in enumerate(qq
): q
.label(qlfmt % i
)
1000 ## Return the number we found.
1001 p
.pop(); return PockPrime(n
, a
, qq
, label
= label
)
1003 ###--------------------------------------------------------------------------
1009 ## Prepare an option parser.
1010 op
= OP
.OptionParser(
1012 pock [-qv] [-s SIEVEBITS] CMD ARGS...
1016 description
= 'Generate or verify certified prime numbers.')
1017 op
.add_option('-v', '--verbose', dest
= 'verbosity',
1018 action
= 'count', default
= 1,
1019 help = 'print mysterious runes while looking for prime numbers')
1020 op
.add_option('-q', '--quiet', dest
= 'quietude',
1021 action
= 'count', default
= 0,
1022 help = 'be quiet while looking for prime numbers')
1023 op
.add_option('-s', '--sievebits', dest
= 'sievebits',
1024 type = 'int', default
= 32,
1025 help = 'size (in bits) of largest small prime')
1026 opts
, argv
= op
.parse_args()
1027 VERBOSITY
= opts
.verbosity
- opts
.quietude
1028 p
= ProgressReporter()
1029 a
= ArgFetcher(argv
, op
.error
)
1031 ## Process arguments and do what the user asked.
1035 ## Generate a prime with no special structure.
1036 initsieve(opts
.sievebits
)
1037 nbits
= a
.int(min = 4)
1039 p
= gen(nbits
, 'p', p
= p
)
1044 ## Generate a Lim--Lee prime.
1045 initsieve(opts
.sievebits
)
1046 nbits
= a
.int(min = 4)
1047 nsubbits
= a
.int(min = 4, max = nbits
)
1049 p
= gen_limlee(nbits
, nsubbits
, 'p', 'q_%d', p
= p
)
1054 ## Check an existing certificate.
1055 fn
= a
.arg(default
= '-', must
= False)
1056 if fn
== '-': f
= stdin
1057 else: f
= open(fn
, 'r')
1062 raise ExpectedError("unknown command `%s'" % w
)
1064 if __name__
== '__main__':
1065 prog
= OS
.path
.basename(argv
[0])
1067 except ExpectedError
: exit('%s: %s' %
(prog
, _excval().message
))
1068 except IOError: exit('%s: %s' %
(prog
, _excval()))
1070 ###----- That's all, folks --------------------------------------------------