diff options
Diffstat (limited to 'theta_lib/utilities')
| -rw-r--r-- | theta_lib/utilities/discrete_log.py | 273 | ||||
| -rw-r--r-- | theta_lib/utilities/fast_sqrt.py | 55 | ||||
| -rw-r--r-- | theta_lib/utilities/order.py | 117 | ||||
| -rw-r--r-- | theta_lib/utilities/strategy.py | 229 | ||||
| -rw-r--r-- | theta_lib/utilities/supersingular.py | 413 |
5 files changed, 1087 insertions, 0 deletions
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 (b<i-1 to avoid doublings on the codomain of the first isogeny) + b, cost = min(((b, C_left[i-b] + C_right[b] + b*mul_c + (i-1-b)*eval_c+first_eval_c) for b in range(1,i-1)), key=lambda t: t[1]) + S_left[i] = [b] + S_left[i-b] + S_right[b] + C_left[i] = cost + + return S_left, S_right, C_left, C_right + +@cached_function +def optimised_strategy_with_first_eval_and_splitting(n,m,mul_c=1,first_eval_c=1): + eval_c = 1.000 + first_eval_c = first_eval_c + mul_c = mul_c + + S_left, S_middle, C_left, C_middle = optimised_strategies_with_first_eval(n-m,mul_c,first_eval_c) + + ## Optimization of the right part (with constraint at the end only) + + # Dictionnary of dictionnaries of translated strategies "on the right". + # trans_S_right[d][i] is an optimal strategy of depth i + # without left edge on the line y=sqrt(3)*(x-(i-1-d)) + trans_S_right={} + trans_C_right={} + + for d in range(1,m+1): + trans_S_right[d]={1:[]} + trans_C_right[d]={1:0} + if d==1: + for i in range(3,n-m+d): + b, cost = min(((b, C_middle[i-b] + trans_C_right[1][b] + b*mul_c + (i-b)*eval_c) for b in [1]+list(range(3,i))), key=lambda t: t[1]) + trans_S_right[1][i] = [b] + S_middle[i-b] + trans_S_right[1][b] + trans_C_right[1][i] = cost + else: + for i in range(2,n-m+d): + if i!=d+1: + b = 1 + cost = trans_C_right[d-b][i-b] + C_middle[b] + b*mul_c + (i-b)*eval_c + for k in range(2,min(i,d)): + cost_k = trans_C_right[d-k][i-k] + C_middle[k] + k*mul_c + (i-k)*eval_c + if cost_k<cost: + b = k + cost = cost_k + # k=d + if i>d: + cost_k = C_middle[i-d] + C_middle[d] + d*mul_c + (i-d)*eval_c + if cost_k<cost: + b = d + cost = cost_k + for k in range(d+2,i): + #print(d,i,k) + cost_k = C_middle[i-k] + trans_C_right[d][k] + k*mul_c + (i-k)*eval_c + if cost_k<cost: + b = k + cost = cost_k + if b<d: + trans_S_right[d][i] = [b] + trans_S_right[d-b][i-b] + S_middle[b] + trans_C_right[d][i] = cost + else: + trans_S_right[d][i] = [b] + S_middle[i-b] + trans_S_right[d][b] + trans_C_right[d][i] = cost + + ## Optimization on the left (last part) taking into account the constraints at the beginning and at the end + for i in range(n-m+1,n+1): + d = i-(n-m) + b = 1 + cost = C_left[i-b] + trans_C_right[d][b] + b*mul_c + (i-1-b)*eval_c + first_eval_c + for k in range(2,i): + if k!=d+1 and k!=n-1: + cost_k = C_left[i-k] + trans_C_right[d][k] + k*mul_c + (i-1-k)*eval_c + first_eval_c + if cost_k<cost: + b = k + cost = cost_k + + S_left[i] = [b] + S_left[i-b] + trans_S_right[d][b] + C_left[i] = cost + + return S_left[n] + +@cached_function +def precompute_strategy_with_first_eval(e,m,M=1,S=0.8,I=100): + r""" + INPUT: + - e: isogeny chain length. + - m: length of the chain in dimension 2 before gluing in dimension 4. + - M: multiplication cost. + - S: squaring cost. + - I: inversion cost. + + OUTPUT: Optimal strategy to compute an isogeny chain without splitting of + length e with m steps in dimension 2 before gluing in dimension 4. + """ + n = e - m + eval_c = 4*(16*M+16*S) + mul_c = 4*(32*M+32*S)+(48*M+I)/log(n*1.0) + first_eval_c = 4*(2*I+(244+6*m)*M+(56+8*m)*S) + + return optimised_strategy_with_first_eval(n, mul_c = mul_c/eval_c, first_eval_c = first_eval_c/eval_c) + +@cached_function +def precompute_strategy_with_first_eval_and_splitting(e,m,M=1,S=0.8,I=100): + r""" + INPUT: + - e: isogeny chain length. + - m: length of the chain in dimension 2 before gluing in dimension 4. + - M: multiplication cost. + - S: squaring cost. + - I: inversion cost. + + OUTPUT: Optimal strategy to compute an isogeny chain of length e + with m steps in dimension 2 before gluing in dimension 4 and + with splitting m steps before the end. + """ + logger.debug(f"Strategy eval split: e={e}, m={m}") + n = e - m + eval_c = 4*(16*M+16*S) + mul_c = 4*(32*M+32*S)+(48*M+I)/log(n*1.0) + first_eval_c = 4*(2*I+(244+6*m)*M+(56+8*m)*S) + + return optimised_strategy_with_first_eval_and_splitting(n, m, mul_c = mul_c/eval_c, first_eval_c = first_eval_c/eval_c) diff --git a/theta_lib/utilities/supersingular.py b/theta_lib/utilities/supersingular.py new file mode 100644 index 0000000..e760c1b --- /dev/null +++ b/theta_lib/utilities/supersingular.py @@ -0,0 +1,413 @@ +""" +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 |