second-round version of SQIsign

Co-authored-by: Marius A. Aardal <marius.andre.aardal@gmail.com>
Co-authored-by: Gora Adj <gora.adj@tii.ae>
Co-authored-by: Diego F. Aranha <dfaranha@cs.au.dk>
Co-authored-by: Andrea Basso <sqisign@andreabasso.com>
Co-authored-by: Isaac Andrés Canales Martínez <icanalesm0500@gmail.com>
Co-authored-by: Jorge Chávez-Saab <jorgechavezsaab@gmail.com>
Co-authored-by: Maria Corte-Real Santos <mariascrsantos98@gmail.com>
Co-authored-by: Luca De Feo <github@defeo.lu>
Co-authored-by: Max Duparc <max.duparc@epfl.ch>
Co-authored-by: Jonathan Komada Eriksen <jonathan.eriksen97@gmail.com>
Co-authored-by: Décio Luiz Gazzoni Filho <decio@decpp.net>
Co-authored-by: Basil Hess <bhe@zurich.ibm.com>
Co-authored-by: Antonin Leroux <antonin.leroux@polytechnique.org>
Co-authored-by: Patrick Longa <plonga@microsoft.com>
Co-authored-by: Luciano Maino <mainoluciano.96@gmail.com>
Co-authored-by: Michael Meyer <michael@random-oracles.org>
Co-authored-by: Hiroshi Onuki <onuki@mist.i.u-tokyo.ac.jp>
Co-authored-by: Lorenz Panny <lorenz@yx7.cc>
Co-authored-by: Giacomo Pope <giacomopope@gmail.com>
Co-authored-by: Krijn Reijnders <reijnderskrijn@gmail.com>
Co-authored-by: Damien Robert <damien.robert@inria.fr>
Co-authored-by: Francisco Rodríguez-Henriquez <francisco.rodriguez@tii.ae>
Co-authored-by: Sina Schaeffler <sschaeffle@student.ethz.ch>
Co-authored-by: Benjamin Wesolowski <benjamin.wesolowski@ens-lyon.fr>

scripts/precomp/cformat.py

@@ -0,0 +1,128 @@
+#!/usr/bin/env python3
+import sys, itertools
+from math import ceil, floor, log
+import sage.all
+
+class Ibz:
+    def __init__(self, v):
+        self.v = int(v)
+    def _literal(self, sz):
+        val = int(self.v)
+        sgn = val < 0
+        num_limbs = (abs(val).bit_length() + sz-1) // sz if val else 0
+        limbs = [(abs(val) >> sz*i) & (2**sz-1) for i in range(num_limbs or 1)]
+        data = {
+                '._mp_alloc': 0,
+                '._mp_size': (-1)**sgn * num_limbs,
+                '._mp_d': '(mp_limb_t[]) {' + ','.join(map(hex,limbs)) + '}',
+            }
+        return '{{' + ', '.join(f'{k} = {v}' for k,v in data.items()) + '}}'
+
+class FpEl:
+    ref_p5248_radix_map  = { 16: 13, 32: 29, 64: 51 }
+    ref_p65376_radix_map = { 16: 13, 32: 28, 64: 55 }
+    ref_p27500_radix_map = { 16: 13, 32: 29, 64: 57 }
+    def __init__(self, n, p, montgomery=True):
+        self.n = n
+        self.p = p
+        self.montgomery = montgomery
+    def __get_radix(self, word_size, arith=None):
+        if arith == "ref" or arith is None:
+            # lvl1
+            if self.p == 0x4ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff:
+                return self.ref_p5248_radix_map[word_size]
+            # lvl3
+            elif self.p == 0x40ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff:
+                return self.ref_p65376_radix_map[word_size]
+            # lvl5
+            elif self.p == 0x1afffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff:
+                return self.ref_p27500_radix_map[word_size]
+            raise ValueError(f'Invalid prime \"{self.p}\"')
+        elif arith == "broadwell":
+            return word_size
+        raise ValueError(f'Invalid arithmetic implementation type \"{arith}\"')
+    def _literal(self, sz, arith=None):
+        radix = self.__get_radix(sz, arith=arith)
+        l = 1 + floor(log(self.p, 2**radix))
+        # If we're using Montgomery representation, we need to multiply
+        # by the Montgomery factor R = 2^nw (n = limb number, w = radix)
+        if self.montgomery:
+            R = 2**(radix * ceil(log(self.p, 2**radix)))
+        else:
+            R = 1
+        el = (self.n * R) % self.p
+        vs = [(int(el) >> radix*i) % 2**radix for i in range(l)]
+        return '{' + ', '.join(map(hex, vs)) + '}'
+
+class Object:
+    def __init__(self, ty, name, obj):
+        if '[' in ty:
+            idx = ty.index('[')
+            depth = ty.count('[]')
+            def rec(os, d):
+                assert d >= 0
+                if not d:
+                    return ()
+                assert isinstance(os,list) or isinstance(os,tuple)
+                r, = {rec(o, d-1) for o in os}
+                return (len(os),) + r
+            dims = rec(obj, depth)
+            self.ty = ty[:idx], ''.join(f'[{d}]' for d in dims)
+        else:
+            self.ty = ty, ''
+        self.name = name
+        self.obj = obj
+
+    def _declaration(self):
+        return f'extern const {self.ty[0]} {self.name}{self.ty[1]};'
+
+    def _literal(self):
+        def rec(obj):
+            if isinstance(obj, int):
+                if obj < 256: return str(obj)
+                else: return hex(obj)
+            if isinstance(obj, sage.all.Integer):
+                if obj < 256: return str(obj)
+                else: return hex(obj)
+            if isinstance(obj, Ibz):
+                literal = "\n#if 0"
+                for sz in (16, 32, 64):
+                    literal += f"\n#elif GMP_LIMB_BITS == {sz}"
+                    literal += f"\n{obj._literal(sz)}"
+                return literal + "\n#endif\n"
+            if isinstance(obj, FpEl):
+                literal = "\n#if 0"
+                for sz in (16, 32, 64):
+                    literal += f"\n#elif RADIX == {sz}"
+                    if sz == 64:
+                        literal += "\n#if defined(SQISIGN_GF_IMPL_BROADWELL)"
+                        literal += f"\n{obj._literal(sz, 'broadwell')}"
+                        literal += "\n#else"
+                        literal += f"\n{obj._literal(sz, 'ref')}"
+                        literal += "\n#endif"
+                    else:
+                        literal += f"\n{obj._literal(sz, 'ref')}"
+                return literal + "\n#endif\n"
+            if isinstance(obj, list) or isinstance(obj, tuple):
+                return '{' + ', '.join(map(rec, obj)) + '}'
+            if isinstance(obj, str):
+                return obj
+            raise NotImplementedError(f'unknown type {type(obj)} in Formatter')
+        return rec(self.obj)
+
+    def _definition(self):
+        return f'const {self.ty[0]} {self.name}{self.ty[1]} = ' + self._literal() + ';'
+
+class ObjectFormatter:
+    def __init__(self, objs):
+        self.objs = objs
+
+    def header(self, file=None):
+        for obj in self.objs:
+            assert isinstance(obj, Object)
+            print(obj._declaration(), file=file)
+
+    def implementation(self, file=None):
+        for obj in self.objs:
+            assert isinstance(obj, Object)
+            print(obj._definition(), file=file)

scripts/precomp/ec_params.sage

@@ -0,0 +1,41 @@
+#!/usr/bin/env python3
+
+from sage.all import *
+from parameters import p, f
+
+if __name__ == '__main__':
+
+    cof = (p+1)//(2**f)
+
+    from cformat import Object, ObjectFormatter
+
+    obj_cof = ObjectFormatter(
+        [
+            Object('digit_t[]', 'p_cofactor_for_2f', [cof]),
+        ]
+    )
+
+    with open("include/ec_params.h", "w") as hfile:
+        with open("ec_params.c", "w") as cfile:
+
+            hfile.write('#ifndef EC_PARAMS_H\n')
+            hfile.write('#define EC_PARAMS_H\n')
+            hfile.write('\n')
+
+            hfile.write('#include <fp.h>\n')
+            cfile.write('#include <ec_params.h>\n')
+            hfile.write('\n')
+
+            hfile.write(f'#define TORSION_EVEN_POWER {f}\n')
+            hfile.write('\n')
+
+            hfile.write('// p+1 divided by the power of 2\n')
+            cfile.write('// p+1 divided by the power of 2\n')
+            obj_cof.header(file=hfile)
+            obj_cof.implementation(file=cfile)
+            hfile.write(f'#define P_COFACTOR_FOR_2F_BITLENGTH {((p+1)//(2**f)).bit_length()}\n')
+            hfile.write('\n')
+            cfile.write('\n')
+
+
+            hfile.write('#endif\n')

scripts/precomp/maxorders.py

@@ -0,0 +1,88 @@
+from sage.all import *
+
+from sage.misc.banner import require_version
+if not require_version(10, 5, print_message=True):
+    exit('')
+
+from parameters import p, num_orders as num
+
+################################################################
+
+# Underlying theory:
+# - Ibukiyama, On maximal orders of division quaternion algebras with certain optimal embeddings
+# - https://ia.cr/2023/106 Lemma 10
+
+from sage.algebras.quatalg.quaternion_algebra import basis_for_quaternion_lattice
+bfql = lambda els: basis_for_quaternion_lattice(els, reverse=True)
+
+Quat1, (i,j,k) = QuaternionAlgebra(-1, -p).objgens()
+assert Quat1.discriminant() == p         # ramifies correctly
+
+O0mat = matrix([list(g) for g in [Quat1(1), i, (i+j)/2, (1+k)/2]])
+O0 = Quat1.quaternion_order(list(O0mat))
+
+orders = [ (1, identity_matrix(QQ,4), O0mat, i, O0mat, vector((1,0,0,0))) ]
+
+q = ZZ(1)
+while len(orders) < num:
+    q = next_prime(q)
+    if q % 4 != 1:  # restricting to q ≡ 1 (mod 4)
+        continue
+
+    Quatq, (ii,jj,kk) = QuaternionAlgebra(-q, -p).objgens()
+    if Quatq.discriminant() != p:       # ramifies incorrectly
+        continue
+
+    x, y = QuadraticForm(QQ, 2, [1,0,p]).solve(q)
+    gamma = x + j*y
+    assert gamma.reduced_norm() == q
+    ims1 = [Quat1(1), i*gamma, j, k*gamma]
+    assert ims1[1]**2 == -q
+    assert ims1[2]**2 == -p
+    assert ims1[1]*ims1[2] == ims1[3]
+    assert ims1[2]*ims1[1] == -ims1[3]
+    # (1,ii,jj,kk)->ims1 is an isomorphism Quatq->Quat1
+    iso1q = ~matrix(map(list, ims1))
+
+    r = min(map(ZZ, Mod(-p, 4*q).sqrt(all=True)))
+
+    basq = [
+            Quatq(1),
+            ii,
+            (1 + jj) / 2,
+            (r + jj) * ii / 2 / q,
+        ]
+
+    Oq = Quatq.quaternion_order(basq)
+    assert Oq.discriminant() == p   # is maximal
+
+    mat1 = matrix(map(list, basq)) * ~iso1q
+    O1 = Quat1.quaternion_order(list(mat1))
+    assert O1.discriminant() == p   # is maximal
+    assert j in O1                  # p-extremal
+
+    # look for an odd connecting ideal
+    I = O0 * O1
+    I *= I.norm().denominator()
+    assert I.is_integral()
+    for v in IntegralLattice(I.gram_matrix()).enumerate_short_vectors():
+        elt = sum(c*g for c,g in zip(v,I.basis()))
+        if ZZ(elt.reduced_norm() / I.norm()) % 2:
+            break
+    I = I * (elt.conjugate() / I.norm())
+    assert I.is_integral()
+    assert I.norm() % 2
+    assert I.left_order() == O0
+
+    O1_ = I.right_order()
+    assert O1_.unit_ideal() == elt * O1 * ~elt
+    idl1 = matrix(map(list, I.basis()))
+
+    # q
+    # isomorphism from (-1,-p) algebra to (-q,-p) algebra
+    # basis of maximal order O₁ in (-1,-p) algebra
+    # element sqrt(-q) in O₁ in (-1,-p) algebra
+    # basis of connecting ideal I from O₀ in (-1,-p) algebra
+    # element γ such that I has right order γ O₁ γ^-1
+    orders.append((q, iso1q, mat1, ims1[1], idl1, vector(elt)))
+

scripts/precomp/parameters.py

@@ -0,0 +1,16 @@
+#!/usr/bin/env python3
+from sage.all import *
+proof.all(False)  # faster
+
+import re
+for l in open('sqisign_parameters.txt'):
+    for k in ('lvl', 'p', 'num_orders'):
+        m = re.search(rf'^\s*{k}\s*=\s*([x0-9a-f]+)', l)
+        if m:
+            v = ZZ(m.groups()[0], 0)
+            globals()[k] = v
+
+f = (p+1).valuation(2)
+
+__all__ = ['lvl', 'p', 'f', 'num_orders']
+

scripts/precomp/precompute_E0_basis.sage

@@ -0,0 +1,38 @@
+#!/usr/bin/env sage
+proof.all(False)  # faster
+
+################################################################
+
+from parameters import p, f
+
+if p % 4 != 3:
+    raise NotImplementedError('requires p ≡ 3 (mod 4)')
+
+assert (1 << f).divides(p + 1)
+Fp2.<i> = GF((p,2), modulus=[1,0,1])
+E0 = EllipticCurve(Fp2, [1, 0])
+
+from torsion_basis import even_torsion_basis_E0
+P, Q = even_torsion_basis_E0(E0, f)
+
+################################################################
+
+from cformat import FpEl, Object, ObjectFormatter
+
+def Fp2_to_list(el):
+    return [FpEl(int(c), p, True) for c in Fp2(el)]
+
+objs = ObjectFormatter([
+        Object('fp2_t', 'BASIS_E0_PX', Fp2_to_list(P.x())),
+        Object('fp2_t', 'BASIS_E0_QX', Fp2_to_list(Q.x())),
+    ])
+
+################################################################
+
+with open('include/e0_basis.h','w') as hfile:
+    with open('e0_basis.c','w') as cfile:
+        print(f'#include <fp2.h>', file=hfile)
+        print(f'#include <e0_basis.h>', file=cfile)
+
+        objs.header(file=hfile)
+        objs.implementation(file=cfile)

scripts/precomp/precompute_endomorphism_action.sage

@@ -0,0 +1,303 @@
+#!/usr/bin/env sage
+proof.all(False)  # faster
+
+from sage.misc.banner import require_version
+if not require_version(10, 0, print_message=True):
+    exit('')
+
+################################################################
+
+from parameters import p, f
+from torsion_basis import even_torsion_basis_E0
+
+################################################################
+
+from sage.groups.generic import order_from_multiple
+pari.allocatemem(1 << 34)  # 16G
+
+if p % 4 != 3:
+    raise NotImplementedError('requires p ≡ 3 (mod 4)')
+
+assert (1 << f).divides(p + 1)
+Fp2.<i> = GF((p,2), modulus=[1,0,1])
+
+sqrtm1 = min(Fp2(-1).sqrt(all=True))
+
+def compute(q, mat, idl, iso1q):
+    print(f'\x1b[33m{q = }\x1b[0m')
+    E0 = EllipticCurve(Fp2, [1,0])
+    E0.set_order((p+1)^2)
+
+    if q == 1:
+        E1 = E0
+        P1, Q1 = even_torsion_basis_E0(E1, f)
+        print(f'E0 = {E1}')
+        print(f'P0 = {P1}')
+        print(f'Q0 = {Q1}')
+
+    else:
+        Quat.<i,j,k> = QuaternionAlgebra(-1, -p)
+        I = Quat.ideal(map(Quat, idl))
+#        print(f'{I = }')
+        O0 = Quat.quaternion_order(list(map(Quat, orders[0][2])))
+#        print(f'{O0 = }')
+        O1 = I.right_order()
+#        print(f'{O1 = }')
+        assert I.left_order() == O0
+        assert O0.is_maximal() and O1.is_maximal()
+        assert I.norm() % 2
+
+        from deuring2d import Deuring2D
+        ctx = Deuring2D(p)
+        assert ctx.O0.order == O0
+        assert ctx.E0 == E0
+        ctx.sqrtm1 = sqrtm1
+
+        P0, Q0 = data[0][1]
+
+        for deg in range(1,10):
+
+            print(f'trying {deg = }...')
+            ctx.e = E0.cardinality(extension_degree=2^deg).sqrt().valuation(2) - 1
+
+            first = True
+            for suitable in ctx.SuitableIdeals(I, attempts=10**6, bound=10**3):
+
+                if first:
+                    Fbig.<U> = Fp2.extension(2^deg)
+                    ctx.E0 = E0.change_ring(Fbig)
+                    ctx.P = P0.change_ring(Fbig)
+                    ctx.Q = Q0.change_ring(Fbig)
+                    assert ctx.e == ctx.E0.order().sqrt().valuation(2) - 1
+                    for _ in range(ctx.e - f):
+                        ctx.P = ctx.P.division_points(2)[0]
+                        ctx.Q = ctx.Q.division_points(2)[0]
+                    ctx.P.set_order(multiple=2^ctx.e)
+                    ctx.Q.set_order(multiple=2^ctx.e)
+                first = False
+
+                try:
+                    E1, P1, Q1 = ctx.IdealToIsogeny(I, suitable=suitable)
+                    break
+                except Deuring2D.Failure:
+                    continue
+
+            else:
+                continue
+            break
+
+        else:
+            raise NotImplementedError('Deuring2D failed')
+
+        E1 = E1.change_ring(Fp2)
+
+        j = GF(p)(E1.j_invariant())
+        X = polygen(GF(p))
+        for A,_ in sorted((256*(X^2-3)^3 - (X^2-4)*j).roots()):
+            E1_ = EllipticCurve(Fp2, [0,A,0,1,0])
+            try:
+                iso = min(E1.isomorphisms(E1_))
+                break
+            except ValueError:
+                pass
+        E1 = iso.codomain()
+        P1 = iso._eval(P1)
+        Q1 = iso._eval(Q1)
+        print(f'{E1 = }')
+
+        P1 *= ctx.P.order() // P0.order()
+        Q1 *= ctx.Q.order() // Q0.order()
+        P1 = P1.change_ring(Fp2)
+        Q1 = Q1.change_ring(Fp2)
+        print(f'{P1 = }')
+        print(f'{Q1 = }')
+        P1.set_order(P0.order())
+        Q1.set_order(Q0.order())
+        assert P0.order() == Q0.order() == P1.order() == Q1.order() == 2^f
+
+        assert P1.weil_pairing(Q1,2^f) == P0.weil_pairing(Q0,2^f)^I.norm()
+
+    if q == 1:
+        endo_i, = (a for a in E1.automorphisms() if a.scaling_factor() == sqrtm1)
+    else:
+        iso = E1.isomorphism(min(Fp2(-q).sqrt(all=True)), is_codomain=True)
+        try:
+            endo_i = iso * E1.isogeny(None, codomain=iso.domain(), degree=q)
+        except ValueError:
+            assert False
+#    assert endo_i^2 == -q
+
+    endo_1 = E1.scalar_multiplication(1)
+    endo_j = E1.frobenius_isogeny()
+    endo_k = endo_i * endo_j
+
+    if __debug__:
+        R = E1.random_point()
+        assert (endo_i^2)(R) == -q*R
+        assert (endo_j^2)(R) == -p*R
+        assert (endo_j*endo_i)(R) == -(endo_i*endo_j)(R)
+
+    denom = mat.denominator()
+    coprime = denom.prime_to_m_part(lcm(P1.order(), Q1.order()))
+    P1d, Q1d = (inverse_mod(coprime, T.order()) * T for T in (P1, Q1))
+
+    denom //= coprime
+
+    extdeg = next(d for d in range(1,denom+1) if ((denom<<f)^2).divides(E1.order(extension_degree=d)))
+    if extdeg == 1:
+        Fbig = Fp2
+    else:
+        Fbig.<U> = Fp2.extension(extdeg)
+
+    P1d, Q1d = (T.change_ring(Fbig) for T in (P1d, Q1d))
+    P1d.set_order(multiple=denom<<f)
+    for l,m in denom.factor():
+        for i in range(m):
+            assert l.divides(P1d.order())
+            P1d = P1d.division_points(l)[0]
+            P1d.set_order(multiple=denom<<f)
+            for Q1d_ in Q1d.division_points(l):
+                o = order_from_multiple(P1d.weil_pairing(Q1d_, P1d.order()), denom<<f, operation='*')
+                if o == P1d.order():
+                    Q1d = Q1d_
+                    break
+            else:
+                assert False
+    assert hasattr(P1d, '_order')
+    Q1d.set_order(multiple=denom<<f)
+
+    denom *= coprime
+
+    PQ1d = P1d, Q1d
+#    mat1 = matrix(Zmod(1<<f), [endo_1._eval(T).log(PQ1d) for T in PQ1d])
+#    assert mat1 == 1            # identity; omit
+    mati = matrix(Zmod(1<<f), [endo_i._eval(T).log(PQ1d) for T in PQ1d])
+    matj = matrix(Zmod(1<<f), [endo_j._eval(T).log(PQ1d) for T in PQ1d])
+#    matk = matrix(Zmod(1<<f), [endo_k._eval(T).log(PQ1d) for T in PQ1d])
+#    assert matk == matj * mati  # redundant; omit
+    matk = matj * mati
+
+    gens = []
+    for row in denom * mat:
+        endo = sum(ZZ(c)*e for c,e in zip(row, (endo_1,endo_i,endo_j,endo_k)))
+        gens.append(endo)
+    gen1, gen2, gen3, gen4 = gens
+
+    assert mat[0] == vector((1,0,0,0))
+#    mat1 = matrix(ZZ, [gen1._eval(T).log(PQ1d) for T in PQ1d]) / denom
+#    assert mat1 == 1            # identity; omit
+    mat2 = matrix(ZZ, [gen2._eval(T).log(PQ1d) for T in PQ1d]) / denom
+    mat3 = matrix(ZZ, [gen3._eval(T).log(PQ1d) for T in PQ1d]) / denom
+    mat4 = matrix(ZZ, [gen4._eval(T).log(PQ1d) for T in PQ1d]) / denom
+    mat2, mat3, mat4 = (M.change_ring(Zmod(1<<f)) for M in (mat2,mat3,mat4))
+
+    A = E1.a2()
+    assert E1.a_invariants() == (0,A,0,1,0)
+
+    return (A, (A+2)/4), (P1, Q1), (mati,matj,matk), (mat2,mat3,mat4)
+
+################################################################
+
+from maxorders import orders
+
+print('qs:', [q for q,_,_,_,_,_ in orders])
+
+todo = [(q, mat*iso1q, idl, iso1q) for q,iso1q,mat,_,idl,_ in orders]
+data = [None] * len(todo)
+
+assert todo[0][0] == 1
+data[0] = compute(*todo[0])  # compute this first; we need it for the others
+print(f'[\x1b[32m+\x1b[0m] finished precomputation for \x1b[36mq = {todo[0][0]}\x1b[0m.')
+
+####XXX
+##for idx,inp in enumerate(todo[1:],1):
+##    data[idx] = compute(*inp)
+##    print(f'[\x1b[32m+\x1b[0m] finished precomputation for \x1b[36mq = {inp[0]}\x1b[0m.')
+##todo = []
+####XXX
+
+for (inp,_),res in parallel(8)(compute)(todo[1:]):
+    q,_,_,_ = inp
+    idx, = (i for i,(qq,_,_,_) in enumerate(todo) if qq == q)
+    assert data[idx] is None
+    data[idx] = res
+    print(f'[\x1b[32m+\x1b[0m] finished precomputation for \x1b[36m{q = }\x1b[0m.')
+
+################################################################
+
+from cformat import FpEl, Ibz, Object, ObjectFormatter
+
+def Fp2_to_list(el):
+    return [FpEl(int(c), p, True) for c in Fp2(el)]
+
+def basis2field(P, Q):
+    vs = [
+            [Fp2_to_list(T[0]), Fp2_to_list(T[2])]
+            for T in (P,Q,P-Q)
+        ]
+    return vs
+
+################################################################
+
+objs = ObjectFormatter([
+        Object('curve_with_endomorphism_ring_t[]', 'CURVES_WITH_ENDOMORPHISMS',
+            [
+                [
+                    [Fp2_to_list(A), Fp2_to_list(1),                    # ec_curve_t A, C
+                     [Fp2_to_list(A24), Fp2_to_list(1)], "true"],       # ec_curve_t A24, is_A24_computed_and_normalized
+                    basis2field(*basis),                                # ec_basis_t
+                    [[Ibz(v) for v in vs] for vs in mati.transpose()],  # ibz_mat_2x2_t
+                    [[Ibz(v) for v in vs] for vs in matj.transpose()],  # ibz_mat_2x2_t
+                    [[Ibz(v) for v in vs] for vs in matk.transpose()],  # ibz_mat_2x2_t
+                    [[Ibz(v) for v in vs] for vs in mat2.transpose()],  # ibz_mat_2x2_t
+                    [[Ibz(v) for v in vs] for vs in mat3.transpose()],  # ibz_mat_2x2_t
+                    [[Ibz(v) for v in vs] for vs in mat4.transpose()],  # ibz_mat_2x2_t
+                ]
+                for (A,A24),basis,(mati,matj,matk),(mat2,mat3,mat4)
+                in data
+            ])
+    ])
+
+with open('include/endomorphism_action.h','w') as hfile:
+    with open('endomorphism_action.c','w') as cfile:
+        print(f'#ifndef ENDOMORPHISM_ACTION_H', file=hfile)
+        print(f'#define ENDOMORPHISM_ACTION_H', file=hfile)
+        print(f'#include <sqisign_namespace.h>', file=hfile)
+        print(f'#include <ec.h>', file=hfile)
+        print(f'#include <quaternion.h>', file=hfile)
+        print(f'#include <stddef.h>', file=cfile)
+        print(f'#include <stdint.h>', file=cfile)
+        print(f'#include <endomorphism_action.h>', file=cfile)
+
+        print('''
+/** Type for precomputed endomorphism rings applied to precomputed torsion bases.
+ *
+ * Precomputed by the precompute scripts.
+ *
+ * @typedef curve_with_endomorphism_ring_t
+ *
+ * @struct curve_with_endomorphism_ring
+ **/
+typedef struct curve_with_endomorphism_ring {
+    ec_curve_t curve;
+    ec_basis_t basis_even;
+    ibz_mat_2x2_t action_i, action_j, action_k;
+    ibz_mat_2x2_t action_gen2, action_gen3, action_gen4;
+} curve_with_endomorphism_ring_t;
+              '''.strip(), file=hfile)
+
+        print(f'#define CURVE_E0 (CURVES_WITH_ENDOMORPHISMS->curve)', file=hfile)
+        print(f'#define BASIS_EVEN (CURVES_WITH_ENDOMORPHISMS->basis_even)', file=hfile)
+        print(f'#define ACTION_I (CURVES_WITH_ENDOMORPHISMS->action_i)', file=hfile)
+        print(f'#define ACTION_J (CURVES_WITH_ENDOMORPHISMS->action_j)', file=hfile)
+        print(f'#define ACTION_K (CURVES_WITH_ENDOMORPHISMS->action_k)', file=hfile)
+        print(f'#define ACTION_GEN2 (CURVES_WITH_ENDOMORPHISMS->action_gen2)', file=hfile)
+        print(f'#define ACTION_GEN3 (CURVES_WITH_ENDOMORPHISMS->action_gen3)', file=hfile)
+        print(f'#define ACTION_GEN4 (CURVES_WITH_ENDOMORPHISMS->action_gen4)', file=hfile)
+        print(f'#define NUM_ALTERNATE_STARTING_CURVES {len(data)-1}', file=hfile)
+        print(f'#define ALTERNATE_STARTING_CURVES (CURVES_WITH_ENDOMORPHISMS+1)', file=hfile)
+
+        objs.header(file=hfile)
+        objs.implementation(file=cfile)
+
+        print(f'#endif', file=hfile)

scripts/precomp/precompute_hd_splitting.sage

@@ -0,0 +1,158 @@
+#!/usr/bin/env sage
+proof.all(False)  # faster
+
+################################################################
+
+from parameters import p
+
+# Field
+Fp2.<i> = GF((p,2), modulus=[1,0,1])
+
+Fp2_constants = [
+    [Fp2(0), Fp2(1), Fp2(i), Fp2(-1), Fp2(-i)],
+    ["FP2_ZERO", "FP2_ONE", "FP2_I", "FP2_MINUS_ONE", "FP2_MINUS_I"]
+]
+
+################################################################
+
+from cformat import FpEl
+
+def Fp2_to_list(el):
+    return [FpEl(int(c), p, True) for c in Fp2(el)]
+
+def Fp2_to_name(el):
+    return Fp2_constants[1][Fp2_constants[0].index(el)]
+
+################################################################
+
+# Splitting Data
+
+chi_eval = [
+    [1,1,1,1],
+    [1,-1,1,-1],
+    [1,1,-1,-1],
+    [1,-1,-1,1]
+]
+
+even_indices = [
+    [0, 0],
+    [0, 1],
+    [0, 2],
+    [0, 3],
+    [1, 0],
+    [1, 2],
+    [2, 0],
+    [2, 1],
+    [3, 0],
+    [3, 3],
+]
+
+splitting_map = {
+    (0, 2): [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, -1, 0]],
+    (3, 3): [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]],
+    (0, 3): [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, -1]],
+    (2, 1): [[1, 1, 1, 1], [1, -1, 1, -1], [1, -1, -1, 1], [1, 1, -1, -1]],
+    (0, 1): [[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, -1, 0, 0]],
+    (1, 2): [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]],
+    (2, 0): [[1, 1, 1, 1], [1, -1, 1, -1], [1, -1, -1, 1], [-1, -1, 1, 1]],
+    (3, 0): [[1, 1, 1, 1], [1, -1, 1, -1], [1, 1, -1, -1], [-1, 1, 1, -1]],
+    (1, 0): [[1, 1, 1, 1], [1, -1, -1, 1], [1, 1, -1, -1], [-1, 1, -1, 1]],
+    (0, 0): [[1, i, 1, i], [1, -i, -1, i], [1, i, -1, -i], [-1, i, -1, i]],
+}
+
+# 24 2x2 maps used for normalization. Applying a uniform matrix to a level 2
+# theta null point ensure that we get a random theta null point in the
+# equivalence class of Γ/Γ(2,4)
+full_normalization_maps_2d = [
+    matrix(2, 2, [1, 0, 0, 1]),
+    matrix(2, 2, [0, 1, 1, 0]),
+    matrix(2, 2, [1, 0, 0, -1]),
+    matrix(2, 2, [0, 1, -1, 0]),
+
+    matrix(2, 2, [1, 1, 1, -1]),
+    matrix(2, 2, [1, -1, 1, 1]),
+    matrix(2, 2, [1, 1, -1, 1]),
+    matrix(2, 2, [-1, 1, 1, 1]),
+
+    matrix(2, 2, [1, 0, 0, i]),
+    matrix(2, 2, [0, i, 1, 0]),
+    matrix(2, 2, [1, 0, 0, -i]),
+    matrix(2, 2, [0, i, -1, 0]),
+
+    matrix(2, 2, [i, 1, 1, i]),
+    matrix(2, 2, [1, i, i, 1]),
+    matrix(2, 2, [i, 1, -1, -i]),
+    matrix(2, 2, [1, i, -i, -1]),
+
+    matrix(2, 2, [1, i, 1, -i]),
+    matrix(2, 2, [1, -i, 1, i]),
+    matrix(2, 2, [1, i, -1, i]),
+    matrix(2, 2, [1, -i, -1, -i]),
+
+    matrix(2, 2, [1, 1, i, -i]),
+    matrix(2, 2, [i, -i, 1, 1]),
+    matrix(2, 2, [1, 1, -i, i]),
+    matrix(2, 2, [i, -i, -1, -1]),
+]
+
+# 6 2x2 maps used for normalisation. A subset of the preceding 24 matrices,
+# that are sufficient to ensure uniform normalisation (under Γ/Γ^0(4))
+# when using the Montgomery model
+normalization_maps_2d = [
+    matrix(2, 2, [1, 0, 0, 1]),
+    matrix(2, 2, [0, 1, 1, 0]),
+    matrix(2, 2, [1, 1, 1, -1]),
+    matrix(2, 2, [-1, 1, 1, 1]),
+    matrix(2, 2, [i, 1, 1, i]),
+    matrix(2, 2, [1, i, i, 1]),
+]
+
+
+# Format from dictionary to list of lists
+splitting_matrices = []
+for ind in even_indices:
+    # Create a list of all ten matrices (represented as 4x4 matrices)
+    splitting_matrices.append([[Fp2_to_name(Fp2(x)) for x in row] for row in splitting_map[tuple(ind)]])
+
+# 6 4x4 maps constructed from the above
+normalization_maps_4d = []
+for m in normalization_maps_2d:
+    M = m.tensor_product(m, subdivide=False).list()
+    matrix_elements_list = list(map(Fp2, M))
+    # Reshape into matrix
+    matrix_elements = [ matrix_elements_list[i:i+4] for i in range(0, len(matrix_elements_list), 4) ]
+    normalization_maps_4d.append([[Fp2_to_name(x) for x in row] for row in matrix_elements])
+
+################################################################
+
+from cformat import Object, ObjectFormatter
+
+objs = ObjectFormatter(
+    [
+        Object('int[][]', 'EVEN_INDEX', even_indices),
+        Object('int[][]', 'CHI_EVAL', chi_eval),
+        Object('fp2_t[]', 'FP2_CONSTANTS', list(map(Fp2_to_list, Fp2_constants[0]))),
+        Object('precomp_basis_change_matrix_t[]', 'SPLITTING_TRANSFORMS', [[x] for x in splitting_matrices]),
+        Object('precomp_basis_change_matrix_t[]', 'NORMALIZATION_TRANSFORMS', [[x] for x in normalization_maps_4d]),
+    ]
+)
+
+with open("include/hd_splitting_transforms.h", "w") as hfile:
+    with open("hd_splitting_transforms.c", "w") as cfile:
+        print("#ifndef HD_SPLITTING_H", file=hfile)
+        print("#define HD_SPLITTING_H", file=hfile)
+        print(f"\n#include <hd.h>", file=hfile)
+        print(f"#include <stdint.h>\n", file=hfile)
+        print("typedef struct precomp_basis_change_matrix {", file=hfile)
+        print("    uint8_t m[4][4];", file=hfile)
+        print("} precomp_basis_change_matrix_t;\n", file=hfile)
+
+        print(f"#include <hd_splitting_transforms.h>\n", file=cfile)
+        for i in range(len(Fp2_constants[1])):
+            print(f"#define {Fp2_constants[1][i]} {i}", file=cfile)
+        print("", file=cfile)
+
+        objs.header(file=hfile)
+        objs.implementation(file=cfile)
+
+        print("\n#endif\n", file=hfile)

scripts/precomp/precompute_quaternion_constants.sage

@@ -0,0 +1,44 @@
+#!/usr/bin/env sage
+proof.all(False)  # faster
+
+################################################################
+
+from parameters import p
+negl = 2**-64
+
+################################################################
+
+logp = ceil(log(p, 2))
+loglogp = ceil(log(logp,2))
+tors2val = (p+1).valuation(2)
+
+defs = dict()
+
+# RepresentInteger data
+small = ceil(log(negl, 2) / -1)
+assert 2**-small <= negl
+
+add_shift = ceil(log(log(negl, 1-1/(64*logp)), 2))
+assert (1 - 1/(64*logp)) ** (2**(add_shift)) <= negl
+
+defs['QUAT_primality_num_iter'] = ceil(-log(negl, 4))
+defs['QUAT_repres_bound_input'] = add_shift
+
+# Equivalent ideal data
+defs['QUAT_equiv_bound_coeff'] = 2**(1 + add_shift//4)
+
+# Find_uv constants
+m = 2 + floor((logp - tors2val) / 4)
+defs['FINDUV_box_size'] = m
+defs['FINDUV_cube_size'] = (2 * m + 1)**4 - 1
+
+################################################################
+
+with open('include/quaternion_constants.h','w') as hfile:
+    print(f'#include <quaternion.h>', file=hfile)
+
+    for k,v in defs.items():
+        v = ZZ(v)
+        print(f'#define {k} {v}', file=hfile)
+
+

scripts/precomp/precompute_quaternion_data.sage

@@ -0,0 +1,91 @@
+#!/usr/bin/env sage
+proof.all(False)  # faster
+
+
+
+from maxorders import p, orders
+
+from cformat import Ibz, Object, ObjectFormatter
+
+# Prime of same size than p for random ideal of fixed norm
+bitlength_p = int(p).bit_length()
+prime_cofactor = next_prime((2^(bitlength_p)))
+
+algobj = [Ibz(p)]
+
+objs = \
+    [
+        [
+            # basis (columns)
+            [
+                Ibz(mat.denominator()),
+                [[Ibz(v) for v in vs]
+                    for vs in mat.transpose()*mat.denominator()],
+            ],
+            # sqrt(-q)
+            [
+                Ibz(mat.denominator()),
+                [Ibz(c) for c in ii*mat.denominator()],
+            ],
+            # sqrt(-p)
+            [
+                Ibz(1),
+                [Ibz(c) for c in (0,0,1,0)]
+            ],
+            q
+        ]
+        for q,_,mat,ii,_,_ in orders
+    ]
+
+idlobjs = \
+    [
+        [
+            # basis (columns)
+            [
+                Ibz(idl.denominator()),
+                [[Ibz(v) for v in vs]
+                    for vs in idl.transpose()*idl.denominator()],
+            ],
+            # norm
+            Ibz(abs(idl.row_space(ZZ).intersection((ZZ^4).submodule([[1,0,0,0]])).basis()[0][0])),
+            # left order
+            '&MAXORD_O0',
+        ]
+        for _,_,mat,_,idl,_ in orders
+    ]
+
+gammaobjs = \
+    [
+        [
+            Ibz(gamma.denominator()),
+            list(map(Ibz, gamma * gamma.denominator())),
+        ]
+        for _,_,_,_,_,gamma in orders
+    ]
+
+objs = ObjectFormatter([
+        Object('ibz_t', 'QUAT_prime_cofactor', Ibz(prime_cofactor)),
+        Object('quat_alg_t', 'QUATALG_PINFTY', algobj),
+        Object('quat_p_extremal_maximal_order_t[]', 'EXTREMAL_ORDERS', objs),
+        Object('quat_left_ideal_t[]', 'CONNECTING_IDEALS', idlobjs),  # ideal corresponding to an isogeny from E0 which acts as identity w.r.t. the basis_even
+        Object('quat_alg_elem_t[]', 'CONJUGATING_ELEMENTS', gammaobjs), # elements γ such that each I has right order γ O₁ γ^-1
+    ])
+
+with open('include/quaternion_data.h','w') as hfile:
+    with open('quaternion_data.c','w') as cfile:
+        print(f'#include <quaternion.h>', file=hfile)
+        print(f'#include <stddef.h>', file=cfile)
+        print(f'#include <stdint.h>', file=cfile)
+        print(f'#include <quaternion_data.h>', file=cfile)
+
+        #FIXME this should eventually go away?
+        print(f'#define MAXORD_O0 (EXTREMAL_ORDERS->order)', file=hfile)
+        print(f'#define STANDARD_EXTREMAL_ORDER (EXTREMAL_ORDERS[0])', file=hfile)
+        print(f'#define NUM_ALTERNATE_EXTREMAL_ORDERS {len(orders)-1}', file=hfile)
+        print(f'#define ALTERNATE_EXTREMAL_ORDERS (EXTREMAL_ORDERS+1)', file=hfile)
+        print(f'#define ALTERNATE_CONNECTING_IDEALS (CONNECTING_IDEALS+1)', file=hfile)
+        print(f'#define ALTERNATE_CONJUGATING_ELEMENTS (CONJUGATING_ELEMENTS+1)', file=hfile)
+
+        objs.header(file=hfile)
+        objs.implementation(file=cfile)
+

scripts/precomp/precompute_sizes.sage

@@ -0,0 +1,93 @@
+#!/usr/bin/env sage
+proof.all(False)  # faster
+
+from sage.misc.banner import require_version
+if not require_version(9, 8, print_message=True):
+    exit('')
+
+################################################################
+
+from parameters import lvl, f, p
+
+################################################################
+
+logp = ceil(log(p, 2))
+tors2val = (p+1).valuation(2)
+tors2part = (p+1).p_primary_part(2)
+tors3part = (p+1).p_primary_part(3)
+
+defs = dict()
+
+TORSION_2POWER_BYTES = (tors2part.bit_length() + 7) // 8
+SECURITY_BITS = round(p.bit_length() / 128) * 64
+RESPONSE_LENGTH = ceil(p.bit_length()/2)
+RESPONSE_BYTES = (RESPONSE_LENGTH + 9) // 8
+
+fpsz = (logp + 63)//64*8
+fp2sz = 2 * fpsz
+defs['SECURITY_BITS'] = SECURITY_BITS
+defs['SQIsign_response_length'] = ceil(logp/2)
+defs['HASH_ITERATIONS'] = 2**(32 * ceil( logp/64 ) - (tors2val - ceil(logp/2)))
+defs['FP_ENCODED_BYTES'] = fpsz
+defs['FP2_ENCODED_BYTES'] = fp2sz
+defs['EC_CURVE_ENCODED_BYTES'] = fp2sz  # just the A
+defs['EC_POINT_ENCODED_BYTES'] = fp2sz  # just the x
+defs['EC_BASIS_ENCODED_BYTES'] = 3 * defs['EC_POINT_ENCODED_BYTES']
+
+defs['PUBLICKEY_BYTES'] = defs['EC_CURVE_ENCODED_BYTES'] + 1  # extra byte for hint
+defs['SECRETKEY_BYTES'] = defs['PUBLICKEY_BYTES'] + 5*defs['FP_ENCODED_BYTES'] + 4*TORSION_2POWER_BYTES
+defs['SIGNATURE_BYTES'] = defs['EC_CURVE_ENCODED_BYTES'] + 2 + 4*RESPONSE_BYTES + (SECURITY_BITS//8) + 1 + 1
+
+size_privkey = defs['SECRETKEY_BYTES']
+size_pubkey = defs['PUBLICKEY_BYTES']
+size_signature = defs['SIGNATURE_BYTES']
+
+algname = f'SQIsign_lvl{lvl}'
+
+################################################################
+
+with open('include/encoded_sizes.h','w') as hfile:
+    for k,v in defs.items():
+        v = ZZ(v)
+        print(f'#define {k} {v}', file=hfile)
+
+################################################################
+
+api = f'''
+// SPDX-License-Identifier: Apache-2.0
+
+#ifndef api_h
+#define api_h
+
+#include <sqisign_namespace.h>
+
+#define CRYPTO_SECRETKEYBYTES {size_privkey}
+#define CRYPTO_PUBLICKEYBYTES {size_pubkey}
+#define CRYPTO_BYTES {size_signature}
+
+#define CRYPTO_ALGNAME "{algname}"
+
+#if defined(ENABLE_SIGN)
+SQISIGN_API
+int
+crypto_sign_keypair(unsigned char *pk, unsigned char *sk);
+
+SQISIGN_API
+int
+crypto_sign(unsigned char *sm, unsigned long long *smlen,
+            const unsigned char *m, unsigned long long mlen,
+            const unsigned char *sk);
+#endif
+
+SQISIGN_API
+int
+crypto_sign_open(unsigned char *m, unsigned long long *mlen,
+                 const unsigned char *sm, unsigned long long smlen,
+                 const unsigned char *pk);
+
+#endif /* api_h */
+'''.strip()
+
+with open(f'../../../nistapi/lvl{lvl}/api.h', 'w') as f:
+    print(api, file=f)
+

scripts/precomp/precompute_torsion_constants.sage

@@ -0,0 +1,40 @@
+#!/usr/bin/env sage
+proof.all(False)  # faster
+
+################################################################
+
+from parameters import p
+
+################################################################
+
+tors2part = (p+1).p_primary_part(2)
+lambda_security = round(p.bit_length() / 128) * 64
+N_sec = next_prime(1 << 4*lambda_security)
+N_com = N_sec
+
+defs = {
+        'TORSION_2POWER_BYTES': (tors2part.bit_length() + 7) // 8,
+    }
+
+from cformat import Ibz, Object, ObjectFormatter
+
+objs = ObjectFormatter([
+        Object('ibz_t', 'TWO_TO_SECURITY_BITS', Ibz(1 << lambda_security)),  # lambda_security = SECURITY_BITS (128, 192, 256)
+        Object('ibz_t', 'TORSION_PLUS_2POWER', Ibz(tors2part)),
+        Object('ibz_t', 'SEC_DEGREE', Ibz(N_sec)),
+        Object('ibz_t', 'COM_DEGREE', Ibz(N_com)),
+    ])
+
+with open('include/torsion_constants.h','w') as hfile:
+    with open('torsion_constants.c','w') as cfile:
+        print(f'#include <quaternion.h>', file=hfile)
+        print(f'#include <stddef.h>', file=cfile)
+        print(f'#include <stdint.h>', file=cfile)
+        print(f'#include <torsion_constants.h>', file=cfile)
+
+        for k,v in defs.items():
+            print(f'#define {k} {v}', file=hfile)
+
+        objs.header(file=hfile)
+        objs.implementation(file=cfile)
+

scripts/precomp/torsion_basis.py

@@ -0,0 +1,75 @@
+from sage.all import ZZ, GF, EllipticCurve, parallel
+
+def even_torsion_basis_E0(E0, f):
+    """
+    For the case when A = 0 we can't use the entangled basis algorithm
+    so we do something "stupid" to simply get something canonical
+    """
+    assert E0.a_invariants() == (0, 0, 0, 1, 0)
+
+    Fp2 = E0.base_ring()
+    p = Fp2.characteristic()
+
+    def points_order_two_f():
+        """
+        Compute a point P of order 2^f with x(P) = 1 + i*x_im
+        """
+        x_im = 0
+        while True:
+            x_im += 1
+            x = Fp2([1, x_im])
+            if not E0.is_x_coord(x):
+                continue
+            # compares a+bi <= c+di iff (a,b) <= (c,d) as tuples, where integers
+            # modulo p are compared via their minimal non-negative representatives
+            P = min(E0.lift_x(x, all=True), key = lambda pt: list(pt.y()))
+            P.set_order(multiple=p+1)
+            if P.order() % (1 << f) == 0:
+                P *= P.order() // (1 << f)
+                P.set_order(1 << f)
+                yield P
+
+    pts = points_order_two_f()
+    P = next(pts)
+    for Q in pts:
+        # Q is picked to be in E[2^f] AND we must ensure that
+        # <P, Q> form a basis, which is the same as e(P, Q) having
+        # full order 1 << f.
+        e = P.weil_pairing(Q, 1 << f)
+        if e ** (1 << f - 1) == -1:
+            break
+
+    # Finally we want to make sure Q is above (0, 0)
+    P2 = (1 << f - 1) * P
+    Q2 = (1 << f - 1) * Q
+    if Q2 == E0(0, 0):
+        pass
+    elif P2 == E0(0, 0):
+        P, Q = Q, P
+    else:
+        Q += P
+
+    assert P.weil_pairing(Q, 1 << f) ** (1 << f - 1) == -1
+    assert (1 << f - 1) * Q == E0(0, 0)
+
+    return P, Q
+
+
+if __name__ == "__main__":
+    # p, f = 5 * 2**248 - 1, 248
+    # p, f = 65 * 2**376 - 1, 376
+    p, f = 27 * 2**500 - 1, 500
+    print(f"p = {ZZ(p+1).factor()} - 1")
+    Fp2 = GF(p**2, modulus=[1, 0, 1], names="i")
+    E = EllipticCurve(Fp2, [1, 0])
+    E.set_order((p + 1) ** 2)
+
+    P, Q = even_torsion_basis_E0(E, f)
+    print(f"{P = }")
+    print(f"{Q = }")
+
+    assert P.order() == 1 << f
+    assert Q.order() == 1 << f
+    e = P.weil_pairing(Q, 1 << f)
+    assert e ** (1 << f - 1) == -1
+    print("all good")