ECIP proof

script/research/zk/ecip/div2.sage

@@ -289,5 +289,8 @@ f_numer = D.func.numerator().mod(eqn)
 f_denom = D.func.denominator().mod(eqn)
 f = f_numer / f_denom
 assert f.denominator() == 1
-print(f.numerator())
+f = f.numerator()
+print(f)
+save([P0.P, P1.P, P2.P, Q.P, f], "div.sobj")
+print("Saved (P₀, P₁, P₂, Q, f)")
 

script/research/zk/ecip/proof.sage

@@ -0,0 +1,93 @@
+# First run div2.sage to generate the function and points which
+# are loaded below.
+
+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")
+M.<z> = L[]
+eqn = y^2 - x^3 - A * x - B
+
+def slope_intercept(P1, P2):
+    P1x, P1y = P1.xy()
+    P2x, P2y = P2.xy()
+    P3x, P3y = (-(P1 + P2)).xy()
+    λ = (P2y - P1y) / (P2x - P1x)
+    μ = P2y - λ*P2x
+    return λ, μ
+
+A0 = E.random_element()
+A1 = E.random_element()
+A2 = -(A0 + A1)
+
+[P0, P1, P2, Q, f] = load("div.sobj")
+
+λ, μ = slope_intercept(A0, A1)
+g = y - λ*x - μ
+
+class Func:
+    def __init__(self, func):
+        self.func = func
+
+    def __call__(self, P):
+        Px, Py = P.xy()
+        return self.func(x=Px, y=Py)
+
+f = Func(f)
+g = Func(g)
+
+# We can actually compute this check in a more efficient way
+assert f(A0) * f(A1) * f(A2) == -g(P0) * g(P1)^2 * g(P2)^3 * g(Q)^5
+
+def dlog(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
+
+# Just confirm results from the paper
+assert λ^2 == A0[0] + A1[0] + A2[0]
+
+dy_dx = (3*x^2 + A)/(2*y)
+dx_dz = 1/(dy_dx - λ)
+
+dx_dz = Func(dx_dz)
+assert dx_dz(A0) + dx_dz(A1) + dx_dz(A2) == 0
+
+# These are the actual checks
+Dlog = dlog(f.func)
+# Another way to calculate dlog:
+#D = f.func
+#a_X = K(D(y=0))
+#b_X = K(D(y=1) - a_X)
+#assert D == a_X + y*b_X
+#diff_f = a_X.differentiate(x) + dy_dx*b_X + y*b_X.differentiate(x)
+#assert diff_f == D.differentiate(x) + D.differentiate(y)*dy_dx
+#diff_f = Func(diff_f)
+
+F = Func(Dlog * dx_dz.func)
+G = Func(-1/g.func)
+
+# The prover constructs a proof of this relation being true
+assert F(A0) + F(A1) + F(A2) == G(P0) + 2*G(P1) + 3*G(P2) + 5*G(Q)
+