ec.py

from sage.all import (
    pari, Zmod, matrix, ZZ, isqrt, gcd, Integers
)

try:
    from . import params
except ImportError:
    import params

def PointCoords(X, P, Q, zeta = None, e = None):
    """
    Given a point X in E[2^e] and a basis (P, Q) of E[2^f] compute the vector
    (a, b) such that X = a*P + b*Q.
    Input:
    - X: a point in E[2^e]
    - P, Q: a basis for E[2^f]
    - zeta: (optional) the scaled tate pairing between P and Q
    - e: (optional) the order of X; if None, params.f is used; if e !=
      params.f, the output (a, b) is such that
        (2**(params.f - e)) * (aP + bQ) = X
    Output:
    - a, b such that X = a*P + b*Q
    """
    if not e:
        e = params.f
    E = X.curve()

    cof = (params.p**2 - 1) // 2**e
    if not zeta:
        zeta = pari.elltatepairing(E, P, Q, 2**params.f) ** cof

    zeta1 = pari.elltatepairing(E, X, Q, 2**e) ** cof
    zeta2 = pari.elltatepairing(E, X, -P, 2**e) ** cof

    a = ZZ(pari.fflog(zeta1, zeta, 2**e))
    b = ZZ(pari.fflog(zeta2, zeta, 2**e))

    assert (2**(params.f - e)) * (a*P + b*Q) == X

    return a, b

def ChangeOfBasis(B1, B2, e=None):
    """
    Compute the matrix M so that M * B1 = B2 where B1 = (P1, P2)^T and B2 =
    (Q1, Q2)^T are two basis of the 2^f torsion [Algorithm 2.5].
    If B1 is None, the standard basis (P0, Q0) on E0 and the precomputed tate
    pairing are used.
    Input:
    - B1 = (P1, P2) a basis for E[2^f]
    - B2 = (Q1, Q2) a basis for E[2^e]
    - e: optional, the order of the points; if None params.f is used
    Output:
    - A matrix M =(x_i) so that Q1 = x1*P1 + x2*P2 and Q2 = x3*P1 + x4*P2
    """
    if not e:
        e = params.f

    cof = (params.p**2 - 1) // 2**e
    if not B1:
        P1, P1 = params.P0, params.Q0
        zeta = params.unred_tate ** cof
    else:
        P1, P2 = B1
        E = P1.curve()
        zeta = pari.elltatepairing(E, P1, P2, 2**params.f) ** cof

    Q1, Q2 = B2

    x1, x2 = PointCoords(Q1, P1, P2, zeta, e)
    x3, x4 = PointCoords(Q2, P1, P2, zeta, e)

    return matrix(Zmod(2**e), 2, [x1, x2, x3, x4])

def SetChangeOfBasisMatrix(E1, E2, P1, Q1, P2, Q2, e):
    """
    Final change of basis in the response [Algorithm 4.8]
    Input:
    - E1, E2: curves
    - P1, Q1: points on E1
    - P2, Q2: points on E2
    - e: scaling factor for the order
    Output:
    - M_2: change of basis matrix on E2
    - B1, B2: basis on E1, E2
    """

    P1_1, Q1_1 = TorsionBasis(E1)
    P2_1, Q2_1 = TorsionBasis(E2)

    # M1 (2^e) * (P1_1, Q1_1) = (P1, Q1)
    M1 = ChangeOfBasis((P1_1, Q1_1), (P1, Q1), e=e)
    # We need M1 * (P1, Q1) = (2^e) * (P1_1, Q1_1)
    M1 = M1**-1
    a, b, c, d = M1.list()
    assert a*P1 + b*Q1 == 2**(params.f - e) * P1_1
    assert c*P1 + d*Q1 == 2**(params.f - e) * Q1_1

    # Apply it to (P2, Q2)
    a,b,c,d = M1.list()
    P2, Q2 = a*P2 + b*Q2, c*P2 + d*Q2

    # Last change of basis
    M2 = ChangeOfBasis((P2_1, Q2_1), (P2, Q2), e=e)
    a,b,c,d = M2.list()
    assert 2**(params.f-e) * (a*P2_1 + b*Q2_1) == P2
    assert 2**(params.f-e) * (c*P2_1 + d*Q2_1) == Q2

    return M2, (P1_1, Q1_1), (P2_1, Q2_1)

def EvalMatrix(M, basis=None):
    """
    Given a matrix M and two points (P, Q) evaluate M * (P, Q)^T
    Input:
    - M: a matrix
    - basis = (P, Q): two points; if None, (P0, Q0) is used
    Output:
    - M * (P, Q)^T
    """
    if basis:
        P, Q = basis
    else:
        P, Q = params.P0, params.Q0

    x1, x2, x3, x4 = M.list()
    return x1*P + x2*Q, x3*P + x4*Q

def EvalBasisPrecomp(theta, basis=None):
    """
    Given a quaternion theta in O0 and a basis (P, Q) of E[2^f] evaluate theta
    on the points using the precomputed action of a basis of O0 on the
    2-torsion.
    Input:
    - theta: an element of O0
    - basis: a basis (P, Q) of E[2^f]; if None, the standard (P0, Q0) is used
    Output:
    - theta(P), theta(Q)
    """
    # Express theta wrt <1, i, (i+j)/2, (1+k)/2>
    a, b, c, d = theta
    x, y, z, w = map(ZZ, (a-d, b-c, 2*c, 2*d))

    M = x*params.mat_1 + y*params.mat_i + z*params.mat_ij2 + w*params.mat_1k2
    return EvalMatrix(M, basis)

def TorsionBasis(E_A, e=None):
    """
    Given a montgomery curve E_A and an integer e <= f, computes a
    deterministic basis (R, S) of E_A[2^e] such that S is above (0, 0).
    Modified version of [Algorithm 2.1] working with (x,y)-points and no hint.
    Input:
    - A: montgomery coefficient
    - e: 2-order of the torsion; if not provided, f is used
    Output:
    - R, S: a basis of E_A[2^e] such that S is above (0, 0)
    """
    A = E_A.a_invariants()[1]
    assert E_A.a_invariants() == (0, A, 0, 1, 0), "not a Montgomery model"

    h = params.Fp(0)
    _i = params.Fp2_i
    if A.is_square():
        while True:
            h += 1
            if (1 + h**2).is_square(): # Norm of (1 + ih)
                continue
            xR = -A / (1 + _i*h) # Now this is non-square
            assert not xR.is_square()
            if E_A.is_x_coord(xR):
                break
    else:
        while True:
            h += 1
            xR = h*A
            if E_A.is_x_coord(xR):
                break

    xRS = -xR - A
    R = E_A.lift_x(xR)
    RS = E_A.lift_x(xRS)
    S = R - RS

    # Clear cofactors
    if not e:
        e = params.f

    cof = (params.p + 1) // 2**e
    R *= cof
    S *= cof

    assert (2**(e-1) * S).xy() == (0, 0)
    return R, S

def IdealToKernel(alpha, e):
    """
    Given a left O0-ideal I = <alpha, 2^e> for e<f, compute (a, b) such that
    kernel(phi_I) is generated by [a2^(f-e)]P0+[b2^(f-e)]Q0 [Algorithm 3.14]
    Input:
    - alpha: generator of I mod 2^e
    - e: integer such that n(I) = 2^e
    Output:
    - (a, b) defined the kernel of phi_I as described above
    """
    # Compute the action of alpha on (P0, Q0)
    a0, a1, a2, a3 = alpha
    x0, x1, x2, x3 = ZZ(a0-a3), ZZ(a1-a2), ZZ(a2*2), ZZ(a3*2)
    M_alpha = x0 * params.mat_1 + x1 * params.mat_i + \
            x2 * params.mat_ij2 + x3 * params.mat_1k2

    M_alpha = M_alpha.change_ring(Integers(2**e)).transpose()

    MK = M_alpha.right_kernel()

    if e == 1:
        # Right kernel here is a vector space
        a, b = [ZZ(i) for i in MK.basis()[0]]
        return a, b

    else:
        # Right kernel is no longer a vector space
        for row in M_alpha.right_kernel().generators_matrix():
            a, b = [ZZ(i) for i in row]
            if a % 2 != 0 or b % 2 != 0:
                # Found full order kernel
                return a, b

    raise ValueError('malformed ideal')

def ComputeEvenNonBacktrackingResponse(E, P, Q, alpha, e1, r):
    """
    Compute the even (non backtracking) part of the response, starting from the
    codomain E of the odd degree response [Algorithm 4.6]
    Input:
    - E: the current response curve
    - P, Q: image points under the odd part of the response (of order 2^e1)
    - alpha: the response quaternion element
    - e1: degree of the 2D iso
    - r: number of even steps
    Output:
    - E1: the actual response / challenge curve
    - P1, Q1: images of P, Q under the even part of the response
    """
    I = params.O0 * alpha + params.O0 * 2**r
    s, t = IdealToKernel(alpha, r)

    E._order = params.p + 1
    K = (2**e1 * s) * P + (2**e1 * t) * Q

    assert K.order() == 2**r
    K._order = ZZ(2**r)
    phi_I = E.isogeny(K)

    return phi_I.codomain(), phi_I(P), phi_I(Q)

def ComputeChallengeIsogeny(E, chl, P, Q, E1, P1, Q1, n):
    """
    Compute the challenge and update the response accordingly.
    Input:
    - E: public key curve
    - chl: challenge integer, so that ker(phi) = P + chl*Q
    - P, Q: a basis of the 2**f torsion on E
    - E1: codomain of the response isogeny
    - P1, Q1: a basis of the 2**f torsion on E1
    - n: number of backtracking steps
    Output:
    - E2: codomain of phi with ker [2^n](P + chl*Q), isomorphic to E1
    - P2, Q2: immages of P1, Q1 under the isomorphism
    """
    E._order = params.p + 1
    K = (2**n) * (P + chl*Q)
    K._order = ZZ(2**(params.f-n))

    phi2 = E.isogeny(K, model="montgomery")
    E2 = phi2.codomain()

    assert phi2.codomain().is_isomorphic(E1)

    eta = E1.isomorphism_to(phi2.codomain())

    return phi2.codomain(), eta(P1), eta(Q1)