Implementation of Galois Rings for SageMath

GaloisRing.sage

# Load using the following command in SageMath:
#
#     load("https://gist.githubusercontent.com/recmo/6393218cab40baa78e44766772d1ec34/raw/GaloisRing.sage")
#
# Get the latest version from https://gist.github.com/recmo/6393218cab40baa78e44766772d1ec34

from sage.structure.richcmp import richcmp

class GaloisRingElement(CommutativeRingElement):

    def __init__(self, parent, p):
        B = parent.base()
        self.p = B(p)
        CommutativeRingElement.__init__(self, parent)

    def __hash__(self):
        return hash(self.p)
    
    def __rshift__(self, other):
        return self.__class__(self.parent(), [c >> other for c in self.list()])

    def _repr_(self):
        str = ""
        if self.p == 0:
            return "0"
        for i, c in enumerate(self.list()):
            if self.parent().repr == "centered":
                c = c.lift_centered()
            elif self.parent().repr == "rational":
                try:
                    c = c.rational_reconstruction()
                except:
                    c = c.lift_centered()
            if c == 0:
                continue
            if str != "":
                str += " + " if c > 0 else " - "
            if self.parent().repr == "centered" or self.parent().repr == "rational":
                c = abs(c)
            if c != 1 or i == 0:
                str += f"{c}"
            if c != 1 and i > 0:
                str += "*"
            if i != 0:
                str += self.parent().variable_name
            if i > 1:
                str += f"^{i}"
        return str
    
    def list(self):
        return self.p.list()
    
    def residue(self):
        '''Returns the residue of this element in the residue field.'''
        F = self.parent().residue_field()
        e = sum([F(c) * F.gen()^i for i, c in enumerate(self.list())])
        return e
    
    def multiplicative_order(self):
        '''Returns the multiplicative order of this element.'''
        if not self.is_unit():
            raise ArithmeticError("multiplicative order of non-units is undefined")
        
        R = self.parent()
        x = self
        # Start with order of the multiplicative group, then remove factors.
        order = R.exponent()
        assert x^order == 1
        for (p, e) in factor(order):
            while order % p == 0:
                if x^(order // p) != 1:
                    break
                order //= p
        return order
    
    def p_adic(self):
        '''Returns the p-adic representation'''
        T = self.parent().teichmuller_set()
        t = dict(zip([t.residue() for t in T] , T))
        r = self
        a = []
        for i in range(self.parent().prime_exponent):
            a.append(t[r.residue()])
            r -= a[-1]
            r >>= 1
        return a
    
    def frobenius(self):
        '''Returns the Frobenius of this element.'''
        R = self.parent()
        return sum([R.prime^i * a^R.prime for i, a in enumerate(self.p_adic())])
    
    def trace(self):
        s = self
        t = s
        for i in range(self.parent().degree - 1):
            s = s.frobenius()
            t += s
        return t
        
    def norm(self):
        s = self
        t = s
        for i in range(self.parent().degree  - 1):
            s = s.frobenius()
            t *= s
        return t

    def _richcmp_(self, other, op):
        return richcmp(self.p, other.p, op)
    
    def _add_(self, other):
        return self.__class__(self.parent(), self.p + other.p)
    
    def _sub_(self, other):
        return self.__class__(self.parent(), self.p - other.p)

    def _mul_(self, other):
        return self.__class__(self.parent(), self.p * other.p)
    
    def _div_(self, other):
        return self * other.inverse()
    
    def is_unit(self):
        return self.residue() != 0
    
    def inverse(self):
        if not self.is_unit():
            raise ZeroDivisionError("division by zero")
        return self^(self.parent().exponent() - 1)

class GaloisRing(UniqueRepresentation, CommutativeRing):
    Element = GaloisRingElement
    
    def __init__(self, p, n, f, names=('x',), repr='centered'):
        # Prime field
        if not p.is_prime():
            raise ValueError("p must be prime")
        self.prime = p
        self.prime_field = GF(p)

        # Coefficient ring
        if n < 1:
            raise ValueError("n must be positive")
        self.prime_exponent = n
        self.coefficient_ring = Integers(self.prime^self.prime_exponent)

        # Galois field
        self.variable_name = names[0]
        P = PolynomialRing(self.prime_field, name=self.variable_name+'bar')
        modulus = f(P.gen())
        if not modulus.is_monic():
            raise ValueError("f must be monic")
        if not modulus.is_irreducible():
            raise ValueError("f must be irreducible")
        self.degree = modulus.degree()
        self.galois_field = GF((p, self.degree), modulus=modulus, name=self.variable_name)

        # Galois ring
        P = PolynomialRing(self.coefficient_ring, name=self.variable_name)
        self.modulus = f(P.gen())
        self.galois_ring = P.quotient(self.modulus, names=(self.variable_name,))

        self.repr = repr

        CommutativeRing.__init__(self, self.galois_ring)

    def _repr_(self):
        return f"Galois Ring in {self.variable_name} of characteristic {self.prime}^{self.prime_exponent} and extension degree {self.degree}"

    def characteristic(self):
        return self.prime^self.prime_exponent
    
    def cardinality(self):
        return self.prime^(self.degree * self.prime_exponent)

    def coefficient_ring(self):
        '''Returns the coefficient ring of the Galois ring, Z/p^dZ.'''
        return self.base().base_ring()
    
    def residue_field(self):
        '''Returns the residue field of this Galois ring, GF(p^d).'''
        return self.galois_field
    
    def gens(self):
        return [self.element_class(self, x) for x in self.base().gens()]
    
    def gen(self):
        '''Generator of the ring.'''
        return self.element_class(self, self.base().gen())
    
    def root(self, order):
        '''Returns a primitive root of given order.'''
        N = self.gen().multiplicative_order()
        if N % order != 0:
            raise ArithmeticError("Primitive root of given order not found.")
        return self.gen()^(N // order)
    
    def teichmuller_set(self):
        n = self.prime^self.degree - 1
        e = self.root(n)
        return [self.zero()] + [e^i for i in range(n)]
    
    def random_element(self):
        return self.element_class(self, self.base().random_element())
    
    def exponent(self):
        '''Returns the number N such that a^N = 1 for all units a.'''
        return (self.prime^self.degree - 1) * self.prime^((self.prime_exponent - 1) * self.degree)
    
    def lenstra_constant(self):
        '''Returns the Lenstra constant for this Galois ring.'''
        return self.prime ** self.degree
    
    def exceptional_sequence(self):
        '''Returns a maximal length exceptional sequence.'''
        return sorted([self.element_class(self, a) for a in self.galois_field])
    
    def _element_constructor_(self, *args, **kwds):
        if len(args)!=1:
            return self.element_class(self, *args, **kwds)
        x = args[0]
        try:
            P = x.parent()
        except AttributeError:
            return self.element_class(self, x, **kwds)
        return self.element_class(self, x, **kwds)

    def _coerce_map_from_(self, S):
        return self.base().has_coerce_map_from(S)