diff --git a/script/research/ec/reduction.sage b/script/research/ec/reduction.sage
new file mode 100644
index 000000000..eee2fc057
--- /dev/null
+++ b/script/research/ec/reduction.sage
@@ -0,0 +1,51 @@
+F = GF(47)
+K.<x, y> = F[]
+EC_A, EC_B = 4, 0
+E = EllipticCurve(F, [EC_A, EC_B])
+C = E.defining_polynomial()
+S = K.quotient(C(x, y, 1)).fraction_field()
+X, Y = S(x), S(y)
+
+inf = E[0]
+
+# If we have points sharing the same x value then construct a vertical
+# line through them to eliminate them.
+
+points = [(34, 30), (44, 14), (7, 18), (28, 31), (27, 45), (12, 15),
+          (43, 22), (11, 23), (38, 9), (0, 1),   (26, 33)]
+
+def lagrange_basis(x_k, domain):
+    assert x_k in domain
+    domain = [x_i for x_i in domain if x_i != x_k]
+    assert x_k not in domain
+
+    l = 1
+    for x_i in domain:
+        l *= (x - x_i)
+    l /= l(x_k, 0)
+
+    # Check everything is correct
+    assert l(x_k, 0) == 1
+    for x_i in domain:
+        assert l(x_i, 0) == 0
+    return l
+
+print(f"P = {points}")
+domain = [Px for Px, _ in points]
+
+f = 0
+for Px, Py in points:
+    f += Py * lagrange_basis(Px, domain)
+# Now make it zero at all the y values
+f = y - f
+# Check polynomial is correct
+for Px, Py in points:
+    assert f(Px, Py) == 0
+
+# Now find remaining points in the support
+I = ideal([C(x, y, 1), f])
+V = [(info[x], info[y]) for info in I.variety()]
+print(f"V(I) = {V}")
+diff = set(P) - set(V)
+print(f"diff = {diff}")
+