divisor construction algos

script/research/zk/ecip/construct.sage

@@ -0,0 +1,40 @@
+load("div.sage")
+
+def construct(points):
+    divs = []
+    for i in range(0, len(points), 2):
+        # Odd last remainder element
+        if i + 1 == len(divs):
+            P = points[i]
+            L = div_line(P, -P)
+            Q = points[i]
+            divs.append((Q, L))
+            break
+
+        P1, P2 = points[i], points[i + 1]
+        L = div_line(P1, P2)
+        Q = P1 + P2
+        divs.append((Q, L))
+
+    # Now apply reduction algorithm repeatedly
+    while len(divs) != 1:
+        divs2 = []
+
+        if len(divs) % 2 == 1:
+            divs2.append(divs[0])
+            divs = divs[1:]
+
+        for i in range(0, len(divs), 2):
+            assert i + 1 < len(divs)
+            Q1, L1 = divs[i]
+            Q2, L2 = divs[i + 1]
+            ℓ = div_line(Q1, Q2)
+            L = ℓ + L1 + L2 - div_line(Q1, -Q1) - div_line(Q2, -Q2)
+            Q = Q1 + Q2
+            divs2.append((Q, L))
+
+        divs = divs2
+
+    assert len(divs) == 1
+    return divs[0][1]
+

script/research/zk/ecip/example-construct-algo.sage

@@ -0,0 +1,63 @@
+load("div.sage")
+
+# 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")
+M.<z> = L[]
+eqn = y^2 - x^3 - A * x - B
+
+P1 = LabelPoint(E.random_element(), {"P₁": 1})
+P2 = LabelPoint(E.random_element(), {"P₂": 1})
+P3 = LabelPoint(E.random_element(), {"P₃": 1})
+P4 = LabelPoint(E.random_element(), {"P₄": 1})
+P5 = LabelPoint(E.random_element(), {"P₅": 1})
+Q = -(P1.P + P2.P + P3.P + P4.P + P5.P)
+Q = LabelPoint(Q, {"Q": 1})
+assert P1.P + P2.P + P3.P + P4.P + P5.P + Q.P == E(0)
+
+L1 = div_line(P1, P2)
+Q1 = P1 + P2
+print(f"L₁ = {L1}")
+print(f"Q₁ = P₁ + P₂")
+L2 = div_line(P3, P4)
+Q2 = P3 + P4
+print(f"L₂ = {L2}")
+print(f"Q₂ = P₃ + P₄")
+L3 = div_line(P5, Q)
+Q3 = P5 + Q
+print(f"L₃ = {L3}")
+print(f"Q₃ = P₅ + Q")
+print()
+
+ℓ4 = div_line(Q1, Q2)
+L4 = ℓ4 + L1 + L2 - div_line(Q1, -Q1) - div_line(Q2, -Q2)
+Q4 = Q1 + Q2
+print(f"ℓ₄ = {ℓ4}")
+print(f"L₄ = ℓ₄ + L₁ + L₂ - div(x - Q₁) - div(x - Q₂)")
+print(f"   = {L4}")
+print(f"Q₄ = Q₁ + Q₂")
+print("Carry L₃ to next level")
+print()
+
+ℓ5 = div_line(Q4, Q3)
+L5 = ℓ5 + L4 + L3 - div_line(Q4, -Q4) - div_line(Q3, -Q3)
+print(f"ℓ₅ = {ℓ5}")
+print(f"L₅ = ℓ₅ + L₄ + L₃ - div(x - Q₄) - div(x - Q₃)")
+print(f"   = {L5}")
+print()
+
+# We should get the same result here:
+load("construct.sage")
+points = [P1, P2, P3, P4, P5, Q]
+L = construct(points)
+print(f"L  = {L}")
+