cuerpos_finitos_robusto.py

import random, re

def primos(n): #devolvemos números primos menores o iguales que n
    if n < 2:
        return []
    criba = [True] * (n + 1) #criba[i] nos dirá si i es primo o no
    criba[0] = criba[1] = False
    for i in range(2, int(n**0.5) + 1):
        if criba[i]:
            for j in range(i*i, n + 1, i):
                criba[j] = False
    return [i for i, es_primo in enumerate(criba) if es_primo]


def descompone(n): #descompone n en factores primos
    factores = []
    lprimos = primos(n)
    i = 0
    while n > 1:
      if(n % lprimos[i] == 0):
        factores.append((lprimos[i], 1))
        n //= lprimos[i]
        while(n % lprimos[i] == 0):
          factores[-1] = (lprimos[i], factores[-1][1] + 1)
          n //= lprimos[i]

      i += 1

    return factores
    
class ElemFp:

  p = 5 #valor por defecto

  def __init__(self, value):
    if not isinstance(value, int):
      raise ValueError(f"El valor introducido: {value}, no es un entero")
    self.value = value % self.__class__.p

  def __repr__(self):
    return f"{self.value!r}"

  def __str__(self):
    return str(self.value)

  @classmethod
  def zero(cls):
    return cls(0)

  @classmethod
  def one(cls):
    return cls(1)

  def __add__(self,other):
    if not isinstance(other, self.__class__):
            raise TypeError("Solo se puede sumar con otro Fp")
    return self.__class__((self.value+other.value)%self.__class__.p)

  def __neg__(self):
    return self.__class__(-self.value % self.__class__.p)

  def __sub__(self, other):
    return self + (-other)

  def __mul__(self, other):
    if not isinstance(other, self.__class__):
        return NotImplemented
    return self.__class__((self.value * other.value) % self.__class__.p)

  def __pow__(self, k):
    if not isinstance(k, int):
        raise TypeError("Exponente debe ser entero")

    p = self.__class__.p

    if self.value == 0:
        if k < 0:
            raise ZeroDivisionError("No se puede calcular el inverso de 0")
        elif k == 0:
            return self.__class__(1) 
        else: # k > 0
            return self.__class__(0) 

    if k == 0:
        return self.__class__(1)
    
    k_eff = k % (p - 1)
    
    return self.__class__(pow(self.value, k_eff, p))

  def __truediv__(self, other):
    if not isinstance(other,self.__class__):
            raise TypeError("Solo se puede dividir entre otro Fp")
    return self*other**-1

  def __eq__(self, other):
    if not isinstance(other, self.__class__):
        return False
    return self.value == other.value
    
class ElemPol:

  clase_cuerpo = ElemFp #por defecto
  var = 'x'

  @classmethod
  def _trim_right_zeros(cls, t):
    t = tuple(t)  # asegurar que es tupla
    i = len(t)
    while i > 0 and t[i-1] == cls.clase_cuerpo.zero():
        i -= 1
    return t[:i]

  def __init__(self, coefs):
    for c in coefs:
      if not isinstance(c, self.__class__.clase_cuerpo):
        raise ValueError(f"El elemento {c} no pertenece al cuerpo {self.__class__.clase_cuerpo}")
    self.coefs = self._trim_right_zeros(coefs)

  def __repr__(self):
    return f"{self.coefs!r}"

  @classmethod
  def pol_to_str(cls, pol, var):
    terms = []
    if issubclass(cls.clase_cuerpo, ElemFp):
      c1 = ''
      c2 = ''
    else:
      c1 = '('
      c2 = ')'

    for i, coef in enumerate(pol):
        if coef == cls.clase_cuerpo.zero():
            continue

        # Decide la parte de la variable
        if i == 0:
            term = f"{coef}"
        elif i == 1:
            if coef == cls.clase_cuerpo.one():
                term = f"{var}"
            else:
                term = f"{c1}{coef}{c2}{var}"
        else:
            if coef == cls.clase_cuerpo.one():
                term = f"{var}^{i}"
            else:
                term = f"{c1}{coef}{c2}{var}^{i}"

        terms.append(term)

    if not terms:
        return "0" #quiza cambiar por zero()

    # Unión de términos con +, excepto el primero si tiene signo
    result = terms[0]
    for t in terms[1:]:
        result += " + " + t
    return result

  def __str__(self):
    return self.__class__.pol_to_str(self.coefs, self.__class__.var)

  @classmethod
  def zero(cls):
    return cls(())

  @classmethod
  def one(cls):
    return cls((cls.clase_cuerpo.one(),))

  def grado(self):
    return len(self.coefs) - 1 #grado de nulo -1

  def __neg__(self):
    return self.__class__(tuple(-c for c in self.coefs))

  def __add__(self, other):
        if not isinstance(other, self.__class__):
            raise TypeError("Solo se puede sumar con otro polinomio en el mismo cuerpo")
        max_len = max(len(self.coefs), len(other.coefs))
        coefs = []
        for i in range(max_len):
            a = self.coefs[i] if i < len(self.coefs) else self.__class__.clase_cuerpo.zero()
            b = other.coefs[i] if i < len(other.coefs) else self.__class__.clase_cuerpo.zero()
            coefs.append(a + b)
        return self.__class__(self._trim_right_zeros(tuple(coefs)))

  def __sub__(self, other):
        return self + (-other)

  def __mul__(self, other):
        if isinstance(other, self.__class__):
          coefs = [self.__class__.clase_cuerpo.zero()] * (len(self.coefs) + len(other.coefs) - 1)
          for i, a in enumerate(self.coefs):
              for j, b in enumerate(other.coefs):
                  coefs[i+j] = coefs[i+j] + (a * b)
          return self.__class__(self._trim_right_zeros(tuple(coefs)))
        elif isinstance(other, self.__class__.clase_cuerpo):
            # pol * escalar
            return other*self
        else:
            raise TypeError("Multiplicación no soportada")

  def __truediv__(self, other):
      if not isinstance(other, self.__class__):
            raise TypeError("Solo se puede dividir con otro polinomio en el mismo cuerpo")
      if other == self.__class__.zero():
            raise ZeroDivisionError("División por el polinomio 0")

      A = list(self.coefs)
      B = list(other.coefs)
      m, n = len(A)-1, len(B)-1

      Q = [self.__class__.clase_cuerpo.zero()] * (m - n + 1 if m >= n else 0)

      while len(A) >= len(B):
          # coeficiente líder
          coef = A[-1] / B[-1]
          d = len(A) - len(B)
          Q[d] = coef

          # resta coef*B*x^d
          for i in range(len(B)):
              A[i+d] = A[i+d] - coef * B[i]

          A = list(self._trim_right_zeros(A))

      return self.__class__(tuple(Q)), self.__class__(tuple(A))

  def __floordiv__(self, other):
        q, _ = self/other
        return q

  def __mod__(self, other):
        _, r = self/other
        return r

  def __pow__(self, k, a = None):
        if not isinstance(k, int):
            raise TypeError("El exponente debe ser un entero")
        if k < 0:
            raise ValueError("No está definida la potencia negativa en polinomios")

        result = self.__class__.one()
        base = self if a is None else self % a
        while k > 0:
          if k & 1:
              result = result * base
              if a is not None:
                  result = result % a
          base = base * base
          if a is not None:
              base = base % a
          k >>= 1
        return result

  def __rmul__(self, other):
    # Permite: escalar * pol
    if isinstance(other, self.__class__.clase_cuerpo):
      if other == self.__class__.clase_cuerpo.zero():
        return self.__class__.zero()
      return self.__class__(self._trim_right_zeros(tuple(c * other for c in self.coefs)))
    raise TypeError("Solo se puede multiplicar por un escalar del mismo cuerpo")

  def __eq__(self, other):
    if not isinstance(other, self.__class__):
        return False
    return self.coefs == other.coefs

  @classmethod
  def gcd_ext(cls, a, b):
      """
      Devuelve (g, x, y) tales que g = a*x + b*y,
      donde g = gcd(a, b) y g es mónico.
      """
      if not isinstance(a, cls) or not isinstance(b, cls):
          raise TypeError("gcd_ext solo está definido entre polinomios del mismo tipo")

      zero = cls.zero()
      one = cls.one()

      # Inicialización
      r0, r1 = a, b
      s0, s1 = one, zero
      t0, t1 = zero, one

      # Algoritmo de Euclides extendido
      while r1 != zero:
          q, r = r0 / r1
          r0, r1 = r1, r
          s0, s1 = s1, s0 - q * s1
          t0, t1 = t1, t0 - q * t1

      # Hacer g mónico
      if r0 != zero:
          lead_inv = r0.coefs[-1] ** -1
          r0 = lead_inv * r0
          s0 = lead_inv * s0
          t0 = lead_inv * t0

      return r0, s0, t0

  @classmethod
  def gcd(cls, a, b):
      g, _, _ = cls.gcd_ext(a, b)
      return g

  @classmethod
  def inv_mod(cls, a, b):
    """
    Devuelve x tal que a*x ≡ 1 (mod b)
    Lanza ValueError si a y b no son coprimos.
    """
    g, x, _ = cls.gcd_ext(a, b)
    one = cls.one()
    if g != one:
        raise ValueError("El inverso modular no existe: gcd(a, b) ≠ 1")
    return x % b

  def derivada(self):
    return self.__class__(tuple(self.clase_cuerpo(i%self.__class__.clase_cuerpo.p)*c for i, c in enumerate(self.coefs[1:], start=1))) #se hace modulo para asi en fq tratar con escalares reales, no polinomios asociados
    
class ElemFq:
    """
    Representa un único elemento del cuerpo F_q.
    """
    clase_pol = None
    g = None
    var = None
    p = None

    def __init__(self, value_pol):
        if isinstance(value_pol, self.__class__.clase_pol):
            self.value = value_pol % self.__class__.g
        else:
            raise TypeError(f"No se acepta el tipo {type(value_pol)} para construir ElemFq")
        """
        elif isinstance(value_pol, int):
            self.value = self.__class__.clase_pol(self.__class__.clase_pol.clase_cuerpo(value_pol))
        """

    def __add__(self, other):
        if not isinstance(other, self.__class__):
          raise TypeError("No se puede operar con un elemento de otro cuerpo")
        resultado_pol = (self.value + other.value) % self.__class__.g
        return self.__class__(resultado_pol)

    def __neg__(self):
        resultado_pol = -self.value % self.__class__.g
        return self.__class__(resultado_pol)

    def __sub__(self, other):
        return self + (-other)

    def __mul__(self, other):
        if not isinstance(other, self.__class__):
            return NotImplemented
        resultado_pol = (self.value * other.value) %  self.__class__.g
        return self.__class__(resultado_pol)

    def __pow__(self, k):
        if k >= 0:
          resultado_pol = pow(self.value, k, self.__class__.g)
        else:
          resultado_pol = pow(self.__class__.clase_pol.inv_mod(self.value, self.__class__.g), -k, self.__class__.g)
        return self.__class__(resultado_pol)

    def __truediv__(self, other):
        inv_other = self.__class__(self.__class__.clase_pol.inv_mod(other.value, self.__class__.g))
        return self * inv_other

    def __eq__(self, other):
        if not isinstance(other, self.__class__):
          return False
        return self.value == other.value #está todo siempre módulo g

    def __repr__(self):
        return f"ElemFq({self.value!r})"

    def __str__(self):
        return self.__class__.clase_pol.pol_to_str(self.value.coefs, self.__class__.var)

    @classmethod
    def zero(cls):
      return cls(cls.clase_pol.zero())

    @classmethod
    def one(cls):
      return cls(cls.clase_pol.one())

class cuerpo_fp:
    def __init__(self, p):  # construye el cuerpo de p elementos Fp = Z/pZ
        if(p not in primos(p)):
            raise ValueError(f"El número {p} no es primo")
        self.p = p

        class ElemFpP(ElemFp): #con esto impedimos que se puedan operar con elementos fabricados por distintos cuerpo_fp
            pass

        ElemFpP.p = p
        ElemFpP.__name__ = f"F{p}"

        self.elem_class = ElemFpP

    def cero(self):        # devuelve el elemento 0
        return self.elem_class.zero()

    def uno(self):         # devuelve el elemento 1
        return self.elem_class.one()

    def elem_de_int(self, n): # fabrica el elemento dado por la clase de n
        return self.elem_class(n)

    def elem_de_str(self, s): # fabrica el elemento a partir de un string (parser)
        try:
            return self.elem_class(int(s.strip()))
        except ValueError:
            raise ValueError(f"No se puede convertir '{s}' en un entero")

    def conv_a_int(self, a):  # devuelve un entero entre 0 y p-1
        if not isinstance(a, self.elem_class):
            raise TypeError(f"No es un F{self.p} fabricado por este cuerpo")
        return a.value

    def conv_a_str(self, a):  # pretty-printer
        if not isinstance(a, self.elem_class):
            raise TypeError(f"No es un F{self.p} fabricado por este cuerpo")
        return str(a)

    def suma(self, a, b):     # a+b
        if not isinstance(a, self.elem_class):
            raise TypeError(f"No es un F{self.p} fabricado por este cuerpo")
        return a + b

    def inv_adit(self, a):    # -a
        if not isinstance(a, self.elem_class):
            raise TypeError(f"No es un F{self.p} fabricado por este cuerpo")
        return -a

    def mult(self, a, b):     # a*b
        if not isinstance(a, self.elem_class):
            raise TypeError(f"No es un F{self.p} fabricado por este cuerpo")
        return a * b

    def pot(self, a, k):      # a^k (k entero)
        if not isinstance(a, self.elem_class):
            raise TypeError(f"No es un F{self.p} fabricado por este cuerpo")
        return a ** k

    def inv_mult(self, a):    # a^(-1)
        if not isinstance(a, self.elem_class):
            raise TypeError(f"No es un F{self.p} fabricado por este cuerpo")
        if a.value == 0:
            raise ZeroDivisionError("No existe el inverso multiplicativo de 0")
        return a ** -1

    def es_cero(self, a):     # a == 0
        if not isinstance(a, self.elem_class):
            raise TypeError(f"No es un F{self.p} fabricado por este cuerpo")
        return a == self.cero()

    def es_uno(self, a):      # a == 1
        if not isinstance(a, self.elem_class):
            raise TypeError(f"No es un F{self.p} fabricado por este cuerpo")
        return a == self.uno()

    def es_igual(self, a, b): # a == b
        if not isinstance(a, self.elem_class):
            raise TypeError(f"No es un F{self.p} fabricado por este cuerpo")
        return a == b

    def aleatorio(self):      # fabrica un elemento aleatorio con prob uniforme
        return self.elem_de_int(random.randint(0, self.p - 1))

    def tabla_suma(self):     # devuelve la matriz de pxp (lista de listas) de la suma
        return [[self.conv_a_int(self.elem_de_int(i) + self.elem_de_int(j)) for j in range(self.p)] for i in range(self.p)]

    def tabla_mult(self):     # devuelve la matriz de pxp (lista de listas) de la mult
        return [[self.conv_a_int(self.elem_de_int(i) * self.elem_de_int(j)) for j in range(self.p)] for i in range(self.p)]

    def tabla_inv_adit(self): # devuelve una lista de long p con los inv_adit
        return [self.conv_a_int(-self.elem_de_int(i)) for i in range(self.p)]

    def tabla_inv_mult(self): # devuelve una lista de long p con los inv_mult (en el índice 0 pone un '*')
        return ['*'] + [self.conv_a_int(self.inv_mult(self.elem_de_int(i))) for i in range(1, self.p)]

    def cuadrado_latino(self, a): # cuadrado latino a*i+j (con a != 0)
        if isinstance(a, int):
          a = self.elem_de_int(a)
        if self.es_cero(a):
            raise ValueError("a debe ser distinto de cero para cuadrado latino")
        return [[self.conv_a_int(a * self.elem_de_int(i) + self.elem_de_int(j)) for j in range(self.p)] for i in range(self.p)]

class anillo_fp_x:
    def __init__(self, fp, var='x'): # construye el anillo Fp[var]
        if isinstance(fp, cuerpo_fp):
            self.fp = fp
        elif isinstance(fp, int):
            #por si usuario pasa como fp el primo a usar
            self.fp = cuerpo_fp(fp)
        else:
            raise TypeError("No se acepta fp")

        if not isinstance(var, str):
            raise TypeError("var debe ser un string")
        self.var = var

        # Creamos una subclase local de ElemPol ligada a este cuerpo
        class Pol(ElemPol):
            clase_cuerpo = self.fp.elem_class
            var = self.var

        Pol.__name__ = f"{self.fp.elem_class}[{self.var}]"

        self.elem_class = Pol
        self.p = self.fp.p
    def cero(self):                 # 0
        return self.elem_class.zero()

    def uno(self):                  # 1
        return self.elem_class.one()

    def elem_de_tuple(self, a):     # fabrica un polinomio a partir de la tupla (a0, a1, ...), ai debe haber sido fabricado por fp!
        return self.elem_class(a)

    def elem_de_int(self, a):       # fabrica un polinomio a partir de los dígitos de a en base p
        if not isinstance(a, int):
            raise ValueError("Debe ser un entero")
        coefs = []
        n = a
        while n > 0:
            coefs.append(self.fp.elem_de_int(n % self.fp.p))
            n //= self.fp.p
        return self.elem_class(tuple(coefs))

    def elem_de_str(self, s):       # fabrica un polinomio a partir de un string (parser), hecho por IA
      terms = re.findall(r'([+-]?[^+-]+)', s.replace(" ",""))
      coeffs_dict = {}
      max_deg = 0
      for term in terms:
          if self.var in term:
              if '*' in term:
                  coef_part, x_part = term.split('*')
              else:
                  match = re.match(r'([+-]?\d*)(%s(\^\d+)?)' % re.escape(self.var), term)
                  if match:
                      coef_part, x_part, _ = match.groups()
                  else:
                      raise ValueError(f"Término inválido: {term}")
              deg = int(x_part.split('^')[1]) if '^' in x_part else 1
              coef = 1 if coef_part in ('','+') else -1 if coef_part=='-' else int(coef_part)
          else:
              deg = 0
              coef = int(term)
          coeffs_dict[deg] = coeffs_dict.get(deg,0) + coef
          max_deg = max(max_deg, deg)
      length = max_deg+1
      coeffs = [0]*length
      for deg, coef in coeffs_dict.items():
          coeffs[deg] = self.fp.elem_de_int(coef)
      return self.elem_de_tuple(coeffs)


    def conv_a_tuple(self, a):      # devuelve la tupla de coeficientes
        if not isinstance(a, self.elem_class):
          raise TypeError("Objeto no fabricado por este anillo")
        return a.coefs

    def conv_a_int(self, a):        # devuelve el entero correspondiente al polinomio
        if not isinstance(a, self.elem_class):
          raise TypeError("Objeto no fabricado por este anillo")
        n = 0
        aux = 1
        for c in a.coefs:
            n += self.fp.conv_a_int(c) * aux
            aux *= self.fp.p
        return n

    def conv_a_str(self, a):        # pretty-printer
        if not isinstance(a, self.elem_class):
          raise TypeError("Objeto no fabricado por este anillo")
        return str(a)

    def suma(self, a, b):           # a+b
        if not isinstance(a, self.elem_class):
          raise TypeError("Objeto no fabricado por este anillo")
        return a + b

    def inv_adit(self, a):          # -a
        if not isinstance(a, self.elem_class):
          raise TypeError("Objeto no fabricado por este anillo")
        return -a

    def mult(self, a, b):           # a*b
        if not isinstance(a, self.elem_class):
          raise TypeError("Objeto no fabricado por este anillo")
        return a * b

    def mult_por_escalar(self, a, e): # a*e (con e en Z/pZ)
        if not isinstance(a, self.elem_class):
          raise TypeError("Objeto no fabricado por este anillo")
        if not isinstance(e, self.fp.elem_class):
          raise TypeError("Escalar no compatible")
        return e * a

    def divmod(self, a, b):        # devuelve q,r tales que a=bq+r y deg(r)<deg(b)
        if not isinstance(a, self.elem_class):
          raise TypeError("Objeto no fabricado por este anillo")
        return a / b

    def div(self, a, b):           # q
        if not isinstance(a, self.elem_class):
          raise TypeError("Objeto no fabricado por este anillo")
        return a // b

    def mod(self, a, b):           # r
        if not isinstance(a, self.elem_class):
          raise TypeError("Objeto no fabricado por este anillo")
        return a % b

    def grado(self, a):            # deg(a)
        if not isinstance(a, self.elem_class):
          raise TypeError("Objeto no fabricado por este anillo")
        return a.grado()

    def gcd(self, a, b):           # devuelve g = gcd(a,b) mónico
        if not isinstance(a, self.elem_class):
          raise TypeError("Objeto no fabricado por este anillo")
        return self.elem_class.gcd(a, b)

    def gcd_ext(self, a, b):       # devuelve g,x,y tales que g=ax+by, g=gcd(a,b) mónico
        if not isinstance(a, self.elem_class):
          raise TypeError("Objeto no fabricado por este anillo")
        return self.elem_class.gcd_ext(a, b)

    def inv_mod(self, a, b):       # devuelve x tal que ax = 1 mod b
        if not isinstance(a, self.elem_class):
          raise TypeError("Objeto no fabricado por este anillo")
        return self.elem_class.inv_mod(a, b)

    def pot_mod(self, a, k, b):    # a^k mod b
        if not isinstance(a, self.elem_class):
          raise TypeError("Objeto no fabricado por este anillo")
        return a.__pow__(k, b)

    def es_cero(self, a):          # a == 0
        if not isinstance(a, self.elem_class):
          raise TypeError("Objeto no fabricado por este anillo")
        return a == self.cero()

    def es_uno(self, a):           # a == 1
        if not isinstance(a, self.elem_class):
          raise TypeError("Objeto no fabricado por este anillo")
        return a == self.uno()

    def es_igual(self, a, b):      # a == b
        if not isinstance(a, self.elem_class):
          raise TypeError("Objeto no fabricado por este anillo")
        return a == b

    def derivada(self, a):
        if not isinstance(a, self.elem_class):
          raise TypeError("Objeto no fabricado por este anillo")
        return a.derivada()
    """
    def sfd(self, a):
        if self.es_cero(a):
          return []

        df = self.derivada(f)
        g = self.gcd(f, df)
        w = f // g
        i = 1
        factores = []

        while not self.es_cero(w):
            y = self.gcd(w, g)
            fi = w // y
            if not self.es_uno(fi):
                factores.append((fi, i))
            w = y
            g = g // y
            i += 1

        return factores
    """
    def sfd(self, f): #para monicos, de grado > 0
      if self.es_cero(f):
          return []

      L = []
      s = 1
      p = self.fp.p

      while True:
          j = 1
          df = self.derivada(f)
          g = self.div(f, self.gcd(f, df))   # g = f / gcd(f, D(f))

          while not self.es_uno(g):
              f = self.div(f, g)
              h = self.gcd(f, g)
              m = self.div(g, h)
              if not self.es_uno(m):
                  L.append((m, j * s))
              g = h
              j += 1

          # Si f aún no es 1, significa que es una p-ésima potencia
          if not self.es_uno(f):
              # Calculamos la raíz p-ésima:
              # basta con tomar los coeficientes de índice múltiplo de p
              tupla_f = self.conv_a_tuple(f)
              coefs = [tupla_f[i] for i in range(0, len(tupla_f), p)]
              f = self.elem_de_tuple(tuple(coefs))
              s *= p
          else:
            break
      return L

    def ddf(self, f): #para monicos, grado > 0, square-free
        if self.es_cero(f):
          return []

        L = []

        pol_x = self.elem_de_int(self.p)
        h = self.mod(pol_x, f)
        p = self.fp.p
        k = 0

        while not self.es_uno(f):
          h = self.pot_mod(h, p, f)
          k += 1
          g = self.gcd(self.suma(h, self.inv_adit(pol_x)), f)
          if not self.es_uno(g):
            L.append((g, k))
            f = self.div(f, g)
            h = self.mod(h, f)

        return L

    def Mk(self, alpha, k, g):
      p = self.fp.p
      if p == 2:
          # caso p=2: suma de potencias de 2
          res = self.cero()
          for j in range(k):
              res = self.suma(res, self.pot_mod(alpha, 2**j, g))
          return res
      else:
          # caso p>2
          exp = (p**k - 1) // 2
          return self.suma(self.pot_mod(alpha, exp, g), self.inv_adit(self.uno()))


    def edf(self, f , k): #para monicos, grado > 0, producto de irreducibles de grado k
      h = [f]
      r = self.grado(f)//k
      while len(h) < r:
        h1 = []
        for g in h:
          alpha = self.elem_de_int(random.randint(0, self.p ** self.grado(g)-1))
          mk_alpha = self.Mk(alpha, k, g)
          d = self.gcd(mk_alpha, g)

          if self.es_uno(d) or self.es_igual(d, g):
                h1.append(g)
          else:
                h1.extend([d, self.div(g, d)])

        h = h1

      return h


    def es_irreducible(self, f):   # test de irreducibilidad de Rabin
        n = self.grado(f)

        if n == 0 or n == -1: #caso polinomio nulo y ctes
          return False

        if n == 1:
          return True

        p = self.fp.p
        pol_x = self.elem_de_int(self.p) #ya esta modulo f, por grado mayor
        pol_x_pn = self.pot_mod(pol_x, p ** n, f)

        if not self.es_igual(pol_x_pn, pol_x):
          return False

        for d,_ in descompone(n):
          q = p ** (n // d)
          pol_x_pq = self.pot_mod(pol_x, q, f)
          g = self.gcd(pol_x_pq - pol_x, f)
          if not self.es_uno(g):
            return False

        return True

    def factorizar(self, f):       # factorización de Cantor-Zassenhaus

          factores_finales = []
          for g, m in self.sfd(f):
              for h, k in self.ddf(g):
                  for t in self.edf(h, k):
                      factores_finales.append((t, m))

          return factores_finales
          
class cuerpo_fq:
    def __init__(self, fp, g, var='a'): # construye el cuerpo Fp[var]/<g(var)>, fp es el cuerpo, g tupla de elementos de fp
        if isinstance(fp, cuerpo_fp):
            self.fp = anillo_fp_x(fp, var)
            self.g = self.fp.elem_de_tuple(g)
        else:
            raise TypeError("No se acepta fp, debe ser un cuerpo")

        if not self.fp.es_irreducible(self.g):
            raise ValueError("g debe ser un polinomio irreducible")

        self.var = var
        self.q = self.fp.p ** self.g.grado()
        self.p = self.fp.p

        class sElemFq(ElemFq):
          clase_pol = self.fp.elem_class
          g = self.g
          var = self.var
          p = self.p

        sElemFq.__name__ = f"F{self.q}"
        self.elem_class = sElemFq

    def cero(self):                # 0
        return self.elem_class.zero()

    def uno(self):                 # 1
        return self.elem_class.one()

    def elem_de_tuple(self, coefs):    # fabrica elemento a partir de tupla de coeficientes de cuerpo_fp
        return self.elem_class(self.fp.elem_de_tuple(coefs))

    def elem_de_int(self, n):      # fabrica elemento a partir de entero
        return self.elem_class(self.fp.elem_de_int(n)) #asumimos n representa polinomio de Fp[x], luego para n > q habrá que reducirlo módulo g. Otra opción es hacer n%=q y entonces representará otra cosa
                                                      #para n < q ambas se comportan igual
    def elem_de_str(self, s):      # fabrica elemento parseando string
        return self.elem_class(self.fp.elem_de_str(s))

    def conv_a_tuple(self, a):     # devuelve tupla de coeficientes sin ceros "extra"
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el cuerpo")
        return self.fp.conv_a_tuple(a.value)

    def conv_a_int(self, a):       # devuelve el entero correspondiente
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el cuerpo")
        return self.fp.conv_a_int(a.value)

    def conv_a_str(self, a):       # pretty-printer
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el cuerpo")
        return str(a)

    def suma(self, a, b):          # a+b
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el cuerpo")
        return a+b

    def inv_adit(self, a):         # -a
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el cuerpo")
        return -a

    def mult(self, a, b):          # a*b
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el cuerpo")
        return a*b

    def pot(self, a, k):           # a^k (k entero)
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el cuerpo")
        return a**k

    def inv_mult(self, a):         # a^(-1)
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el cuerpo")
        return a**(-1)

    def es_cero(self, a):          # a == 0
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el cuerpo")
        return a == self.cero()

    def es_uno(self, a):           # a == 1
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el cuerpo")
        return a == self.uno()

    def es_igual(self, a, b):      # a == b, en todos los métodos elementos que se pasan ya están módulo g, han tenido que ser fabricados.
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el cuerpo")
        return a == b

    def aleatorio(self):           # devuelve un elemento aleatorio con prob uniforme
        return self.elem_de_int(random.randint(0, self.q - 1))

    def tabla_suma(self):          # matriz de qxq correspondiente a la suma (con la notación int)
        elems = [self.elem_de_int(i) for i in range(self.q)]
        return [[self.conv_a_int(self.suma(a,b)) for b in elems] for a in elems]

    def tabla_mult(self):          # matriz de qxq correspondiente a la mult (con la notación int)
        elems = [self.elem_de_int(i) for i in range(self.q)]
        return [[self.conv_a_int(self.mult(a,b)) for b in elems] for a in elems]

    def tabla_inv_adit(self):      # lista de inv_adit (con la notación int)
        elems = [self.elem_de_int(i) for i in range(self.q)]
        return [self.conv_a_int(self.inv_adit(a)) for a in elems]

    def tabla_inv_mult(self):      # lista de inv_mult (con la notación int)
        elems = [self.elem_de_int(i) for i in range(self.q)]
        res = ['*']  # índice 0 no tiene inverso
        for a in elems[1:]:
            try:
                res.append(self.conv_a_int(self.inv_mult(a)))
            except ValueError:
                res.append(None)
        return res

    def cuadrado_latino(self, a):  # cuadrado latino para a != 0 (con notación int)
        if isinstance(a, int):
          a = self.elem_de_int(a)
        if self.es_cero(a):
            raise ValueError("a debe ser distinto de cero")
        elems = [self.elem_de_int(i) for i in range(self.q)]
        return [[self.conv_a_int(self.suma(self.mult(a, elems[i]), elems[j])) for j in range(self.q)] for i in range(self.q)]


class anillo_fq_x: #copia de anillo_fp_x
    def __init__(self, fq, var='x'): # # Fq[var], var debe ser distinta que la de fq
        if isinstance(fq, cuerpo_fq):
            self.fq = fq
        else:
            raise TypeError("No se acepta fq")

        if var == fq.var:
            raise ValueError("var debe ser distinta de la de fq")
        self.var = var

        # Creamos una subclase local de ElemPol ligada a este cuerpo
        class Pol(ElemPol):
            clase_cuerpo = fq.elem_class
            var = self.var

        Pol.__name__ = f"{fq.elem_class}[{self.var}]"
        self.elem_class = Pol
        self.q = self.fq.q
        self.p = self.fq.p
    def cero(self):                 # 0
        return self.elem_class.zero()

    def uno(self):                  # 1
        return self.elem_class.one()

    def elem_de_tuple(self, a):     # fabrica un polinomio a partir de la tupla (a0, a1, ...)
        return self.elem_class(a)

    def elem_de_int(self, a):       # fabrica un polinomio a partir de los dígitos de a en base q
        if not isinstance(a, int):
            raise ValueError("Debe ser un entero")
        coefs = []
        n = a
        while n > 0:
            coefs.append(self.fq.elem_de_int(n % self.q))
            n //= self.q
        return self.elem_class(tuple(coefs))


    def _eval_coef_expr(self, expr_str):
        """
        Función de arranque para el nuevo parser.
        Crea una instancia de _CoefParser y ejecuta el parseo.
        """
        # 1. Tokenizer: Convierte el string en una lista de tokens.
        # Ej: "a^2 + 3" -> ['a', '^', '2', '+', '3']
        tokens = re.findall(r'(\d+|[a-zA-Z]+|[+\-*/^()])', expr_str)

        # 2. Creamos y ejecutamos el parser.
        parser = self._CoefParser(tokens, self.fq)
        return parser.parse()



    class _CoefParser:
        """
        Clase interna que implementa un parser recursivo descendente
        para las expresiones de los coeficientes.
        """
        def __init__(self, tokens, fq_factory):
            self.tokens = tokens
            self.fq = fq_factory
            self.pos = 0
            self.var_obj = self.fq.elem_de_tuple((self.fq.fp.fp.cero(), self.fq.fp.fp.uno()))
            self.var_str = self.fq.var

        def parse(self):
            if not self.tokens:
                return self.fq.cero()
            result = self.parse_expression()
            if self.pos < len(self.tokens):
                raise ValueError(f"Token inesperado al final de la expresión: '{self.peek()}'")
            return result

        def peek(self):
            return self.tokens[self.pos] if self.pos < len(self.tokens) else None

        def consume(self, expected_token=None):
            if self.pos >= len(self.tokens):
                if expected_token:
                    raise ValueError(f"Final inesperado de la expresión, se esperaba '{expected_token}'")
                return None

            token = self.tokens[self.pos]
            if expected_token and token != expected_token:
                raise ValueError(f"Se esperaba el token '{expected_token}' pero se encontró '{token}'")

            self.pos += 1
            return token

        def parse_atom(self):
            token = self.peek()
            if token is None:
                raise ValueError("Final inesperado de la expresión")

            if token.isdigit():
                return self.fq.elem_de_int(int(self.consume()))
            elif token == self.var_str:
                self.consume()
                return self.var_obj
            elif token == '(':
                self.consume('(')
                result = self.parse_expression()
                self.consume(')')
                return result
            elif token == '-':
                self.consume('-')
                return -self.parse_atom()
            else:
                raise ValueError(f"Token inesperado: '{token}'")

        # --- MÉTODO CORREGIDO ---
        def parse_factor(self): # Maneja potencias (^)
            left = self.parse_atom()
            while self.peek() == '^':
                self.consume('^') # Consume el operador '^'

                # El exponente DEBE ser un entero. No hacemos llamada recursiva.
                exp_token = self.peek()
                if exp_token is None or not exp_token.isdigit():
                    raise ValueError(f"Se esperaba un exponente entero después de '^', pero se encontró '{exp_token}'")

                exponent = int(self.consume())

                # Ahora la operación es ElemFq ** int, que sí funciona.
                left = left ** exponent
            return left

        def parse_term(self): # Maneja multiplicaciones (*) y divisiones (/)
            left = self.parse_factor()
            while self.peek() in ('*', '/'):
                op = self.consume()
                right = self.parse_factor()
                if op == '*':
                    left = left * right
                elif op == '/':
                    left = left / right
            return left

        def parse_expression(self): # Maneja sumas (+) y restas (-)
            left = self.parse_term()
            while self.peek() in ('+', '-'):
                op = self.consume()
                right = self.parse_term()
                if op == '+':
                    left = left + right
                elif op == '-':
                    left = left - right
            return left

    def elem_de_str(self, s): #hecho por IA
        # Paso 1: Normalizar el string para manejar la resta de forma consistente.
        # Reemplazamos ' - ' por ' + -' solo si está espaciado, para no tocar cosas como 'a-1'.
        s = s.replace(" - ", " + -")

        # Paso 2: Dividir el string en términos, respetando los paréntesis.
        terms_str = []
        current_term_start = 0
        paren_level = 0
        for i, char in enumerate(s):
            if char == '(':
                paren_level += 1
            elif char == ')':
                paren_level -= 1
            elif char == '+' and paren_level == 0:
                # Encontramos un '+' que separa términos (no está dentro de un paréntesis).
                terms_str.append(s[current_term_start:i].strip())
                current_term_start = i + 1

        # Añadir el último término después del último '+'
        terms_str.append(s[current_term_start:].strip())

        if not terms_str or terms_str == ['']: return self.cero()

        coeffs_dict = {}
        max_deg = 0

        for term in terms_str:
            term = term.strip()
            if not term: continue
            if term.startswith('+'): term = term[1:].strip()

            coef_str = "1"
            var_part = ""

            # El resto de la lógica para procesar un término individual funciona bien.
            if term.startswith('('):
                match_paren = term.rfind(')')
                if match_paren == -1:
                    raise ValueError(f"Paréntesis no balanceado en el término: '{term}'")

                coef_str = term[1:match_paren]
                var_part = term[match_paren+1:].strip()

                if var_part.startswith('*'):
                    var_part = var_part[1:].strip()
            else:
                try:
                    var_index = term.index(self.var)
                    coef_str = term[:var_index].strip()
                    var_part = term[var_index:].strip()

                    if coef_str.endswith('*'):
                        coef_str = coef_str[:-1].strip()

                    if not coef_str or coef_str == '+': coef_str = "1"
                    if coef_str == '-': coef_str = "-1"

                except ValueError:
                    coef_str = term
                    var_part = ""

            if not var_part:
                deg = 0
            elif '^' in var_part:
                try:
                    deg = int(var_part.split('^')[1])
                except (ValueError, IndexError):
                    raise ValueError(f"Exponente mal formado en el término: '{term}'")
            else:
                deg = 1

            coef_obj = self._eval_coef_expr(coef_str)

            current_coef = coeffs_dict.get(deg, self.fq.cero())
            coeffs_dict[deg] = current_coef + coef_obj
            max_deg = max(max_deg, deg)

        final_coeffs = []
        for i in range(max_deg + 1):
            final_coeffs.append(coeffs_dict.get(i, self.fq.cero()))

        return self.elem_de_tuple(tuple(final_coeffs))

    def conv_a_tuple(self, a):      # devuelve la tupla de coeficientes
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el anillo")
        return a.coefs

    def conv_a_int(self, a):        # devuelve el entero correspondiente al polinomio
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el anillo")
        n = 0
        aux = 1
        for c in a.coefs:
            n += self.fq.conv_a_int(c) * aux
            aux *= self.q
        return n

    def conv_a_str(self, a):        # pretty-printer
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el anillo")
        return str(a)

    def suma(self, a, b):           # a+b
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el anillo")
        return a + b

    def inv_adit(self, a):          # -a
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el anillo")
        return -a

    def mult(self, a, b):           # a*b
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el anillo")
        return a * b

    def mult_por_escalar(self, a, e): # a*e (con e en Fq)
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el anillo")
        return e * a

    def divmod(self, a, b):        # devuelve q,r tales que a=bq+r y deg(r)<deg(b)
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el anillo")
        return a / b

    def div(self, a, b):           # q
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el anillo")
        return a // b

    def mod(self, a, b):           # r
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el anillo")
        return a % b

    def grado(self, a):            # deg(a)
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el anillo")
        return a.grado()

    def gcd(self, a, b):           # devuelve g = gcd(a,b) mónico
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el anillo")
        return self.elem_class.gcd(a, b)

    def gcd_ext(self, a, b):       # devuelve g,x,y tales que g=ax+by, g=gcd(a,b) mónico
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el anillo")
        return self.elem_class.gcd_ext(a, b)

    def inv_mod(self, a, b):       # devuelve x tal que ax = 1 mod b
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el anillo")
        return self.elem_class.inv_mod(a, b)

    def pot_mod(self, a, k, b):    # a^k mod b
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el anillo")
        return a.__pow__(k, b)

    def es_cero(self, a):          # a == 0
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el anillo")
        return a == self.cero()

    def es_uno(self, a):           # a == 1
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el anillo")
        return a == self.uno()

    def es_igual(self, a, b):      # a == b
        if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el anillo")
        return a == b

    def derivada(self, a):
          if not isinstance(a, self.elem_class):
            raise TypeError("Objeto no fabricado por el anillo")

          nuevos_coefs = []
          for i, c in enumerate(a.coefs[1:], start=1):
              escalar_i = self.fq.elem_de_int(i%self.p) #modulo p para no multiplicar por polinomio de Fq

              nuevo_coef = escalar_i * c
              nuevos_coefs.append(nuevo_coef)

          return self.elem_de_tuple(tuple(nuevos_coefs))
    """
    def sfd(self, a):
        if self.es_cero(a):
          return []

        df = self.derivada(f)
        g = self.gcd(f, df)
        w = f // g
        i = 1
        factores = []

        while not self.es_cero(w):
            y = self.gcd(w, g)
            fi = w // y
            if not self.es_uno(fi):
                factores.append((fi, i))
            w = y
            g = g // y
            i += 1

        return factores
    """
    def sfd(self, f): #para monicos, de grado > 0
      if self.es_cero(f):
          return []

      L = []
      s = 1
      q = self.fq.q
      p = self.p
      while True:
          j = 1
          df = self.derivada(f)
          g = self.div(f, self.gcd(f, df))   # g = f / gcd(f, D(f))

          while not self.es_uno(g):
              f = self.div(f, g)
              h = self.gcd(f, g)
              m = self.div(g, h)
              if not self.es_uno(m):
                  L.append((m, j * s))
              g = h
              j += 1

          # Si f aún no es 1, significa que es una p-ésima potencia
          if not self.es_uno(f):
              # Calculamos la raíz p-ésima:
              # basta con tomar los coeficientes de índice múltiplo de p
              tupla_f = self.conv_a_tuple(f)
              coefs = [tupla_f[i] for i in range(0, len(tupla_f), p)]
              f = self.elem_de_tuple(tuple(coefs))
              s *= p
          else:
            break
      return L

    def ddf(self, f): #para monicos, grado > 0, square-free
        if self.es_cero(f):
          return []

        L = []

        pol_x = self.elem_de_int(self.q)
        h = self.mod(pol_x, f)
        q = self.q
        k = 0

        while not self.es_uno(f):
          h = self.pot_mod(h, q, f)
          k += 1
          g = self.gcd(self.suma(h, self.inv_adit(pol_x)), f)
          if not self.es_uno(g):
            L.append((g, k))
            f = self.div(f, g)
            h = self.mod(h, f)

        return L

    def Mk(self, alpha, k, g):
      q = self.q
      if self.p == 2:
          # caso p=2: suma de potencias de 2
          res = self.cero()
          for j in range(k):
              res = self.suma(res, self.pot_mod(alpha, 2**j, g))
          return res
      else:
          # caso p>2
          exp = (q**k - 1) // 2
          return self.suma(self.pot_mod(alpha, exp, g), self.inv_adit(self.uno()))


    def edf(self, f , k): #para monicos, grado > 0, producto de irreducibles de grado k
      h = [f]
      r = self.grado(f)//k
      while len(h) < r:
        h1 = []
        for g in h:
          alpha = self.elem_de_tuple(tuple(self.fq.aleatorio() for _ in range(self.grado(g))))
          mk_alpha = self.Mk(alpha, k, g)
          d = self.gcd(mk_alpha, g)

          if self.es_uno(d) or self.es_igual(d, g):
                h1.append(g)
          else:
                h1.extend([d, self.div(g, d)])

        h = h1

      return h


    def es_irreducible(self, f):   # test de irreducibilidad de Rabin
        n = self.grado(f)

        if n == 0 or n == -1: #caso polinomio nulo y ctes
          return False

        if n == 1:
          return True

        q = self.q
        pol_x = self.elem_de_int(q) #ya esta modulo f, por grado mayor
        pol_x_pn = self.pot_mod(pol_x, q ** n, f)

        if not self.es_igual(pol_x_pn, pol_x):
          return False

        for d,_ in descompone(n):
          q2 = q ** (n // d)
          pol_x_pq = self.pot_mod(pol_x, q2, f)
          g = self.gcd(pol_x_pq - pol_x, f)
          if not self.es_uno(g):
            return False

        return True

    def factorizar(self, f):       # factorización de Cantor-Zassenhaus

          factores_finales = []
          for g, m in self.sfd(f):
              for h, k in self.ddf(g):
                  for t in self.edf(h, k):
                      factores_finales.append((t, m))

          return factores_finales