dlog calculations

script/research/ec/dlog.sage

@@ -0,0 +1,280 @@
+# Initialize an elliptic curve
+p = 115792089237316195423570985008687907853269984665640564039457584007908834671663
+r = 115792089237316195423570985008687907852837564279074904382605163141518161494337
+Fp = GF(p)  # Base Field
+Fr = GF(r)  # Scalar Field
+A = 0
+B = 7
+E = EllipticCurve(GF(p), [A,B])
+assert(E.cardinality() == r)
+
+K.<x> = PolynomialRing(Fp, implementation="generic")
+L.<y> = PolynomialRing(K, implementation="generic")
+eqn = y^2 - x^3 - A * x - B
+
+# Returns line passing through points, works for all points and returns 1 for O + O = O
+def line(A, B):
+    if A == 0 and B == 0:
+        return 1
+    else:
+        [a, b, c] = Matrix([A, B, -(A+B)]).transpose().kernel().basis()[0]
+        return a*x + b*y + c
+
+def dlog(D):
+    # Derivative via partials
+    Dx = D.differentiate(x)
+    Dy = D.differentiate(y)
+    Dz = Dx + Dy * ((3*x^2 + A) / (2*y))
+    assert D != 0
+    return Dz/D
+
+def dlog_alt(D):
+    Dx = D.differentiate(x)
+    Dy = D.differentiate(y)
+    #Dz = Dx + Dy * ((3*x^2 + A) / (2*y))
+
+    # Normally we calculate:
+    #   Dz/D
+    # Due to a bug in sage, we will make the denominator D
+    # solely an equation in x by taking its norm.
+
+    # Denominator = V · V'
+    V = 2*y * D
+    # 2y Dz
+    Dz_numer = (2*y*Dx + Dy * (3*x^2 + A) * V(y=-y)).mod(eqn)
+    # Change denominator to the norm
+    D_denom = (V * V(y=-y)).mod(eqn)
+
+    return Dz_numer / D_denom
+
+P0 = E.random_element()
+P1 = E.random_element()
+P2 = E.random_element()
+Q = -int(Fr(5)^-1) * (P0 + 2*P1 + 3*P2)
+assert P0 + 2*P1 + 3*P2 + 5*Q == 0
+
+def div_add(div_f, div_g):
+    div = div_f.copy()
+    for P, n in div_g.items():
+        if P in div:
+            div[P] += n
+        else:
+            div[P] = n
+    div = dict((P, n) for P, n in div.items() if n != 0)
+    return div
+
+def div_invert(div):
+    return dict((P, -n) for P, n in div.items())
+
+def div_sub(div_f, div_g):
+    inv_div_g = div_invert(div_g)
+    return div_add(div_f, inv_div_g)
+
+# 2[P₂] + [-2P₂] - 3[∞]
+f1 = line(P2, P2)
+D1 = {"P2": 2, "-2P2": 1, "∞": -3}
+# 2[P₁] + [-2P₁] - 3[∞]
+f2 = line(P1, P1)
+D2 = {"P1": 2, "-2P1": 1, "∞": -3}
+# [P₂] + [-P₂] - 2[∞]
+f3 = line(P2, -P2)
+D3 = {"P2": 1, "-P2": 1, "∞": -2}
+# [P₀] + [-P₀] - 2[∞]
+f4 = line(P0, -P0)
+D4 = {"P0": 1, "-P0": 1, "∞": -2}
+
+# (2[P₂] + [-2P₂] - 3[∞]
+#  + 2[P₁] + [-2P₁] - 3[∞]
+#  + [P₂] + [-P₂] - 2[∞]
+#  + [P₀] + [-P₀] - 2[∞])
+# =
+# [P₀] + 2[P₁] + 3[P₂] + [-P₀] + [-2P₁] + [-2P₂] + [-P₂] - 10[∞]
+f5 = f1*f2*f3*f4
+D5 = div_add(div_add(D1, D2), div_add(D3, D4))
+assert D5 == {
+    "P0":   1,
+    "P1":   2,
+    "P2":   3,
+    "-P0":  1,
+    "-2P1": 1,
+    "-2P2": 1,
+    "-P2":  1,
+    "∞":  -10
+}
+
+# [-2P₂] + [-2P₁] + [2(P₁ + P₂)] - 3[∞]
+f6 = line(-2*P2, -2*P1)
+D6 = {"-2P2": 1, "-2P1": 1, "2P1 + 2P2": 1, "∞": -3}
+# [-P₂] + [-P₀] + [P₀ + P₂] - 3[∞]
+f7 = line(-P2, -P0)
+D7 = {"-P2": 1, "-P0": 1, "P0 + P2": 1, "∞": -3}
+
+# ([P₀] + 2[P₁] + 3[P₂] + [-P₀] + [-2P₁] + [-2P₂] + [-P₂] - 10[∞]
+#  - [-2P₂] - [-2P₁] - [2(P₁ + P₂)] + 3[∞]
+#  - [-P₂] - [-P₀] - [P₀ + P₂] + 3[∞])
+# =
+# [P₀] + 2[P₁] + 3[P₂] - [2(P₁ + P₂)] - [P₀ + P₂] - 4[∞]
+f8 = f5/(f6*f7)
+D8 = div_sub(D5, div_add(D6, D7))
+assert D8 == {
+    "P0":         1,
+    "P1":         2,
+    "P2":         3,
+    "2P1 + 2P2": -1,
+    "P0 + P2":   -1,
+    "∞":         -4
+}
+
+# [P₀ + P₂] + [2(P₁ + P₂)] + [-(P₀ + 2P₁ + 3P₂)] - 3[∞]
+f9 = line(P0 + P2, 2*(P1 + P2))
+D9 = {"P0 + P2": 1, "2P1 + 2P2": 1, "5Q": 1, "∞": -3}
+# ([P₀] + 2[P₁] + 3[P₂] - [2(P₁ + P₂)] - [P₀ + P₂] - 4[∞]
+#  + [P₀ + P₂] + [2(P₁ + P₂)] + [-(P₀ + 2P₁ + 3P₂)] - 3[∞])
+# =
+# [P₀] + 2[P₁] + 3[P₂] + [-(P₀ + 2P₁ + 3P₂)] - 7[∞]
+# = [P₀] + 2[P₁] + 3[P₂] + [5Q] - 7[∞]
+f10 = f8*f9
+D10 = div_add(D8, D9)
+assert D10 == {
+    "P0": 1,
+    "P1": 2,
+    "P2": 3,
+    "5Q": 1,
+    "∞": -7
+}
+
+# Now construct 5[Q]
+
+# 2[Q] + [-2Q] - 3[∞]
+f11 = line(Q, Q)
+D11 = {"Q": 2, "-2Q": 1, "∞": -3}
+
+# [-2Q] + [2Q] - 2[∞]
+f12 = line(-2*Q, 2*Q)
+D12 = {"-2Q": 1, "2Q": 1, "∞": -2}
+
+# (2[Q] + [-2Q] - 3[∞]) - ([-2Q] + [2Q] - 2[∞])
+# ==
+# 2[Q] - [2Q] - [∞]
+f13 = f11/f12
+D13 = div_sub(D11, D12)
+assert D13 == {
+    "Q": 2,
+    "2Q": -1,
+    "∞": -1
+}
+
+# multiply by 3
+# 6[Q] - 3[2Q] - 3[∞]
+f14 = f13*f13*f13
+D14 = div_add(div_add(D13, D13), D13)
+assert D14 == {
+    "Q": 6,
+    "2Q": -3,
+    "∞": -3
+}
+
+# 2[2Q] + [-4Q] - 3[∞]
+f15 = line(2*Q, 2*Q)
+D15 = {"2Q": 2, "-4Q": 1, "∞": -3}
+
+# (6[Q] - 3[2Q] - 3[∞]) + (2[2Q] + [-4Q] - 3[∞])
+# ==
+# 6[Q] - [2Q] + [-4Q] - 6[∞]
+f16 = f14*f15
+D16 = div_add(D14, D15)
+assert D16 == {
+    "Q": 6,
+    "2Q": -1,
+    "-4Q": 1,
+    "∞": -6
+}
+
+# [2Q] + [-2Q] - 2[∞]
+f17 = line(2*Q, -2*Q)
+D17 = {"2Q": 1, "-2Q": 1, "∞": -2}
+
+# (6[Q] - [2Q] + [-4Q] - 6[∞]) + ([2Q] + [-2Q] - 2[∞])
+# ==
+# 6[Q] + [-2Q] + [-4Q] - 8[∞]
+f18 = f16*f17
+D18 = div_add(D16, D17)
+assert D18 == {
+    "Q": 6,
+    "-2Q": 1,
+    "-4Q": 1,
+    "∞": -8
+}
+
+# [-2Q] + [-4Q] + [6Q] - 3[∞]
+f19 = line(-2*Q, -4*Q)
+D19 = {"-2Q": 1, "-4Q": 1, "6Q": 1, "∞": -3}
+
+# (6[Q] + [-2Q] + [-4Q] - 8[∞]) - ([-2Q] + [-4Q] + [6Q] - 3[∞])
+# ==
+# 6[Q] - [6Q] - 5[∞]
+f20 = f18/f19
+D20 = div_sub(D18, D19)
+assert D20 == {
+    "Q": 6,
+    "6Q": -1,
+    "∞": -5
+}
+
+# [6Q] + [-6Q] - 2[∞]
+f21 = line(6*Q, -6*Q)
+D21 = {"6Q": 1, "-6Q": 1, "∞": -2}
+
+# (6[Q] - [6Q] - 5[∞]) + ([6Q] + [-6Q] - 2[∞])
+# ==
+# 6[Q] + [-6Q] - 7[∞]
+f22 = f20*f21
+D22 = div_add(D20, D21)
+assert D22 == {"Q": 6, "-6Q": 1, "∞": -7}
+
+# [Q] + [-6Q] + [5Q] - 3[∞]
+f23 = line(Q, -6*Q)
+D23 = {"Q": 1, "-6Q": 1, "5Q": 1, "∞": -3}
+
+# (6[Q] + [-6Q] - 7[∞]) - ([Q] + [-6Q] + [5Q] - 3[∞])
+# ==
+# 5[Q] - [5Q] - 4[∞]
+f24 = f22/f23
+D24 = div_sub(D22, D23)
+assert D24 == {"Q": 5, "5Q": -1, "∞": -4}
+
+# Now combine the result
+f = f10*f24
+D = div_add(D10, D24)
+assert D == {
+    "P0": 1,
+    "P1": 2,
+    "P2": 3,
+    "Q":  5,
+    "∞": -11
+}
+
+f_numer = f.numerator().mod(eqn)
+f_denom = f.denominator().mod(eqn)
+# ZeroDivisionError
+#DLog = dlog(f_numer)
+
+assert f(x=P0[0], y=P0[1]) == 0
+assert f(x=P1[0], y=P1[1]) == 0
+assert f(x=P2[0], y=P2[1]) == 0
+# Need to modify f because this is 0/0
+#assert f(x=Q[0], y=Q[1]) == 0
+
+Ps = [P0] + 2*[P1] + 3*[P2] + 5*[Q]
+D = construct_function(Ps)
+
+assert D(x=P0[0], y=P0[1]) == 0
+assert D(x=P1[0], y=P1[1]) == 0
+assert D(x=P2[0], y=P2[1]) == 0
+assert D(x=Q[0], y=Q[1]) == 0
+
+# This will fail due to a bug in sage:
+#DLog = dlog(D)
+# ZeroDivisionError
+DLog = dlog_alt(D)
+