modified tate lichenbaum pairing

script/research/ec/pairing-modified-tate.sage

@@ -0,0 +1,110 @@
+# for more info check Washington example 11.7
+# https://github.com/narodnik/elliptic-curves-washington-solutions
+p = 11
+K = GF(p)
+E = EllipticCurve(K, [-1, 1])
+R.<x, y> = PolynomialRing(K)
+
+n = 5
+P = E(3, 6)
+assert P.order() == 5
+Q = P
+inf = E(0)
+
+# We are computing <P, P>_5
+#   v_5 = f_5(P) / f_5(P)
+# where
+#   div(f_5) = n[P + R] - n[R]
+
+D_P = ((1, P), (-1, inf))
+
+# Random point for D_Q
+R = E(0, 1)
+D_Q = ((1, Q + R), (-1, R))
+
+# n = 5 = 1 + 4
+# so we perform one addition, two doublings and another addition
+
+# Step 1
+i = n
+j = 0
+k = 1
+fj = fk = 1
+# Compute f1 such that
+#   div(f1) = [P + R] - [P] - [R] + [∞]
+# But D_P = [(3, 6)] - [∞]
+# so R = ∞, and hence D_1 = [P + ∞] - [P] - [∞] + [∞]
+#                         = 0
+# hence f1 = 1
+
+# First valuations of f0, f1
+(vj, vk) = (1, 1)
+
+dydx = lambda x, y: (3*x^2 - 1) / (2*y)
+
+print(f"i = {i}, j = {j}, k = {k}")
+while i != 0:
+    if i % 2 == 0:
+        i /= 2
+
+        # We are computing div(l) which is
+        # the line between kP and kP
+        kP = k*P
+        m = dydx(kP[0], kP[1])
+        c = kP[1] - m*kP[0]
+        l = y - m*x - c
+        # And now the vertical line through 2kP
+        _2kP = 2*kP
+        v = x - _2kP[0]
+
+        f = l / v
+
+        f_valuation = 1
+        for d, X in D_Q:
+            f_valuation *= f(X[0], X[1])^d
+        vk = vk^2 * f_valuation
+        print(f"  f = {f}")
+        print(f"  vk = {vk}")
+
+        k *= 2
+    else:
+        i -= 1
+
+        if j + k == 1:
+            assert k == 1
+            # fj = fk
+            # vj = vk
+        else:
+            # Interpolate jP and kP
+            assert j != k
+            jP = j*P
+            kP = k*P
+            print(f"{jP}, {kP}")
+            if jP[0] != kP[0]:
+                m = (jP[1] - kP[1]) / (jP[0] - kP[0])
+                c = jP[1] - m*jP[0]
+                l = y - m*x - c
+            else:
+                l = x - jP[0]
+            # Vertical line through (j + k)P
+            jkP = (j + k)*P
+            if jkP != inf:
+                v = x - jkP[0]
+            else:
+                v = K(1)
+
+            f = l / v
+            print(f"  f = {f}")
+            
+            f_valuation = 1
+            for d, X in D_Q:
+                f_valuation *= f(X[0], X[1])^d
+            vj = vj * vk * f_valuation
+            print(f"  vj = {vj}")
+
+        j += k
+
+    print(f"i = {i}, j = {j}, k = {k}")
+
+modified_tate = vj^((p - 1) / n)
+print(f"result = {modified_tate}")