Initial commitHEADmain

1 files changed, 1267 insertions, 0 deletions

diff --git a/utils.py b/utils.py
new file mode 100644
index 0000000..cdaec4b
--- /dev/null
+++ b/utils.py
@@ -0,0 +1,1267 @@
+#!/usr/bin/env python
+
+# Silence warnings for now
+# https://stackoverflow.com/a/15778297
+from warnings import simplefilter
+
+import sage.all
+from sage.arith.misc import (
+    divisors,
+    factor,
+    gcd,
+    is_prime_power,
+    is_square,
+    jacobi_symbol,
+    kronecker_symbol,
+    xgcd,
+)
+from math import prod
+from sage.all import pari
+from sage.arith.misc import is_power_of_two
+from sage.functions.other import ceil, floor
+from sage.matrix.constructor import Matrix
+from sage.misc.functional import sqrt
+from sage.modules.free_module_element import vector
+from sage.rings.finite_rings.finite_field_constructor import GF
+from sage.rings.integer import Integer
+from sage.rings.integer_ring import ZZ
+from sage.rings.number_field.number_field import QuadraticField
+from sage.schemes.elliptic_curves.cm import hilbert_class_polynomial
+from sage.schemes.elliptic_curves.constructor import EllipticCurve
+from sage.sets.primes import Primes
+from sage.structure.factorization_integer import IntegerFactorization
+
+from math import log
+
+from random import randint
+from itertools import product, cycle, islice
+
+# Silence warnings about slow implementations in the pre-processing step e.g.
+#     verbose 0 (4176: multi_polynomial_ideal.py, groebner_basis)
+#     Warning: falling back to very slow toy implementation.
+try:
+    from sage.misc.verbose import set_verbose
+    set_verbose(-1)
+except ImportError:
+    try:
+        from sage.misc.misc import set_verbose
+        set_verbose(-1)
+    except ImportError:
+        pass
+
+simplefilter(action='ignore', category=FutureWarning)
+simplefilter(action='ignore', category=DeprecationWarning)
+
+
+def padicval(p, N):
+    if N == 0:
+        return 0
+
+    if N == 1:
+        return 0
+
+    k = 0
+    while N % p**k == 0:
+        k += 1
+    return k - 1
+
+
+class Cideal():
+    def __init__(self, O, a, b, m):
+        self.O = O
+        self.a = a
+        self.b = b % a
+        self.m = m
+        self.f = O.conductor()
+        self.D = O.ambient().discriminant()
+
+        if (self.b**2 + self.f**2*(self.D**2 - self.D)//4 + self.f*self.b*self.D) % self.a != 0:
+            print(f"{a}, {b} not good choices")
+            raise ValueError
+
+    def norm(self):
+        return self.m**2 * self.a
+
+    def discriminant(self):
+        return self.D
+
+    def gens(self):
+        return [self.m*self.a, self.m*(self.b % self.a)]
+
+    def __mul__(self, other):
+        if isinstance(other, int) or isinstance(other, sage.rings.integer.Integer):
+            return Cideal(self.O, self.a, self.b, other*self.m)
+
+        if isinstance(other, list) or isinstance(other, tuple) or isinstance(other, sage.modules.vector_integer_dense.Vector_integer_dense):
+            assert len(other) == 2
+
+            # We expect other to be in the [1, fw] basis!
+
+            if other[1] == 0:
+                # In other words, we are multiplying by a + 0*f*w
+                return self.__mul__(other[0])
+
+            assert other in Cideal(self.O, 1, 0, 1)
+
+            f, D = self.f, self.D
+            K = (D**2 - D)//4
+
+            a, b = self.a, self.b
+            c, d = other
+
+            f, D = self.f, self.D
+            K = (D**2 - D)//4
+
+            # Generators of product written as rows in the basis [fw, 1]
+            # i.e. a + b*fw is written as [b, a]
+            # The hermite row reduced matrix with no zero rows will have the
+            # form
+            #     [ C, B ]
+            #     [ 0, A ]
+            # This is because Sage's hermite form does row reduction
+            G = Matrix(ZZ, [
+                            [a*d,             a*c],
+                            [c + b*d + d*D*f, c*b - d*f**2*K],
+                        ])
+
+            G = G.hermite_form(include_zero_rows=False)
+            # print()
+            # print(G)
+            # print(f"Multiplying ideal I = [{a}, {b} + f*w] by ({c}, {d//f}*f*w)")
+            A = G[1, 1]
+            B = G[0, 1]
+            C = G[0, 0]
+            # print(G)
+            M = self.m
+
+            # alpha = a*c
+            # beta = a*d
+            # gamma = c*b - d*f**2*K
+            # delta = b*d + c + d*D*f
+
+            # C, u, v = xgcd(beta, delta)
+            # A = gcd(alpha*delta, beta*gamma)/C
+            # B = u*alpha + v*gamma
+            # M = self.m
+
+            if C != 1:
+                if A % C != 0 or B % C != 0:
+                    print("There has been an error in the representation of the ideal")
+                    exit(148)
+                A = A//C
+                B = B//C
+                M = M*C
+
+            assert G[1, 0] == 0
+
+            return Cideal(self.O, A, B, M)
+
+        if isinstance(other, Cideal):
+            assert self.O == other.O
+
+            a1, b1 = self.a, self.b
+            a2, b2 = other.a, other.b
+            f, D = self.f, self.D
+            K = (D**2 - D)//4
+
+            # Generators of product written as rows in the basis [fw, 1]
+            # i.e. a + b*fw is written as [b, a]
+            # This is because Sage's hermite form does row reduction
+            # (we want columns, so transposing the bases gives the "right"
+            # result)
+            G = Matrix(ZZ, [
+                            [0,             a1*a2         ],
+                            [a1,            a1*b2         ],
+                            [a2,            a2*b1         ],
+                            [b1 + b2 + f*D, b1*b2 - f**2*K],
+                        ])
+
+            G = G.hermite_form(include_zero_rows=False)
+
+            if G.dimensions() != (2, 2):
+                print("Ideal multiplication failed")
+
+            C = G[0, 0]
+            B = G[0, 1]
+            A = G[1, 1]
+            M = self.m * other.m
+
+            if C != 1:
+                if A % C != 0 or B % C != 0:
+                    print("There has been an error in the representation of the ideal")
+                A = A//C
+                B = B//C
+                M = M*C
+
+            return Cideal(self.O, A, B, M)
+
+    __rmul__ = __mul__
+
+
+    def __pow__(self, other):
+        if isinstance(other, int) or isinstance(other, sage.rings.integer.Integer):
+            if other == 0:
+                return Cideal(self.O, 1, 0, 1)
+            # Use double and add method
+            result = Cideal(self.O, 1, 0, 1)
+            power = self
+            while other >= 1:
+                # print(other)
+                if other %2 == 1:
+                    result = result * power
+                    other = other - 1
+                power = power*power
+                other = other //2
+            return result
+
+    def __repr__(self):
+        # b_str = f"{self.b} + " if self.b != 0 else ""
+        # m_str = f"{self.m}*" if self.m != 1 else ""
+        # return f"{m_str}[{self.a}, {b_str}{self.f}*w]"
+
+        # True repr
+        def _factor(n): return 0 if n == 0 else factor(n)
+        return f"Cideal(QuadraticField({self.D}).order_of_conductor({_factor(self.f)}), {_factor(self.a)}, {_factor(self.b)}, {_factor(self.m)})"
+
+    def __shortrepr__(self):
+        b_str = f"{self.b} + " if self.b != 0 else ""
+        m_str = f"{self.m}*" if self.m != 1 else ""
+        return f"{m_str}[{self.a}, {b_str}{self.f}*w]"
+
+    def __eq__(self, other):
+        return self.a == other.a \
+                and self.b == other.b \
+                and self.m == other.m \
+                and self.O == other.O
+
+    def __hash__(self):
+        return hash((self.O, self.a, self.b, self.m))
+
+    def __contains__(self, other):
+
+        # Cast lists, tuples, vectors all as lists
+        other = list(other)
+
+        if len(other) == 2:
+            # In the [1, fw] basis!
+            c, d = other
+            # So mathematically
+            #   other = c + dfw
+            # This lies in O = m[a, b + fw] = [ma, m(b + fw)] if and only if
+            #   c + dfw = mxa + myb + myfw
+            # for some integers x, y. In other words,
+            #   (i)  d = my for some integer y
+            #        <=> d \in m\Z
+            #   (ii) c = mxa + myb for some  integer x (and same y from (i))
+            #        <=> c = mxa + db
+            #        <=> c - db = mxa
+            #        <=> c - db \in ma\Z
+            return (d % self.m == 0) and (c - self.b*d) % (self.m*self.a) == 0
+
+    def new_basis(self, other):
+
+        # Cast lists, tuples, vectors all as lists
+        other = list(other)
+
+        if len(other) == 2:
+            # In the [1, fw] basis!
+            c, d = other
+
+            # So mathematically
+            #   other = c + dfw
+            #
+            #   this lies in m[a, b + fw] if and only if c and d are multiples
+            #   of m and (c/m) + (d/m) fw lies in [a, b + fw]
+
+            assert c % self.m == 0  and d % self.m == 0
+
+            c = c // self.m
+            d = d // self.m
+
+            # Now we check that c + d fw lies in [a, b + fw]
+            #
+            # c + d fw = c + d (b + fw) - db = (c - db)/a * a + d (b + fw)
+            assert (c - self.b*d) % self.a == 0
+
+            alpha = (c - self.b*d) // self.a
+            beta = d
+
+            print(f"{c} + {d} fw = {self.m} * [{alpha} * {self.a} + {beta} * ({self.b} + fw)]")
+
+    def __floordiv__(self, other):
+        if isinstance(other, int) or isinstance(other, sage.rings.integer.Integer):
+            if self.m % other == 0:
+                return Cideal(self.O, self.a, self.b, self.m // other)
+            print("Could not divide exactly!")
+            exit(123)
+
+    def conjugate(self):
+        return Cideal(self.O, self.a, -(self.b + self.f*self.D), self.m)
+
+    def scalar_product(self, v, u):
+        # In the [1, fw] basis!
+        a, b = v
+        c, d = u
+        return ((2*a + b*self.f*self.D)*(2*c + d*self.f*self.D) - b*d*self.f**2*self.D)/4
+
+    def norm_sq(self, v):
+        x, y = v
+        return x**2 + x*y*self.f*self.D + y**2 * self.f**2 * (self.D**2-self.D)//4
+        # return self.scalar_product(v, v)
+
+    def reduced_basis(self):
+
+        def scalar_product(v, w): return self.scalar_product(v, w)
+        def norm_sq(v): return self.norm_sq(v)
+
+        def rround(n):
+            return -1*round(-n) if n < 0 else round(n)
+
+        # Using Galbraith's Algorithm 23 to compute a minimal basis
+        # In the [1, fw_D] basis!
+        A = self.m*vector(ZZ, [self.a, 0])
+        B = self.m*vector(ZZ, [self.b, 1])
+
+        b_1 = A
+        b_2 = B
+
+        B_1 = norm_sq(b_1)
+        mu = scalar_product(b_1, b_2)/B_1
+        b_2 = b_2 - rround(mu)*b_1
+        B_2 = norm_sq(b_2)
+
+        while B_2 < B_1:
+            b_1, b_2 = b_2, b_1
+            B_1 = norm_sq(b_1)
+            mu = scalar_product(b_1, b_2) / B_1
+            b_2 = b_2 - rround(mu)*b_1
+            B_2 = norm_sq(b_2)
+
+        # # Wikipedia's variant
+        # if norm_sq(b_1) <= norm_sq(b_2):
+        #     b_1, b_2 = b_2, b_1
+
+        # # b_2 has the smaller norm here
+        # while norm_sq(b_2) <= norm_sq(b_1):
+        #     # Get the projection factor of b_2 on b_1
+        #     mu = rround(scalar_product(b_2, b_1)/norm_sq(b_2))
+        #     r = b_1 - mu*b_2
+        #     b_1 = b_2
+        #     b_2 = r
+
+        if b_1 not in self:
+            print(self.a, self.b, self.m)
+            print(b_1)
+            assert b_1 in self
+
+        if b_2 not in self:
+            print(self.a, self.b, self.m)
+            print(b_2)
+            assert b_2 in self
+
+        if b_1 != b_2:
+            # Some easy sanity checks
+            assert norm_sq(b_1) <= norm_sq(b_2)
+            assert norm_sq(b_1) <= norm_sq(A)
+            assert norm_sq(b_1) <= norm_sq(B)
+            assert norm_sq(b_1 - b_2) >= norm_sq(b_1)
+            assert norm_sq(b_1 + b_2) >= norm_sq(b_1)
+
+        return b_1, b_2
+
+    def is_principal(self):
+        # if self.a == 1 and self.b == 0:
+        #     # This is an integer multiple of the whole ring
+        #     return True
+
+        b_1, b_2 = self.reduced_basis()
+
+        # print(f"Shortest vector in {self}: {b_1}, self.m = {self.m}")
+        # Now b_1 is the element of the lattice with the smallest norm
+        # Now we check whether this element generates the whole ring
+
+        return self == b_1 * Cideal(self.O, 1, 0, 1)
+
+    def is_equivalent(self, other):
+        return (self.conjugate()*other).is_principal()
+
+
+def find_ideals(O, norm, equivalent_to=None):
+    f = O.conductor()
+    D = O.ambient().discriminant()
+    n_w_D = (D**2 - D) // 4
+    for a in divisors(norm):
+        if Integer(norm//a).is_square():
+            for b in range(a):
+                if (b**2 + b*f*D + f**2*n_w_D) % a == 0:
+                    m = Integer(sqrt(norm//a))
+                    I = Cideal(O, a, b, m)
+                    if not equivalent_to:
+                        yield I
+                    else:
+                        if I.is_equivalent(equivalent_to):
+                            yield I
+
+
+def find_ideal(O, norm, equivalent_to=None):
+    ideals_iterator = find_ideals(O, norm, equivalent_to=equivalent_to)
+    return next(ideals_iterator)
+
+
+def qfbsolve(a, b, c, enn, nonzero=False, not2power=False, flag=3):
+    if isinstance(enn, list):
+        enn = IntegerFactorization(enn)
+    # From Pari's documentation (Page 218)
+    # https://pari.math.u-bordeaux.fr/pub/pari/manuals/2.15.5/users.pdf
+    # The flag is a bitmask
+    # Bit 1: decides if all solutions (1) or only one solution (0) are returned
+    # Bit 2: decides if the returned solutions are primitive (0) or not (1)
+    # That is
+    # flag=0, Return only one primitive solution
+    # flag=1, Return all primitive solutions
+    # flag=2, Return only one (possibly imprimitive) solution
+    # flag=3, Return all primitive and imprimitive solutions
+    #
+    # qfbsolve returns objects of type
+    #   <class 'cypari2.gen.Gen'> <class 'cypari2.gen.Gen'>
+    # ...
+
+    if flag == 0 or flag == 2:
+        # Pari returns a single pair in this case (_not_ a 1 element list)
+        _solutions = pari.qfbsolve(pari.Qfb(a, b, c), enn, flag=flag)
+        if _solutions:
+            solutions = [_solutions]
+        else:
+            solutions = []
+    else:
+        # Here pari returns a list of pairs
+        solutions = pari.qfbsolve(pari.Qfb(a, b, c), enn, flag=flag)
+
+    # ... we want them as integers
+
+    if b == 0:
+        # Some normalisation
+        solutions = [(abs(int(x)), abs(int(y))) for (x, y) in solutions]
+    else:
+        solutions = [(int(x), int(y)) for (x, y) in solutions]
+
+    if nonzero:
+        solutions = [(x, y) for (x, y) in solutions if x != 0 and y != 0]
+
+    if not2power:
+        solutions = [(x, y) for (x, y) in solutions if not is_power_of_two(x) and not is_power_of_two(y)]
+
+    # Stupid dedupe
+    solutions = list(set(solutions))
+
+    return solutions
+
+
+def has_torsion(E, n=2, card=False):
+    # Try to avoid counting points multiple times by calling abelian_group once
+    Eab = E.abelian_group()
+    invs = Eab.invariants()
+    if len(invs) == 2:
+        m1, m2 = invs
+        if card:
+            cardinality = Eab.cardinality()
+            return min(padicval(n, m1), padicval(n, m2)), cardinality
+        return min(padicval(n, m1), padicval(n, m2))
+    else:
+        # E is cyclic
+        if card:
+            return 0, Eab.cardinality()
+        return 0
+
+
+def get_point(E, order):
+    while True:
+        P = E.random_point()
+        if P == E.zero:
+            continue
+        if P.order() % order == 0:
+            yield (P.order()//order)*P
+
+
+def info(E, D=None, N=None, T=None, name=None, tries=None, indent=-1):
+    """Pretty prints information about Elliptic Curves"""
+    # N is the target torsion
+    # T is the actual torsion
+    # D is the discriminant (of End(E))
+
+    ret_str = []
+
+    C = E.cardinality()
+    p = E.base_field().characteristic()
+    q = E.base_field().cardinality()
+    k = padicval(p, q)
+
+    if not name:
+        name = 'E'
+
+    if N:
+        (n, l) = is_prime_power(N, get_data=True)
+
+    q_str = f"{p}^{k}" if k != 1 else f"{p}"
+
+    E_str = f"{name} = EllipticCurve(GF({q_str}), {E.a_invariants()})"
+    ret_str += [E_str]
+
+    C_str = f"card({name}) = {factor(C)}"
+    ret_str += [C_str]
+
+    if N:
+        if not T:
+            t = has_torsion(E, n)
+        else:
+            t = padicval(n, T)
+
+        tors_str = f"{n}-Torsion: {n}^{t} (Target {n}^{l})"
+        ret_str += [tors_str]
+
+    if tries:
+        ret_str += [f"Iterations to find {name}: {tries}"]
+
+    if N and D:
+        v = int(sqrt(((q+1-C)**2 - 4*q)/D))
+        v_val = padicval(-D, v)
+        two_val = padicval(2, v)
+        # if two_val < l:
+        #     print("Found one!")
+        #     exit(3)
+        v_str = f"v = {factor(v)}, Torsion losses: v_{-D}(v) = {v_val}, v_2(v) = {two_val}"
+
+        # cls_no = Integer(D).class_number()
+        # ret_str += [f"{D = }, h(D) = {cls_no}"]
+
+        ret_str += [v_str]
+
+    if indent == -1:
+        return " | ".join(ret_str)
+    else:
+        return " "*indent + ("\n"+" "*indent).join(ret_str)
+
+
+def find_E(D, N, cls_no=None, cls_poly=None, max_multiple=None, rand=False):
+    """Finds ordinary curve E whose N-torsion is rational
+
+    Assumes N % 4 == 0
+
+    Since the N-torsion is a subgroup of order N^2, our general strategy is to
+    first find curves with cardinality m*N^2 and then check whether this
+    constitutes the N-torsion
+
+    Passing cls_no and cls_poly is useful when finding many curves with a fixed
+    discriminant.
+
+    Passing rand may be beneficial when the method gets "stuck" on smaller
+    cardinality attempts.
+    """
+
+    if not max_multiple:
+        max_multiple = N**4
+
+    if not cls_no:
+        # Must cast as sage integer for the class_number() method to exist
+        cls_no = Integer(D).class_number()
+
+    if not cls_poly:
+        cls_poly = hilbert_class_polynomial(D)
+
+    if rand:
+        # Generator expression for lazy sampling
+        # There is probably a better inbuilt solution for this
+        ms = (randint(1, max_multiple) for x in range(1, max_multiple))
+    else:
+        ms = range(1, max_multiple)
+
+    for m in ms:
+        # If the N-torsion is to be rational, the curve must have
+        # F_q-cardinality that is a multiple of N^2
+        C = m*N**2
+
+        # Now searching for prime powers q that are of the form
+        #     q = C + 1 \pm sqrt(4C + v^2 D)
+        # We do this by first finding v's so that 4C - v^2 D is a square
+        # In other words,
+        #     x^2 = 4C + v^2 D    <=>    4C = x^2 - v^2 D
+        for x, v in qfbsolve(1, 0, -D, 4*C):
+            for sign in [-1, 1]:
+                q = C + 1 + sign * x
+                # Verify that this is a prime power
+                (p, k) = is_prime_power(q, get_data=True)
+                if k == 0:
+                    continue
+
+                # Ensure the resulting curve is ordinary
+                if jacobi_symbol(D, p) != 1:
+                    continue
+
+                # Get roots over F_q (listed as pairs (root, multiplicity))
+                for j, _ in cls_poly.roots(GF(q)):
+                    E = EllipticCurve(GF(q), j=j)
+
+                    assert E.is_ordinary()
+
+                    M = E.cardinality()
+
+                    if M % N**2 != 0:
+                        # Got the wrong twist! Choose the other...
+                        # (assume j != 0, 1728)
+                        E = E.quadratic_twist()
+
+                    if not has_torsion(E, N):
+                        continue
+
+                    # Because there is a point of order 4, there should be a
+                    # montgomery model for this curve
+                    E = E.montgomery_model()
+                    return E
+
+
+def positive_bezout_coefficients(x, y, N, squares=False, allowzero=False):
+    """Compute all positive integer pairs (a, b) so that ax + by = N
+
+    Assumes that x, y, N > 0
+    """
+
+    if x == 0 and y == 0:
+        return ()
+
+    if y == 0:
+        if allowzero:
+            return ((N // x, 0)) if N % x == 0 else ()
+        return ()
+
+    if x == 0:
+        if allowzero:
+            return ((0, N // y)) if N % y == 0 else ()
+        return ()
+
+    (g, a, b) = xgcd(x, y)
+
+    if N % g != 0:
+        return ()
+
+    _N = N // g
+    _y = y // g
+    _x = x // g
+
+    coeffs = ((N*a + k*y, N*b - k*x) for k in range(ceil(float(-N*a/y)), floor(float(N*b/x))+1))
+
+    for A, B in coeffs:
+
+        if squares and not is_square(A) and not is_square(B):
+            continue
+
+        if not allowzero and (A == 0 or B == 0):
+            continue
+
+        yield A, B
+
+
+def compute_good_norms(O, T, prime_at_least=2):
+    f = O.conductor()
+    D = O.discriminant()//f**2
+
+    P = Primes()
+    good_primes = []
+    p = P.first()  # p = 2
+    # Not allowing even numbers for now, getting next prime
+    p = P.next(p)  # p = 3
+
+    while p < T:
+        if f % 2 == 0:
+            if kronecker_symbol(f**2*D//4, p) != -1:
+                good_primes += [p]
+        else:
+            if kronecker_symbol(f**2*D, 4*p) != -1:
+                good_primes += [p]
+
+        p = P.next(p)
+
+    print(f"There are {len(good_primes)} good primes to choose from")
+    # This could be done way, way faster using some "reverse sieving" technique
+    # For now, our numbers are small enough that it doesn't make a big
+    # difference
+    good_primes = set(good_primes)
+    good_norms = reverse_sieve(good_primes, T, prime_at_least=prime_at_least)
+    good_norms = sorted(good_norms)
+    print(f"This gives us {len(good_norms)} many possible norms")
+
+    return good_norms
+
+
+def get_interesting_ideals(O, T, at_most_ideals=0, prime_at_least=2):
+    """Return pairs of coprime integers (m, j) for which there exist (x, y) such
+    that m*x^2 + j*y^2 = T
+    """
+    good_norms = compute_good_norms(O, T)
+
+    ideal_combinations = []
+
+    for i, m in enumerate(good_norms):
+        if i % 50 == 0:
+            print(f"\rFinding an ideal for every possible norm up to {T}: {i/(len(good_norms))*100:.2f}% Done", end='')
+
+        ideals_norm_m = list(find_ideals(O, m))
+
+        if len(ideals_norm_m) == 0:
+            print("No ideals found, even though {m = } was a good norm!")
+            continue
+
+        for j in good_norms:
+            if j >= m:
+                break
+
+            if gcd(m, j) != 1:
+                continue
+
+            sols = qfbsolve(m, 0, j, T, nonzero=True, not2power=True)
+
+            if not sols:
+                continue
+
+            ideals_norm_j = list(find_ideals(O, j))
+
+            if len(ideals_norm_m) == 0:
+                print("No ideals found, even thought {j = } was a good norm!")
+                continue
+
+            success = False
+            for I, J in product(ideals_norm_m, ideals_norm_j):
+                assert I.norm() == m and J.norm() == j
+                if I.is_equivalent(J):
+                    success = True
+                    break
+
+            if not success:
+                continue
+
+            print("Success!")
+
+            # print(f"norm(I) = {I.norm()}, I = {I}")
+            # print(f"norm(J) = {J.norm()}, J = {J}")
+            # print(f"I.conjugate()*J = {I.conjugate()*J}")
+            # print(f"Are equivalent = {I.is_equivalent(J)}")
+            # print(f"Smallest elements of I.conjugate()*J = {(I.conjugate()*J).reduced_basis()}")
+
+            for x, y in sols:
+                ideal_combinations += [(I, x, J, y)]
+                if at_most_ideals != 0 and len(ideal_combinations) >= at_most_ideals:
+                    print()
+                    return ideal_combinations
+
+    print()
+    return ideal_combinations
+
+
+
+# def sample(a, b, c, D, f):
+
+
+# def get_interesting_ideals(O, T, at_most_ideals=0, prime_at_least=1):
+#     """Return pairs of coprime integers (m, j) for which there exist (x, y) such
+#     that m*x^2 + j*y^2 = T
+#     """
+#     f = O.conductor()
+#     D = O.discriminant()//f**2
+
+#     P = Primes()
+#     good_primes = []
+#     p = P.first()  # p = 2
+#     # Not allowing even numbers for now, getting next prime
+#     p = P.next(p)  # p = 3
+
+#     while p < T:
+#         if f % 2 == 0:
+#             if kronecker_symbol(f**2*D//4, p) != -1:
+#                 good_primes += [p]
+#         else:
+#             if kronecker_symbol(f**2*D, 4*p) != -1:
+#                 good_primes += [p]
+
+#         p = P.next(p)
+
+#     print(f"There are {len(good_primes)} good primes to choose from")
+#     good_primes = set(good_primes)
+#     good_primes = sorted(list(good_primes))
+#     print(good_primes[-1])
+#     good_norms = reverse_sieve(good_primes, 10*T, prime_at_least=prime_at_least)
+#     good_norms = sorted(good_norms)
+#     print(good_norms[-1])
+#     exit(123)
+#     print(f"This gives us {len(good_norms)} many possible norms")
+
+#     for _m in range(len(good_norms)):
+#         for _j in range(_m):
+#             m = good_norms[_m]
+#             j = good_norms[_j]
+
+#             if not is_square(4*m*j + f**2*D):
+#                 print(f"{4*m*j + f**2*D} is not a square")
+#                 continue
+
+#             # Since f is even, b is even
+#             b = int(sqrt(4*m*j + f**2*D))
+
+#             I = Cideal(O, a, -(b + f**2*D)//2, 1)
+#             J = Cideal(O, c, (b - f**2*D)//2, 1)
+
+#             print("The ideals are equivalent:", I.is_equivalent(J))
+#             print("{I.norm()}, {J.norm()}")
+#             assert I.is_equivalent(J)
+
+#     ideal_combinations = []
+
+#     return ideal_combinations
+
+
+def reverse_sieve(primes, bound, prime_at_least=1):
+    if prime_at_least >= bound:
+        print("Starting point is greater than the upper bound")
+        raise ValueError
+    # This is very dumb
+    # Could be done nicer with some constructive (recusive even) algorithm to
+    # multiply numbers together from the primes
+    # The hard part is to know when to terminate
+    return [m for m in range(prime_at_least, bound) if set([p for p, _ in factor(m)]).issubset(set(primes))]
+
+
+def prod(it):
+    result = 1
+    for i in it:
+        result = result * i
+    return result
+
+
+def wD(E, D):
+    """Computing the isogeny corresponding to w_D, i.e. the generator of the
+    endomorphism ring with conductor f
+    """
+    q = E.base_field().cardinality()
+    C = E.cardinality()
+    t = q + 1 - C
+    # Throws error if not coercable to int
+    v = ZZ(sqrt((t**2 - 4*q)/D))
+
+    assert (v*D - t) % 2 == 0
+    assert gcd(v, D) == 1
+
+    deg_wD = (D**2 - D)//4
+
+    # Try smaller extensions first, before going to the full extension
+    for extension in [deg_wD - 1, deg_wD + 1, deg_wD**2 - 1]:
+        base_field = E.base_field()
+        irr_element = base_field['x'].irreducible_element(extension)
+        extension_field = base_field.extension(modulus=irr_element, names='v')
+
+        R = E.change_ring(extension_field)
+
+        pi_E = E.frobenius_endomorphism()
+        pi_R = R.frobenius_endomorphism()
+
+        wD_kernel = []
+        for P in R.zero().division_points(deg_wD):
+            if P.order() == deg_wD and (v*D - t)//2 * P + pi_R(P) == R.zero():
+                wD_kernel = [P]
+                wD_R = R.isogeny(wD_kernel)
+                if wD_R.codomain().j_invariant() == E.j_invariant():
+                    break
+
+        if len(wD_kernel) == 0:
+            continue
+
+        wD_R = wD_R.codomain().isomorphism_to(R) * wD_R
+
+        if wD_R.degree() != deg_wD:
+            continue
+
+        # Cheap method to restrict wD_R to E
+        def wD(P):
+            # Assumes P lies in E
+            return E(wD_R(R(P)))
+
+        # Ensure the map we created is the correct one. Mathematically
+        #   wD = ((vD - t)/2 + \pi)/v
+        if not all([v*wD(P) == (v*D - t)//2 * P + pi_E(P) for P in [E.random_point() for _ in range(100)]]):
+            continue
+
+        # w_D satisfies w_D^2 - D w_D + (D^2 - D)/4 = 0
+        if not all([wD(wD(P)) - D * wD(P) + (D**2 - D)//4 * P == E.zero() for P in [E.random_point() for _ in range(100)]]):
+            continue
+
+        return wD
+
+
+def good_coprime_norm_pairs(n, D, f):
+    for norm_I in range(2, f**2 * (int(1 + sqrt(-D)))):
+        if norm_I % D == 0 or norm_I % f == 0:
+            continue
+        for norm_J in range(2, min(norm_I, 2**n - norm_I)):
+            if norm_J % D == 0 or norm_J % f == 0:
+                continue
+            if gcd(norm_I, norm_J) != 1:
+                continue
+
+            # If
+            #   norm_I*a^2 + norm_J*b^2 = 2^m for some m <= n with m + 2k = n
+            # then
+            #   norm_I*(2^k a)^2 + norm_J*(2^k b)^2 = 2^m * 2^(2k) = 2^n
+            # Conversely
+            #   norm_I*a^2 + norm_J*b^2 = 2^m for some m <= n with m + 2k + 1 = n
+            # then
+            #   norm_I*(2^k a)^2 + norm_J*(2^k b)^2 = 2^m * 2^(2k) = 2^{n - 1}
+            #
+            # So by finding a solution for
+            #   norm_I*a^2 + norm_J*b^2 = 2^n
+            # or
+            #   norm_I*a^2 + norm_J*b^2 = 2^{n-1}
+            # we have verified that there exists a solution to (*)
+
+            # Use flag=2 to return after finding one solution
+            # We only care that there is one solution
+
+            sols = qfbsolve(norm_I, 0, norm_J, [(2, n)], flag=2)
+
+            if not sols:
+                sols = qfbsolve(norm_I, 0, norm_J, [(2, n - 1)], flag=2)
+
+            if sols:
+                a, b = sols[0]
+                yield norm_I, a, norm_J, b
+
+
+def find_class_group(d, f):
+    D = d if d % 4 == 1 else 4 * d
+    K = QuadraticField(D, names="i")
+    O = K.order_of_conductor(f)
+    (_f, e_f) = is_prime_power(f, get_data=True)
+    class_number = f - (f // _f) * kronecker_symbol(D, _f)
+    assert class_number == Integer(D * f**2).class_number()
+
+    # int() rounds down, which is what we want here
+    maxnorm = f * int(sqrt(-D))
+    class_group = [Cideal(O, 1, 0, 1)]
+    for norm in range(1, maxnorm):
+        for I in find_ideals(O, norm):
+            # For simplicity, only find ideals with m = 1
+            # (Means coprimality conditions are easier to meet later on)
+            if I.m > 1:
+                continue
+
+            alpha = I.a
+            beta = 2 * I.b + I.f * I.D
+            gamma = (I.b**2 + I.b * I.f * I.D + I.f**2 * (I.D**2 - I.D) // 4) // I.a
+
+            if gcd([alpha, beta, gamma]) > 1:
+                continue
+
+            if not any([I.is_equivalent(_) for _ in class_group]):
+                class_group += [I]
+                print(f"Found ideal: {len(class_group)}/{class_number} norm = {I.norm()} (maxnorm: {maxnorm})")
+                if len(class_group) >= class_number:
+                    break
+            if len(class_group) >= class_number:
+                break
+        if len(class_group) >= class_number:
+            break
+
+    assert len(class_group) == class_number
+
+    print(f"Computed class group of size {class_number}")
+    return class_group
+
+
+def roundrobin(*iterables):
+    # https://docs.python.org/3/library/itertools.html#recipes
+    "Visit input iterables in a cycle until each is exhausted."
+    # roundrobin('ABC', 'D', 'EF') → A D E B F C
+    # Algorithm credited to George Sakkis
+    iterators = map(iter, iterables)
+    for num_active in range(len(iterables), 0, -1):
+        iterators = cycle(islice(iterators, num_active))
+        yield from map(next, iterators)
+
+
+def find_pairs_squares(d, n, f):
+    D = d if d % 4 == 1 else 4*d
+    K = QuadraticField(D, names='i')
+    (i, ) = K._first_ngens(1)
+    O = K.order_of_conductor(f)
+
+    found = []
+    count_pairs = 0
+    (_f, e_f) = is_prime_power(f, get_data=True)
+    class_number = f - (f//_f)*kronecker_symbol(D, _f)
+
+    for norm_I, a, norm_J, b in good_coprime_norm_pairs(n, D, f):
+        count_pairs += 1
+        if len(found) >= class_number:
+            return found
+
+        # Iter over the ideals with norm_I and find J with norm(J) = norm_J
+        # Since norm(m[a, b + fw]) = m^2a == norm_I, there are only
+        # finitely many such ideals
+        for I in find_ideals(O, norm_I):
+            # We are only interested in finding classes of ideals
+            if any([I.is_equivalent(_) for _ in found]):
+                continue
+
+            alpha = I.a
+            beta = 2*I.b + f*D
+            gamma = (I.b**2 + I.b*I.f*I.D + I.f**2*(I.D**2 - I.D)//4) // I.a
+
+            # This returns the coefficients of beta, if they exist
+            # We only care if there is one solution (flag=2)
+            equiv_sol = qfbsolve(alpha, beta, gamma, norm_J, flag=2)
+
+            if equiv_sol:
+                (x, y) = equiv_sol[0]
+
+                # Recall that conjugate(w) = D - w. Hence
+                #   conjugate(x + yfw) = x + yfD - yfw
+                # so in the [1, fw] basis  we get
+                #   conjugate([x, y]) = [x + yfD, -y]
+                beta = I.m*vector(ZZ, [x*I.a + y*I.b, y])
+                beta_conj = vector(ZZ, [beta[0] + beta[1]*I.f*I.D, -beta[1]])
+                J = (beta_conj*I) // I.norm()
+
+                # I or J, doesn't matter. We are interested in the class
+                found += [I]
+
+                kani_endo = (J.conjugate()*I).reduced_basis()[0]
+                vstring = []
+                vstring += [f"(v1)"]
+                vstring += [f"Norm pairs tried: {count_pairs:>10}"]
+                vstring += [f"Distinct classes: {len(found)}/{class_number}"]
+                vstring += [f"2^{padicval(2, a**2*norm_I + b**2*norm_J)} == {a}^2*{norm_I} + {b}^2*{norm_J}"]
+                vstring += [f"kani_endo == {kani_endo[0]} + {kani_endo[1]} fw"]
+                print(" | ".join(vstring))
+
+    print("Could not find pairs for the whole class group")
+    return found
+
+def find_pairs_full(class_group, d, n, f):
+
+    D = d if d % 4 == 1 else 4 * d
+
+    found = []
+    for n_H, H in enumerate(class_group):
+        tries = 0
+
+        # Norm form of equivalent ideals to H
+        alpha = H.a
+        beta = 2 * H.b + H.f * H.D
+        gamma = (H.b**2 + H.b * H.f * H.D + H.f**2 * (H.D**2 - H.D) // 4) // H.a
+
+        # Should already be true, because of how the class group was computed
+        assert gcd([alpha, beta, gamma]) == 1
+
+        # Norm form of O = [1, fw]
+        delta = 1
+        epsilon = f * D
+        phi = f**2 * (D**2 - D) // 4
+
+        finish = False
+
+        while True:
+
+            if finish:
+                break
+
+            upper_bound = 2 ** (n // 6)
+
+            x_I = randint(-upper_bound, upper_bound)
+            y_I = randint(-upper_bound, upper_bound)
+            norm_I = alpha * x_I**2 + beta * x_I * y_I + gamma * y_I**2
+
+            x_J = randint(-upper_bound, upper_bound)
+            y_J = randint(-upper_bound, upper_bound)
+            norm_J = alpha * x_J**2 + beta * x_J * y_J + gamma * y_J**2
+
+            if gcd(norm_I, norm_J) > 1:
+                continue
+
+            tries += 1
+            # print(f"{tries} Coprime pair")
+
+            if tries > 2**12:
+                # Don't get stuck in infinite loop
+                print(f"Class {n_H+1:>2}/{len(class_group)} ({len(found)/(n_H+1):.1%}): Giving up...")
+                finish = True
+                break
+
+            bezout_coefficients = 0
+            for a, b in roundrobin(
+                positive_bezout_coefficients(norm_I, norm_J, 2 ** (n - 1)),
+                positive_bezout_coefficients(norm_I, norm_J, 2**n),
+            ):
+                bezout_coefficients += 1
+
+                # g = gcd(a, b)
+                # a, b = a//gcd(a, b), b//gcd(a, b)
+
+                if finish:
+                    break
+
+                if bezout_coefficients > 2**10:
+                    break
+
+                if gcd(a*norm_I, b*norm_J) > 1:
+                    continue
+
+                if is_square(a):
+                    end_a = (int(sqrt(a)), 0)
+                    end_a_str = f"{end_a[0]}^2"
+                elif sols := qfbsolve(delta, epsilon, phi, a, flag=2):
+                    end_a = sols[0]
+                    end_a_str = f"deg({end_a[0]} + {end_a[1]} fw)"
+                else:
+                    continue
+
+                # Assert that the degree of the endmorphism is correct
+                assert delta*end_a[0]**2 + epsilon*end_a[0]*end_a[1] + phi*end_a[1]**2 == a
+
+                if is_square(b):
+                    end_b = (int(sqrt(b)), 0)
+                    end_b_str = f"{int(sqrt(b))}^2"
+                elif sols := qfbsolve(delta, epsilon, phi, b, flag=2):
+                    end_b = sols[0]
+                    end_b_str = f"deg({end_b[0]} + {end_b[1]} fw)"
+                else:
+                    continue
+
+                # Assert that the degree of the endmorphism is correct
+                assert delta*end_b[0]**2 + epsilon*end_b[0]*end_b[1] + phi*end_b[1]**2 == b
+
+                finish = True
+
+                # Recover the ideals
+                beta_I_conj = H.m * vector(
+                    ZZ, [x_I * H.a + y_I * (H.b + H.f * H.D), -y_I]
+                )
+                I = beta_I_conj * H
+                I = I // H.norm()
+
+                assert norm_I == I.norm()
+
+                # beta_J = H.m*vector(ZZ, [x_J*H.a + y_J*H.b, y_J])
+                beta_J_conj = H.m * vector(
+                    ZZ, [x_J * H.a + y_J * (H.b + H.f * H.D), -y_J]
+                )
+                J = beta_J_conj * H
+                J = J // H.norm()
+
+                assert norm_J == J.norm()
+
+                # assert gcd(a * norm_I, b * norm_J) == 1
+                assert I.is_equivalent(H)
+                assert J.is_equivalent(H)
+                assert I.is_equivalent(J)
+                assert J.is_equivalent(I)
+
+                print()
+
+                kani = (J.conjugate() * I).reduced_basis()[0]
+
+                assert delta*kani[0]**2 + epsilon*kani[0]*kani[1] + phi*kani[1]**2 == J.norm()*I.norm()
+
+                found += [(H, a, I, b, J, end_a, end_b, kani)]
+
+                verb_str = []
+                verb_str += [f"Class {n_H+1:>2}/{len(class_group)}"]
+                verb_str += [f"({len(found)/(n_H+1):.1%})"]
+                verb_str += [f"(tries: {tries:>4}, b_coefs: {bezout_coefficients:>4})"]
+                # Class group computed so that H.m = 1
+                verb_str += [f"Found for class [H] = [{H.a}, {H.b} + fw]"]
+                verb_str += [f"{norm_I}*{end_a_str} + {norm_J}*{end_b_str} = 2^{padicval(2, norm_I*a + norm_J*b)}"]
+                verb_str += [f"(Kani end = {kani[0]} + {kani[1]} fw)"]
+                print(" ".join(verb_str))
+
+    return found
+
+
+def volcano_walk(d, n, _f, e_f):
+    # # Chosen parameters
+    # # Discriminant of the Elliptic curve we find
+    # d = -7
+    # # Available 2 torsion on middle curve
+    # n = 32
+    # # Conductor of endomorphism ring of middle curve
+    # _f = 3
+    # e_f = 5
+
+    # Derived parameters
+    # Discriminant, called \Delta_D in the thesis
+    D = d if d % 4 == 1 else 4*d
+    K = QuadraticField(D, names='i')
+    (i, ) = K._first_ngens(1)
+    f = _f**e_f
+    N = f*2**n
+
+    # print(f"Pre-computing class polynomial and class number of {D}...")
+    cls_no = Integer(D).class_number(proof=False)
+    cls_poly = hilbert_class_polynomial(D)
+    # print("Done")
+    # print()
+
+    sim_cls_no = f - (f//_f)*kronecker_symbol(D, _f)
+
+    cases = {1: 'Elkies', -1: 'Atikin', 0: 'Ramified'}
+
+    print("Setup:")
+    print(f"    {d = }, ({log(-d, 2):.2f} bits)")
+    print(f"    {D = }, Discriminant of {K = }")
+    print(f"    h(D) = {cls_no} ({log(cls_no, 2):.2f} bits)")
+    print(f"    h(f**2D) = {sim_cls_no} ({log(sim_cls_no, 2):.2f} bits)")
+    print(f"    N = {factor(N)}, Target rational torsion of curve")
+    print(f"    In the {cases[kronecker_symbol(D, 2)]} case")
+    print()
+
+    print("Looking for curve...")
+    print()
+    result = find_E(D, N, cls_no=cls_no, cls_poly=cls_poly)
+    if result:
+        E = result
+        print(info(E))
+    else:
+        print()
+        print("No curve found!")
+        exit(1)
+
+    q = E.base_field().cardinality()
+    p, k = q.is_prime_power(get_data=True)
+    q_str = f"{p}^{k}" if k != 1 else f"{p}"
+
+    assert has_torsion(E, 2) >= n
+    assert has_torsion(E, _f) >= e_f
+
+    # Walk down the _f-isogeny volcano to the floor.
+
+    P, Q = E.torsion_basis(f)
+    _P, _Q = (_f)**(e_f-1)*P, (_f)**(e_f-1)*Q
+
+    one_away_success = False
+    for kernel in [_P + k*_Q for kernel in range(_f)] + [_Q]:
+        if E.isogeny(kernel).codomain().j_invariant() != E.j_invariant():
+            one_away_success = True
+            correct_torsion_point = P + k*Q
+            break
+
+    assert one_away_success
+
+    phi = E.isogeny(correct_torsion_point, algorithm='factored')
+    F = phi.codomain()
+
+    assert has_torsion(F, _f) == 0
+
+    assert E == phi.domain()
+    assert F == phi.codomain()
+
+    print(f"Found curve with conductor {f = } and rational 2^{has_torsion(F, 2)}-torsion")
+    print(f"F = EllipticCurve(GF({q_str}), {F.a_invariants()})")
+
+    assert has_torsion(F, 2) >= n
+
+    return E, F, phi