1  DIV_POINT = 1
  2  DIV_FUNC  = 2
  3  
  4  class Divisor:
  5  
  6      def __init__(self, field):
  7          self.field = field
  8          self._div = []
  9  
 10      def __call__(self, Px, Py, Pz=1):
 11          K = self.field
 12          # Convert to base field
 13          Px, Py, Pz = K(Px), K(Py), K(Pz)
 14          # Normalize coordinates
 15          if Pz > 0:
 16              Px /= Pz
 17              Py /= Pz
 18              Pz = 1
 19          P = (Px, Py, Pz)
 20  
 21          D = Divisor(K)
 22          D._div += [(DIV_POINT, 1, P)]
 23          return D
 24  
 25      def div(self, f):
 26          K = self.field
 27          D = Divisor(K)
 28          D._div += [(DIV_FUNC, 1, f)]
 29          return D
 30  
 31      def _clean(self):
 32          self._div = [(type_id, self._deg_obj(P), P)
 33                       for type_id, P in self._objs()]
 34          self._div = [(type_id, n, P) for type_id, n, P in self._div if n != 0]
 35      def _deg_obj(self, P):
 36          return sum(n for type_id, n, Q in self._div if P == Q)
 37      def _objs(self):
 38          return set((type_id, P) for type_id, _, P in self._div)
 39  
 40      def __add__(self, other):
 41          K = self.field
 42          D = Divisor(K)
 43          D._div = self._div[:] + other._div[:]
 44          D._clean()
 45          return D
 46  
 47      def __sub__(self, other):
 48          K = self.field
 49          D = self + -1*other
 50          return D
 51  
 52      def __mul__(self, n):
 53          K = self.field
 54          D = Divisor(K)
 55          D._div = [(type_id, n*m, P) for type_id, m, P in self._div]
 56          return D
 57  
 58      __rmul__ = __mul__
 59  
 60      def __repr__(self):
 61          out = ""
 62          if not self._div:
 63              out += "0"
 64          for i, (type_id, n, obj) in enumerate(self._div):
 65              assert n != 0
 66              if i > 0:
 67                  if n > 0:
 68                      out += " + "
 69                  else:
 70                      out += " - "
 71              else:
 72                  if n < 0:
 73                      out += "-"
 74              assert type_id in (DIV_POINT, DIV_FUNC)
 75              if abs(n) > 1:
 76                  out += f"{abs(n)}"
 77              if type_id == DIV_POINT:
 78                  out += self._format_point(obj)
 79              elif type_id == DIV_FUNC:
 80                  out += f" div({obj})"
 81          return out
 82  
 83      def _format_point(self, P):
 84          Px, Py, Pz = P
 85          assert Pz in (0, 1)
 86          if Pz == 0:
 87              return f"[∞]"
 88          return f"[({Px}, {Py})]"
 89  
 90      def deg(self):
 91          return sum(n for _, n, _ in self._div)
 92  
 93      def supp(self):
 94          return set(P for _, _, P in self._div)
 95  
 96  K.<x, y> = GF(47)[]
 97  D = Divisor(K)
 98  D = 4*D(2, 3) + D(2, 4) + 2*D(0, 1, 0) + 4*D.div(x^2 + y)
 99  E = 6*D(4, 2) - 6*D(6, 3) + 6*D(6, 3)
100  #D -= 4*D.div(x^2 + y)
101  #E -= 6*D(4, 2)
102  print(f"D = {D}")
103  print(f"E = {E}")
104  print(f"D + E = {D + E}")
105  print(f"6E = {6*E}")
106  print(f"deg(D) = {D.deg()}")
107  print(f"supp(D) = {D.supp()}")
108  print(f"deg(E) = {E.deg()}")
109  print(f"supp(E) = {E.supp()}")
110