ideals.py

"""
Helper functions for various computations associated to
the quaternion algebra, and ideals and orders of the 
quaternion algebra.

Some of these functions could be ported to SageMath. It's a TODO
for when SQISign is being less actively worked on.
"""

# Sage Imports
from sage.all import (
    randint,
    ZZ,
    ceil,
    log,
    gcd,
    Matrix,
    prod
)

# Local imports
from setup import p, O0, ω

# ================================================ #
#  Helpers for elements of the quaternion algebra  #
# ================================================ #

def quadratic_norm(x, y):
    """
    Given two integers x,y, which correspond 
    to the element x + ωy ∈ R[ω], returns
    Nrd(x + ωy)
    

    Note: This function implements the norm
    function f(x,y) in the SQISign papers.

    For SQISign, we have ω = i and i^2 = -1
    so f(x,y) = x**2 + y**2 which is more
    efficient
    """
    if ω**2 != -1:
        raise ValueError(f"quadratic_norm() requires Z[ω] with ω^2 = -1. {ω = }")
    return ZZ(x) ** 2 + ZZ(y) ** 2


def is_integral(I):
    """
    Checks whether the input ideal is integral.
    """
    return all([b in I.left_order() for b in I.basis()])


def ideal_basis_gcd(I):
    """
    Computes the gcd of the coefficients of
    the ideal written as a linear combination
    of the basis of its left order.
    """
    I_basis = I.basis_matrix()
    O_basis = I.left_order().unit_ideal().basis_matrix()

    # Write I in the basis of its left order
    M = I_basis * O_basis.inverse()
    return gcd((gcd(M_row) for M_row in M))


def is_cyclic(I):
    """
    Computes whether the input ideal is cyclic,
    all the work is done by the helper function
    `ideal_basis_gcd()`.
    """
    return ideal_basis_gcd(I) == 1


def reduced_basis(I, check=False):
    """
    Computes the Minkowski reduced basis of the
    input ideal. Note: this produces the same 
    result for all ideals in the equivalence class
    so corresponds to the reduced basis of the
    smallest equivalent ideal to I
    
    Input: an ideal
    Output: A Minkowski-reduced basis

    Optional: when check is True, the basis is
              checked against the Minkowski bounds
    """

    def _matrix_to_gens(M, B):
        """
        Converts from a matrix to generators in the quat.
        algebra
        """
        return [sum(c * g for c, g in zip(row, B)) for row in M]

    B = I.basis()
    G = I.gram_matrix()
    U = G.LLL_gram().transpose()

    reduced_basis_elements = _matrix_to_gens(U, B)

    if check:
        norm_product = 16 * prod([x.reduced_norm() for x in reduced_basis_elements])
        tmp = p**2 * I.norm() ** 4
        assert norm_product <= 4 * tmp, "Minkowski reduced basis is too large"
        assert norm_product >= tmp, "Minkowski reduced basis is too small"

    return reduced_basis_elements

def left_isomorphism(I, J):
    """
    Given two isomorphic left ideals I, J computes
    α such that J = I*α
    """
    B = I.quaternion_algebra()

    if B != J.quaternion_algebra():
        raise ValueError("Arguments must be ideals in the same algebra.")

    if I.left_order() != J.left_order():
        raise ValueError("Arguments must have the same left order.")

    IJ = I.conjugate() * J
    L = reduced_basis(IJ)
    for t in L:
        α = t / I.norm()
        if J == I * α:
            return α

    raise ValueError("Could not find a left isomorphism...")


def ideal_generator(I, coprime_factor=1):
    """
    Given an ideal I of norm D, finds a generator
    α such that I = O(α,D) = Oα + OD

    Optional: Enure the norm of the generator is coprime 
    to the integer coprime_factor
    """
    OI = I.left_order()
    D = ZZ(I.norm())
    bound = ceil(4 * log(p))

    gcd_norm = coprime_factor * D**2

    # Stop infinite loops.
    for _ in range(1000):
        α = sum([b * randint(-bound, bound) for b in I.basis()])
        if gcd(ZZ(α.reduced_norm()), gcd_norm) == D:
            assert I == OI * α + OI * D
            return α
    raise ValueError(f"Cannot find a good α for D = {D}, I = {I}, n(I) = {D}")