diff --git a/script/research/zk/circle-stark.sage b/script/research/zk/circle-stark.sage
new file mode 100644
index 000000000..3df76f5ed
--- /dev/null
+++ b/script/research/zk/circle-stark.sage
@@ -0,0 +1,91 @@
+def random_mersenne_prime():
+    while True:
+        p = random_prime(100, 200)
+        m = 2^p - 1
+        if is_prime(m):
+            return m
+
+#p = random_mersenne_prime()
+p = 8191
+
+# -1 should not be a quadratic residue modulo p
+assert legendre_symbol(-1, p) == -1
+assert p % 4 == 3
+
+F = GF(p)
+R.<x> = F[]
+K.<i> = F.extension(x^2 + 1)
+
+def get_point(t):
+    x = (1 - t^2)/(1 + t^2)
+    y = 2*t / (1 + t^2)
+    return (x, y)
+
+(p1_x, p1_y) = get_point(K(3))
+assert p1_x^2 + p1_y^2 == 1
+z = p1_x + i*p1_y
+
+(p2_x, p2_y) = get_point(K(7))
+assert p2_x^2 + p2_y^2 == 1
+w = p2_x + i*p2_y
+
+def abs(z):
+    return z[0]^2 + z[1]^2
+
+# z and w are now elements of F_p(i) which is a field
+assert abs(z) == 1
+assert abs(w) == 1
+assert abs(z * w) == 1
+
+# We can also construct the inverse
+def conjugate(z):
+    return z[0] - i*z[1]
+# Remember that (x + iy)(x - iy) = x^2 - i^2 y^2 = x^2 + y^2
+assert z * conjugate(z) in F
+# Now we find the multiplicative inverse
+z_inv = conjugate(z) / (z * conjugate(z))
+assert z * z_inv == 1
+
+# Size of K is p^2
+assert len(K) == p^2
+
+# Because in sage we cannot construct a homomorphism to GF(p^2) directly
+# we instead construct the isomorphic field extension, and use that instead.
+Fp2.<a> = GF(p^2)
+conway_fp2 = a.minimal_polynomial()
+L.<j> = F.extension(conway_fp2(x=x))
+
+y = L.multiplicative_generator()
+phi = K.hom([y^(y.multiplicative_order()/i.multiplicative_order())])
+# K ≌ GF(p^2)
+assert phi.is_injective() and phi.is_surjective()
+
+assert a.multiplicative_order() == p^2 - 1
+g_K = phi.inverse()(y)
+g1 = g_K^int((p^2 - 1)/(p + 1))
+assert abs(g1) == 1
+
+g2 = K.multiplicative_generator()
+g2 = g2^(p - 1)
+assert abs(g2) == 1
+
+# bug in sage where .unit_group is missing
+# https://ask.sagemath.org/question/62822/make-morphism-from-gfp2s-multiplicative-group-to-gfps-multiplicative-group/
+C = AbelianGroup([p + 1])
+g, = C.gens()
+while True:
+    print(f"Group of order = {g.order()}")
+    for C2 in C.subgroups():
+        print(f"  {C2}")
+    print()
+
+    if g.order() == 1:
+        break
+
+    # For some annoying reason, I cannot iterate on subgroups
+    # C = C.subgroup([g^2])
+    # Trying to get the subgroup of this will give me some error.
+    # Lets just construct it manually.
+    C = AbelianGroup([(g^2).order()])
+    g, = C.gens()
+