diff --git a/src/circle_fft.ipynb b/src/circle_fft.ipynb
index 341ed16..d70fa9c 100644
--- a/src/circle_fft.ipynb
+++ b/src/circle_fft.ipynb
@@ -1 +1 @@
-{"cells":[{"cell_type":"markdown","id":"8816ed9e","metadata":{},"source":["## Circle Group\n","\n","Before we start building circle FFT, let's first build the circle group.\n","\n","We choose a basic field GF(31), which has a order of 31. And then we extend it to a 2-degree extension field GF(31^2), to get a multiplicative subgroup of order 32."]},{"cell_type":"code","execution_count":1,"id":"2c77b9d4-0cab-4100-8eae-bfded3cf543d","metadata":{},"outputs":[{"data":{"text/plain":["Univariate Polynomial Ring in X over Finite Field of size 31"]},"execution_count":1,"metadata":{},"output_type":"execute_result"}],"source":["F31 = GF(31)\n","R31.<X> = F31[]\n","R31"]},{"cell_type":"markdown","id":"a64c4622","metadata":{},"source":["Now we extend F31 to a 2-degree extension field GF(31^2)."]},{"cell_type":"code","execution_count":2,"id":"f25afd7d-90a4-4c10-bcf9-ba392e61d69d","metadata":{},"outputs":[{"data":{"text/plain":["Finite Field in i of size 31^2"]},"execution_count":2,"metadata":{},"output_type":"execute_result"}],"source":["C31.<i> = F31.extension(X^2 + 1)\n","C31"]},{"cell_type":"markdown","id":"b9bc3979-ea87-4c31-9a4e-5b5a2d55f252","metadata":{},"source":["We want to find a generator of order 32, so we take the generator of GF(31^2) and get 30th power of it.\n","\n","**Fact**: For any prime field $\\mathbb{F}_p$, where $p=2^n - 1$, the 2-degree extension over $\\mathbb{F}_p$ is \n","$K=\\mathbb{F}_{p^2}$, then $K$ has a multiplicative subgroup with order $2^n$.\n","\n","Proof:\n","\n","$$\n","p^2 - 1 = (2^n - 1) * (2^n - 1) - 1 = 2^{2n} - 2^{n+1} = 2^n(2^n - 2)\n","$$"]},{"cell_type":"code","execution_count":3,"id":"f8f575a0-4a95-43e7-b0f7-ad2b7cd42103","metadata":{},"outputs":[{"data":{"text/plain":["11*i + 29"]},"execution_count":3,"metadata":{},"output_type":"execute_result"}],"source":["\n","g = C31.multiplicative_generator()^(F31.order() - 1)\n","g"]},{"cell_type":"markdown","id":"ca647bd5","metadata":{},"source":["Define a function to calculate the order of an element in the circle group."]},{"cell_type":"code","execution_count":4,"id":"0089303e-1e4c-440d-a7c2-fc791b4d23cf","metadata":{},"outputs":[],"source":["def group_ord(elem):\n","    K = elem.parent()\n","\n","    for i in range(1, K.order()):\n","        if elem**i == K.one():\n","            return i\n","    return \"False\""]},{"cell_type":"markdown","id":"609163d2","metadata":{},"source":["Check that g is of order 32."]},{"cell_type":"code","execution_count":5,"id":"6de6738e-e132-4825-a0ff-3185089d6a17","metadata":{},"outputs":[{"data":{"text/plain":["32"]},"execution_count":5,"metadata":{},"output_type":"execute_result"}],"source":["group_ord(g)"]},{"cell_type":"markdown","id":"743f8247","metadata":{},"source":["Define the group multiplication, which is similar to complex number multiplication."]},{"cell_type":"code","execution_count":6,"id":"98f66494-5eca-42f9-9150-32d83cd8ce0f","metadata":{},"outputs":[],"source":["def group_mul(g1, g2):\n","    x1, y1 = g1\n","    x2, y2 = g2\n","    return (x1 * x2 - y1 * y2) + (x1 * y2 + x2 * y1) * i "]},{"cell_type":"markdown","id":"0ecf4e47","metadata":{},"source":["Define the group inverse, which is also the J mapping."]},{"cell_type":"code","execution_count":7,"id":"10780cbd-3786-4642-8851-13642de9f72e","metadata":{},"outputs":[],"source":["def group_inv(g1):\n","    x1, y1 = g1\n","    return (x1 - y1 * i)"]},{"cell_type":"markdown","id":"03439c18","metadata":{},"source":["Define show function to print group."]},{"cell_type":"code","execution_count":8,"id":"98046a5b-d75a-40fe-90e9-c4a9e3940512","metadata":{},"outputs":[],"source":["def show(GG): \n","    for j in range(0, len(GG)):\n","        t = GG[j]\n","        t0 = t[0]\n","        t1 = t[1]\n","        if t1 > F31(15) and t1 <= F31(30):\n","            t1_str = \"-\" + str(31-t1)\n","        else:\n","            t1_str = str(t1)\n","        print(\"i={}, log={}, ord={}, ({},{})\".format(j, log_gen(g, t), group_ord(t), t0, t1_str))"]},{"cell_type":"markdown","id":"6249f760","metadata":{},"source":["Define log_gen function to calculate the power of an element in the circle group."]},{"cell_type":"code","execution_count":9,"id":"e3aa4193-e683-4834-a18c-c835b718a878","metadata":{},"outputs":[],"source":["def log_gen(gen, t):\n","    K = gen.parent()\n","    for i in range(1, K.order()):\n","        if gen^i == t:\n","            return i\n","    return \"Fail\""]},{"cell_type":"markdown","id":"40185187-0f6f-40f5-9869-5611576e5ce5","metadata":{},"source":["## Jump over the circle\n","\n","Here we define the generator for `Circle Group` as chosen by Vitalik https://vitalik.eth.limo/general/2024/07/23/circlestarks.html"]},{"cell_type":"code","execution_count":12,"id":"e96e1eba-4ef3-4f3a-9003-43d090044a96","metadata":{},"outputs":[{"data":{"text/plain":["32"]},"execution_count":12,"metadata":{},"output_type":"execute_result"}],"source":["# optional\n","g = 10 + 5*i\n","group_ord(g)"]},{"cell_type":"code","execution_count":13,"id":"bfd12d98-21e7-4b48-9f6c-a1a7d7e4aee2","metadata":{},"outputs":[{"data":{"text/plain":["7*i + 13"]},"execution_count":13,"metadata":{},"output_type":"execute_result"}],"source":["g * g"]},{"cell_type":"code","execution_count":14,"id":"77e0b387-b325-4cf4-8c82-cc07368b21b3","metadata":{},"outputs":[{"data":{"text/plain":["[1,\n"," 5*i + 10,\n"," 7*i + 13,\n"," 11*i + 2,\n"," 27*i + 27,\n"," 2*i + 11,\n"," 13*i + 7,\n"," 10*i + 5,\n"," i,\n"," 10*i + 26,\n"," 13*i + 24,\n"," 2*i + 20,\n"," 27*i + 4,\n"," 11*i + 29,\n"," 7*i + 18,\n"," 5*i + 21,\n"," 30,\n"," 26*i + 21,\n"," 24*i + 18,\n"," 20*i + 29,\n"," 4*i + 4,\n"," 29*i + 20,\n"," 18*i + 24,\n"," 21*i + 26,\n"," 30*i,\n"," 21*i + 5,\n"," 18*i + 7,\n"," 29*i + 11,\n"," 4*i + 27,\n"," 20*i + 2,\n"," 24*i + 13,\n"," 26*i + 10]"]},"execution_count":14,"metadata":{},"output_type":"execute_result"}],"source":["G5 = [g^i for i in range(0, 2^5)]\n","# G5's elements must be on the circle x^2 + y^2 = 1\n","# due to:\n","#   a^(p + 1) == 1, a = x + y * i\n","#   a * a ^ p == (x + y * i) * (x - y * i) == x^2 + y^2 == 1\n","for t in G5:\n","    x, y = t\n","    assert x^2 + y^2 == 1\n","G5"]},{"cell_type":"code","execution_count":15,"id":"d138513c-a1a8-48fa-a93f-d50b131c2a2b","metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["i=0, log=32, ord=1, (1,0)\n","i=1, log=1, ord=32, (10,5)\n","i=2, log=2, ord=16, (13,7)\n","i=3, log=3, ord=32, (2,11)\n","i=4, log=4, ord=8, (27,27)\n","i=5, log=5, ord=32, (11,2)\n","i=6, log=6, ord=16, (7,13)\n","i=7, log=7, ord=32, (5,10)\n","i=8, log=8, ord=4, (0,1)\n","i=9, log=9, ord=32, (26,10)\n","i=10, log=10, ord=16, (24,13)\n","i=11, log=11, ord=32, (20,2)\n","i=12, log=12, ord=8, (4,27)\n","i=13, log=13, ord=32, (29,11)\n","i=14, log=14, ord=16, (18,7)\n","i=15, log=15, ord=32, (21,5)\n","i=16, log=16, ord=2, (30,0)\n","i=17, log=17, ord=32, (21,26)\n","i=18, log=18, ord=16, (18,24)\n","i=19, log=19, ord=32, (29,20)\n","i=20, log=20, ord=8, (4,4)\n","i=21, log=21, ord=32, (20,29)\n","i=22, log=22, ord=16, (24,18)\n","i=23, log=23, ord=32, (26,21)\n","i=24, log=24, ord=4, (0,30)\n","i=25, log=25, ord=32, (5,21)\n","i=26, log=26, ord=16, (7,18)\n","i=27, log=27, ord=32, (11,29)\n","i=28, log=28, ord=8, (27,4)\n","i=29, log=29, ord=32, (2,20)\n","i=30, log=30, ord=16, (13,24)\n","i=31, log=31, ord=32, (10,26)\n"]}],"source":["for j in range(0, 32):\n","    t = g^j\n","    print(\"i={}, log={}, ord={}, ({},{})\".format(j, log_gen(g, t), group_ord(t), t[0], t[1]))"]},{"cell_type":"code","execution_count":16,"id":"331e689e-b199-4717-ac08-266f1a52e82c","metadata":{},"outputs":[{"data":{"text/plain":["True"]},"execution_count":16,"metadata":{},"output_type":"execute_result"}],"source":["G5[1] * G5[1] == G5[2]"]},{"cell_type":"code","execution_count":17,"id":"e91a4d9e-3d07-480b-9505-6af79f03ad0b","metadata":{},"outputs":[{"data":{"text/plain":["True"]},"execution_count":17,"metadata":{},"output_type":"execute_result"}],"source":["G5[3] * G5[4] == G5[7]"]},{"cell_type":"code","execution_count":18,"id":"e5695c4e-4e25-491e-80fb-cda2061dc62f","metadata":{},"outputs":[{"data":{"text/plain":["(10, 5)"]},"execution_count":18,"metadata":{},"output_type":"execute_result"}],"source":["x1, y1 = G5[1]\n","x1, y1"]},{"cell_type":"code","execution_count":19,"id":"091fc1f5-9ea7-4855-a3c0-945ea384650f","metadata":{},"outputs":[{"data":{"text/plain":["True"]},"execution_count":19,"metadata":{},"output_type":"execute_result"}],"source":["group_mul(G5[1], G5[2]) == G5[3]"]},{"cell_type":"code","execution_count":20,"id":"6967f7e3-1adb-4498-96ae-6969ce0cfcf5","metadata":{},"outputs":[{"data":{"text/plain":["(True, 1)"]},"execution_count":20,"metadata":{},"output_type":"execute_result"}],"source":["group_inv(G5[1]) ==  G5[31], G5[1] * group_inv(G5[1])"]},{"cell_type":"code","execution_count":21,"id":"8ed48611-82dc-47c6-8655-2fd1fdf4e074","metadata":{},"outputs":[{"data":{"text/plain":["[1,\n"," 7*i + 13,\n"," 27*i + 27,\n"," 13*i + 7,\n"," i,\n"," 13*i + 24,\n"," 27*i + 4,\n"," 7*i + 18,\n"," 30,\n"," 24*i + 18,\n"," 4*i + 4,\n"," 18*i + 24,\n"," 30*i,\n"," 18*i + 7,\n"," 4*i + 27,\n"," 24*i + 13]"]},"execution_count":21,"metadata":{},"output_type":"execute_result"}],"source":["G4 = [t^2 for t in G5[:16]]\n","G4"]},{"cell_type":"code","execution_count":22,"id":"673480c6-0531-46d7-a57f-950f5b46edd8","metadata":{},"outputs":[{"data":{"text/plain":["[1, 27*i + 27, i, 27*i + 4, 30, 4*i + 4, 30*i, 4*i + 27]"]},"execution_count":22,"metadata":{},"output_type":"execute_result"}],"source":["G3 = [t^2 for t in G4[:8]]\n","G3"]},{"cell_type":"markdown","id":"ce8e103e-2e1d-495e-b129-7e441455a3ed","metadata":{},"source":["## Standard Position Coset\n","\n","$D$ is a standard position coset of size $2^n$\n","\n","$$\n","D = Q\\cdot G_n\n","$$\n","\n","where $ord(Q)=2^{n+1}$.\n","\n","And $D$ is also a twin-coset of size $2^n$ defined with $G_{n-1}$:\n","\n","$$\n","D = Q\\cdot G_{n-1} \\uplus Q^{-1}\\cdot G_{n-1}\n","$$\n","\n","The rotation $Q^2$ is exactly the generator of $G_n$."]},{"cell_type":"code","execution_count":23,"id":"396deee6-21a2-4f38-8969-f620664bdd6b","metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["i=0, log=1, ord=32, (10,5)\n","i=1, log=3, ord=32, (2,11)\n","i=2, log=5, ord=32, (11,2)\n","i=3, log=7, ord=32, (5,10)\n","i=4, log=9, ord=32, (26,10)\n","i=5, log=11, ord=32, (20,2)\n","i=6, log=13, ord=32, (29,11)\n","i=7, log=15, ord=32, (21,5)\n","i=8, log=17, ord=32, (21,-5)\n","i=9, log=19, ord=32, (29,-11)\n","i=10, log=21, ord=32, (20,-2)\n","i=11, log=23, ord=32, (26,-10)\n","i=12, log=25, ord=32, (5,-10)\n","i=13, log=27, ord=32, (11,-2)\n","i=14, log=29, ord=32, (2,-11)\n","i=15, log=31, ord=32, (10,-5)\n"]}],"source":["# Define a coset or G4 from a rotation Q, ord(Q) = 32\n","#\n","# Let Q = g (generator of G, whose order is 32)\n","\n","D_standard_coset_4 = [g * t for t in G4]; show(D_standard_coset_4)"]},{"cell_type":"code","execution_count":24,"id":"f24c075e-89eb-4db5-a533-c331b8f91a4c","metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["i=0, log=1, ord=32, (10,5)\n","i=1, log=3, ord=32, (2,11)\n","i=2, log=5, ord=32, (11,2)\n","i=3, log=7, ord=32, (5,10)\n","i=4, log=9, ord=32, (26,10)\n","i=5, log=11, ord=32, (20,2)\n","i=6, log=13, ord=32, (29,11)\n","i=7, log=15, ord=32, (21,5)\n","i=8, log=17, ord=32, (21,-5)\n","i=9, log=19, ord=32, (29,-11)\n","i=10, log=21, ord=32, (20,-2)\n","i=11, log=23, ord=32, (26,-10)\n","i=12, log=25, ord=32, (5,-10)\n","i=13, log=27, ord=32, (11,-2)\n","i=14, log=29, ord=32, (2,-11)\n","i=15, log=31, ord=32, (10,-5)\n"]}],"source":["# D_twin_coset_4 is equal to D_standard_coset_4 in the set-theoretic sense, with different order\n","\n","D_twin_coset_4 = [g * t for t in G3] + [group_inv(g) * t for t in G3]; show(D_standard_coset_4)\n","\n","# `log` specifies the order of an element in the circle"]},{"cell_type":"code","execution_count":25,"id":"12dc1f2a-d03c-46ea-bedc-217e53c9cb9b","metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["i=0, log=3, ord=32, (2,11)\n","i=1, log=7, ord=32, (5,10)\n","i=2, log=11, ord=32, (20,2)\n","i=3, log=15, ord=32, (21,5)\n","i=4, log=19, ord=32, (29,-11)\n","i=5, log=23, ord=32, (26,-10)\n","i=6, log=27, ord=32, (11,-2)\n","i=7, log=31, ord=32, (10,-5)\n","i=8, log=29, ord=32, (2,-11)\n","i=9, log=1, ord=32, (10,5)\n","i=10, log=5, ord=32, (11,2)\n","i=11, log=9, ord=32, (26,10)\n","i=12, log=13, ord=32, (29,11)\n","i=13, log=17, ord=32, (21,-5)\n","i=14, log=21, ord=32, (20,-2)\n","i=15, log=25, ord=32, (5,-10)\n"]}],"source":["D_twin_coset_4_1 = [G5[3] * t for t in G3] + [group_inv(G5[3]) * t for t in G3]; show(D_twin_coset_4_1)"]},{"cell_type":"code","execution_count":26,"id":"b08f58a7-52ce-48c0-91a1-6e692344c627","metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["i=0, log=5, ord=32, (11,2)\n","i=1, log=9, ord=32, (26,10)\n","i=2, log=13, ord=32, (29,11)\n","i=3, log=17, ord=32, (21,-5)\n","i=4, log=21, ord=32, (20,-2)\n","i=5, log=25, ord=32, (5,-10)\n","i=6, log=29, ord=32, (2,-11)\n","i=7, log=1, ord=32, (10,5)\n","i=8, log=27, ord=32, (11,-2)\n","i=9, log=31, ord=32, (10,-5)\n","i=10, log=3, ord=32, (2,11)\n","i=11, log=7, ord=32, (5,10)\n","i=12, log=11, ord=32, (20,2)\n","i=13, log=15, ord=32, (21,5)\n","i=14, log=19, ord=32, (29,-11)\n","i=15, log=23, ord=32, (26,-10)\n"]}],"source":["D_twin_coset_4_2 = [G5[5] * t for t in G3] + [group_inv(G5[5]) * t for t in G3]; show(D_twin_coset_4_2)"]},{"cell_type":"code","execution_count":27,"id":"17b78d53-4cb2-446b-ac00-1006df2568d9","metadata":{},"outputs":[{"data":{"text/plain":["True"]},"execution_count":27,"metadata":{},"output_type":"execute_result"}],"source":["# Fact:\n","#   g^2 is the generator of G4\n","#   g * G4 is the standard position coset\n","\n","G4 == [(g^2)^i for i in range(0, 16)]"]},{"cell_type":"code","execution_count":28,"id":"3976036e-e166-4084-adb6-28594cc4d1f5","metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["i=0, log=2, ord=16, (13,7)\n","i=1, log=6, ord=16, (7,13)\n","i=2, log=10, ord=16, (24,13)\n","i=3, log=14, ord=16, (18,7)\n","i=4, log=18, ord=16, (18,-7)\n","i=5, log=22, ord=16, (24,-13)\n","i=6, log=26, ord=16, (7,-13)\n","i=7, log=30, ord=16, (13,-7)\n"]}],"source":["# Define a coset or G3 from a rotation Q, ord(Q) = 16\n","#\n","# Let Q = g^2 (g is the generator of G5, whose order is 32)\n","\n","D_standard_coset_3 = [g^2 * t for t in G3]; show(D_standard_coset_3)"]},{"cell_type":"markdown","id":"261d3556-cc74-4961-b271-88b03fe3045c","metadata":{},"source":["## Squaring map\n","\n","**Lemma**: For any twin-coset $D$ of size $N=2^n$, the image of $\\pi(D)$ must be a twin-coset of size $2^{n-1}$.\n","\n","$$\n","D = Q\\cdot G_{n-1} \\uplus Q^{-1}\\cdot G_{n-1}\n","$$\n","\n","$$\n","\\pi(D) = \\pi(Q)\\cdot G_{n-2} \\uplus \\pi(Q)^{-1}\\cdot G_{n-2}\n","$$\n","\n","**Fact**: If $D$ is a standard position coset, so is $\\pi(D)$."]},{"cell_type":"code","execution_count":29,"id":"afd781ca-8a1c-4a2a-af97-cfc8a2b63179","metadata":{},"outputs":[],"source":["def sq(D):\n","    rs = []\n","    for t in D:\n","        if t^2 not in rs:\n","            rs += [t^2]  \n","    return rs"]},{"cell_type":"code","execution_count":30,"id":"99e26070-ecd6-4795-a54d-1721e9319209","metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["i=0, log=2, ord=16, (13,7)\n","i=1, log=6, ord=16, (7,13)\n","i=2, log=10, ord=16, (24,13)\n","i=3, log=14, ord=16, (18,7)\n","i=4, log=18, ord=16, (18,-7)\n","i=5, log=22, ord=16, (24,-13)\n","i=6, log=26, ord=16, (7,-13)\n","i=7, log=30, ord=16, (13,-7)\n"]}],"source":["D_standard_coset_3 = sq(D_standard_coset_4); show(D_standard_coset_3)"]},{"cell_type":"markdown","id":"78a82ea8-1720-4a3d-9a81-5dedbf918502","metadata":{},"source":["## Vanishing Polynomial"]},{"cell_type":"code","execution_count":31,"id":"bdeed1fc-f78c-4315-a8ba-62869b20137a","metadata":{},"outputs":[{"data":{"text/plain":["[10, 13, 27, 0, 30]"]},"execution_count":31,"metadata":{},"output_type":"execute_result"}],"source":["v = [x]\n","for idx in range(1, 5):\n","    v.append(2 * v[idx-1] ^ 2 - 1)\n","\n","v"]},{"cell_type":"markdown","id":"89ebe6b9","metadata":{},"source":["## Circle FFT"]},{"cell_type":"code","execution_count":32,"id":"1414592f","metadata":{},"outputs":[{"data":{"text/plain":["{5*i + 10: 3,\n"," 11*i + 2: 28,\n"," 2*i + 11: 15,\n"," 10*i + 5: 7,\n"," 10*i + 26: 29,\n"," 2*i + 20: 0,\n"," 11*i + 29: 2,\n"," 5*i + 21: 5,\n"," 26*i + 21: 20,\n"," 20*i + 29: 15,\n"," 29*i + 20: 0,\n"," 21*i + 26: 6,\n"," 21*i + 5: 12,\n"," 29*i + 11: 16,\n"," 20*i + 2: 6,\n"," 26*i + 10: 2}"]},"execution_count":32,"metadata":{},"output_type":"execute_result"}],"source":["f = {}\n","from random import randint\n","\n","for idx, t in enumerate(D_standard_coset_4):\n","    f[t] = C31(randint(0, 100))\n","\n","f"]},{"cell_type":"code","execution_count":33,"id":"0f1f54b5","metadata":{},"outputs":[],"source":["def pi_value(t):\n","    # x^2 - y^2 == 2 * x^2 - 1 (x^2 + y^2 = 1)\n","    return C31(2 * t^2 - 1)"]},{"cell_type":"code","execution_count":34,"id":"2a87e2ae","metadata":{},"outputs":[],"source":["def ifft_first_step(f):\n","    f0 = {}\n","    f1 = {}\n","    for t in f:\n","        x, y = t\n","\n","        f0[x] = (f[t] + f[group_inv(t)]) / C31(2)\n","        f1[x] = (f[t] - f[group_inv(t)]) / (C31(2) * y)\n","\n","        assert f[t] == f0[x] + y * f1[x]\n","\n","        # print(\"{}: {} = {} + {} * {}\".format(t, f[t], f0[x], f1[x], y))\n","\n","    return f0, f1"]},{"cell_type":"code","execution_count":35,"id":"936da672","metadata":{},"outputs":[],"source":["def ifft_normal_step(f, debug_output=None):\n","    if debug_output:\n","        debug_output.append(f)\n","\n","    if len(f) == 1:\n","        res = []\n","        for x in f:\n","            res.append(f[x])\n","        return res\n","\n","    f0 = {}\n","    f1 = {}\n","\n","    for x in f:\n","        assert x != 0, \"f should be on coset\"\n","        f0[pi_value(x)] = (f[x] + f[-x]) / C31(2)\n","        f1[pi_value(x)] = (f[x] - f[-x]) / (C31(2) * x)\n","        # print('f[{}] - f[{}] = {} - {}'.format(x, -x, f[x], f[-x]))\n","\n","        assert f[x] == f0[pi_value(x)] + x * f1[pi_value(x)]\n","\n","        # print(\"x={}: {} = {} + {} * {}\".format(x, f[x], f0[pi(x)], f1[pi(x)], x))\n","        # print(\"pi(x)={}\".format(pi(x)))\n","\n","    return ifft_normal_step(f0, debug_output=debug_output) + ifft_normal_step(f1, debug_output=debug_output)"]},{"cell_type":"code","execution_count":36,"id":"b8b62ba5","metadata":{},"outputs":[],"source":["def ifft(f, debug_output=None):\n","    f0, f1 = ifft_first_step(f)\n","    if debug_output:\n","        debug_output.append(f0)\n","        debug_output.append(f1)\n","    f0 = ifft_normal_step(f0, debug_output=debug_output)\n","    f1 = ifft_normal_step(f1, debug_output=debug_output)\n","\n","    return f0 + f1"]},{"cell_type":"code","execution_count":37,"id":"f7349218","metadata":{},"outputs":[{"data":{"text/plain":["[22, 30, 1, 17, 19, 27, 2, 14, 24, 7, 30, 10, 22, 20, 26, 4]"]},"execution_count":37,"metadata":{},"output_type":"execute_result"}],"source":["coeffs = ifft(f)\n","coeffs"]},{"cell_type":"code","execution_count":38,"id":"2ea6a0e0","metadata":{},"outputs":[],"source":["def pie_group(D):\n","    D_new = []\n","    for x in D:\n","        x_new = C31(2 * x^2 - 1)\n","        if x_new not in D_new:\n","            D_new.append(x_new)\n","\n","    assert len(D_new) * 2 == len(D), \"len(D_new) * 2 != len(D), {} * 2 != {}, D_new={}, D={}\".format(len(D_new), len(D), D_new, D)\n","    \n","    return D_new"]},{"cell_type":"code","execution_count":39,"id":"d509cd22","metadata":{},"outputs":[],"source":["def fft_first_step(f, D):\n","    assert len(f) == len(D), \"len(f) != len(D), {} != {}, f={}, D={}\".format(len(f), len(D), f, D)\n","\n","    len_f = len(f)\n","    f0 = f[:len_f//2]\n","    f1 = f[len_f//2:]\n","\n","    D_new = []\n","    for t in D:\n","        x, _ = t\n","        if x not in D_new:\n","            D_new.append(x)\n","\n","    assert len(D_new) * 2 == len(D), \"len(D_new) * 2 != len(D), {} * 2 != {}, D_new={}, D={}\".format(len(D_new), len(D), D_new, D)\n","\n","    return f0, f1, D_new"]},{"cell_type":"code","execution_count":40,"id":"4413ffd8","metadata":{},"outputs":[],"source":["def fft_normal_step(f, D, debug_output=None):\n","    if debug_output:\n","        debug_output.append(f)\n","\n","    if len(f) == 1:\n","        return {D[0]: f[0]}\n","    \n","    next_domain = pie_group(D)\n","    assert len(next_domain) * 2 == len(D), \"len(next_domain) * 2 != len(D), {} * 2 != {}, next_domain={}, D={}\".format(len(next_domain), len(D), next_domain, D)\n","    assert len(f) == len(D), \"len(f) != len(D), {} != {}, f={}, D={}\".format(len(f), len(D), f, D)\n","    f0 = fft_normal_step(f[:len(f)//2], next_domain, debug_output)\n","    f1 = fft_normal_step(f[len(f)//2:], next_domain, debug_output)\n","\n","    f_new = {}\n","    for x in D:\n","        f_new[x] = f0[pi_value(x)] + f1[pi_value(x)] * x\n","\n","    for x in D:\n","        if x != 0:\n","            assert f0[pi_value(x)] == (f_new[x] + f_new[-x]) / C31(2), \"f0[pi(x)] = {}\".format(f0[pi_value(x)])\n","            assert f1[pi_value(x)] == (f_new[x] - f_new[-x]) / (C31(2) * x), \"f1[pi(x)] = {}\".format(f1[pi_value(x)])\n","        else:\n","            assert f0[pi_value(x)] == f_new[x], \"f0[pi(x)] = {}\".format(f0[pi_value(x)])\n","\n","    assert len(f) == len(f_new), \"len(f) != len(f_new), {} != {}, f={}, f_new={}, D={}\".format(len(f), len(f_new), f, f_new, D)\n","    assert ifft_normal_step(f_new) == f, \"ifft(f_new) != f, {} != {}\".format(ifft_normal_step(f_new), f)\n","        \n","    return f_new"]},{"cell_type":"code","execution_count":41,"id":"2b5ce548","metadata":{},"outputs":[],"source":["def fft(f, D, debug_output=None):\n","\n","    assert len(f) == len(D), \"len(f) != len(D), {} != {}, f={}, D={}\".format(len(f), len(D), f, D)\n","\n","    D_copy = D[:]\n","    f0, f1, D = fft_first_step(f, D)\n","\n","    assert len(f0) == len(D), \"len(f0) != len(D), {} != {}, f0={}, D={}\".format(len(f0), len(D), f0, D)\n","    assert len(f1) == len(D), \"len(f1) != len(D), {} != {}, f1={}, D={}\".format(len(f1), len(D), f1, D)\n","\n","    if debug_output:\n","        debug_output.append(f0)\n","        debug_output.append(f1)\n","\n","    f0 = fft_normal_step(f0, D, debug_output)\n","    f1 = fft_normal_step(f1, D, debug_output)\n","\n","    f = {}\n","    for t in D_copy:\n","        x, y = t\n","        f[t] = f0[x] + f1[x] * y\n","\n","    return f\n"]},{"cell_type":"code","execution_count":42,"id":"89d2a2c7","metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["f[5*i + 10]=3, f_prime[5*i + 10]=3\n","f[11*i + 2]=28, f_prime[11*i + 2]=28\n","f[2*i + 11]=15, f_prime[2*i + 11]=15\n","f[10*i + 5]=7, f_prime[10*i + 5]=7\n","f[10*i + 26]=29, f_prime[10*i + 26]=29\n","f[2*i + 20]=0, f_prime[2*i + 20]=0\n","f[11*i + 29]=2, f_prime[11*i + 29]=2\n","f[5*i + 21]=5, f_prime[5*i + 21]=5\n","f[26*i + 21]=20, f_prime[26*i + 21]=20\n","f[20*i + 29]=15, f_prime[20*i + 29]=15\n","f[29*i + 20]=0, f_prime[29*i + 20]=0\n","f[21*i + 26]=6, f_prime[21*i + 26]=6\n","f[21*i + 5]=12, f_prime[21*i + 5]=12\n","f[29*i + 11]=16, f_prime[29*i + 11]=16\n","f[20*i + 2]=6, f_prime[20*i + 2]=6\n","f[26*i + 10]=2, f_prime[26*i + 10]=2\n"]}],"source":["f_prime = fft(coeffs, D_standard_coset_4)\n","\n","for t, s in zip(f, f_prime):\n","    assert t == s\n","\n","for t in f:\n","    print(\"f[{}]={}, f_prime[{}]={}\".format(t, f[t], t, f_prime[t]))\n","    assert f[t] == f_prime[t], \"f[{}] != f_prime[{}], {} != {}\".format(t, t, f[t], f_prime[t])"]}],"metadata":{"kernelspec":{"display_name":"SageMath 9.8","language":"sage","name":"SageMath-9.8"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.11.1"}},"nbformat":4,"nbformat_minor":5}
+{"cells":[{"attachments":{},"cell_type":"markdown","id":"8816ed9e","metadata":{},"source":["## Circle Group\n","\n","Before we start building circle FFT, let's first build the circle group, which will be used as the domain of circle fft.\n","\n","We choose a basic field GF(31), which has a order of 31. And then we extend it to a 2-degree extension field GF($31^2$), to get a multiplicative subgroup of order 32."]},{"cell_type":"code","execution_count":1,"id":"2c77b9d4-0cab-4100-8eae-bfded3cf543d","metadata":{},"outputs":[{"data":{"text/plain":["Univariate Polynomial Ring in X over Finite Field of size 31"]},"execution_count":1,"metadata":{},"output_type":"execute_result"}],"source":["F31 = GF(31)\n","R31.<X> = F31[]\n","R31"]},{"attachments":{},"cell_type":"markdown","id":"a64c4622","metadata":{},"source":["Now we extend F31 to a 2-degree extension field GF($31^2$)."]},{"cell_type":"code","execution_count":2,"id":"f25afd7d-90a4-4c10-bcf9-ba392e61d69d","metadata":{},"outputs":[{"data":{"text/plain":["Finite Field in i of size 31^2"]},"execution_count":2,"metadata":{},"output_type":"execute_result"}],"source":["C31.<i> = F31.extension(X^2 + 1)\n","C31"]},{"attachments":{},"cell_type":"markdown","id":"b9bc3979-ea87-4c31-9a4e-5b5a2d55f252","metadata":{},"source":["We want to find a group element of order 32, so we take the generator of GF($31^2$) and get 30th power of it.\n","\n","**Fact**: For any prime field $\\mathbb{F}_p$, where $p=2^n - 1$, the 2-degree extension over $\\mathbb{F}_p$ is \n","$K=\\mathbb{F}_{p^2}$, then $K$ has a multiplicative subgroup with order $2^n$.\n","\n","Proof:\n","\n","$$\n","p^2 - 1 = (2^n - 1) * (2^n - 1) - 1 = 2^{2n} - 2^{n+1} = 2^n(2^n - 2)\n","$$\n","\n","In our case, F31.order() - 1 is $2^n - 2$."]},{"cell_type":"code","execution_count":3,"id":"f8f575a0-4a95-43e7-b0f7-ad2b7cd42103","metadata":{},"outputs":[{"data":{"text/plain":["11*i + 29"]},"execution_count":3,"metadata":{},"output_type":"execute_result"}],"source":["\n","g = C31.multiplicative_generator()^(F31.order() - 1)\n","g"]},{"attachments":{},"cell_type":"markdown","id":"ca647bd5","metadata":{},"source":["Define a function to calculate the order of an element in the circle group."]},{"cell_type":"code","execution_count":4,"id":"0089303e-1e4c-440d-a7c2-fc791b4d23cf","metadata":{},"outputs":[],"source":["def group_ord(elem):\n","    K = elem.parent()\n","\n","    for i in range(1, K.order()):\n","        # An order of an element is defined by how much power we need to raise the element to get the identity element\n","        if elem**i == K.one():\n","            return i\n","    return \"False\""]},{"attachments":{},"cell_type":"markdown","id":"609163d2","metadata":{},"source":["Check that g is of order 32."]},{"cell_type":"code","execution_count":5,"id":"6de6738e-e132-4825-a0ff-3185089d6a17","metadata":{},"outputs":[{"data":{"text/plain":["32"]},"execution_count":5,"metadata":{},"output_type":"execute_result"}],"source":["group_ord(g)"]},{"attachments":{},"cell_type":"markdown","id":"743f8247","metadata":{},"source":["Define the group multiplication, which is similar to complex number multiplication."]},{"cell_type":"code","execution_count":6,"id":"98f66494-5eca-42f9-9150-32d83cd8ce0f","metadata":{},"outputs":[],"source":["def group_mul(g1, g2):\n","    x1, y1 = g1\n","    x2, y2 = g2\n","    return (x1 * x2 - y1 * y2) + (x1 * y2 + x2 * y1) * i "]},{"attachments":{},"cell_type":"markdown","id":"0ecf4e47","metadata":{},"source":["Define the group inverse, which is also the J mapping in the paper."]},{"cell_type":"code","execution_count":7,"id":"10780cbd-3786-4642-8851-13642de9f72e","metadata":{},"outputs":[],"source":["def group_inv(g1):\n","    x1, y1 = g1\n","    return (x1 - y1 * i)"]},{"attachments":{},"cell_type":"markdown","id":"03439c18","metadata":{},"source":["Define show function to print group."]},{"cell_type":"code","execution_count":8,"id":"98046a5b-d75a-40fe-90e9-c4a9e3940512","metadata":{},"outputs":[],"source":["def show(GG): \n","    for j in range(0, len(GG)):\n","        t = GG[j]\n","        t0 = t[0]\n","        t1 = t[1]\n","        if t1 > F31(15) and t1 <= F31(30):\n","            t1_str = \"-\" + str(31-t1)\n","        else:\n","            t1_str = str(t1)\n","        print(\"i={}, log={}, ord={}, ({},{})\".format(j, log_gen(g, t), group_ord(t), t0, t1_str))"]},{"attachments":{},"cell_type":"markdown","id":"6249f760","metadata":{},"source":["Define log_gen function to calculate the power of an element in the circle group."]},{"cell_type":"code","execution_count":9,"id":"e3aa4193-e683-4834-a18c-c835b718a878","metadata":{},"outputs":[],"source":["def log_gen(gen, t):\n","    K = gen.parent()\n","    for i in range(1, K.order()):\n","        if gen^i == t:\n","            return i\n","    return \"Fail\""]},{"attachments":{},"cell_type":"markdown","id":"40185187-0f6f-40f5-9869-5611576e5ce5","metadata":{},"source":["## Jump over the circle\n","\n","Here we define the generator for `Circle Group` as chosen by Vitalik https://vitalik.eth.limo/general/2024/07/23/circlestarks.html"]},{"cell_type":"code","execution_count":10,"id":"e96e1eba-4ef3-4f3a-9003-43d090044a96","metadata":{},"outputs":[{"data":{"text/plain":["32"]},"execution_count":10,"metadata":{},"output_type":"execute_result"}],"source":["# optional\n","g = 10 + 5*i\n","group_ord(g)"]},{"cell_type":"markdown","id":"d7555c5e","metadata":{},"source":["Then we generate G5, which is a normal multiplicative group of order $2^5$.\n","\n","G5's elements must be on the circle $x^2 + y^2 = 1$, due to:\n","\n","$$\n","a^{p + 1} = 1, a = x + y * i\n","$$\n","\n","$$  \n","a * a ^ p = (x + y * i) * (x - y * i) = x^2 + y^2 = 1\n","$$\n","\n","Since $a^p$ is a Frobenius map, $a^p$ equals to the conjugate of $a$, which is $x - y * i$."]},{"cell_type":"code","execution_count":11,"id":"77e0b387-b325-4cf4-8c82-cc07368b21b3","metadata":{},"outputs":[{"data":{"text/plain":["[1,\n"," 5*i + 10,\n"," 7*i + 13,\n"," 11*i + 2,\n"," 27*i + 27,\n"," 2*i + 11,\n"," 13*i + 7,\n"," 10*i + 5,\n"," i,\n"," 10*i + 26,\n"," 13*i + 24,\n"," 2*i + 20,\n"," 27*i + 4,\n"," 11*i + 29,\n"," 7*i + 18,\n"," 5*i + 21,\n"," 30,\n"," 26*i + 21,\n"," 24*i + 18,\n"," 20*i + 29,\n"," 4*i + 4,\n"," 29*i + 20,\n"," 18*i + 24,\n"," 21*i + 26,\n"," 30*i,\n"," 21*i + 5,\n"," 18*i + 7,\n"," 29*i + 11,\n"," 4*i + 27,\n"," 20*i + 2,\n"," 24*i + 13,\n"," 26*i + 10]"]},"execution_count":11,"metadata":{},"output_type":"execute_result"}],"source":["G5 = [g^i for i in range(0, 2^5)]\n","# check that G5's elements are on the circle x^2 + y^2 = 1\n","for t in G5:\n","    x, y = t\n","    assert x^2 + y^2 == 1\n","G5"]},{"cell_type":"markdown","id":"ecca3bbe","metadata":{},"source":["Print information of G5's elements."]},{"cell_type":"code","execution_count":12,"id":"d138513c-a1a8-48fa-a93f-d50b131c2a2b","metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["i=0, log=32, ord=1, (1,0)\n","i=1, log=1, ord=32, (10,5)\n","i=2, log=2, ord=16, (13,7)\n","i=3, log=3, ord=32, (2,11)\n","i=4, log=4, ord=8, (27,27)\n","i=5, log=5, ord=32, (11,2)\n","i=6, log=6, ord=16, (7,13)\n","i=7, log=7, ord=32, (5,10)\n","i=8, log=8, ord=4, (0,1)\n","i=9, log=9, ord=32, (26,10)\n","i=10, log=10, ord=16, (24,13)\n","i=11, log=11, ord=32, (20,2)\n","i=12, log=12, ord=8, (4,27)\n","i=13, log=13, ord=32, (29,11)\n","i=14, log=14, ord=16, (18,7)\n","i=15, log=15, ord=32, (21,5)\n","i=16, log=16, ord=2, (30,0)\n","i=17, log=17, ord=32, (21,26)\n","i=18, log=18, ord=16, (18,24)\n","i=19, log=19, ord=32, (29,20)\n","i=20, log=20, ord=8, (4,4)\n","i=21, log=21, ord=32, (20,29)\n","i=22, log=22, ord=16, (24,18)\n","i=23, log=23, ord=32, (26,21)\n","i=24, log=24, ord=4, (0,30)\n","i=25, log=25, ord=32, (5,21)\n","i=26, log=26, ord=16, (7,18)\n","i=27, log=27, ord=32, (11,29)\n","i=28, log=28, ord=8, (27,4)\n","i=29, log=29, ord=32, (2,20)\n","i=30, log=30, ord=16, (13,24)\n","i=31, log=31, ord=32, (10,26)\n"]}],"source":["for j in range(0, 32):\n","    t = g^j\n","    print(\"i={}, log={}, ord={}, ({},{})\".format(j, log_gen(g, t), group_ord(t), t[0], t[1]))"]},{"cell_type":"markdown","id":"09f52870","metadata":{},"source":["Check elements' operation on G5."]},{"cell_type":"code","execution_count":13,"id":"331e689e-b199-4717-ac08-266f1a52e82c","metadata":{},"outputs":[],"source":["assert G5[1] * G5[1] == G5[2]\n","assert G5[3] * G5[4] == G5[7]\n","assert group_mul(G5[1], G5[2]) == G5[3]\n","assert group_inv(G5[1]) ==  G5[31]\n","assert G5[1] * group_inv(G5[1]) == 1"]},{"cell_type":"markdown","id":"0100cc02","metadata":{},"source":["Generate G4 from G5 by squaring each element."]},{"cell_type":"code","execution_count":14,"id":"8ed48611-82dc-47c6-8655-2fd1fdf4e074","metadata":{},"outputs":[{"data":{"text/plain":["[1,\n"," 7*i + 13,\n"," 27*i + 27,\n"," 13*i + 7,\n"," i,\n"," 13*i + 24,\n"," 27*i + 4,\n"," 7*i + 18,\n"," 30,\n"," 24*i + 18,\n"," 4*i + 4,\n"," 18*i + 24,\n"," 30*i,\n"," 18*i + 7,\n"," 4*i + 27,\n"," 24*i + 13]"]},"execution_count":14,"metadata":{},"output_type":"execute_result"}],"source":["G4 = [t^2 for t in G5[:16]]\n","G4"]},{"cell_type":"markdown","id":"97fb1c1c","metadata":{},"source":["Generate G3 from G4 by squaring each element."]},{"cell_type":"code","execution_count":15,"id":"673480c6-0531-46d7-a57f-950f5b46edd8","metadata":{},"outputs":[{"data":{"text/plain":["[1, 27*i + 27, i, 27*i + 4, 30, 4*i + 4, 30*i, 4*i + 27]"]},"execution_count":15,"metadata":{},"output_type":"execute_result"}],"source":["G3 = [t^2 for t in G4[:8]]\n","G3"]},{"attachments":{},"cell_type":"markdown","id":"ce8e103e-2e1d-495e-b129-7e441455a3ed","metadata":{},"source":["## Standard Position Coset\n","\n","Let $Q$ be the generator of $G_{n+1}$, $D$ is a `standard position coset` of size $2^n$\n","\n","$$\n","D = Q\\cdot G_n\n","$$\n","\n","where $ord(Q)=2^{n+1}$.\n","\n","And $D$ is also a `twin-coset` of size $2^n$ defined with $G_{n-1}$:\n","\n","$$\n","D = Q\\cdot G_{n-1} \\uplus Q^{-1}\\cdot G_{n-1}\n","$$\n","\n","The rotation $Q^2$ is exactly the generator of $G_n$.\n","\n","Here we construct a series of cosets, including standard position cosets and twin cosets, where we can randomly find a coset for circle fft."]},{"cell_type":"code","execution_count":16,"id":"396deee6-21a2-4f38-8969-f620664bdd6b","metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["i=0, log=1, ord=32, (10,5)\n","i=1, log=3, ord=32, (2,11)\n","i=2, log=5, ord=32, (11,2)\n","i=3, log=7, ord=32, (5,10)\n","i=4, log=9, ord=32, (26,10)\n","i=5, log=11, ord=32, (20,2)\n","i=6, log=13, ord=32, (29,11)\n","i=7, log=15, ord=32, (21,5)\n","i=8, log=17, ord=32, (21,-5)\n","i=9, log=19, ord=32, (29,-11)\n","i=10, log=21, ord=32, (20,-2)\n","i=11, log=23, ord=32, (26,-10)\n","i=12, log=25, ord=32, (5,-10)\n","i=13, log=27, ord=32, (11,-2)\n","i=14, log=29, ord=32, (2,-11)\n","i=15, log=31, ord=32, (10,-5)\n"]}],"source":["# Define a coset or G4 from a rotation Q, ord(Q) = 32\n","#\n","# Let Q = g (generator of G, whose order is 32)\n","\n","D_standard_coset_4 = [g * t for t in G4]; show(D_standard_coset_4)"]},{"cell_type":"code","execution_count":17,"id":"f24c075e-89eb-4db5-a533-c331b8f91a4c","metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["i=0, log=1, ord=32, (10,5)\n","i=1, log=3, ord=32, (2,11)\n","i=2, log=5, ord=32, (11,2)\n","i=3, log=7, ord=32, (5,10)\n","i=4, log=9, ord=32, (26,10)\n","i=5, log=11, ord=32, (20,2)\n","i=6, log=13, ord=32, (29,11)\n","i=7, log=15, ord=32, (21,5)\n","i=8, log=17, ord=32, (21,-5)\n","i=9, log=19, ord=32, (29,-11)\n","i=10, log=21, ord=32, (20,-2)\n","i=11, log=23, ord=32, (26,-10)\n","i=12, log=25, ord=32, (5,-10)\n","i=13, log=27, ord=32, (11,-2)\n","i=14, log=29, ord=32, (2,-11)\n","i=15, log=31, ord=32, (10,-5)\n"]}],"source":["# D_twin_coset_4 is equal to D_standard_coset_4 in the set-theoretic sense, with different order\n","\n","D_twin_coset_4 = [g * t for t in G3] + [group_inv(g) * t for t in G3]; show(D_standard_coset_4)\n","\n","# `log` specifies the order of an element in the circle"]},{"cell_type":"code","execution_count":18,"id":"12dc1f2a-d03c-46ea-bedc-217e53c9cb9b","metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["i=0, log=3, ord=32, (2,11)\n","i=1, log=7, ord=32, (5,10)\n","i=2, log=11, ord=32, (20,2)\n","i=3, log=15, ord=32, (21,5)\n","i=4, log=19, ord=32, (29,-11)\n","i=5, log=23, ord=32, (26,-10)\n","i=6, log=27, ord=32, (11,-2)\n","i=7, log=31, ord=32, (10,-5)\n","i=8, log=29, ord=32, (2,-11)\n","i=9, log=1, ord=32, (10,5)\n","i=10, log=5, ord=32, (11,2)\n","i=11, log=9, ord=32, (26,10)\n","i=12, log=13, ord=32, (29,11)\n","i=13, log=17, ord=32, (21,-5)\n","i=14, log=21, ord=32, (20,-2)\n","i=15, log=25, ord=32, (5,-10)\n"]}],"source":["# This is also a twin-coset of size 16\n","\n","D_twin_coset_4_1 = [G5[3] * t for t in G3] + [group_inv(G5[3]) * t for t in G3]; show(D_twin_coset_4_1)"]},{"cell_type":"code","execution_count":19,"id":"b08f58a7-52ce-48c0-91a1-6e692344c627","metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["i=0, log=5, ord=32, (11,2)\n","i=1, log=9, ord=32, (26,10)\n","i=2, log=13, ord=32, (29,11)\n","i=3, log=17, ord=32, (21,-5)\n","i=4, log=21, ord=32, (20,-2)\n","i=5, log=25, ord=32, (5,-10)\n","i=6, log=29, ord=32, (2,-11)\n","i=7, log=1, ord=32, (10,5)\n","i=8, log=27, ord=32, (11,-2)\n","i=9, log=31, ord=32, (10,-5)\n","i=10, log=3, ord=32, (2,11)\n","i=11, log=7, ord=32, (5,10)\n","i=12, log=11, ord=32, (20,2)\n","i=13, log=15, ord=32, (21,5)\n","i=14, log=19, ord=32, (29,-11)\n","i=15, log=23, ord=32, (26,-10)\n"]}],"source":["# Again\n","\n","D_twin_coset_4_2 = [G5[5] * t for t in G3] + [group_inv(G5[5]) * t for t in G3]; show(D_twin_coset_4_2)"]},{"cell_type":"code","execution_count":20,"id":"17b78d53-4cb2-446b-ac00-1006df2568d9","metadata":{},"outputs":[{"data":{"text/plain":["True"]},"execution_count":20,"metadata":{},"output_type":"execute_result"}],"source":["# Fact:\n","#   g^2 is the generator of G4\n","#   g * G4 is the standard position coset\n","\n","G4 == [(g^2)^i for i in range(0, 16)]"]},{"cell_type":"code","execution_count":21,"id":"3976036e-e166-4084-adb6-28594cc4d1f5","metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["i=0, log=2, ord=16, (13,7)\n","i=1, log=6, ord=16, (7,13)\n","i=2, log=10, ord=16, (24,13)\n","i=3, log=14, ord=16, (18,7)\n","i=4, log=18, ord=16, (18,-7)\n","i=5, log=22, ord=16, (24,-13)\n","i=6, log=26, ord=16, (7,-13)\n","i=7, log=30, ord=16, (13,-7)\n"]}],"source":["# Define a coset or G3 from a rotation Q, ord(Q) = 16\n","#\n","# Let Q = g^2 (g is the generator of G5, whose order is 32)\n","\n","D_standard_coset_3 = [g^2 * t for t in G3]; show(D_standard_coset_3)"]},{"attachments":{},"cell_type":"markdown","id":"261d3556-cc74-4961-b271-88b03fe3045c","metadata":{},"source":["## Squaring map\n","\n","**Lemma**: For any twin-coset $D$ of size $N=2^n$, the image of $\\pi(D)$ must be a twin-coset of size $2^{n-1}$.\n","\n","$$\n","D = Q\\cdot G_{n-1} \\uplus Q^{-1}\\cdot G_{n-1}\n","$$\n","\n","$$\n","\\pi(D) = \\pi(Q)\\cdot G_{n-2} \\uplus \\pi(Q)^{-1}\\cdot G_{n-2}\n","$$\n","\n","**Fact**: If $D$ is a standard position coset, so is $\\pi(D)$."]},{"cell_type":"code","execution_count":22,"id":"afd781ca-8a1c-4a2a-af97-cfc8a2b63179","metadata":{},"outputs":[],"source":["def sq(D):\n","    rs = []\n","    for t in D:\n","        if t^2 not in rs:\n","            # x' == 2 * x^2 - 1\n","            rs += [t^2]\n","    return rs"]},{"cell_type":"code","execution_count":23,"id":"99e26070-ecd6-4795-a54d-1721e9319209","metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["i=0, log=2, ord=16, (13,7)\n","i=1, log=6, ord=16, (7,13)\n","i=2, log=10, ord=16, (24,13)\n","i=3, log=14, ord=16, (18,7)\n","i=4, log=18, ord=16, (18,-7)\n","i=5, log=22, ord=16, (24,-13)\n","i=6, log=26, ord=16, (7,-13)\n","i=7, log=30, ord=16, (13,-7)\n"]}],"source":["# This is exactly the same coset we generated from G4 before\n","\n","D_standard_coset_3 = sq(D_standard_coset_4); show(D_standard_coset_3)"]},{"attachments":{},"cell_type":"markdown","id":"89ebe6b9","metadata":{},"source":["## Circle FFT"]},{"cell_type":"markdown","id":"0cd1bb39","metadata":{},"source":["To apply circle fft, let's generate a random polynomial's evaluation based on D_twin_coset_4_1 (or any other twin coset you prefer)."]},{"cell_type":"code","execution_count":24,"id":"1414592f","metadata":{},"outputs":[{"data":{"text/plain":["{11*i + 2: 23,\n"," 10*i + 5: 8,\n"," 2*i + 20: 24,\n"," 5*i + 21: 2,\n"," 20*i + 29: 0,\n"," 21*i + 26: 25,\n"," 29*i + 11: 10,\n"," 26*i + 10: 0,\n"," 20*i + 2: 22,\n"," 5*i + 10: 4,\n"," 2*i + 11: 4,\n"," 10*i + 26: 16,\n"," 11*i + 29: 27,\n"," 26*i + 21: 0,\n"," 29*i + 20: 21,\n"," 21*i + 5: 2}"]},"execution_count":24,"metadata":{},"output_type":"execute_result"}],"source":["f = {}\n","from random import randint\n","\n","for t in D_twin_coset_4_1:\n","    f[t] = C31(randint(0, 100))\n","\n","f"]},{"cell_type":"markdown","id":"a93fac2b","metadata":{},"source":["Define the `pi` mapping applying to point's x coordinate.\n"]},{"cell_type":"code","execution_count":25,"id":"0f1f54b5","metadata":{},"outputs":[],"source":["def pi(t):\n","    # x^2 - y^2 == 2 * x^2 - 1 (x^2 + y^2 = 1)\n","    return C31(2 * t^2 - 1)"]},{"cell_type":"markdown","id":"f2d3d49e","metadata":{},"source":["In the first step of ifft, we try to eliminate y coordinate to simplifier our computations in later steps.\n","\n","We extract y from the polynomial by dividing polynomial in to 2 parts, f0 and f1, which correspond to the equation f[t] == f0[x] + y * f1[x].\n","\n","t is a group element, x and y are the x and y coordinate of t respectively.\n","\n","So that f0 and f1 are only functions of x coordinate."]},{"cell_type":"code","execution_count":26,"id":"2a87e2ae","metadata":{},"outputs":[],"source":["def ifft_first_step(f):\n","    f0 = {}\n","    f1 = {}\n","    for t in f:\n","        x, y = t\n","\n","        f0[x] = (f[t] + f[group_inv(t)]) / C31(2)\n","        f1[x] = (f[t] - f[group_inv(t)]) / (C31(2) * y)\n","\n","        # Check that f is divided into 2 parts correctly\n","        assert f[t] == f0[x] + y * f1[x]\n","\n","        # print(\"{}: {} = {} + {} * {}\".format(t, f[t], f0[x], f1[x], y))\n","\n","    return f0, f1"]},{"cell_type":"markdown","id":"8ce68c05","metadata":{},"source":["In normal steps, we deal with the polynomial exactly the same as the classical polynomial ifft, deeming the polynomial as a univariate polynomial, except that x is applied by the `pi` mapping.\n","\n","Basically, we divide the polynomial into 2 parts, $f_0$ and $f_1$, which correspond to the equation:\n","\n","$$\n","f[x] = f_0[x] + x * f_1[x]\n","$$\n","\n","Thus\n","\n","$$\n","f_0[x] = \\frac{f[x] + f[-x]}{2}\n","$$\n","$$\n","f_1[x] = \\frac{f[x] - f[-x]}{2x}\n","$$\n"]},{"cell_type":"code","execution_count":27,"id":"936da672","metadata":{},"outputs":[],"source":["def ifft_normal_step(f):\n","\n","    if len(f) == 1:\n","        res = []\n","        for x in f:\n","            res.append(f[x])\n","        return res\n","\n","    f0 = {}\n","    f1 = {}\n","\n","    for x in f:\n","        assert x != 0, \"f should be on coset\"\n","        f0[pi(x)] = (f[x] + f[-x]) / C31(2)\n","        f1[pi(x)] = (f[x] - f[-x]) / (C31(2) * x)\n","\n","        # Check that f is divided into 2 parts correctly\n","        assert f[x] == f0[pi(x)] + x * f1[pi(x)]\n","\n","    return ifft_normal_step(f0) + ifft_normal_step(f1)"]},{"cell_type":"markdown","id":"5af274b7","metadata":{},"source":["Now we have ifft_first_step and ifft_normal_step, we can define the ifft function."]},{"cell_type":"code","execution_count":28,"id":"b8b62ba5","metadata":{},"outputs":[],"source":["def ifft(f):\n","    f0, f1 = ifft_first_step(f)\n","    f0 = ifft_normal_step(f0)\n","    f1 = ifft_normal_step(f1)\n","\n","    return f0 + f1"]},{"cell_type":"markdown","id":"0fe0200e","metadata":{},"source":["Apply ifft to our random polynomial and see what we get."]},{"cell_type":"code","execution_count":29,"id":"f7349218","metadata":{},"outputs":[{"data":{"text/plain":["[4, 6, 22, 3, 6, 23, 18, 30, 27, 26, 13, 24, 17, 16, 29, 2]"]},"execution_count":29,"metadata":{},"output_type":"execute_result"}],"source":["coeffs = ifft(f)\n","coeffs"]},{"cell_type":"markdown","id":"35516858","metadata":{},"source":["Now we try to define the fft function.\n","\n","Before that, we define `pie_group` function to generate next domain of a domain, which is equivalent to `sq` function in the x dimension (you can check that by focusing on the x coordinate in `sq` operation).\n"]},{"cell_type":"code","execution_count":30,"id":"2ea6a0e0","metadata":{},"outputs":[],"source":["def pie_group(D):\n","    D_new = []\n","    for x in D:\n","        x_new = pi(x)\n","        if x_new not in D_new:\n","            D_new.append(x_new)\n","\n","    # Check that the new domain is exactly half size of the old domain\n","    assert len(D_new) * 2 == len(D), \"len(D_new) * 2 != len(D), {} * 2 != {}, D_new={}, D={}\".format(len(D_new), len(D), D_new, D)\n","    \n","    return D_new"]},{"cell_type":"markdown","id":"55e8d597","metadata":{},"source":["Define the first step of fft, which contains the same logic as ifft_first_step.\n","\n","The basis of coefficients is like:\n","\n","$$\n","1, x, x^2, x^3, \\cdots, x^{\\frac{N}{2}}, y, yx, yx^2, yx^3, \\cdots, yx^{\\frac{N}{2}} \\tag{1}\n","$$\n","\n","So to divide the polynomial into 2 parts, we just need to divide the coefficients into 2 parts, it's basis are also divided into 2 parts as follows:\n","\n","$$\n","1, x, x^2, x^3, \\cdots, x^{\\frac{N}{2}} \\tag{2}\n","$$\n","\n","$$\n","y, yx, yx^2, yx^3, \\cdots, yx^{\\frac{N}{2}} \\tag{3}\n","$$\n","\n","To eliminate y in normal steps' computation, we just consider basis of the second part of coefficients (which is actually basis (3)) as same as basis (2), and supply the y at the end of fft.\n","\n","First, we define the first step of fft."]},{"cell_type":"code","execution_count":31,"id":"d509cd22","metadata":{},"outputs":[],"source":["def fft_first_step(f, D):\n","    # Check that the polynomial and the domain have the same length\n","    assert len(f) == len(D), \"len(f) != len(D), {} != {}, f={}, D={}\".format(len(f), len(D), f, D)\n","\n","    # divide the polynomial into 2 parts\n","    len_f = len(f)\n","    f0 = f[:len_f//2]\n","    f1 = f[len_f//2:]\n","\n","    # halve the domain by simply removing the y coordinate\n","    D_new = []\n","    for t in D:\n","        x, _ = t\n","        if x not in D_new:\n","            D_new.append(x)\n","\n","    # Check that the new domain is exactly half size of the old domain\n","    assert len(D_new) * 2 == len(D), \"len(D_new) * 2 != len(D), {} * 2 != {}, D_new={}, D={}\".format(len(D_new), len(D), D_new, D)\n","\n","    return f0, f1, D_new"]},{"cell_type":"markdown","id":"2c8684be","metadata":{},"source":["In normal steps, we deal with the polynomial exactly the same as the classical polynomial fft, deeming the polynomial as a univariate polynomial, except that $x$ is applied by the `pi` mapping.\n","\n","Divide the polynomial into 2 parts, $f_0$ and $f_1$, which correspond to the equation:\n","\n","$$\n","f[x] = f_0[x] + x * f_1[x]\n","$$\n","\n","Merge $f_0$ and $f_1$ using $x$ from the domain, which corresponds to the equation above too."]},{"cell_type":"code","execution_count":32,"id":"4413ffd8","metadata":{},"outputs":[],"source":["def fft_normal_step(f, D):\n","\n","    if len(f) == 1:\n","        return {D[0]: f[0]}\n","    \n","    next_domain = pie_group(D)\n","\n","    # Check that the new domain is exactly half size of the old domain\n","    assert len(next_domain) * 2 == len(D), \"len(next_domain) * 2 != len(D), {} * 2 != {}, next_domain={}, D={}\".format(len(next_domain), len(D), next_domain, D)\n","    # Check that the polynomial and the domain have the same length\n","    assert len(f) == len(D), \"len(f) != len(D), {} != {}, f={}, D={}\".format(len(f), len(D), f, D)\n","\n","    f0 = fft_normal_step(f[:len(f)//2], next_domain)\n","    f1 = fft_normal_step(f[len(f)//2:], next_domain)\n","\n","    f_new = {}\n","    for x in D:\n","        f_new[x] = f0[pi(x)] + f1[pi(x)] * x\n","\n","    # Check that f is divided into 2 parts correctly\n","    for x in D:\n","        if x != 0:\n","            assert f0[pi(x)] == (f_new[x] + f_new[-x]) / C31(2), \"f0[pi(x)] = {}\".format(f0[pi(x)])\n","            assert f1[pi(x)] == (f_new[x] - f_new[-x]) / (C31(2) * x), \"f1[pi(x)] = {}\".format(f1[pi(x)])\n","        else:\n","            assert f0[pi(x)] == f_new[x], \"f0[pi(x)] = {}\".format(f0[pi(x)])\n","\n","    # Check that the polynomial and the domain have the same length\n","    assert len(f) == len(f_new), \"len(f) != len(f_new), {} != {}, f={}, f_new={}, D={}\".format(len(f), len(f_new), f, f_new, D)\n","\n","    # Check that ifft and fft are correct inverse operations\n","    assert ifft_normal_step(f_new) == f, \"ifft(f_new) != f, {} != {}\".format(ifft_normal_step(f_new), f)\n","        \n","    return f_new"]},{"cell_type":"markdown","id":"72f35757","metadata":{},"source":["We have fft_first_step and fft_normal_step, we can define the fft function."]},{"cell_type":"code","execution_count":33,"id":"2b5ce548","metadata":{},"outputs":[],"source":["def fft(f, D):\n","\n","    # Check that the polynomial and the domain have the same length\n","    assert len(f) == len(D), \"len(f) != len(D), {} != {}, f={}, D={}\".format(len(f), len(D), f, D)\n","\n","    D_copy = D[:]\n","    f0, f1, D = fft_first_step(f, D)\n","\n","    # Check that the polynomial and the domain have the same length\n","    assert len(f0) == len(D), \"len(f0) != len(D), {} != {}, f0={}, D={}\".format(len(f0), len(D), f0, D)\n","    assert len(f1) == len(D), \"len(f1) != len(D), {} != {}, f1={}, D={}\".format(len(f1), len(D), f1, D)\n","\n","    f0 = fft_normal_step(f0, D)\n","    f1 = fft_normal_step(f1, D)\n","\n","    f = {}\n","    # supply y to the polynomial\n","    for t in D_copy:\n","        x, y = t\n","        f[t] = f0[x] + f1[x] * y\n","\n","    return f\n"]},{"cell_type":"markdown","id":"1a6d1dad","metadata":{},"source":["Check that ifft and fft are correct inverse operations."]},{"cell_type":"code","execution_count":34,"id":"89d2a2c7","metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["f[11*i + 2]=23, f_prime[11*i + 2]=23\n","f[10*i + 5]=8, f_prime[10*i + 5]=8\n","f[2*i + 20]=24, f_prime[2*i + 20]=24\n","f[5*i + 21]=2, f_prime[5*i + 21]=2\n","f[20*i + 29]=0, f_prime[20*i + 29]=0\n","f[21*i + 26]=25, f_prime[21*i + 26]=25\n","f[29*i + 11]=10, f_prime[29*i + 11]=10\n","f[26*i + 10]=0, f_prime[26*i + 10]=0\n","f[20*i + 2]=22, f_prime[20*i + 2]=22\n","f[5*i + 10]=4, f_prime[5*i + 10]=4\n","f[2*i + 11]=4, f_prime[2*i + 11]=4\n","f[10*i + 26]=16, f_prime[10*i + 26]=16\n","f[11*i + 29]=27, f_prime[11*i + 29]=27\n","f[26*i + 21]=0, f_prime[26*i + 21]=0\n","f[29*i + 20]=21, f_prime[29*i + 20]=21\n","f[21*i + 5]=2, f_prime[21*i + 5]=2\n"]}],"source":["f_prime = fft(coeffs, D_twin_coset_4_1)\n","\n","for t, s in zip(f, f_prime):\n","    assert t == s\n","\n","for t in f:\n","    print(\"f[{}]={}, f_prime[{}]={}\".format(t, f[t], t, f_prime[t]))\n","    assert f[t] == f_prime[t], \"f[{}] != f_prime[{}], {} != {}\".format(t, t, f[t], f_prime[t])"]}],"metadata":{"kernelspec":{"display_name":"SageMath 9.8","language":"sage","name":"SageMath-9.8"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.11.1"}},"nbformat":4,"nbformat_minor":5}