Clean up PEGASIS submodule inclusion

1 files changed, 0 insertions, 413 deletions

diff --git a/theta_lib/utilities/supersingular.py b/theta_lib/utilities/supersingular.py
deleted file mode 100644
index e760c1b..0000000
--- a/theta_lib/utilities/supersingular.py
+++ /dev/null
@@ -1,413 +0,0 @@
-"""
-Helper functions for the supersingular elliptic curve computations in FESTA
-"""
-
-# =========================================== #
-# Compute points of order D and Torsion Bases #
-# =========================================== #
-
-"""
-Most of this code has been taken from:
-https://github.com/FESTA-PKE/FESTA-SageMath
-
-Copyright (c) 2023 Andrea Basso, Luciano Maino and Giacomo Pope.
-
-Functions with another Copyright mention are not from the above authors.
-"""
-
-# Sage Imports
-from sage.all import ZZ, GF, Integers
-
-# Local imports
-from .order import has_order_D
-from .discrete_log import weil_pairing_pari, ell_discrete_log_pari, discrete_log_pari
-from .fast_sqrt import sqrt_Fp2, sqrt_Fp
-
-
-
-
-def random_point(E):
-    """
-    Returns a random point on the elliptic curve E
-    assumed to be in Montgomery form with a base 
-    field which characteristic p = 3 mod 4
-    """
-    A = E.a_invariants()[1]
-    if E.a_invariants() != (0,A,0,1,0):
-        raise ValueError("Function `generate_point` assumes the curve E is in the Montgomery model")
-
-    # Try 10000 times then give up, just protection
-    # for infinite loops
-    F = E.base_ring()
-    for _ in range(10000):
-        x = F.random_element()
-        y2 = x*(x**2 + A*x + 1)
-        if y2.is_square():
-            y = sqrt_Fp2(y2)
-            return E(x, y)
-
-    raise ValueError("Generated 10000 points, something is probably going wrong somewhere.")
-
-def generate_point(E, x_start=0):
-    """
-    Generate points on a curve E with x-coordinate 
-    i + x for x in Fp and i is the generator of Fp^2
-    such that i^2 = -1.
-    """
-    F = E.base_ring()
-    one = F.one()
-
-    if x_start:
-        x = x_start + one
-    else:
-        x = F.gen() + one
-
-    A = E.a_invariants()[1]
-    if E.a_invariants() != (0,A,0,1,0):
-        raise ValueError("Function `generate_point` assumes the curve E is in the Montgomery model")
-
-    # Try 10000 times then give up, just protection
-    # for infinite loops
-    for _ in range(10000):
-        y2 = x*(x**2 + A*x + 1)
-        if y2.is_square():
-            y = sqrt_Fp2(y2)
-            yield E(x, y)
-        x += one
-
-    raise ValueError("Generated 10000 points, something is probably going wrong somewhere.")
-
-def generate_Fp_point(E):
-    """
-    Returns a random Fp-rational point on the 
-    elliptic curve E assumed to be in Montgomery 
-    form with a base field which characteristic 
-    p = 3 mod 4 and coefficients defined over Fp.
-
-    Copyright (c) Pierrick Dartois 2025.
-    """
-
-    A = E.a_invariants()[1]
-    if E.a_invariants() != (0,A,0,1,0):
-        raise ValueError("Function `generate_Fp_point` assumes the curve E is in the Montgomery model")
-    if A[1] != 0:
-        raise ValueError("The curve E should be Fp-rational")
-
-    p=E.base_field().characteristic()
-    Fp=GF(p,proof=False)
-    one=Fp(1)
-    x=one
-
-    # Try 10000 times then give up, just protection
-    # for infinite loops
-    for _ in range(10000):
-        y2 = x*(x**2 + A*x + 1)
-        if y2.is_square():
-            y = sqrt_Fp(y2)
-            yield E(x, y)
-        x += one
-
-    raise ValueError("Generated 10000 points, something is probably going wrong somewhere.")
-
-def generate_point_order_D(E, D, x_start=0):
-    """
-    Input:  An elliptic curve E / Fp2
-            An integer D dividing (p +1)
-    Output: A point P of order D.
-    """
-    p = E.base().characteristic()
-    n = (p + 1) // D
-
-    Ps = generate_point(E, x_start=x_start)
-    for G in Ps:
-        P = n * G
-
-        # Case when we randomly picked
-        # a point in the n-torsion
-        if P.is_zero():
-            continue
-
-        # Check that P has order exactly D
-        if has_order_D(P, D):
-            P._order = ZZ(D)
-            yield P
-
-    raise ValueError(f"Never found a point P of order D.")
-
-def generate_Fp_point_order_D(E, D):
-    """
-    Input:  An elliptic curve E / Fp2
-            An integer D dividing (p +1)
-    Output: A point P of order D.
-
-    Copyright (c) Pierrick Dartois 2025.
-    """
-    p = E.base().characteristic()
-    n = (p + 1) // D
-    if D%2==0:
-        n=n//2
-
-    Ps = generate_Fp_point(E)
-    for G in Ps:
-        P = n * G
-
-        # Case when we randomly picked
-        # a point in the n-torsion
-        if P.is_zero():
-            continue
-
-        # Check that P has order exactly D
-        if has_order_D(P, D):
-            P._order = ZZ(D)
-            yield P
-
-    raise ValueError(f"Never found a point P of order D.")
-
-def compute_point_order_D(E, D, x_start=0):
-    """
-    Wrapper function around a generator which returns the first
-    point of order D
-    """
-    return generate_point_order_D(E, D, x_start=x_start).__next__()
-
-def compute_Fp_point_order_D(E, D):
-    """
-    Wrapper function around a generator which returns the first
-    point of order D
-
-    Copyright (c) Pierrick Dartois 2025.
-    """
-    return generate_Fp_point_order_D(E, D).__next__()
-
-def compute_linearly_independent_point_with_pairing(E, P, D, x_start=0):
-    """
-    Input:  An elliptic curve E / Fp2
-            A point P ∈ E[D]
-            An integer D dividing (p +1)
-    Output: A point Q such that E[D] = <P, Q>
-            The Weil pairing e(P,Q)
-    """
-    Qs = generate_point_order_D(E, D, x_start=x_start)
-    for Q in Qs:
-        # Make sure the point is linearly independent
-        pair = weil_pairing_pari(P, Q, D)
-        if has_order_D(pair, D, multiplicative=True):
-            Q._order = ZZ(D)
-            return Q, pair
-    raise ValueError("Never found a point Q linearly independent to P")
-
-def compute_linearly_independent_point(E, P, D, x_start=0):
-    """
-    Wrapper function around `compute_linearly_independent_point_with_pairing`
-    which only returns a linearly independent point
-    """
-    Q, _ = compute_linearly_independent_point_with_pairing(E, P, D, x_start=x_start)
-    return Q
-
-def torsion_basis_with_pairing(E, D):
-    """
-    Generate basis of E(Fp^2)[D] of supersingular curve
-
-    While computing E[D] = <P, Q> we naturally compute the 
-    Weil pairing e(P,Q), which we also return as in some cases
-    the Weil pairing is then used when solving the BiDLP
-    """
-    p = E.base().characteristic()
-
-    # Ensure D divides the curve's order
-    if (p + 1) % D != 0:
-        print(f"{ZZ(D).factor() = }")
-        print(f"{ZZ(p+1).factor() = }")
-        raise ValueError(f"D must divide the point's order")
-
-    P = compute_point_order_D(E, D)
-    Q, ePQ = compute_linearly_independent_point_with_pairing(E, P, D, x_start=P[0])
-
-    return P, Q, ePQ
-
-def torsion_basis(E, D):
-    """
-    Wrapper function around torsion_basis_with_pairing which only
-    returns the torsion basis <P,Q> = E[D]
-    """
-    P, Q, _ = torsion_basis_with_pairing(E, D)
-    return P, Q
-
-def torsion_basis_2f_frobenius(E,f):
-    """
-    Returns a basis (P, Q) of E[2^f] with pi_p(P)=P and pi_p(Q)=-Q,
-    pi_p being the p-th Frobenius endomorphism of E.
-
-    Copyright (c) Pierrick Dartois 2025.
-    """
-
-    P = compute_Fp_point_order_D(E, 2**f)
-
-    i = E.base_field().gen()
-    p = E.base_field().characteristic() 
-
-    R = compute_linearly_independent_point(E, P, 2**f)
-    piR = E(R[0]**p,R[1]**p)
-    piR_p_R = piR + R
-
-    two_lamb = ell_discrete_log_pari(E,piR_p_R,P,2**f)
-    lamb = two_lamb//2
-    Q = R-lamb*P
-
-    piQ = E(Q[0]**p,Q[1]**p)
-
-    assert piQ == -Q
-
-    return P, Q 
-
-def torsion_basis_to_Fp_rational_point(E,P,Q,D):
-    """
-    Returns an Fp-rational point R=[lamb]P+[mu]Q
-    of D-torsion given a D-torsion basis (P,Q) of E.
-
-    Copyright (c) Pierrick Dartois 2025.
-    """
-
-    p=E.base_field().characteristic()
-    piP = E(P[0]**p,P[1]**p)
-    piQ = E(Q[0]**p,Q[1]**p)
-
-    e_PQ = weil_pairing_pari(P,Q,D)
-    e_PpiP = weil_pairing_pari(P, piP, D)
-    e_QpiP = weil_pairing_pari(Q, piP, D)
-    e_PpiQ = weil_pairing_pari(P, piQ, D)
-    e_QpiQ = weil_pairing_pari(Q, piQ, D)
-
-    # pi(P)=[a]P+[b]Q, pi(Q)=[c]P+[d]Q
-    a = -discrete_log_pari(e_QpiP,e_PQ,D)
-    b = discrete_log_pari(e_PpiP,e_PQ,D)
-    c = -discrete_log_pari(e_QpiQ,e_PQ,D)
-    d = discrete_log_pari(e_PpiQ,e_PQ,D)
-
-    ZD = Integers(D)
-
-    if a%D==1:
-        mu = ZD(b/d)
-        R=P+mu*Q
-    elif c==0:
-        # c=0 and a!=1 => d=1 so pi(Q)=Q
-        R=Q
-    else:
-        # c!=0 and a!=1
-        mu = ZD((1-a)/c)
-        R=P+mu*Q
-
-    piR = E(R[0]**p,R[1]**p)
-
-    assert piR==R
-
-    return R
-
-# =========================================== #
-#   Entangled torsion basis for fast E[2^k]   #
-# =========================================== #
-
-def precompute_elligator_tables(F):
-    """
-    Precomputes tables of quadratic non-residue or 
-    quadratic residue in Fp2. Used to compute entangled
-    torsion bases following https://ia.cr/2017/1143
-    """
-    u = 2*F.gen()
-
-    T1 = dict()
-    T2 = dict()
-    # TODO: estimate how large r should be
-    for r in range(1, 30):
-        v = 1 / (1 + u*r**2)
-        if v.is_square():
-            T2[r] = v
-        else:
-            T1[r] = v
-    return T1, T2
-
-def entangled_torsion_basis(E, elligator_tables, cofactor):
-    """
-    Optimised algorithm following https://ia.cr/2017/1143
-    which modifies the elligator method of hashing to points
-    to find points P,Q of order k*2^b. Clearing the cofactor
-    gives the torsion basis without checking the order or
-    computing a Weil pairing. 
-
-    To do this, we need tables TQNR, TQR of pairs of values
-    (r, v) where r is an integer and v = 1/(1 + ur^2) where
-    v is either a quadratic non-residue or quadratic residue 
-    in Fp2 and u = 2i = 2*sqrt(-1).
-    """
-    F = E.base_ring()
-    p = F.characteristic()
-    p_sqrt = (p+1)//4
-
-    i = F.gen()
-    u =  2 * i
-    u0 = 1 + i
-
-    TQNR, TQR = elligator_tables
-    
-    # Pick the look up table depending on whether
-    # A = a + ib is a QR or NQR
-    A = E.a_invariants()[1]
-    if (0,A,0,1,0) != E.a_invariants():
-        raise ValueError("The elliptic curve E must be in Montgomery form")
-    if A.is_square():
-        T = TQNR
-    else:
-        T = TQR
-
-    # Look through the table to find point with 
-    # rational (x,y)
-    y = None
-    for r, v in T.items():
-        x = -A * v
-
-        t = x * (x**2 + A*x + 1)
-        
-        # Break when we find rational y: t = y^2 
-        c, d = t.list()
-        z = c**2 + d**2
-        s = z**p_sqrt
-        if s**2 == z:
-            y = sqrt_Fp2(t)
-            break
-    
-    if y is None:
-        raise ValueError("Never found a y-coordinate, increase the lookup table size")
-
-    z = (c + s) // 2
-    alpha = z**p_sqrt
-    beta  = d / (2*alpha)
-
-    if alpha**2 == z:
-        y =  F([alpha, beta])
-    else:
-        y = -F([beta, alpha])
-
-    S1 = E([x, y])
-    S2 = E([u*r**2*x, u0*r*y])
-
-    return cofactor*S1, cofactor*S2
-
-# =============================================== #
-#  Ensure Basis <P,Q> of E[2^k] has (0,0) under Q #
-# =============================================== #
-
-def fix_torsion_basis_renes(P, Q, k):
-    """
-    Set the torsion basis P,Q such that
-    2^(k-1)Q = (0,0) to ensure that (0,0)
-    is never a kernel of a two isogeny
-    """
-    cofactor = 2**(k-1)
-
-    R = cofactor*P
-    if R[0] == 0:
-        return Q, P
-    R = cofactor*Q
-    if R[0] == 0:
-        return P, Q
-    return P, P + Q