try to break algo

script/research/zk/fft/fft4.sage

@@ -1,11 +1,5 @@
 import itertools
 
-p = 199
-#n = 16
-n = 8
-
-assert p.is_prime()
-
 def find_ext_order(p, n):
     N = 1
     while True:
@@ -17,7 +11,7 @@ def find_ext_order(p, n):
 
         N += 1
 
-def find_nth_root_unity(K, n):
+def find_nth_root_unity(K, p, N, n):
     # It cannot be a quadratic residue if n is odd
     #assert n % 2 == 1
 
@@ -31,25 +25,7 @@ def find_nth_root_unity(K, n):
 
     return ω
 
-N = find_ext_order(p, n)
-print(f"p = {p}")
-print(f"n = {n}")
-print(f"N = {N}")
-print(f"p^N = {p^N}")
-K.<a> = GF(p^N, repr="int")
-ω = find_nth_root_unity(K, n)
-print(f"ω = {ω}")
-print()
-
-L.<X> = K[]
-
-#f = 9*X^7 + 45*X^6 + 33*X^5 + 7*X^3 + X^2 + 110*X + 4
-f = 7*X^3 + X^2 + 110*X + 4
-assert f.degree() < n/2
-print(f"f = {f}")
-print()
-
-def vectorify(f):
+def vectorify(X, f, n):
     assert f.degree() < n
     fT = vector([f[i] for i in range(f.degree() + 1)] +
                 # Zero padding
@@ -67,7 +43,7 @@ def dot(a, b):
 def alternate(list1, list2):
     return itertools.chain(*zip(list1, list2))
 
-def calc_dft(ω_powers, f):
+def calc_dft(n, ω_powers, f):
     m = len(f)
     indent = " " * (n - m)
     print(f"{indent}calc_dft({ω_powers}, {f})")
@@ -86,18 +62,43 @@ def calc_dft(ω_powers, f):
     print()
 
     ω_powers = vector(ω_i for ω_i in ω_powers[::2])
-    rT = calc_dft(ω_powers, r)
-    sT = calc_dft(ω_powers, s)
+    rT = calc_dft(n, ω_powers, r)
+    sT = calc_dft(n, ω_powers, s)
 
     result = list(alternate(rT, sT))
     print(f"{indent}return {result}")
     return result
 
-ω_powers = vector(ω^i for i in range(n/2))
-fT = vectorify(f)
-dft = calc_dft(ω_powers, fT)
-print()
-print(f"DFT(f) = {dft}")
-f_evals = [f(X=ω^i) for i in range(n)]
-print(f"f(ω^i) = {f_evals}")
+def test():
+    p = 199
+    #n = 16
+    n = 8
+    assert p.is_prime()
+    N = find_ext_order(p, n)
+    print(f"p = {p}")
+    print(f"n = {n}")
+    print(f"N = {N}")
+    print(f"p^N = {p^N}")
+    K.<a> = GF(p^N, repr="int")
+    ω = find_nth_root_unity(K, p, N, n)
+    print(f"ω = {ω}")
+    print()
+
+    L.<X> = K[]
+
+    #f = 9*X^7 + 45*X^6 + 33*X^5 + 7*X^3 + X^2 + 110*X + 4
+    f = 7*X^3 + X^2 + 110*X + 4
+    assert f.degree() < n/2
+    print(f"f = {f}")
+    print()
+
+    ω_powers = vector(ω^i for i in range(n/2))
+    fT = vectorify(X, f, n)
+    dft = calc_dft(n, ω_powers, fT)
+    print()
+    print(f"DFT(f) = {dft}")
+    f_evals = [f(X=ω^i) for i in range(n)]
+    print(f"f(ω^i) = {f_evals}")
+
+test()
 

script/research/zk/fft/fft5.sage

@@ -0,0 +1,128 @@
+import itertools
+
+def find_ext_order(p, n):
+    N = 1
+    while True:
+        pNx_order = p^N - 1
+
+        # Does n divide the group order 𝔽_{p^N}^×?
+        if pNx_order % n == 0:
+            return N
+
+        N += 1
+
+def find_nth_root_unity(K, p, N, n):
+    # It cannot be a quadratic residue if n is odd
+    #assert n % 2 == 1
+
+    # So there is an nth root of unity in p^N. Now we have to find it.
+    pNx_order = p^N - 1
+
+    ω = K.gens()[0]
+    ω = ω^(pNx_order/n)
+    assert ω^n == 1
+    assert ω^(n - 1) != 1
+
+    return ω
+
+def vectorify(X, f, n):
+    assert f.degree() < n
+    fT = vector([f[i] for i in range(f.degree() + 1)] +
+                # Zero padding
+                [0 for _ in range(n - f.degree() - 1)])
+    assert len(fT) == n
+    # Just check decomposed polynomial is in the correct order
+    assert sum([fT[i]*X^i for i in range(n)]) == f
+    return fT
+
+def dot(a, b):
+    assert len(a) == len(b)
+    return [a_i*b_i for a_i, b_i in zip(a, b)]
+
+# ABC, DEF -> ADBECF
+def alternate(list1, list2):
+    return itertools.chain(*zip(list1, list2))
+
+def calc_dft(n, ω_powers, f):
+    m = len(f)
+    if m == 1:
+        return f
+    g, h = vector(f[:m/2]), vector(f[m/2:])
+
+    r = g + h
+    s = dot(g - h, ω_powers)
+
+    ω_powers = vector(ω_i for ω_i in ω_powers[::2])
+    rT = calc_dft(n, ω_powers, r)
+    sT = calc_dft(n, ω_powers, s)
+
+    result = list(alternate(rT, sT))
+    return result
+
+def test1():
+    p = 199
+    #n = 16
+    n = 8
+    assert p.is_prime()
+    N = find_ext_order(p, n)
+    print(f"p = {p}")
+    print(f"n = {n}")
+    print(f"N = {N}")
+    print(f"p^N = {p^N}")
+    K.<a> = GF(p^N, repr="int")
+    ω = find_nth_root_unity(K, p, N, n)
+    print(f"ω = {ω}")
+    print()
+
+    L.<X> = K[]
+
+    #f = 9*X^7 + 45*X^6 + 33*X^5 + 7*X^3 + X^2 + 110*X + 4
+    f = 7*X^3 + X^2 + 110*X + 4
+    assert f.degree() < n/2
+    print(f"f = {f}")
+    print()
+
+    ω_powers = vector(ω^i for i in range(n/2))
+    fT = vectorify(X, f, n)
+    dft = calc_dft(n, ω_powers, fT)
+    print()
+    print(f"DFT(f) = {dft}")
+    f_evals = [f(X=ω^i) for i in range(n)]
+    print(f"f(ω^i) = {f_evals}")
+
+def random_test():
+    p = random_prime(1000)
+    #n = 16
+    n = int(2^ZZ.random_element(2, 10))
+    assert p.is_prime()
+    N = find_ext_order(p, n)
+    print(f"p = {p}")
+    print(f"n = {n}")
+    print(f"N = {N}")
+    print(f"p^N = {p^N}")
+    K.<a> = GF(p^N, repr="int")
+    ω = find_nth_root_unity(K, p, N, n)
+    print(f"ω = {ω}")
+    print()
+
+    L.<X> = K[]
+
+    #f = 9*X^7 + 45*X^6 + 33*X^5 + 7*X^3 + X^2 + 110*X + 4
+    f = 0
+    for i in range(n/2):
+        f += ZZ.random_element(0, 200) * X^i
+    assert f.degree() < n/2
+    print(f"f = {f}")
+    print()
+
+    ω_powers = vector(ω^i for i in range(n/2))
+    fT = vectorify(X, f, n)
+    dft = calc_dft(n, ω_powers, fT)
+    print()
+    print(f"DFT(f) = {dft}")
+    f_evals = [f(X=ω^i) for i in range(n)]
+    print(f"f(ω^i) = {f_evals}")
+
+#test1()
+random_test()
+