--- /dev/null
+#! /usr/bin/python
+### -*- mode: python, coding: utf-8 -*-
+###
+### Tool for generating and verifying primality certificates
+###
+### (c) 2017 Straylight/Edgeware
+###
+
+###----- Licensing notice ---------------------------------------------------
+###
+### This file is part of the Python interface to Catacomb.
+###
+### Catacomb/Python is free software; you can redistribute it and/or modify
+### it under the terms of the GNU General Public License as published by
+### the Free Software Foundation; either version 2 of the License, or
+### (at your option) any later version.
+###
+### Catacomb/Python is distributed in the hope that it will be useful,
+### but WITHOUT ANY WARRANTY; without even the implied warranty of
+### MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+### GNU General Public License for more details.
+###
+### You should have received a copy of the GNU General Public License
+### along with Catacomb/Python; if not, write to the Free Software Foundation,
+### Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+
+###--------------------------------------------------------------------------
+### Imported modules.
+
+from sys import argv, stdin, stdout, stderr
+import os as OS
+import itertools as I
+import math as M
+import optparse as OP
+
+import catacomb as C
+
+###--------------------------------------------------------------------------
+### Utilities.
+
+class ExpectedError (Exception):
+ """
+ I represent an expected error, which should be reported in a friendly way.
+ """
+ pass
+
+def prod(ff, one = 1):
+ """
+ Return ONE times the product of the elements of FF.
+
+ This is not done very efficiently.
+ """
+ return reduce(lambda prod, f: prod*f, ff, one)
+
+def parse_label(line):
+ """
+ Split LINE at an `=' character and return the left and right halves.
+
+ The returned pieces have leading and trailing whitespace trimmed.
+ """
+ eq = line.find('=')
+ if eq < 0: raise ExpectedError('expected `LABEL = ...\'')
+ return line[:eq].strip(), line[eq + 1:].strip()
+
+def parse_list(s, n):
+ l = s.split(',', n - 1)
+ if n is not None and len(l) != n:
+ raise ExpectedError('expected `,\'-separated list of %d items' % n)
+ return l
+
+def conv_int(s):
+ """Convert S to a integer."""
+ try: return C.MP(s, 0)
+ except TypeError: raise ExpectedError('invalid integer `%s\'' % s)
+
+VERBOSITY = 1
+
+class ProgressReporter (object):
+ """
+ I keep users amused while the program looks for large prime numbers.
+
+ My main strategy is the printing of incomprehensible runes. I can be
+ muffled by lowering by verbosity level.
+
+ Prime searches are recursive in nature. When a new recursive level is
+ started, call `push'; and call `pop' when the level is finished. This must
+ be done around the top level too.
+ """
+ def __init__(me):
+ """Initialize the ProgressReporter."""
+ me._level = -1
+ me._update()
+ def _update(me):
+ """
+ Update our idea of whether we're active.
+
+ We don't write inscrutable runes when inactive. The current policy is to
+ write nothing if verbosity is zero, to write runes for the top level only
+ if verbosity is 1, and to write runes always if verbosity is higher than
+ that.
+ """
+ me._active = VERBOSITY >= 2 or (VERBOSITY == 1 and me._level == 0)
+ def push(me):
+ """Push a new search level."""
+ me._level += 1
+ me._update()
+ if me._level > 0: me.p('[')
+ else: me.p(';; ')
+ def pop(me):
+ """Pop a search level."""
+ if me._level > 0: me.p(']')
+ else: me.p('\n')
+ me._level -= 1
+ me._update()
+ def p(me, ch):
+ """Print CH as a progress rune."""
+ if me._active: stderr.write(ch); stderr.flush()
+
+def combinations(r, v):
+ """
+ Return an iterator which yields all combinations of R elements from V.
+
+ V must be an indexable sequence. The each combination is returned as a
+ list, containing elements from V in their original order.
+ """
+
+ ## Set up the selection vector. C will contain the indices of the items of
+ ## V we've selected for the current combination. At all times, C contains
+ ## a strictly increasing sequence of integers in the interval [0, N).
+ n = len(v)
+ c = range(r)
+
+ while True:
+
+ ## Yield up the current combination.
+ vv = [v[i] for i in c]
+ yield vv
+
+ ## Now advance to the next one. Find the last index in C which we can
+ ## increment subject to the rules. As we iterate downwards, i will
+ ## contain the index into C, and j will be the maximum acceptable value
+ ## for the corresponding item. We'll step the last index until it
+ ## reaches the limit, and then step the next one down, resetting the last
+ ## index, and so on.
+ i, j = r, n
+ while True:
+
+ ## If i is zero here, then we've advanced everything as far as it will
+ ## go. We're done.
+ if i == 0: return
+
+ ## Move down to the next index.
+ i -= 1; j -= 1
+
+ ## If this index isn't at its maximum value, then we've found the place
+ ## to step.
+ if c[i] != j: break
+
+ ## Step this index on by one, and set the following indices to the
+ ## immediately following values.
+ j = c[i] + 1
+ while i < r: c[i] = j; i += 1; j += 1
+
+class ArgFetcher (object):
+ """
+ I return arguments from a list, reporting problems when they occur.
+ """
+ def __init__(me, argv, errfn):
+ """
+ Initialize, returning successive arguments from ARGV.
+
+ Errors are reported to ERRFN.
+ """
+ me._argv = argv
+ me._argc = len(argv)
+ me._errfn = errfn
+ me._i = 0
+ def arg(me, default = None, must = True):
+ """
+ Return the next argument.
+
+ If MUST is true, then report an error (to the ERRFN) if there are no more
+ arguments; otherwise, return the DEFAULT.
+ """
+ if me._i >= me._argc:
+ if must: me._errfn('missing argument')
+ return default
+ arg = me._argv[me._i]; me._i += 1
+ return arg
+ def int(me, default = None, must = True, min = None, max = None):
+ """
+ Return the next argument converted to an integer.
+
+ If MUST is true, then report an error (to the ERRFN) if there are no more
+ arguments; otherwise return the DEFAULT. Report an error if the next
+ argument is not a valid integer, or if the integer is beyond the MIN and
+ MAX bounds.
+ """
+ arg = me.arg(default = None, must = must)
+ if arg is None: return default
+ try: arg = int(arg)
+ except ValueError: me._errfn('bad integer')
+ if (min is not None and arg < min) or (max is not None and arg > max):
+ me._errfn('out of range')
+ return arg
+
+###--------------------------------------------------------------------------
+### Sieving for small primes.
+
+class Sieve (object):
+ """
+ I represent a collection of small primes, up to some chosen limit.
+
+ The limit is available as the `limit' attribute. Let L be this limit;
+ then, if N < L^2 is some composite, then N has at least one prime factor
+ less than L.
+ """
+
+ ## Figure out the number of bits in a (nonnegative) primitive `int'. We'll
+ ## use a list of these as our sieve.
+ _NBIT = 15
+ while type(1 << (_NBIT + 1)) == int: _NBIT += 1
+
+ def __init__(me, limit):
+ """
+ Initialize a sieve holding all primes below LIMIT.
+ """
+
+ ## The sieve is maintained in the `_bits' attribute. This is a list of
+ ## integers, used as a bitmask: let 2 < n < L be an odd integer; then bit
+ ## (n - 3)/2 will be clear iff n is prime. Let W be the value of
+ ## `_NBIT', above; then bit W i + j in the sieve is stored in bit j of
+ ## `_bits[i]'.
+
+ ## Store the limit for later inspection.
+ me.limit = limit
+
+ ## Calculate the size of sieve we'll need and initialize the bit list.
+ n = (limit - 2)/2
+ sievesz = (n + me._NBIT - 1)/me._NBIT
+ me._sievemax = sievesz*me._NBIT
+ me._bits = n*[0]
+
+ ## This is standard Sieve of Eratosthenes. For each index i: if
+ ## bit i is clear, then p = 2 i + 3 is prime, so set the bits
+ ## corresponding to each multiple of p, i.e., bits (k p - 3)/2 =
+ ## (2 k i + 3 - 3)/2 = k i for k > 1.
+ for i in xrange(me._sievemax):
+ if me._bitp(i): i += 1; continue
+ p = 2*i + 3
+ if p >= limit: break
+ for j in xrange(i + p, me._sievemax, p): me._setbit(j)
+ i += 1
+
+ def _bitp(me, i): i, j = divmod(i, me._NBIT); return (me._bits[i] >> j)&1
+ def _setbit(me, i): i, j = divmod(i, me._NBIT); me._bits[i] |= 1 << j
+
+ def smallprimes(me):
+ """
+ Return an iterator over the known small primes.
+ """
+ yield 2
+ n = 3
+ for b in me._bits:
+ for j in xrange(me._NBIT):
+ if not (b&1): yield n
+ b >>= 1; n += 2
+
+## We generate the sieve on demand.
+SIEVE = None
+
+def initsieve(sievebits):
+ """
+ Generate the sieve.
+
+ Ensure that it can be used to check the primality of numbers up to (but not
+ including) 2^SIEVEBITS.
+ """
+ global SIEVE
+ if SIEVE is not None: raise ValueError('sieve already defined')
+ if sievebits < 6: sievebits = 6
+ SIEVE = Sieve(1 << (sievebits + 1)/2)
+
+###--------------------------------------------------------------------------
+### Primality checking.
+
+def small_test(p):
+ """
+ Check that P is a small prime.
+
+ If not, raise an `ExpectedError'. The `SIEVE' variable must have been
+ initialized.
+ """
+ if p < 2: raise ExpectedError('%d too small' % p)
+ if SIEVE.limit*SIEVE.limit < p:
+ raise ExpectedError('%d too large for small prime' % p)
+ for q in SIEVE.smallprimes():
+ if q*q > p: return
+ if p%q == 0: raise ExpectedError('%d divides %d' % (q, p))
+
+def pock_test(p, a, qq):
+ """
+ Check that P is prime using Pocklington's criterion.
+
+ If not, raise an `ExpectedError'.
+
+ Let Q be the product of the elements of the sequence QQ. The test works as
+ follows. Suppose p is the smallest prime factor of P. If A^{P-1} /== 1
+ (mod P) then P is certainly composite (Fermat's test); otherwise, we have
+ established that the order of A in (Z/pZ)^* divides P - 1. Next, let t =
+ A^{(P-1)/q} for some prime factor q of Q, and let g = gcd(t - 1, P). If g
+ = P then the proof is inconclusive; if 1 < g < P then g is a nontrivial
+ factor of P, so P is composite; otherwise, t has order q in (Z/pZ)^*, so
+ (Z/pZ)^* contains a subgroup of size q, and therefore q divides p - 1. If
+ QQ is a sequence of distinct primes, and the preceding criterion holds for
+ all q in QQ, then Q divides p - 1. If Q^2 < P then the proof is
+ inconclusive; otherwise, let p' be any prime dividing P/p. Then p' >= p >
+ Q, so p p' > Q^2 > P; but p p' divides P, so this is a contradiction.
+ Therefore P/p has no prime factors, and P is prime.
+ """
+
+ ## We don't actually need the distinctness criterion. Suppose that q^e
+ ## divides Q. Then gcd(t - 1, P) = 1 implies that A^{(P-1)/q^{e-1}} has
+ ## order q^e in (Z/pZ)^*, which accounts for the multiplicity.
+
+ Q = prod(qq)
+ if p < 2: raise ExpectedError('%d too small' % p)
+ if Q*Q <= p:
+ raise ExpectedError('too few Pocklington factors for %d' % p)
+ if pow(a, p - 1, p) != 1:
+ raise ExpectedError('%d is Fermat witness for %d' % (a, p))
+ for q in qq:
+ if Q%(q*q) == 0:
+ raise ExpectedError('duplicate Pocklington factor %d for %d' % (q, p))
+ g = p.gcd(pow(a, (p - 1)/q, p) - 1)
+ if g == p:
+ raise ExpectedError('%d order not multiple of %d mod %d' % (a, q, p))
+ elif g != 1:
+ raise ExpectedError('%d divides %d' % (g, p))
+
+def ecpp_test(p, a, b, x, y, qq):
+ """
+ Check that P is prime using Goldwasser and Kilian's ECPP method.
+
+ If not, raise an `ExpectedError'.
+
+ Let Q be the product of the elements of the sequence QQ. Suppose p is the
+ smallest prime factor of P. Let g = gcd(4 A^3 + 27 B^2, P). If g = P then
+ the test is inconclusive; otherwise, if g /= 1 then g is a nontrivial
+ factor of P. Define E(GF(p)) = { (x, y) | y^2 = x^3 + A x + B } U { inf }
+ to be the elliptic curve over p with short-Weierstraß coefficients A and B;
+ we have just checked that this curve is not singular. If R = (X, Y) is not
+ a point on this curve, then the test is inconclusive. If Q R is not the
+ point at infinity, then the test fails; otherwise we deduce that P has
+ Q-torsion in E. Let S = (Q/q) R for some prime factor q of Q. If S is the
+ point at infinity then the test is inconclusive; otherwise, q divides the
+ order of S in E. If QQ is a sequence of distinct primes, and the preceding
+ criterion holds for all q in QQ, then Q divides the order of S. Therefore
+ #E(p) >= Q. If Q <= (qrrt(P) + 1)^2 then the test is inconclusive.
+ Otherwise, Hasse's theorem tells us that |p + 1 - #E(p)| <= 2 sqrt(p);
+ hence we must have p + 1 + 2 sqrt(p) = (sqrt(p) + 1)^2 >= #E(p) >= Q >
+ (qrrt(P) + 1)^2; so sqrt(p) + 1 > qrrt(P) + 1, i.e., p^2 > P. As for
+ Pocklington above, if p' is any prime factor of P/p, then p p' >= p^2 > P,
+ which is a contradiction, and we conclude that P is prime.
+ """
+
+ ## This isn't going to work if gcd(P, 6) /= 1: we're going to use the
+ ## large-characteristic addition formulae.
+ g = p.gcd(6)
+ if g != 1: raise ExpectedError('%d divides %d' % (g, p))
+
+ ## We want to check that Q > (qrrt(P) + 1)^2 iff sqrt(Q) > qrrt(P) + 1; but
+ ## calculating square roots is not enjoyable (partly because we have to
+ ## deal with the imprecision). Fortunately, some algebra will help: the
+ ## condition holds iff qrrt(P) < sqrt(Q) - 1 iff P < Q^2 - 4 Q sqrt(Q) +
+ ## 6 Q - 4 sqrt(Q) + 1 = Q (Q + 6) + 1 - 4 sqrt(Q) (Q + 1) iff Q (Q + 6) -
+ ## P + 1 > 4 sqrt(Q) (Q + 1) iff (Q (Q + 6) - P + 1)^2 > 16 Q (Q + 1)^2
+ Q = prod(qq)
+ t, u = Q*(Q + 6) - p + 1, 4*(Q + 1)
+ if t*t <= Q*u*u: raise ExpectedError('too few subgroups for ECPP')
+
+ ## Construct the curve.
+ E = C.PrimeField(p).ec(a, b) # careful: may not be a prime!
+
+ ## Find the base point.
+ R = E(x, y)
+ if not R.oncurvep():
+ raise ExpectedError('(%d, %d) is not on the curve' % (x, y))
+
+ ## Check that it has Q-torsion.
+ if Q*R: raise ExpectedError('(%d, %d) not a %d-torsion point' % (x, y, Q))
+
+ ## Now check the individual factors.
+ for q in qq:
+ if Q%(q*q) == 0:
+ raise ExpectedError('duplicate ECPP factor %d for %d' % (q, p))
+ S = (Q/q)*R
+ if not S:
+ raise ExpectedError('(%d, %d) order not a multiple of %d' % (x, y, q))
+ g = p.gcd(S._z)
+ if g != 1:
+ raise ExpectedError('%d divides %d' % (g, p))
+
+###--------------------------------------------------------------------------
+### Proof steps and proofs.
+
+class BaseStep (object):
+ """
+ I'm a step in a primality proof.
+
+ I assert that a particular number is prime, and can check this.
+
+ This class provides basic protocol for proof steps, mostly to do with
+ handling labels.
+
+ The step's label is kept in its `label' attribute. It can be set by the
+ constructor, and is `None' by default. Users can modify this attribute if
+ they like. Labels beginning `$' are assumed to be internal and
+ uninteresting; other labels cause `check' lines to be written to the output
+ listing the actual number of interest.
+
+ Protocol that proof steps should provide:
+
+ label A string labelling the proof step and the associated prime
+ number.
+
+ p The prime number which this step proves to be prime.
+
+ check() Check that the proof step is actually correct, assuming that
+ any previous steps have already been verified.
+
+ out(FILE) Write an appropriate encoding of the proof step to the output
+ FILE.
+ """
+ def __init__(me, label = None, *arg, **kw):
+ """Initialize a proof step, setting a default label if necessary."""
+ super(BaseStep, me).__init__(*arg, **kw)
+ me.label = label
+ def out(me, file):
+ """
+ Write the proof step to an output FILE.
+
+ Subclasses must implement a method `_out' which actually does the work.
+ Here, we write a `check' line to verify that the proof actually applies
+ to the number we wanted, if the label says that this is an interesting
+ step.
+ """
+ me._out(file)
+ if me.label is not None and not me.label.startswith('$'):
+ file.write('check %s, %d, %d\n' % (me.label, me.p.nbits, me.p))
+
+class SmallStep (BaseStep):
+ """
+ I represent a claim that a number is a small prime.
+
+ Such claims act as the base cases in a complicated primality proof. When
+ verifying, the claim is checked by trial division using a collection of
+ known small primes.
+ """
+ def __init__(me, pp, p, *arg, **kw):
+ """
+ Initialize a small-prime step.
+
+ PP is the overall PrimeProof object of which this is a step; P is the
+ small number whose primality is asserted.
+ """
+ super(SmallStep, me).__init__(*arg, **kw)
+ me.p = p
+ def check(me):
+ """Check that the number is indeed a small prime."""
+ return small_test(me.p)
+ def _out(me, file):
+ """Write a small-prime step to the FILE."""
+ file.write('small %s = %d\n' % (me.label, me.p))
+ def __repr__(me): return 'SmallStep(%d)' % (me.p)
+ @classmethod
+ def parse(cls, pp, line):
+ """
+ Parse a small-prime step from a LINE in a proof file.
+
+ SMALL-STEP ::= `small' LABEL `=' P
+
+ PP is a PrimeProof object holding the results from the previous steps.
+ """
+ if SIEVE is None: raise ExpectedError('missing `sievebits\' line')
+ label, p = parse_label(line)
+ return cls(pp, conv_int(p), label = label)
+
+class PockStep (BaseStep):
+ """
+ I represent a Pocklington certificate for a number.
+
+ The number is not explicitly represented in a proof file. See `pock_test'
+ for the underlying mathematics.
+ """
+ def __init__(me, pp, a, R, qqi, *arg, **kw):
+ """
+ Inititialize a Pocklington step.
+
+ PP is the overall PrimeProof object of which this is a step; A is the
+ generator of a substantial subgroup of units; R is a cofactor; and QQI is
+ a sequence of labels for previous proof steps. If Q is the product of
+ the primes listed in QQI, then the number whose primality is asserted is
+ 2 Q R + 1.
+ """
+ super(PockStep, me).__init__(*arg, **kw)
+ me._a = a
+ me._R = R
+ me._qqi = qqi
+ me._qq = [pp.get_step(qi).p for qi in qqi]
+ me.p = prod(me._qq, 2*R) + 1
+ def check(me):
+ """Verify a proof step based on Pocklington's theorem."""
+ return pock_test(me.p, me._a, me._qq)
+ def _out(me, file):
+ """Write a Pocklington step to the FILE."""
+ file.write('pock %s = %d, %d, [%s]\n' % \
+ (me.label, me._a,
+ me._R, ', '.join('%s' % qi for qi in me._qqi)))
+ def __repr__(me): return 'PockStep(%d, %d, %s)' % (me._a, me._R, me._qqi)
+ @classmethod
+ def parse(cls, pp, line):
+ """
+ Parse a Pocklington step from a LINE in a proof file.
+
+ POCK-STEP ::= `pock' LABEL `=' A `,' R `,' `[' Q-LIST `]'
+ Q-LIST ::= Q [`,' Q-LIST]
+
+ PP is a PrimeProof object holding the results from the previous steps.
+ """
+ label, rest = parse_label(line)
+ a, R, qq = parse_list(rest, 3)
+ qq = qq.strip()
+ if not qq.startswith('[') or not qq.endswith(']'):
+ raise ExpectedError('missing `[...]\' around Pocklington factors')
+ return cls(pp, conv_int(a), conv_int(R),
+ [q.strip() for q in qq[1:-1].split(',')], label = label)
+
+class ECPPStep (BaseStep):
+ """
+ I represent a Goldwasser--Kilian ECPP certificate for a number.
+ """
+ def __init__(me, pp, p, a, b, x, y, qqi, *arg, **kw):
+ """
+ Inititialize an ECPP step.
+
+ PP is the overall PrimeProof object of which this is a step; P is the
+ number whose primality is asserted; A and B are the short Weierstraß
+ curve coefficients; X and Y are the base point coordinates; and QQI is a
+ sequence of labels for previous proof steps.
+ """
+ super(ECPPStep, me).__init__(*arg, **kw)
+ me._a, me._b = a, b
+ me._x, me._y = x, y
+ me._qqi = qqi
+ me._qq = [pp.get_step(qi).p for qi in qqi]
+ me.p = p
+ def check(me):
+ """Verify a proof step based on Goldwasser and Kilian's theorem."""
+ return ecpp_test(me.p, me._a, me._b, me._x, me._y, me._qq)
+ def _out(me, file):
+ """Write an ECPP step to the FILE."""
+ file.write('ecpp %s = %d, %d, %d, %d, %d, [%s]\n' % \
+ (me.label, me.p, me._a, me._b, me._x, me._y,
+ ', '.join('%s' % qi for qi in me._qqi)))
+ def __repr__(me):
+ return 'ECPPstep(%d, %d, %d, %d, %d, %s)' % \
+ (me.p, me._a, me._b, me._x, me._y, me._qqi)
+ @classmethod
+ def parse(cls, pp, line):
+ """
+ Parse an ECPP step from a LINE in a proof file.
+
+ ECPP-STEP ::= `ecpp' LABEL `=' P `,' A `,' B `,' X `,' Y `,'
+ `[' Q-LIST `]'
+ Q-LIST ::= Q [`,' Q-LIST]
+
+ PP is a PrimeProof object holding the results from the previous steps.
+ """
+ label, rest = parse_label(line)
+ p, a, b, x, y, qq = parse_list(rest, 6)
+ qq = qq.strip()
+ if not qq.startswith('[') or not qq.endswith(']'):
+ raise ExpectedError('missing `[...]\' around ECPP factors')
+ return cls(pp, conv_int(p), conv_int(a), conv_int(b),
+ conv_int(x), conv_int(y),
+ [q.strip() for q in qq[1:-1].split(',')], label = label)
+
+def check(pp, line):
+ """
+ Handle a `check' line in a proof file.
+
+ CHECK ::= `check' LABEL, B, N
+
+ Verify that the proof step with the given LABEL asserts the primality of
+ the integer N, and that 2^{B-1} <= N < 2^B.
+ """
+ label, nb, p = parse_list(line, 3)
+ label, nb, p = label.strip(), conv_int(nb), conv_int(p)
+ pi = pp.get_step(label).p
+ if pi != p:
+ raise ExpectedError('check failed: %s = %d /= %d' % (label, pi, p))
+ if p.nbits != nb:
+ raise ExpectedError('check failed: nbits(%s) = %d /= %d' % \
+ (label, p.nbits, nb))
+ if VERBOSITY: print ';; %s = %d [%d]' % (label, p, nb)
+
+def setsievebits(pp, line):
+ """
+ Handle a `sievebits' line in a proof file.
+
+ SIEVEBITS ::= `sievebits' N
+
+ Ensure that the verifier is willing to accept small primes up to 2^N.
+ """
+ initsieve(int(line))
+
+class PrimeProof (object):
+ """
+ I represent a proof of primality for one or more numbers.
+
+ I can encode my proof as a line-oriented text file, in a simple format, and
+ read such a proof back to check it.
+ """
+
+ ## A table to dispatch on keywords read from a file.
+ STEPMAP = { 'small': SmallStep.parse,
+ 'pock': PockStep.parse,
+ 'ecpp': ECPPStep.parse,
+ 'sievebits': setsievebits,
+ 'check': check }
+
+ def __init__(me):
+ """
+ Initialize a proof object.
+ """
+ me._steps = {} # Maps labels to steps.
+ me._stepseq = [] # Sequence of labels, in order.
+ me._pmap = {} # Maps primes to steps.
+ me._i = 0
+
+ def addstep(me, step):
+ """
+ Add a new STEP to the proof.
+
+ The STEP may have a label already. If not, a new internal label is
+ chosen. The proof step is checked before being added to the proof. The
+ label is returned.
+ """
+
+ ## If there's already a step for this prime, and the new step doesn't
+ ## have a label, then return the old one instead.
+ if step.label is None:
+ try: return me._pmap[step.p]
+ except KeyError: pass
+
+ ## Make sure the step is actually correct.
+ step.check()
+
+ ## Generate a label if the step doesn't have one already.
+ if step.label is None: step.label = '$t%d' % me._i; me._i += 1
+
+ ## If the label is already taken then we have a problem.
+ if step.label in me._steps:
+ raise ValueError('duplicate label `%s\'' % step.label)
+
+ ## Store the proof step.
+ me._pmap[step.p] = step.label
+ me._steps[step.label] = step
+ me._stepseq.append(step.label)
+ return step.label
+
+ def get_step(me, label):
+ """
+ Check that LABEL labels a known step, and return that step.
+ """
+ try: return me._steps[label]
+ except KeyError: raise ExpectedError('unknown label `%s\'' % label)
+
+ def write(me, file):
+ """
+ Write the proof to the given FILE.
+ """
+
+ ## Prefix the main steps with a `sievebits' line.
+ file.write('sievebits %d\n' % (2*(SIEVE.limit.bit_length() - 1)))
+
+ ## Write the steps out one by one.
+ for label in me._stepseq: me._steps[label].out(file)
+
+ def read(me, file):
+ """
+ Read a proof from a given FILE.
+
+ FILE ::= {STEP | CHECK | SIEVEBITS} [FILE]
+ STEP ::= SMALL-STEP | POCK-STEP
+
+ Comments (beginning `;') and blank lines are ignored. Other lines begin
+ with a keyword.
+ """
+ lastp = None
+ for lno, line in enumerate(file, 1):
+ line = line.strip()
+ if line.startswith(';'): continue
+ ww = line.split(None, 1)
+ if not ww: continue
+ w = ww[0]
+ if len(ww) > 1: tail = ww[1]
+ else: tail = ''
+ try:
+ try: op = me.STEPMAP[w]
+ except KeyError:
+ raise ExpectedError('unrecognized keyword `%s\'' % w)
+ step = op(me, tail)
+ if step is not None:
+ me.addstep(step)
+ lastp = step.p
+ except ExpectedError, e:
+ raise ExpectedError('%s:%d: %s' % (file.name, lno, e.message))
+ return lastp
+
+###--------------------------------------------------------------------------
+### Finding provable primes.
+
+class BasePrime (object):
+ """
+ I represent a prime number which has been found and can be proven.
+
+ This object can eventually be turned into a sequence of proof steps and
+ added to a PrimeProof. This isn't done immediately, because some
+ prime-search strategies want to build a pool of provable primes and will
+ then select some subset of them to actually construct the number of final
+ interest. This way, we avoid cluttering the output proof with proofs of
+ uninteresting numbers.
+
+ Protocol required.
+
+ p The prime number in question.
+
+ label(LABEL) Associate LABEL with this prime, and the corresponding proof
+ step. A label can be set in the constructor, or later using
+ this method.
+
+ register(PP) Register the prime with a PrimeProof, adding any necessary
+ proof steps. Returns the label of the proof step for this
+ number.
+
+ _mkstep(PP, **KW)
+ Return a proof step for this prime.
+ """
+ def __init__(me, label = None, *args, **kw):
+ """Initialize a provable prime number object."""
+ super(BasePrime, me).__init__(*args, **kw)
+ me._index = me._pp = None
+ me._label = label
+ def label(me, label):
+ """Set this number's LABEL."""
+ me._label = label
+ def register(me, pp):
+ """
+ Register the prime's proof steps with PrimeProof PP.
+
+ Return the final step's label.
+ """
+ if me._pp is not None:
+ assert me._pp == pp
+ else:
+ me._pp = pp
+ me._index = pp.addstep(me._mkstep(pp, label = me._label))
+ ##try: me._index = pp.addstep(me._mkstep(pp, label = me._label))
+ ##except: raise RuntimeError('generated proof failed sanity check')
+ return me._index
+
+class SmallPrime (BasePrime):
+ """I represent a prime small enough to be checked in isolation."""
+ def __init__(me, p, *args, **kw):
+ super(SmallPrime, me).__init__(*args, **kw)
+ me.p = p
+ def _mkstep(me, pp, **kw):
+ return SmallStep(pp, me.p, **kw)
+
+class PockPrime (BasePrime):
+ """I represent a prime proven using Pocklington's theorem."""
+ def __init__(me, p, a, qq, *args, **kw):
+ super(PockPrime, me).__init__(*args, **kw)
+ me.p = p
+ me._a = a
+ me._qq = qq
+ def _mkstep(me, pp, **kw):
+ return PockStep(pp, me._a, (me.p - 1)/prod((q.p for q in me._qq), 2),
+ [q.register(pp) for q in me._qq], **kw)
+
+def gen_small(nbits, label = None, p = None):
+ """
+ Return a new small prime.
+
+ The prime will be exactly NBITS bits long. The proof step will have the
+ given LABEL attached. Report progress to the ProgressReporter P.
+ """
+ while True:
+
+ ## Pick a random NBITS-bit number.
+ n = C.rand.mp(nbits, 1)
+ assert n.nbits == nbits
+
+ ## If it's probably prime, then check it against the small primes we
+ ## know. If it passes then we're done. Otherwise, try again.
+ if n.primep():
+ for q in SIEVE.smallprimes():
+ if q*q > n: return SmallPrime(n, label = label)
+ if n%q == 0: break
+
+def gen_pock(nbits, nsubbits = 0, label = None, p = ProgressReporter()):
+ """
+ Return a new prime provable using Pocklington's theorem.
+
+ The prime N will be exactly NBITS long, of the form N = 2 Q R + 1. If
+ NSUBBITS is nonzero, then each prime factor of Q will be NSUBBITS bits
+ long; otherwise a suitable default will be chosen. The proof step will
+ have the given LABEL attached. Report progress to the ProgressReporter P.
+
+ The prime numbers this function returns are a long way from being uniformly
+ distributed.
+ """
+
+ ## Pick a suitable value for NSUBBITS if we don't have one.
+ if not nsubbits:
+
+ ## This is remarkably tricky. Picking about 1/3 sqrt(NBITS) factors
+ ## seems about right for large numbers, but there's serious trouble
+ ## lurking for small sizes.
+ nsubbits = int(3*M.sqrt(nbits))
+ if nbits < nsubbits + 3: nsubbits = nbits//2 + 1
+ if nbits == 2*nsubbits: nsubbits += 1
+
+ ## Figure out how many subgroups we'll need.
+ npiece = ((nbits + 1)//2 + nsubbits - 1)//nsubbits
+ p.push()
+
+ ## Keep searching...
+ while True:
+
+ ## Come up with a collection of known prime factors.
+ p.p('!'); qq = [gen(nsubbits, p = p) for i in xrange(npiece)]
+ Q = prod(q.p for q in qq)
+
+ ## Come up with bounds on the cofactor. If we're to have N = 2 Q R + 1,
+ ## and 2^{B-1} <= N < 2^B, then we must have 2^{B-2}/Q <= R < 2^{B-1}/Q.
+ Rbase = (C.MP(0).setbit(nbits - 2) + Q - 1)//Q
+ Rwd = C.MP(0).setbit(nbits - 2)//Q
+
+ ## Probe the available space of cofactors. If the space is kind of
+ ## narrow, then we want to give up quickly if we're not finding anything
+ ## suitable.
+ step = 0
+ while step < Rwd:
+ step += 1
+
+ ## Pick a random cofactor and examine the number we ended up with.
+ ## Make sure it really does have the length we expect.
+ R = C.rand.range(Rwd) + Rbase
+ n = 2*Q*R + 1
+ assert n.nbits == nbits
+
+ ## As a complication, if NPIECE is 1, it's just about possible that Q^2
+ ## <= n, in which case this isn't going to work.
+ if Q*Q < n: continue
+
+ ## If n has small factors, then pick another cofactor.
+ if C.PrimeFilter.smallfactor(n) == C.PGEN_FAIL: continue
+
+ ## Work through the small primes to find a suitable generator. The
+ ## value 2 is almost always acceptable, so don't try too hard here.
+ for a in I.islice(SIEVE.smallprimes(), 16):
+
+ ## First, try the Fermat test. If that fails, then n is definitely
+ ## composite.
+ if pow(a, n - 1, n) != 1: p.p('.'); break
+ p.p('*')
+
+ ## Work through the subgroup orders, checking that suitable powers of
+ ## a generate the necessary subgroups.
+ for q in qq:
+ if n.gcd(pow(a, (n - 1)/q.p, n) - 1) != 1:
+ p.p('@'); ok = False; break
+ else:
+ ok = True
+
+ ## we're all good.
+ if ok: p.pop(); return PockPrime(n, a, qq, label = label)
+
+def gen(nbits, label = None, p = ProgressReporter()):
+ """
+ Generate a prime number with NBITS bits.
+
+ Give it the LABEL, and report progress to P.
+ """
+ if SIEVE.limit >> (nbits + 1)/2: g = gen_small
+ else: g = gen_pock
+ return g(nbits, label = label, p = p)
+
+def gen_limlee(nbits, nsubbits,
+ label = None, qlfmt = None, p = ProgressReporter()):
+ """
+ Generate a Lim--Lee prime with NBITS bits.
+
+ Let p be the prime. Then we'll have p = 2 q_0 q_1 ... q_k, with all q_i at
+ least NSUBBITS bits long, and all but q_k exactly that long.
+
+ The prime will be given the LABEL; progress is reported to P. The factors
+ q_i will be labelled by filling in the `printf'-style format string QLFMT
+ with the argument i.
+ """
+
+ ## Figure out how many factors (p - 1)/2 will have.
+ npiece = nbits//nsubbits
+ if npiece < 2: raise ExpectedError('too few pieces')
+
+ ## Decide how big to make the pool of factors.
+ poolsz = max(3*npiece + 5, 25) # Heuristic from GnuPG
+
+ ## Prepare for the main loop.
+ disp = nstep = 0
+ qbig = None
+ p.push()
+
+ ## Try to make a prime.
+ while True:
+ p.p('!')
+
+ ## Construct a pool of NSUBBITS-size primes. There's a problem with very
+ ## small sizes: we might not be able to build a pool of distinct primes.
+ pool = []; qmap = {}
+ for i in xrange(poolsz):
+ for j in xrange(64):
+ q = gen(nsubbits, p = p)
+ if q.p not in qmap: break
+ else:
+ raise ExpectedError('insufficient diversity')
+ qmap[q.p] = q
+ pool.append(q)
+
+ ## Work through combinations of factors from the pool.
+ for qq in combinations(npiece - 1, pool):
+
+ ## Construct the product of the selected factors.
+ qsmall = prod(q.p for q in qq)
+
+ ## Maybe we'll need to replace the large factor. Try not to do this
+ ## too often. DISP measures the large factor's performance at
+ ## producing candidates with the right length. If it looks bad then
+ ## we'll have to replace it.
+ if 3*disp*disp > nstep*nstep:
+ qbig = None
+ if disp < 0: p.p('<')
+ else: p.p('>')
+
+ ## If we don't have a large factor, then make one.
+ if qbig is None:
+ qbig = gen(nbits - qsmall.nbits, p = p)
+ disp = 0; nstep = 0
+
+ ## We have a candidate. Calculate it and make sure it has the right
+ ## length.
+ n = 2*qsmall*qbig.p + 1
+ nstep += 1
+ if n.nbits < nbits: disp -= 1
+ elif n.nbits > nbits: disp += 1
+ elif C.PrimeFilter.smallfactor(n) == C.PGEN_FAIL: pass
+ else:
+
+ ## The candidate has passed the small-primes test. Now check it
+ ## against Pocklington.
+ for a in I.islice(SIEVE.smallprimes(), 16):
+
+ ## Fermat test.
+ if pow(a, n - 1, n) != 1: p.p('.'); break
+ p.p('*')
+
+ ## Find a generator of a sufficiently large subgroup.
+ if n.gcd(pow(a, (n - 1)/qbig.p, n) - 1) != 1: p.p('@'); continue
+ ok = True
+ for q in qq:
+ if n.gcd(pow(a, (n - 1)/q.p, n) - 1) != 1:
+ p.p('@'); ok = False; break
+
+ ## We're done.
+ if ok:
+
+ ## Label the factors.
+ qq.append(qbig)
+ if qlfmt:
+ for i, q in enumerate(qq): q.label(qlfmt % i)
+
+ ## Return the number we found.
+ p.pop(); return PockPrime(n, a, qq, label = label)
+
+###--------------------------------------------------------------------------
+### Main program.
+
+def __main__():
+ global VERBOSITY
+
+ ## Prepare an option parser.
+ op = OP.OptionParser(
+ usage = '''\
+pock [-qv] CMD ARGS...
+ gen NBITS
+ ll NBITS NSUBBITS
+ [check] [FILE]''',
+ description = 'Generate or verify certified prime numbers.')
+ op.add_option('-v', '--verbose', dest = 'verbosity',
+ action = 'count', default = 1,
+ help = 'Print mysterious runes while looking for prime numbers.')
+ op.add_option('-q', '--quiet', dest = 'quietude',
+ action = 'count', default = 0,
+ help = 'be quiet while looking for prime numbers.')
+ op.add_option('-s', '--sievebits', dest = 'sievebits',
+ type = 'int', default = 32,
+ help = 'Size (in bits) of largest small prime.')
+ opts, argv = op.parse_args()
+ VERBOSITY = opts.verbosity - opts.quietude
+ p = ProgressReporter()
+ a = ArgFetcher(argv, op.error)
+
+ ## Process arguments and do what the user asked.
+ w = a.arg()
+
+ if w == 'gen':
+ ## Generate a prime with no special structure.
+ initsieve(opts.sievebits)
+ nbits = a.int(min = 4)
+ pp = PrimeProof()
+ p = gen(nbits, 'p', p = p)
+ p.register(pp)
+ pp.write(stdout)
+
+ elif w == 'll':
+ ## Generate a Lim--Lee prime.
+ initsieve(opts.sievebits)
+ nbits = a.int(min = 4)
+ nsubbits = a.int(min = 4, max = nbits)
+ pp = PrimeProof()
+ p = gen_limlee(nbits, nsubbits, 'p', 'q_%d', p = p)
+ p.register(pp)
+ pp.write(stdout)
+
+ elif w == 'check':
+ ## Check an existing certificate.
+ fn = a.arg(default = '-', must = False)
+ if fn == '-': f = stdin
+ else: f = open(fn, 'r')
+ pp = PrimeProof()
+ p = pp.read(f)
+
+ else:
+ raise ExpectedError("unknown command `%s'" % w)
+
+if __name__ == '__main__':
+ prog = OS.path.basename(argv[0])
+ try: __main__()
+ except ExpectedError, e: exit('%s: %s' % (prog, e.message))
+ except IOError, e: exit('%s: %s' % (prog, e))
+
+###----- That's all, folks --------------------------------------------------
~mdw / catacomb-python / blobdiff