From d40de259097c5e8d8fd35539560ca7c3d47523e7 Mon Sep 17 00:00:00 2001 From: Ryan Rueger Date: Sat, 1 Mar 2025 20:25:41 +0100 Subject: Initial Commit MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Damien Robert Co-Authored-By: Frederik Vercauteren Co-Authored-By: Jonathan Komada Eriksen Co-Authored-By: Pierrick Dartois Co-Authored-By: Riccardo Invernizzi Co-Authored-By: Ryan Rueger [0.01s] Co-Authored-By: Benjamin Wesolowski Co-Authored-By: Arthur Herlédan Le Merdy Co-Authored-By: Boris Fouotsa --- theta_lib/utilities/discrete_log.py | 273 +++++++++++++++++++++++ theta_lib/utilities/fast_sqrt.py | 55 +++++ theta_lib/utilities/order.py | 117 ++++++++++ theta_lib/utilities/strategy.py | 229 +++++++++++++++++++ theta_lib/utilities/supersingular.py | 413 +++++++++++++++++++++++++++++++++++ 5 files changed, 1087 insertions(+) create mode 100644 theta_lib/utilities/discrete_log.py create mode 100644 theta_lib/utilities/fast_sqrt.py create mode 100644 theta_lib/utilities/order.py create mode 100644 theta_lib/utilities/strategy.py create mode 100644 theta_lib/utilities/supersingular.py (limited to 'theta_lib/utilities') diff --git a/theta_lib/utilities/discrete_log.py b/theta_lib/utilities/discrete_log.py new file mode 100644 index 0000000..dff04dc --- /dev/null +++ b/theta_lib/utilities/discrete_log.py @@ -0,0 +1,273 @@ +# ===================================== # +# Fast DLP solving using Weil pairing # +# ===================================== # + +""" +This code has been taken from: +https://github.com/FESTA-PKE/FESTA-SageMath + +Copyright (c) 2023 Andrea Basso, Luciano Maino and Giacomo Pope. +""" + + + +# Sage imports +from sage.all import ( + ZZ +) + +# import pari for fast dlog +import cypari2 + +# Make instance of Pari +pari = cypari2.Pari() + +def discrete_log_pari(a, base, order): + """ + Wrapper around pari discrete log. Works like a.log(b), + but allows us to use the optional argument order. This + is important as we skip Pari attempting to factor the + full order of Fp^k, which is slow. + """ + x = pari.fflog(a, base, order) + return ZZ(x) + +def ell_discrete_log_pari(E,P,Base,order): + """ + Wrapper around pari elllog function to compute discrete + logarithms in elliptic curves over finite fields. + Works like P.log(Base). + """ + x = pari.elllog(E,P,Base,order) + return ZZ(x) + +def weil_pairing_pari(P, Q, D, check=False): + """ + Wrapper around Pari's implementation of the Weil pairing + Allows the check of whether P,Q are in E[D] to be optional + """ + if check: + nP, nQ = D*P, D*Q + if nP.is_zero() or nQ.is_zero(): + raise ValueError("points must both be n-torsion") + + return pari.ellweilpairing(P.curve(), P, Q, D) + +def tate_pairing_pari(P, Q, D): + """ + Wrapper around Pari's implementation of the Tate pairing + NOTE: this is not the reduced Tate pairing, so to make this + match with SageMath you need + + P.tate_pairing(Q, D, k) == pari.elltatepairing(E, P, Q, D)**((p^k - 1) / D) + """ + E = P.curve() + return pari.elltatepairing(E, P, Q, D) + +def _precompute_baby_steps(base, step, e): + """ + Helper function to compute the baby steps for + pohlig_hellman_base and windowed_pohlig_hellman. + """ + baby_steps = [base] + for _ in range(e): + base = base**step + baby_steps.append(base) + return baby_steps + +def pohlig_hellman_base(a, base, e): + """ + Solve the discrete log for a = base^x for + elements base,a of order 2^e using the + Pohlig-Hellman algorithm. + """ + baby_steps = _precompute_baby_steps(base, 2, e) + + dlog = 0 + exp = 2**(e-1) + + # Solve the discrete log mod 2, only 2 choices + # for each digit! + for i in range(e): + if a**exp != 1: + a /= baby_steps[i] + dlog += (2**i) + + if a == 1: + break + + exp //= 2 + + return dlog + +def windowed_pohlig_hellman(a, base, e, window): + """ + Solve the discrete log for a = base^x for + elements base,a of order 2^e using the + windowed Pohlig-Hellman algorithm following + https://ia.cr/2016/963. + + Algorithm runs recursively, computing in windows + l^wi for window=[w1, w2, w3, ...]. + + Runs the base case when window = [] + """ + # Base case when all windows have been used + if not window: + return pohlig_hellman_base(a, base, e) + + # Collect the next window + w, window = window[0], window[1:] + step = 2**w + + # When the window is not a divisor of e, we compute + # e mod w and solve for both the lower e - e mod w + # bits and then the e mod w bits at the end. + e_div_w, e_rem = divmod(e, w) + e_prime = e - e_rem + + # First force elements to have order e - e_rem + a_prime = a**(2**e_rem) + base_prime = base**(2**e_rem) + + # Compute base^(2^w*i) for i in (0, ..., e/w-1) + baby_steps = _precompute_baby_steps(base_prime, step, e_div_w) + + # Initialise some pieces + dlog = 0 + if e_prime: + exp = 2**(e_prime - w) + + # Work in subgroup of size 2^w + s = base_prime**(exp) + + # Windowed Pohlig-Hellman to recover dlog as + # alpha = sum l^(i*w) * alpha_i + for i in range(e_div_w): + # Solve the dlog in 2^w + ri = a_prime**exp + alpha_i = windowed_pohlig_hellman(ri, s, w, window) + + # Update a value and dlog computation + a_prime /= baby_steps[i]**(alpha_i) + dlog += alpha_i * step**i + + if a_prime == 1: + break + + exp //= step + + # TODO: + # I don't know if this is a nice way to do + # this last step... Works well enough but could + # be improved I imagine... + exp = 2**e_prime + if e_rem: + base_last = base**exp + a_last = a / base**dlog + dlog_last = pohlig_hellman_base(a_last, base_last, e_rem) + dlog += exp * dlog_last + + return dlog + +def BiDLP(R, P, Q, D, ePQ=None): + """ + Given a basis P,Q of E[D] finds + a,b such that R = [a]P + [b]Q. + + Uses the fact that + e([a]P + [b]Q, [c]P + [d]Q) = e(P,Q)^(ad-bc) + + Optional: include the pairing e(P,Q) which can be precomputed + which is helpful when running multiple BiDLP problems with P,Q + as input. This happens, for example, during compression. + """ + # e(P,Q) + if ePQ: + pair_PQ = ePQ + else: + pair_PQ = weil_pairing_pari(P, Q, D) + + # Write R = aP + bQ for unknown a,b + # e(R, Q) = e(P, Q)^a + pair_a = weil_pairing_pari(R, Q, D) + + # e(R,-P) = e(P, Q)^b + pair_b = weil_pairing_pari(R, -P, D) + + # Now solve the dlog in Fq + a = discrete_log_pari(pair_a, pair_PQ, D) + b = discrete_log_pari(pair_b, pair_PQ, D) + + return a, b + +def BiDLP_power_two(R, P, Q, e, window, ePQ=None): + r""" + Same as the above, but uses optimisations using that + D = 2^e. + + First, rather than use the Weil pairing, we can use + the Tate pairing which is approx 2x faster. + + Secondly, rather than solve the discrete log naively, + we use an optimised windowed Pohlig-Hellman. + + NOTE: this could be optimised further, following the + SIKE key-compression algorithms and for the Tate pairing + we could reuse Miller loops to save even more time, but + that seemed a little overkill for a SageMath PoC + + Finally, as the Tate pairing produces elements in \mu_n + we also have fast inversion from conjugation, but SageMath + has slow conjugation, so this doesn't help for now. + """ + p = R.curve().base_ring().characteristic() + D = 2**e + exp = (p**2 - 1) // D + + # e(P,Q) + if ePQ: + pair_PQ = ePQ + else: + pair_PQ = tate_pairing_pari(P, Q, D)**exp + + # Write R = aP + bQ for unknown a,b + # e(R, Q) = e(P, Q)^a + pair_a = tate_pairing_pari(Q, -R, D)**exp + + # e(R,-P) = e(P, Q)^b + pair_b = tate_pairing_pari(P, R, D)**exp + + # Now solve the dlog in Fq + a = windowed_pohlig_hellman(pair_a, pair_PQ, e, window) + b = windowed_pohlig_hellman(pair_b, pair_PQ, e, window) + + return a, b + +def DLP_power_two(R, P, Q, e, window, ePQ=None, first=True): + r""" + This is the same as BiDLP but it only returns either a or b + depending on whether first is true or false. + This is used in compression, where we only send 3 of the 4 + scalars from BiDLP + """ + p = R.curve().base_ring().characteristic() + D = 2**e + exp = (p**2 - 1) // D + + # e(P,Q) + if ePQ: + pair_PQ = ePQ + else: + pair_PQ = tate_pairing_pari(P, Q, D)**exp + + if first: + pair_a = tate_pairing_pari(Q, -R, D)**exp + x = windowed_pohlig_hellman(pair_a, pair_PQ, e, window) + + else: + pair_b = tate_pairing_pari(P, R, D)**exp + x = windowed_pohlig_hellman(pair_b, pair_PQ, e, window) + + return x + diff --git a/theta_lib/utilities/fast_sqrt.py b/theta_lib/utilities/fast_sqrt.py new file mode 100644 index 0000000..0609fe5 --- /dev/null +++ b/theta_lib/utilities/fast_sqrt.py @@ -0,0 +1,55 @@ + +# ============================================ # +# Fast square root and quadratic roots # +# ============================================ # + +""" +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. +""" + +def sqrt_Fp2(a): + """ + Efficiently computes the sqrt + of an element in Fp2 using that + we always have a prime p such that + p ≡ 3 mod 4. + """ + Fp2 = a.parent() + p = Fp2.characteristic() + i = Fp2.gen() # i = √-1 + + a1 = a ** ((p - 3) // 4) + x0 = a1 * a + alpha = a1 * x0 + + if alpha == -1: + x = i * x0 + else: + b = (1 + alpha) ** ((p - 1) // 2) + x = b * x0 + + return x + +def n_sqrt(a, n): + for _ in range(n): + a = sqrt_Fp2(a) + return a + +def sqrt_Fp(a): + """ + Efficiently computes the sqrt + of an element in Fp using that + we always have a prime p such that + p ≡ 3 mod 4. + + Copyright (c) Pierrick Dartois 2025. + """ + Fp = a.parent() + p = Fp.characteristic() + + return a**((p+1)//4) diff --git a/theta_lib/utilities/order.py b/theta_lib/utilities/order.py new file mode 100644 index 0000000..223f568 --- /dev/null +++ b/theta_lib/utilities/order.py @@ -0,0 +1,117 @@ +# ================================================== # +# Code to check whether a group element has order D # +# ================================================== # + +""" +This code has been taken from: +https://github.com/FESTA-PKE/FESTA-SageMath + +Copyright (c) 2023 Andrea Basso, Luciano Maino and Giacomo Pope. +""" + +from sage.all import( + ZZ, + prod, + cached_function +) + +def batch_cofactor_mul_generic(G_list, pis, group_action, lower, upper): + """ + Input: A list of elements `G_list`, such that + G is the first entry and the rest is empty + in the sublist G_list[lower:upper] + A list `pis` of primes p such that + their product is D + The `group_action` of the group + Indices lower and upper + Output: None + + NOTE: G_list is created in place + """ + + # check that indices are valid + if lower > upper: + raise ValueError(f"Wrong input to cofactor_multiples()") + + # last recursion step does not need any further splitting + if upper - lower == 1: + return + + # Split list in two parts, + # multiply to get new start points for the two sublists, + # and call the function recursively for both sublists. + mid = lower + (upper - lower + 1) // 2 + cl, cu = 1, 1 + for i in range(lower, mid): + cu = cu * pis[i] + for i in range(mid, upper): + cl = cl * pis[i] + # cl = prod(pis[lower:mid]) + # cu = prod(pis[mid:upper]) + + G_list[mid] = group_action(G_list[lower], cu) + G_list[lower] = group_action(G_list[lower], cl) + + batch_cofactor_mul_generic(G_list, pis, group_action, lower, mid) + batch_cofactor_mul_generic(G_list, pis, group_action, mid, upper) + + +@cached_function +def has_order_constants(D): + """ + Helper function, finds constants to + help with has_order_D + """ + D = ZZ(D) + pis = [p for p, _ in D.factor()] + D_radical = prod(pis) + Dtop = D // D_radical + return Dtop, pis + + +def has_order_D(G, D, multiplicative=False): + """ + Given an element G in a group, checks if the + element has order exactly D. This is much faster + than determining its order, and is enough for + many checks we need when computing the torsion + basis. + + We allow both additive and multiplicative groups + which means we can use this when computing the order + of points and elements in Fp^k when checking the + multiplicative order of the Weil pairing output + """ + # For the case when we work with elements of Fp^k + if multiplicative: + group_action = lambda a, k: a**k + is_identity = lambda a: a == 1 + identity = 1 + # For the case when we work with elements of E / Fp^k + else: + group_action = lambda a, k: k * a + is_identity = lambda a: a.is_zero() + identity = G.curve()(0) + + if is_identity(G): + return False + + D_top, pis = has_order_constants(D) + + # If G is the identity after clearing the top + # factors, we can abort early + Gtop = group_action(G, D_top) + if is_identity(Gtop): + return False + + G_list = [identity for _ in range(len(pis))] + G_list[0] = Gtop + + # Lastly we have to determine whether removing any prime + # factors of the order gives the identity of the group + if len(pis) > 1: + batch_cofactor_mul_generic(G_list, pis, group_action, 0, len(pis)) + if not all([not is_identity(G) for G in G_list]): + return False + + return True \ No newline at end of file diff --git a/theta_lib/utilities/strategy.py b/theta_lib/utilities/strategy.py new file mode 100644 index 0000000..550ef09 --- /dev/null +++ b/theta_lib/utilities/strategy.py @@ -0,0 +1,229 @@ +# ============================================================================ # +# Compute optimised strategy for 2-isogeny chains (in dimensions 2 and 4) # +# ============================================================================ # + +""" +The function optimised_strategy has been taken from: +https://github.com/FESTA-PKE/FESTA-SageMath + +Copyright (c) 2023 Andrea Basso, Luciano Maino and Giacomo Pope. + +Other functions are original work. +""" + +from sage.all import log, cached_function +import logging +logger = logging.getLogger(__name__) +#logger.setLevel("DEBUG") + +@cached_function +def optimised_strategy(n, mul_c=1): + """ + Algorithm 60: https://sike.org/files/SIDH-spec.pdf + Shown to be appropriate for (l,l)-chains in + https://ia.cr/2023/508 + + Note: the costs we consider are: + eval_c: the cost of one isogeny evaluation + mul_c: the cost of one element doubling + """ + + eval_c = 1.000 + mul_c = mul_c + + S = {1:[]} + C = {1:0 } + for i in range(2, n+1): + b, cost = min(((b, C[i-b] + C[b] + b*mul_c + (i-b)*eval_c) for b in range(1,i)), key=lambda t: t[1]) + S[i] = [b] + S[i-b] + S[b] + C[i] = cost + + return S[n] + +@cached_function +def optimised_strategy_with_first_eval(n,mul_c=1,first_eval_c=1): + r""" + Adapted from optimised_strategy when the fist isogeny evaluation is more costly. + This is well suited to gluing comptations. Computes optimal strategies with constraint + at the beginning. This takes into account the fact that doublings on the codomain of + the first isogeny are impossible (because of zero dual theta constants). + + CAUTION: When splittings are involved, do not use this function. Use + optimised_strategy_with_first_eval_and_splitting instead. + + INPUT: + - n: number of leaves of the strategy (length of the isogeny). + - mul_c: relative cost of one doubling compared to one generic 2-isogeny evaluation. + - first_eval_c: relative cost of an evaluation of the first 2-isogeny (gluing) + compared to one generic 2-isogeny evaluation. + + OUTPUT: + - S_left[n]: an optimal strategy of depth n with constraint at the beginning + represented as a sequence [s_0,...,s_{t-2}], where there is an index i for every + internal node of the strategy, where indices are ordered depth-first left-first + (as the way we move on the strategy) and s_i is the number of leaves to the right + of internal node i (see https://sike.org/files/SIDH-spec.pdf, pp. 16-17). + """ + + S_left, _, _, _ = optimised_strategies_with_first_eval(n,mul_c,first_eval_c) + + return S_left[n] + +@cached_function +def optimised_strategies_with_first_eval(n,mul_c=1,first_eval_c=1): + r""" + Adapted from optimised_strategy when the fist isogeny evaluation is more costly. + This is well suited to gluing comptations. Computes optimal strategies with constraint + at the beginning. This takes into account the fact that doublings on the codomain of + the first isogeny are impossible (because of zero dual theta constants). + + CAUTION: When splittings are involved, do not use this function. Use + optimised_strategy_with_first_eval_and_splitting instead. + + INPUT: + - n: number of leaves of the strategy (length of the isogeny). + - mul_c: relative cost of one doubling compared to one generic 2-isogeny evaluation. + - first_eval_c: relative cost of an evaluation of the first 2-isogeny (gluing) + compared to one generic 2-isogeny evaluation. + + OUTPUT: + - S_left: Optimal strategies "on the left"/with constraint at the beginning i.e. meeting the + first left edge that do not contain any left edge on the line y=sqrt(3)*(x-1). + - S_right: Optimal strategies "on the right" i.e. not meeting the fisrt left edge (no constraint). + """ + + # print(f"Strategy eval: n={n}, mul_c={mul_c}, first_eval_c={first_eval_c}") + + eval_c = 1.000 + first_eval_c = first_eval_c + mul_c = mul_c + + S_left = {1:[], 2:[1]} # Optimal strategies "on the left" i.e. meeting the first left edge + S_right = {1:[]} # Optimal strategies "on the right" i.e. not meeting the first left edge + C_left = {1:0, 2:mul_c+first_eval_c } # Cost of strategies on the left + C_right = {1:0 } # Cost of strategies on the right + for i in range(2, n+1): + # Optimisation on the right + b, cost = min(((b, C_right[i-b] + C_right[b] + b*mul_c + (i-b)*eval_c) for b in range(1,i)), key=lambda t: t[1]) + S_right[i] = [b] + S_right[i-b] + S_right[b] + C_right[i] = cost + + for i in range(3,n+1): + # Optimisation on the left (bd: + cost_k = C_middle[i-d] + C_middle[d] + d*mul_c + (i-d)*eval_c + if cost_k + 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] = 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 = 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 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 -- cgit v1.2.3-70-g09d2