gen_base_table.py

#!/usr/bin/env python3
"""Build the table s2n-bignum's p256_scalarmulbase reads.

    ./gen_base_table.py 5 p256_base_table.c

Its own words for the layout: for each i, j with blocksize*i <= 256 and
1 <= j <= B, where B = 2^(blocksize-1), the multiple 2^(blocksize*i) * j * P
goes at tab + 64*(B*i + (j-1)), as a Montgomery-affine pair, four
little-endian 64-bit words each.

Five is the blocksize crypton ships: on an Apple M4 it is 6.40 microseconds
against 7.15 for four and 6.10 for six, and the table is 52 KiB against 33
and 88.  Nothing outside this file knows the number -- it is written into
the generated file and read from there.

The curve constants come from `openssl ecparam`, not from memory, and the
generated table is checked against crypton's own base-point multiplication
by the test suite, so a mistake here does not pass quietly.
"""
import subprocess, re, sys

def curve():
    out = subprocess.run(["openssl","ecparam","-name","prime256v1",
                          "-param_enc","explicit","-text","-noout"],
                         capture_output=True, text=True).stdout
    f, cur = {}, None
    for line in out.splitlines():
        m = re.match(r'^(Prime|A|B|Generator \(uncompressed\)|Order):?\s*$', line.strip())
        if m:
            cur = m.group(1); f[cur] = ""; continue
        if cur and re.match(r'^\s+[0-9a-f]{2}[:0-9a-f]*:?\s*$', line):
            f[cur] += line.strip()
        elif cur and line and not line.startswith(' '):
            cur = None
    def num(s):
        return int.from_bytes(bytes(int(v,16) for v in s.strip(':').split(':') if v), 'big')
    g = bytes(int(v,16) for v in f['Generator (uncompressed)'].strip(':').split(':') if v)
    assert g[0] == 4 and len(g) == 65
    return (num(f['Prime']), num(f['A']), num(f['B']),
            int.from_bytes(g[1:33],'big'), int.from_bytes(g[33:],'big'), num(f['Order']))

P, A, B, GX, GY, N = curve()
assert (GY*GY - GX**3 - A*GX - B) % P == 0, "the generator is not on the curve"

def add(p, q):
    if p is None: return q
    if q is None: return p
    (x1,y1),(x2,y2) = p,q
    if x1 == x2:
        if (y1 + y2) % P == 0: return None
        l = (3*x1*x1 + A) * pow(2*y1, -1, P) % P
    else:
        l = (y2-y1) * pow(x2-x1, -1, P) % P
    x3 = (l*l - x1 - x2) % P
    return (x3, (l*(x1-x3) - y1) % P)

R = 1 << 256
def mont(v):  return v * R % P
def words(v): return [(v >> (64*k)) & 0xFFFFFFFFFFFFFFFF for k in range(4)]

def table(blocksize):
    Bn = 1 << (blocksize - 1)
    blocks = 256 // blocksize + 1
    out = []
    base = (GX, GY)                      # 2^(blocksize*i) * P
    for i in range(blocks):
        acc = None
        for j in range(1, Bn + 1):
            acc = add(acc, base)         # j * base
            out.append(acc)
        for _ in range(blocksize):
            base = add(base, base)
    return blocks, Bn, out

def emit(blocksize, path):
    blocks, Bn, pts = table(blocksize)
    with open(path, 'w') as fh:
        fh.write("/* Generated by cbits/p256/gen_base_table.py; do not hand-edit.\n"
                 " * The multiples of the base point that s2n-bignum's\n"
                 " * p256_scalarmulbase reads, in Montgomery-affine form. */\n")
        fh.write("#include <stdint.h>\n\n")
        fh.write(f"const uint64_t crypton_p256_s2n_base_blocksize = {blocksize};\n\n")
        fh.write(f"const uint64_t crypton_p256_s2n_base_table[{len(pts)*8}] = {{\n")
        for (x, y) in pts:
            ws = words(mont(x)) + words(mont(y))
            fh.write("\t" + ",".join(f"0x{w:016x}ULL" for w in ws) + ",\n")
        fh.write("};\n")
    return blocks, Bn, len(pts)*64

if __name__ == "__main__":
    b = int(sys.argv[1]); path = sys.argv[2]
    blocks, Bn, size = emit(b, path)
    print(f"blocksize {b}: {blocks} blocks x {Bn} = {size} bytes")