algorithms.c: Add bindings for STROBE.

[catacomb-python] / pock

#! /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.

import sys as SYS; 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

if SYS.version_info >= (3,):
  xrange = range
  range = lambda *args: list(xrange(*args))

###--------------------------------------------------------------------------
### Utilities.

def _excval():
  """Return the most recent exception object."""
  return SYS.exc_info()[1]

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."""
  return C.MPMul().factor(one).factor(ff).done()

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.
  try: _MAX = SYS.maxint
  except AttributeError: _MAX = SYS.maxsize
  _NBIT = 15
  while 1 << (_NBIT + 1) < _MAX: _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 = sievesz*[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 ExpectedError('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:
        raise ExpectedError('%s:%d: %s' %
                            (file.name, lno, _excval().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] [-s SIEVEBITS] 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: exit('%s: %s' % (prog, _excval().message))
  except IOError: exit('%s: %s' % (prog, _excval()))

###----- That's all, folks --------------------------------------------------