New script contrib/proveprime.py.

This generates primality certificates for numbers, in the form of
Python / testcrypt code that calls Pockle methods. It factors p-1 by
calling out to the 'yafu' utility, which is a moderately sophisticated
integer factoring tool (including ECC and quadratic sieve methods)
that runs as a standalone command-line program.

Also added a Pockle test generated as output from this script, which
verifies the primality of the three NIST curves' moduli and their
generators' orders. I already had Pockle certificates for the moduli
and orders used in EdDSA, so this completes the set, and it does it
without me having had to do a lot of manual work.

contrib/proveprime.py

@@ -0,0 +1,162 @@
+#!/usr/bin/env python3
+
+import argparse
+import functools
+import math
+import os
+import re
+import subprocess
+import sys
+import itertools
+
+def gen_names():
+    for i in itertools.count():
+        name = "p{:d}".format(i)
+        if name not in nameset:
+            yield name
+nameset=set()
+names = gen_names()
+
+class YafuError(Exception):
+    pass
+
+verbose = False
+def diag(*args):
+    if verbose:
+        print(*args, file=sys.stderr)
+
+factorcache = set()
+factorcachefile = None
+def cache_factor(f):
+    if f not in factorcache:
+        factorcache.add(f)
+    if factorcachefile is not None:
+        factorcachefile.write("{:d}\n".format(f))
+        factorcachefile.flush()
+
+yafu = None
+yafu_pattern = re.compile(rb"^P\d+ = (\d+)$")
+def call_yafu(n):
+    n_orig = n
+    diag("starting yafu", n_orig)
+    p = subprocess.Popen([yafu, "-v", "-v"], stdin=subprocess.PIPE,
+                         stdout=subprocess.PIPE)
+    p.stdin.write("{:d}\n".format(n).encode("ASCII"))
+    p.stdin.close()
+    factors = []
+    for line in iter(p.stdout.readline, b''):
+        line = line.rstrip(b"\r\n")
+        diag("yafu output:", line.decode())
+        m = yafu_pattern.match(line)
+        if m is not None:
+            f = int(m.group(1))
+            if n % f != 0:
+                raise YafuError("bad yafu factor {:d}".format(f))
+            factors.append(f)
+            if f >> 64:
+                cache_factor(f)
+            n //= f
+    p.wait()
+    diag("done yafu", n_orig)
+    return factors, n
+
+def factorise(n):
+    allfactors = []
+    for f in factorcache:
+        if n % f == 0:
+            n //= f
+            allfactors.append(f)
+    while n > 1:
+        factors, n = call_yafu(n)
+        allfactors.extend(factors)
+    return sorted(allfactors)
+
+def product(ns):
+    return functools.reduce(lambda a,b: a*b, ns, 1)
+
+smallprimes = set()
+commands = {}
+
+def proveprime(p, name=None):
+    if p >> 32 == 0:
+        smallprimes.add(p)
+        return "{:d}".format(p)
+
+    if name is None:
+        name = next(names)
+    print("{} = {:d}".format(name, p))
+
+    fs = factorise(p-1)
+    fs.reverse()
+    prod = product(fs)
+    qs = []
+    for q in fs:
+        newprod = prod // q
+        if newprod * newprod * newprod > p:
+            prod = newprod
+        else:
+            qs.append(q)
+    assert prod == product(qs)
+    assert prod * prod * prod > p
+    qset = set(qs)
+    qnamedict = {q: proveprime(q) for q in qset}
+    qnames = [qnamedict[q] for q in qs]
+    for w in itertools.count(2):
+        assert pow(w, p-1, p) == 1, "{}={:d} is not prime!".format(name, p)
+        diag("trying witness", w, "for", p)
+        for q in qset:
+            wpower = pow(w, (p-1) // q, p) - 1
+            if math.gcd(wpower, p) != 1:
+                break
+        else:
+            diag("found witness", w, "for", p)
+            break
+    commands[p]= (name, w, qnames)
+    return name
+
+def main():
+    parser = argparse.ArgumentParser(description='')
+    parser.add_argument("prime", nargs="+",
+                        help="Number to prove prime. Can be prefixed by a "
+                        "variable name and '=', e.g. 'x=9999999967'.")
+    parser.add_argument("--cryptsuite", action="store_true",
+                        help="Generate abbreviated Pockle calls suitable "
+                        "for the tests in cryptsuite.py.")
+    parser.add_argument("--yafu", default="yafu",
+                        help="yafu binary to help with factoring.")
+    parser.add_argument("-v", "--verbose", action="store_true",
+                        help="Write diagnostics to standard error.")
+    parser.add_argument("--cache", help="Cache of useful factors of things.")
+    args = parser.parse_args()
+
+    global verbose, yafu
+    verbose = args.verbose
+    yafu = args.yafu
+
+    if args.cache is not None:
+        with open(args.cache, "r") as fh:
+            for line in iter(fh.readline, ""):
+                factorcache.add(int(line.rstrip("\r\n")))
+        global factorcachefile
+        factorcachefile = open(args.cache, "a")
+
+    for ps in args.prime:
+        name, value = (ps.split("=", 1) if "=" in ps
+                       else (None, ps))
+        proveprime(int(value, 0), name)
+
+    print("po = pockle_new()")
+    if len(smallprimes) > 0:
+        if args.cryptsuite:
+            print("add_small(po, {})".format(
+                ", ".join("{:d}".format(q) for q in sorted(smallprimes))))
+        else:
+            for q in sorted(smallprimes):
+                print("pockle_add_small_prime(po, {:d})".format(q))
+    for p, (name, w, qnames) in sorted(commands.items()):
+        print("{cmd}(po, {name}, [{qs}], {w:d})".format(
+            cmd = "add" if args.cryptsuite else "pockle_add_prime",
+            name=name, w=w, qs=", ".join(qnames)))
+
+if __name__ == '__main__':
+    main()

test/cryptsuite.py

@@ -1079,6 +1079,49 @@ class keygen(MyTestBase):
         add(po, p6, [p3,p4], 2)
         add(po, p, [p2,p5,p6], 2)
 
+        # Combined certificate for the moduli and generator orders of
+        # the three NIST curves, generated by contrib/proveprime.py
+        # (with some cosmetic tidying)
+        p256 = 2**256 - 2**224 + 2**192 + 2**96 - 1
+        p384 = 2**384 - 2**128 - 2**96 + 2**32 - 1
+        p521 = 2**521 - 1
+        order256 = p256 - 0x4319055358e8617b0c46353d039cdaae
+        order384 = p384 - 0x389cb27e0bc8d21fa7e5f24cb74f58851313e696333ad68c
+        t = 0x5ae79787c40d069948033feb708f65a2fc44a36477663b851449048e16ec79bf6
+        order521 = p521 - t
+        p0 = order384 // 12895580879789762060783039592702
+        p1 = 1059392654943455286185473617842338478315215895509773412096307
+        p2 = 55942463741690639
+        p3 = 37344768852931
+        p4 = order521 // 1898873518475180724503002533770555108536
+        p5 = p4 // 994165722
+        p6 = 144471089338257942164514676806340723
+        p7 = p384 // 2054993070433694
+        p8 = 1357291859799823621
+        po = pockle_new()
+        add_small(po, 2, 3, 5, 11, 17, 19, 31, 41, 53, 67, 71, 109, 131, 149,
+                  157, 257, 521, 641, 1613, 2731, 3407, 6317, 8191, 8389,
+                  14461, 17449, 38189, 38557, 42641, 51481, 61681, 65537,
+                  133279, 248431, 312289, 409891, 490463, 858001, 6700417,
+                  187019741)
+        add(po, p3, [149, 11, 5, 3, 2], 3)
+        add(po, p2, [p3], 2)
+        add(po, p8, [6317, 67, 2, 2], 2)
+        add(po, p6, [133279, 14461, 109, 3], 7)
+        add(po, p1, [p2, 248431], 2)
+        add(po, order256, [187019741, 38189, 17449, 3407, 131, 71, 2, 2, 2, 2],
+            7)
+        add(po, p256, [6700417, 490463, 65537, 641, 257, 17, 5, 5, 3, 2], 6)
+        add(po, p0, [p1], 2)
+        add(po, p7, [p8, 312289, 38557, 8389, 11, 2], 3)
+        add(po, p5, [p6, 19], 2)
+        add(po, order384, [p0], 2)
+        add(po, p384, [p7], 2)
+        add(po, p4, [p5], 2)
+        add(po, order521, [p4], 2)
+        add(po, p521, [858001, 409891, 61681, 51481, 42641, 8191, 2731, 1613,
+                       521, 157, 131, 53, 41, 31, 17, 11, 5, 5, 3, 2], 3)
+
     def testPockleNegative(self):
         def add_small(po, p):
             self.assertEqual(pockle_add_small_prime(po, p), 'POCKLE_OK')