compute-orders-deterministic-eisenstein.py

from sage.all import random_prime, GF, EllipticCurve, GF, is_prime
import math
ZETA = lambda n,x,p: pow(x, ((p-1)//n), p)
MODSQRT = lambda n, p: pow(n % p, (p + 1) // 4, p)

def cornacchia(d, p):
    """
    Standard Cornacchia algorithm for x^2 + d*y^2 = p
    """
    assert p % 12 == 7
    if d <= 0 or d >= p:
        raise ValueError("invalid input")
    # Check if -d is a quadratic residue mod p
    if pow((-d) % p, (p - 1) // 2, p) != 1:
        raise ValueError("no solution")
    # Find sqrt(-d) mod p
    assert p%4==3
    x0 = MODSQRT(-d, p)
    # Choose the larger square root
    if x0 < p // 2:
        x0 = p - x0

    # Extended Euclidean algorithm
    a, b = p, x0
    limit = int(p**0.5)
    while b > limit:
        a, b = b, a % b
        assert a > 0
        assert b > 0

    # Check if we can complete the solution
    remainder = p - b * b
    if remainder % d != 0:
        raise ValueError("no solution 1")

    c = remainder // d
    t = math.isqrt(c) # Python 3.8+
    if t * t != c:
        raise ValueError("no solution 2")

    return b, t

def sample_primes(bitlen):
    while True:
        p = random_prime(2**bitlen, lbound=2**(bitlen-1))
        if p % 12 == 7:
            yield p

def make_norms_cd(c,d,p):
    return [p + 1 + (x_0*c) + (x_1*d) for x_0, x_1 in [
        (-2,1), (1,-2), (-1,-1), (2,-1), (-1,2), (1,1)
    ]]

ALL_SUB_PATTERNS = [
    (0, 2, 1, 3, 5, 4),
    (0, 4, 5, 3, 1, 2),
    (3, 1, 2, 0, 4, 5),
    (3, 5, 4, 0, 2, 1),
]

def classify_triple_xor(arr):
    # Given four binary flags, return 0..3
    # found using machine learning, scikit, xgvboost etc.
    return arr[1] + ((arr[0] ^ arr[2] ^ arr[3]) * (3 - 2*arr[1]))

def calculate_curve_orders(p):
    assert p % 12 == 7
    F = GF(p)
    g = F.multiplicative_generator()
    a, b = cornacchia(3,p)
    c, d = a + b, 2 * b
    assert c**2 - c*d + d**2 == p
    assert p + 1 + a - (3*b) == p + 1 + c - (2*d)
    seq = [
        int(ZETA(2,d,p)==1),  # zeta_2(c) is anti-correlated,
        int((ZETA(3,g,p) * c + d) == 0),
        int(ZETA(2,a,p) == 1),
        int(ZETA(2,b,p) == 1),
    ]
    idx = classify_triple_xor(seq)
    result = [0] * 6
    norms_cd = make_norms_cd(c,d,p)
    for i,j in enumerate(ALL_SUB_PATTERNS[idx]):
        result[i] = norms_cd[j]
    return result

def find_generator(g,p):
    g = int(g)
    x = 1
    while True:
        yy = (pow(x,3,p) + g) % p
        y = MODSQRT(yy, p)
        if (y*y) == yy:
            if y & 1:
                y = p - y
            E = EllipticCurve(GF(p), [0, g])
            if E.point((x,y)).order() == E.order():
                return x,y
        x += 1

def example_secp256k1():
    p = 115792089237316195423570985008687907853269984665640564039457584007908834671663
    F = GF(p)
    g = F.multiplicative_generator()
    q = 115792089237316195423570985008687907852837564279074904382605163141518161494337
    G = find_generator(g, p)
    orders = calculate_curve_orders(p)
    assert orders.index(q) == 5
    print('secpk256k1')
    print('         p', p, hex(p))
    print('         g', g)
    print('         G', G)
    print('    orders')
    for i, n in enumerate(orders):
        print(f'         {i}', n, hex(n), 'prime' if is_prime(n) else '')
    print('       idx', orders.index(q))
    print()

def main():
    example_secp256k1()
    print('Running Tests', end='')
    for bitlen in range(32,256):
        i = 0
        for p in sample_primes(bitlen):
            i += 1
            if i > 50:
                break
            F = GF(p)
            g = F.multiplicative_generator()
            curves = [EllipticCurve(F, [0, g**i]) for i in range(6)]
            orders = [_.order() for _ in curves]
            predicted_orders = calculate_curve_orders(p)
            assert predicted_orders == orders
        print('.', end='', flush=True)
    print('TESTS OK ^_^')

if __name__ == "__main__":
    main()

mapping-lemma.md

Deterministic Curve Order Mapping via Eisenstein Integers in $\mathbb{Q}(\sqrt{-3})$

See: lemma/2-eisenstein-mapping-magic.py

The history of elliptic curve point counting has rich and deep roots which have been refined over the decades, we can trace a lineage of sorts which leads to the most useful results for us:

• In 1986: Lenstra [@lenstra1986elliptic, 11, §4] provides an intuitive formula for j-invariant 0 curves over $\mathbb{Q}(\sqrt{-3})$, which comes surprisingly close to Wu & Xu's [@GuangwuXu_2020_1136] conclusions in 2020.
• In 1995: Schoof [@schoof1995counting§4] expains how to count the number of points when the endomorphism ring of E is known.
• In 2005: Nogami and Morikawa [@Nogami2005] propose a method to obtain the six orders of these curves by counting the order of only one curve.
• In 2006: Hess, Smart, and Vercauteren [@Hess_Smart_Vercauteren_2006_110] propose similar methods for twists $\phi_d : E' \mapsto E$ over extension fields $\mathbb{F}_{q^d}$ where $d = 6$ if $j(E) = 0$.
• In 2018: Kim, Bahr, Neyman and Taylor [@kim2018orders§2.1] then completely characterize, by j-invariant, the number of orders of elliptic curves over all finite fields $F_{p^r}$ using combinatorial arguments and elementary number theory.

We present a simple, explicit and computationally efficient method for a deterministic mapping between the isomorphisms $E_{p,j} : y^2=x^3+g^j$ and their orders using properties of Eisenstein integers and primitive characters:

• Any prime $p \equiv 7 \pmod{12}$ can be represented as $p = c^2 - cd + d^2 = a^2 + 3b^2$.

• This corresponds to a prime ideal factorization $(p) = \pi\bar{\pi}$ in the ring of Eisenstein integers $\mathbb{Z}[\omega]$, where $\pi = c + d\omega$ and $\omega = e^{2\pi i/3}$.

• The orders of the six curves $E_j: y^2 = x^3 + g^j$ correspond to the norms of the six elements $\pi \pm{{1,\omega,\omega^2}} \in \mathbb{Z}[\omega]$, where $p + 1 + \text{Trace}(u) = N(\pi + u) = \#E_j$, which we define as traces =

  • traces[0] = $\text{Trace}(1) = - 2a = - 2c + d$
  • traces[0] = $\text{Trace}(\omega^2) = - a - 3b = - c - d$
  • traces[0] = $\text{Trace}(\omega) = a - 3b = c - 2d$
  • traces[0] = $\text{Trace}(- 1) = 2a = 2c - d$
  • traces[0] = $\text{Trace}(- \omega^2) = a + 3b = c + d$
  • traces[0] = $\text{Trace}(- \omega) = -a + 3b = - c + 2d$

  A lookup table for the coefficients $(t_0 c) + (t_1 d)$ of traces[i] can be stated simply as:

        t = [(-2,1), (-1,-1), (1,-2), (2,-1), (1,1), (-1,2)][i%6]

We have empirically verified that the set of curve orders matches the set of traces. However, the specific mapping between the generator index $j$ of $g^j$ and the index $i$ of traces[i] depends on the chirality and orientation w.r.t the Eisenstein lattice. We solve this mapping in a novel way:

• Find the canonical $(a,b)$ to represent $p = a^2 + 3b^2$ (e.g. using Cornacchia's algorithm)
  • Where $a$ is even, $b$ is odd, and both are positive integers
• Derive $c = a + b$ and $d = 2b$

We can then utilize a lookup table, indexed by:

• $f: {0,1}^2 \to {0,1,2,3}$
  • $f(u_0,u_1) \to 2 u_0 + u_1$
  • $u_0: \zeta_2(a b d) = 1$
  • $u_1: \zeta_3(g) c + d = 0$

When $n | (p-1)$, $\zeta_n(x)$ denotes $x$'s character in the primitive $n$-th roots of unity via $x^{(p-1)/n} \bmod p$. The resulting $f$ provides a bijection between the traces and the power of the generator $g$:

• $f(u_0,u_1) = 0$: $\text{traces}[i] + p + 1 \leftrightarrow \text{\#}E_j : y^2 = x^3 + g^{(0, 1, 2, 3, 4, 5)[i]}$
• $f(u_0,u_1) = 1$: $\text{traces}[i] + p + 1 \leftrightarrow \text{\#}E_j : y^2 = x^3 + g^{(0, 5, 4, 3, 2, 1)[i]}$
• $f(u_0,u_1) = 2$: $\text{traces}[i] + p + 1 \leftrightarrow \text{\#}E_j : y^2 = x^3 + g^{(3, 4, 5, 0, 1, 2)[i]}$
• $f(u_0,u_1) = 3$: $\text{traces}[i] + p + 1 \leftrightarrow \text{\#}E_j : y^2 = x^3 + g^{(3, 2, 1, 0, 5, 4)[i]}$

While the specific order of our traces in relation to the mapping indices are somewhat arbitrary, they can be thought of as permutations of the Dirichlet characters $\chi_9$ or $\chi_7$. This method provides a complete and deterministic characterization of the distinct six isomorphism classes of j-invariant 0 curves over $\mathrm{GF}(p)$ when $p \equiv 7 \pmod{12}$.

This aformentioned method requires 3 modulo exponentiations, which includes finding $\sqrt{-3}$ for Cornacchia's algorithm, and two integer square roots. The remaining operations are either table lookups, or simple binary and integer arithmetic without negative intermediates. This computes all six curve orders simultaneously, compared to approaches requiring separate computations for each curve.

safe-prime-glv-compatible-curves.py

from sage.all import random_prime, is_prime, GF, EllipticCurve, factor
from sage.rings.factorint import factor_trial_division
from math import gcd, isqrt

ZETA = lambda n,x,p: pow(x%p, ((p-1)//n), p)
MODSQRT = lambda n, p: pow(n%p, (p + 1) // 4, p)

MAX_C = 2**31
bitlen = 64

def check_glv_endomorphism2(curve, p, n, g_base, generator=None):
    # Use provided generator or get default
    if generator is None:
        try:
            generator = curve.gens()[0]
        except:
            return {"supports_glv": False, "reason": "No generator specified and unable to find default generator"}

    """
    # Find a cube root of unity in the base field
    beta = ZETA(3, g_base, p) #g_base**((p-1)//3)
    # Verify beta is a non-trivial cube root of unity
    if beta == 1 or beta**3 != 1:
        assert False
    # Verify beta satisfies the minimal polynomial
    if beta**2 + beta + 1 != 0:
        assert False
    betas = [beta,beta**2]
    """

    beta_val_6th = ZETA(6, g_base, p) # g_base**((n-1)//6)
    beta_val_3rd = ZETA(3, g_base, p) # g_base**((n-1)//3)
    beta_vals = [pow(beta_val_3rd,i,p) for i in range(1, 3)] # + [pow(beta_val_6th, i, p) for i in range(1, 6)]

    Fq = GF(n)
    Fq_star = Fq.unit_group()
    g_scalar = Fq(Fq_star.gen(0))
    lambda_val_6th = g_scalar**((n-1)//6)
    lambda_val_3rd = g_scalar**((n-1)//3)
    lambda_vals = [lambda_val_3rd**i for i in range(1, 3)] # [lambda_val_6th**i for i in range(1, 6)] +
    #print('lambda 6', lambda_val_6th, [lambda_val_6th**i for i in range(1, 6)])
    #print('lambda 3', lambda_val_3rd, [lambda_val_3rd**i for i in range(1, 3)])
    # Test both values to see if either works
    #lambda_vals = [lambda_val_6th, lambda_val_3rd]

    for beta_i, beta in enumerate(beta_vals):
        try:
            for lambda_i,lambda_val in enumerate(lambda_vals):
                # Verify lambda^6 = 1
                if lambda_val**6 != 1:
                    continue

                # Check if lambda satisfies the minimal polynomial
                if lambda_val**2 + lambda_val + 1 != 0:
                    continue

                # The critical test: check if lambda*G = (beta*G_x, G_y)
                try:
                    endo_point = curve(beta * generator[0], generator[1])
                    scalar_point = lambda_val * generator
                    if endo_point == scalar_point:
                        return {
                            "supports_glv": True,
                            "beta": beta,
                            "lambda": lambda_val,
                            "verification": "lambda*G == E(beta*G[0], G[1])",
                            "lambda_i": lambda_i,
                            "beta_i": beta_i
                        }
                except Exception:
                    continue

        except Exception as e:
            return {"supports_glv": False, "reason": f"Error finding lambda value: {e}"}

    return {
        "supports_glv": False,
        "reason": "No suitable lambda value found that satisfies the endomorphism relation",
        "beta_candidates": beta_vals,
        "lambda_candidates": lambda_vals
    }


def cornacchia(d, p):
    """Standard Cornacchia algorithm for x^2 + d*y^2 = p"""
    assert p % 12 == 7
    if d <= 0 or d >= p:
        raise ValueError("invalid input")
    if ZETA(2, -d, p) != 1: raise ValueError("no solution")
    x0 = MODSQRT(-d, p)
    # Choose the larger square root
    if x0 < p // 2:
        x0 = p - x0
    # Extended Euclidean algorithm
    a, b = p, x0
    limit = isqrt(p) # Python 3.8+
    while b > limit:
        a, b = b, a % b
        assert a > 0 and b > 0 # guaranteed positive integers
    remainder = p - b * b
    assert remainder % d == 0  # guaranteed congruence
    c = remainder // d
    t = isqrt(c) # Python 3.8+
    assert t * t == c  # guaranteed exact squares
    return b, t

def sample_primes(bitlen):
    while True:
        p = random_prime(2**bitlen, lbound=2**(bitlen-1))
        if p % 12 == 7:
            yield p

def make_terms_cd(c,d):
    return [(x_0*c) + (x_1*d) for x_0, x_1 in [
        (-2,1), (-1,-1), (1,-2), (2,-1), (1,1), (-1,2)
    ]]

def make_norms_cd(c,d,p):
    return [p + 1 + _ for _ in make_terms_cd(c,d)]

ALL_SUB_PATTERNS = [
    (0, 1, 2, 3, 4, 5),
    (0, 5, 4, 3, 2, 1),
    (3, 4, 5, 0, 1, 2),
    (3, 2, 1, 0, 5, 4),
]

def calculate_curve_orders(p, g):
    assert p % 12 == 7
    #F = GF(p)
    #g = F.multiplicative_generator()
    a, b = cornacchia(3,p)
    assert a%2 == 0 and b%2 == 1
    c, d = a + b, 2 * b
    assert c**2 - c*d + d**2 == p
    assert p + 1 + a - (3*b) == p + 1 + c - (2*d)
    u0 = ZETA(2,a*b*d,p)
    u1 = int((ZETA(3,g,p) * c + d) == 0)
    idx = (int(u0 == 1) * 2) + u1
    result = [0] * 6
    norms_cd = make_norms_cd(c,d,p)
    for i,j in enumerate(ALL_SUB_PATTERNS[idx]):
        result[i] = norms_cd[j]
    return result

def analyze_large_prime_compatibility(k, mod, target_residue, min_log2_q=250):
    """Analyze if k can produce large primes, not just any primes"""
    target_kq_mod = (target_residue - 1) % mod
    # Find all q residues that work
    valid_q_mods = []
    for q_mod in range(mod):
        if (k * q_mod) % mod == target_kq_mod:
            valid_q_mods.append(q_mod)
    # Analyze each residue class for large prime compatibility
    large_prime_compatible = []
    small_prime_only = []
    for q_mod in valid_q_mods:
        if q_mod == 0:
            continue  # q=0 is not prime
        elif q_mod == 1:
            large_prime_compatible.append(q_mod)  # Many large primes ≡ 1 (mod m)
        elif gcd(q_mod, mod) == 1:
            large_prime_compatible.append(q_mod)  # Coprime to mod, can contain large primes
        else:
            # q_mod shares a factor with mod
            # Check if this residue class contains any prime > mod
            #g = gcd(q_mod, mod) #?
            if q_mod <= mod and q_mod in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31]:
                # This residue class only contains one small prime
                if 2**min_log2_q > mod:  # We need much larger primes
                    small_prime_only.append(q_mod)
                else:
                    large_prime_compatible.append(q_mod)
            else:
                small_prime_only.append(q_mod)
    if not large_prime_compatible:
        return False

    # Calculate effective density for large primes
    # This is approximate - we're estimating based on coprimality
    total_coprime_residues = sum(1 for r in range(mod) if gcd(r, mod) == 1)
    large_prime_density = len(large_prime_compatible) / total_coprime_residues

    if large_prime_density == 1:
        print(f"\nAnalyzing k={k} for p ≡ {target_residue} (mod {mod})")
        print(f"Need: k×q ≡ {target_kq_mod} (mod {mod})")
        print(f"All solutions for q (mod {mod}): {valid_q_mods}")
        print(f"Large-prime-compatible q residues: {large_prime_compatible}")
        print(f"Small-prime-only q residues: {small_prime_only}")
        print(f"Large prime density: {len(large_prime_compatible)}/{total_coprime_residues} = {large_prime_density:.3f}")
        return True

def get_prime_range(k, min_log2_q):
    p_min = (2**bitlen) - (2**32) - MAX_C
    p_max = (2**bitlen) - (2**32)
    q_min = max((p_min - 1) // k, int(2**min_log2_q))
    q_max = (p_max - 1) // k
    if q_min < q_max:
        return p_min, p_max, q_min, q_max

def test_k_practically(k, mod, target_residue, min_log2_q):
    """Actually test a k value to see if it finds primes quickly"""
    target_kq_mod = (target_residue - 1) % mod
    possible_q_mods = []
    for q_mod in range(mod):
        if (k * q_mod) % mod == target_kq_mod:
            possible_q_mods.append(q_mod)
    p_min, p_max,q_min,q_max = get_prime_range(k, min_log2_q)
    assert q_min < q_max

    while True:
        q = random_prime(q_max, False, lbound=q_min)

        if q % mod not in possible_q_mods:
            continue

        p = k * q + 1

        if not (p_min <= p <= p_max) or p % mod != target_residue:
            continue

        if is_prime(p):
            #print("\tfound\tq", hex(q), log2(q))
            #print("\t\tp", hex(p), log2(p))
            #print("\t\tp-1 factors", list(factor(p-1)))
            yield p,q

from collections import defaultdict

def main():
    viable_k_values = defaultdict(set)
    min_log2_q = bitlen-8
    for need_mod, need_residue in [(12,7),(4,3)]:
        for k in range(1,100):
            range_ok = get_prime_range(k, min_log2_q) is not None
            if not range_ok:
                break
            if analyze_large_prime_compatibility(k, need_mod, need_residue):
                viable_k_values[(need_mod,need_residue)].add(k)
    print(viable_k_values)

    ok = defaultdict(int)
    fail = defaultdict(int)
    dorp = 0
    for curve_order_k in viable_k_values[(12,7)]:
        # First, find safe primes for the order of our curve
        for curve_order, curve_order_bigfactor in test_k_practically(curve_order_k, 12, 7, min_log2_q):
            #assert (curve_order_k*curve_order_bigfactor)+1 == curve_order
            uprime,vprime = cornacchia(3, curve_order)
            u, v = 2*uprime, 2*vprime
            possible_values = [
                (v + 1, (u + (v + 1) + 1) // 2),
                (v - 1, (u + (v - 1) - 1) // 2),
                (v + 1, (u + v) // 2),
                (v - 1, (u + v) // 2),
                (v + 1, (u + (v + 1)) // 2),
                (v - 1, (u + (v - 1)) // 2),
            ]
            for i,(d,c) in enumerate(possible_values):
                p = c**2 - c*d + d**2
                if p%12==7 and is_prime(p):
                    # We only care if p-1 has small factors below 2**16, or one big factor!
                    p_minus_one_factors = factor_trial_division(p-1, 2**16)
                    if not all([is_prime(_[0]**[1][0]) for _ in p_minus_one_factors]):
                        continue
                    blah = make_norms_cd(c,d,p)
                    assert curve_order in blah
                    print('i',i, 'c', c, 'd', d, 'p', p, 'is_prime(p)', is_prime(p), 'curve valid?', curve_order in blah)
                    F = GF(p)
                    g_base = F.multiplicative_generator()
                    new_orders = calculate_curve_orders(p, g_base)
                    assert curve_order in new_orders
                    print(' factors', factor(p-1), ",", factor(curve_order-1))
                    E = EllipticCurve(F, [0, g_base])
                    # We get about 50/50% chance of it supporting the GLV endomorphism
                    # Either both work or neither works
                    glv_result = check_glv_endomorphism2(E, p, curve_order, g_base)
                    if glv_result['supports_glv']:
                        ok[i] += 1
                    else:
                        fail[i] += 1
                    dorp += 1
                    print(' GLV? 1', glv_result)
                    print(' GLV? 5', check_glv_endomorphism2(E, p, curve_order, g_base**5))
                    print()
                    if dorp % 20 == 0:
                        for morp in ok.keys():
                            print(" c d pair", morp, 'success', ok[morp], 'fail', fail[morp])

if __name__ == "__main__":
    main()

wip-glv.py

#!/usr/bin/env python3
import json
from math import log2
from random import randint
from time import perf_counter
from collections import defaultdict
from contextlib import contextmanager
from sage.all import random_prime, is_prime, GF, EllipticCurve, sqrt as sage_sqrt, primitive_root

STATS = defaultdict(lambda: defaultdict(int))
COUNT = 25
BITLENS = range(30,150)

@contextmanager
def measure(key, action):
    start = perf_counter()
    try:
        yield
    finally:
        STATS[key][action] += perf_counter() - start

def check_glv_endomorphism(curve, generator=None):
    """
    Check if an elliptic curve supports the specific endomorphism needed for GLV method.

    Parameters:
        curve: An elliptic curve defined over a finite field
        generator: Optional generator point. If None, will use the default generator

    Returns:
        A dictionary containing the test results and relevant values
    """
    # Get the curve parameters
    F = curve.base_field()
    p = F.order()
    n = curve.order()

    # Check if curve has the correct form for j-invariant 0
    if curve.a_invariants() != (0, 0, 0, 0, curve.a_invariants()[4]):
        return {"supports_glv": False, "reason": "Curve is not in the form y^2 = x^3 + b"}

    # Check field characteristics
    if p % 3 != 1:
        return {"supports_glv": False, "reason": "Field characteristic p must satisfy p ≡ 1 (mod 3)"}

    # Check if n is prime
    if not n.is_prime():
        return {"supports_glv": False, "reason": "Curve order must be prime"}

    # Use provided generator or get default
    if generator is None:
        try:
            generator = curve.gens()[0]
        except:
            return {"supports_glv": False, "reason": "No generator specified and unable to find default generator"}

    # Find a cube root of unity in the base field
    try:
        # primitive_root() is usually faster than via F.unit_group(), e.g. F_star = F.unit_group(), g_base = F(F_star.gen(0))
        g_base = F(primitive_root(F.order()))
        beta = g_base**((p-1)//3)

        # Verify beta is a non-trivial cube root of unity
        if beta == 1 or beta**3 != 1:
            return {"supports_glv": False, "reason": "Failed to find valid cube root of unity in base field"}

        # Verify beta satisfies the minimal polynomial
        if beta**2 + beta + 1 != 0:
            return {"supports_glv": False, "reason": "Cube root of unity does not satisfy minimal polynomial"}
    except Exception as e:
        return {"supports_glv": False, "reason": f"Error finding cube root of unity: {e}"}

    # Find primitive 6th root of unity in the scalar field
    try:
        Fq = GF(n)
        g_scalar = Fq(primitive_root(n))

        # Try both (n-1)/6 and (n-1)/3 for lambda
        lambda_val_1 = g_scalar**((n-1)//6)
        lambda_val_2 = g_scalar**((n-1)//3)

        # Test both values to see if either works
        lambda_vals = [lambda_val_1, lambda_val_2]
        working_lambda = None

        for lambda_val in lambda_vals:
            # Verify lambda^6 = 1
            if lambda_val**6 != 1:
                continue

            # Check if lambda satisfies the minimal polynomial
            if lambda_val**2 + lambda_val + 1 != 0:
                continue

            # The critical test: check if lambda*G = (beta*G_x, G_y)
            try:
                endo_point = curve(beta * generator[0], generator[1])
                scalar_point = lambda_val * generator

                if endo_point == scalar_point:
                    working_lambda = lambda_val
                    break
            except Exception:
                continue

        if working_lambda is None:
            return {
                "supports_glv": False,
                "reason": "No suitable lambda value found that satisfies the endomorphism relation",
                "beta": beta,
                "lambda_candidates": lambda_vals,
                "test_points": {
                    "endomorphism": endo_point,
                    "scalar_mult": scalar_point
                }
            }

    except Exception as e:
        return {"supports_glv": False, "reason": f"Error finding lambda value: {e}"}

    # If we get here, the curve supports GLV
    return {
        "supports_glv": True,
        "beta": beta,
        "lambda": working_lambda,
        "verification": "lambda*G == E(beta*G[0], G[1])"
    }

def find_generator(E:EllipticCurve):
    x = E.base_field()(1)
    while True:
        try:
            G = E.lift_x(x)
        except ValueError:
            x += 1
            continue
        if G.order() == E.order():
            # Generator must be even?
            if int(G[1]) & 1:
                G = -G
            return G
        x += 1

def find_split_constants_explicit_tof(p:int, E:EllipticCurve):
    # from libsecp256k1/sage
    assert p % 3 == 1
    assert E.j_invariant() == 0
    t = int(E.trace_of_frobenius())
    c = int(sage_sqrt((4*p - t**2)//3))
    b1, b2 = c, (1 - (t - c)//2) % E.order()
    return b1, int(b2)

def decompose_scalar(k, n, g1, g2, b1, b2, lambda_val, bitsize2=384):
    """Decompose scalar k into k1 and k2 using GLV method."""
    k = k % n

    # Compute c1 and c2 using precomputed constants and bit operations
    # This replaces division with multiplication and right shift
    c1 = (k * g1) >> bitsize2
    c2 = (k * g2) >> bitsize2

    # Handle final bit for rounding as per the paper
    if (k * g1) & (1 << (bitsize2-1)):
        c1 += 1
    if (k * g2) & (1 << (bitsize2-1)):
        c2 += 1

    neg_b1 = -b1 % n
    neg_b2 = -b2 % n

    # Compute k2 = -c1·(-b1) - c2·(-b2)
    k2 = (c1 * neg_b1 + c2 * neg_b2) % n

    # Compute k1 = k - k2·λ mod n
    k1 = (k - (k2 * lambda_val) % n) % n

    # Ensure k1 and k2 are properly minimized
    if k1 > (n >> 1):
        k1 = k1 - n
    if k2 > (n >> 1):
        k2 = k2 - n

    return k1, k2

def find_g1_g2(b1, b2, n, bitsize=256, bitsize2=384):
    # Python's round() is off by 1
    quotient_g1 = (2**bitsize2)*(-b2)//n
    remainder_g1 = (2**bitsize2)*(-b2)%n
    g1 = (quotient_g1 + (1 if remainder_g1 >= n//2 else 0)) % (2**bitsize)
    # Python's round() is off by one
    quotient_g2 = (2**bitsize2)*(b1)//n
    remainder_g2 = (2**bitsize2)*(b1)%n
    g2 = (quotient_g2 + (1 if remainder_g2 >= n//2 else 0)) % (2**bitsize)
    return g1, g2

def sample_primes(bitlen):
    while True:
        with measure(bitlen, '1_time_random_prime'):
            p = random_prime(2**bitlen, lbound=2**(bitlen-1))
        if p % 12 == 7:
            yield p
            STATS[bitlen]['1_prime_congruent'] += 1
            continue
        STATS[bitlen]['1_prime_not_congruent'] += 1

def sample_curves(bitlen):
    for p in sample_primes(bitlen):
        with measure(bitlen, '2_time_primitive_root'):
            F = GF(p)
            g = F(primitive_root(p))
        n_found = 0
        # 6th primitive roots g^1 and g^-1 are only ones which can be prime
        for b in [int(F(g)), int(F(g**5))]:
            with measure(bitlen, '2_time_check_curve'):
                E = EllipticCurve(F, [0, b])
                E_q = E.order()
                E_q_prime = is_prime(E_q)
            if E_q_prime:
                n_found += 1
                yield p, b, E
        if n_found == 0:
            STATS[bitlen]['2_no_prime_curve'] += 1
        else:
            STATS[bitlen]['2_prime_curve'] += 1
            if n_found == 2:
                STATS[bitlen]['2_1_prime_curve_both'] += 1
            else:
                STATS[bitlen]['2_1_prime_curve_one'] += 1

def try_decompose(E_q, bitlen, b1, b2, lambda_val):
    g1, g2 = find_g1_g2(b1, b2, E_q, bitlen, bitlen//2 + bitlen)
    count = 10
    k1_log2_total = 0
    k2_log2_total = 0
    for _ in range(count):
        k = randint(2**510, 2**512) % E_q
        k1, k2 = decompose_scalar(k, E_q, g1, g2, b1, b2, lambda_val, (bitlen//2)+bitlen)
        recomposed = (k1 + (k2 * lambda_val) % E_q) % E_q
        if k != recomposed:
            return None
        k1_log2 = log2(abs(int(k1))) if k1 != 0 else 0
        k2_log2 = log2(abs(int(k2))) if k2 != 0 else 0
        k1_log2_total += k1_log2
        k2_log2_total += k2_log2
    k1_log2 = round(k1_log2_total / count,2)
    k2_log2 = round(k2_log2_total / count,2)
    return k1_log2, k2_log2, k1_log2 <= bitlen//2 and k2_log2 <= bitlen//2, int(lambda_val) & 1

with open('graphs/3-glv-decomp.jsonl', 'w') as handle:
    for bitlen in BITLENS:
        FOUND = 0
        bitlen_start = perf_counter()
        with measure(bitlen, 'time_total'):
            for p, b, E in sample_curves(bitlen):
                if FOUND >= COUNT:
                    break
                with measure(bitlen, '2_time_basic_field_checks'):
                    b_log2 = log2(b)
                    E_q = E.order()
                    E_q_prime = is_prime(E_q)

                STATS[bitlen]['3_found_curves'] += 1

                with measure(bitlen, '3_time_find_generator'):
                    G = find_generator(E)
                if log2(G[1]) <= 8:
                    STATS[bitlen]['3_1_G_y_is_under_8bit'] += 1
                else:
                    STATS[bitlen]['3_1_G_y_is_over_8bit'] += 1
                with measure(bitlen, '3_time_check_glv_endomorphism'):
                    glv_info = check_glv_endomorphism(E, G)
                supports_glv = glv_info['supports_glv']
                if not supports_glv:
                    STATS[bitlen]['3_not_glv'] += 1
                    continue
                STATS[bitlen]['3_glv'] += 1

                lambda_val = glv_info['lambda']
                neg_lambda = -int(lambda_val) % E_q

                with measure(bitlen, '4_time_find_split_constants_explicit_tof'):
                    b1, b2 = find_split_constants_explicit_tof(p, E)

                ok = False
                ok_i = None
                with measure(bitlen, '4_time_try_decompose'):
                    checks = [
                        try_decompose(E_q, bitlen, b1, b2, lambda_val),
                        try_decompose(E_q, bitlen, -b1 % E_q, b2, lambda_val),
                        try_decompose(E_q, bitlen, b1, -b2 % E_q, lambda_val),
                        try_decompose(E_q, bitlen, -b1 % E_q, -b2 % E_q, lambda_val),
                        try_decompose(E_q, bitlen, b1, b2, neg_lambda),
                        try_decompose(E_q, bitlen, -b1 % E_q, b2, neg_lambda),
                        try_decompose(E_q, bitlen, b1, -b2 % E_q, neg_lambda),
                        try_decompose(E_q, bitlen, -b1 % E_q, -b2 % E_q, neg_lambda)
                    ]
                for i, x in enumerate(checks):
                    if x[2]:
                        ok = True
                        ok_i = i

                if ok:
                    STATS[bitlen]['4_glv_efficient_decompose'] += 1
                    print()
                    print(f"{bitlen}-bit Curve with b = {b} (log2={round(b_log2,3)},{round(b_log2/log2(p),3)}), j-invariant = {E.j_invariant()}, order = {E_q} (prime? {E_q_prime})")
                    print(' - G', G)
                    print(' - GLV', glv_info)
                    for i, x in enumerate(checks):
                        print(i, x)
                    print(ok)
                    print()
                    print()
                    FOUND += 1
                else:
                    STATS[bitlen]['4_glv_inefficient_decompose'] += 1

        for b, x in STATS.items():
            for y in sorted(x.keys()):
                print(b, y, x[y])
        print()

        handle.write(json.dumps((bitlen,STATS[bitlen])) + "\n")
        handle.flush()