From 4a2dbeca6d382efbd2cc4c5ffb928465b158fa61 Mon Sep 17 00:00:00 2001 From: Neil Date: Fri, 25 Oct 2024 18:19:01 +0700 Subject: [PATCH] feat(circle): add circle FRI Add circle FRI to circle.ipynb --- src/circle.ipynb | 1 + src/circle_fft.ipynb | 1 - 2 files changed, 1 insertion(+), 1 deletion(-) create mode 100644 src/circle.ipynb delete mode 100644 src/circle_fft.ipynb diff --git a/src/circle.ipynb b/src/circle.ipynb new file mode 100644 index 0000000..9420ad1 --- /dev/null +++ b/src/circle.ipynb @@ -0,0 +1 @@ +{"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":3003,"id":"2c77b9d4-0cab-4100-8eae-bfded3cf543d","metadata":{},"outputs":[{"data":{"text/plain":["Univariate Polynomial Ring in X over Finite Field of size 31"]},"execution_count":3003,"metadata":{},"output_type":"execute_result"}],"source":["F31 = GF(31)\n","R31. = 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":3004,"id":"f25afd7d-90a4-4c10-bcf9-ba392e61d69d","metadata":{},"outputs":[{"data":{"text/plain":["Finite Field in i of size 31^2"]},"execution_count":3004,"metadata":{},"output_type":"execute_result"}],"source":["C31. = 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":3005,"id":"f8f575a0-4a95-43e7-b0f7-ad2b7cd42103","metadata":{},"outputs":[{"data":{"text/plain":["11*i + 29"]},"execution_count":3005,"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":3006,"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":3007,"id":"6de6738e-e132-4825-a0ff-3185089d6a17","metadata":{},"outputs":[{"data":{"text/plain":["32"]},"execution_count":3007,"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":3008,"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":3009,"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":3010,"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":3011,"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":3012,"id":"e96e1eba-4ef3-4f3a-9003-43d090044a96","metadata":{},"outputs":[{"data":{"text/plain":["32"]},"execution_count":3012,"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":3013,"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":3013,"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":3014,"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":3015,"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":3016,"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":3016,"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":3017,"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":3017,"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":3018,"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":3019,"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":3020,"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":3021,"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":3022,"id":"17b78d53-4cb2-446b-ac00-1006df2568d9","metadata":{},"outputs":[{"data":{"text/plain":["True"]},"execution_count":3022,"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":3023,"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":3024,"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":3025,"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":3026,"id":"1414592f","metadata":{},"outputs":[{"data":{"text/plain":["{11*i + 2: 23,\n"," 10*i + 5: 5,\n"," 2*i + 20: 8,\n"," 5*i + 21: 8,\n"," 20*i + 29: 8,\n"," 21*i + 26: 19,\n"," 29*i + 11: 6,\n"," 26*i + 10: 27,\n"," 20*i + 2: 1,\n"," 5*i + 10: 0,\n"," 2*i + 11: 11,\n"," 10*i + 26: 10,\n"," 11*i + 29: 30,\n"," 26*i + 21: 8,\n"," 29*i + 20: 10,\n"," 21*i + 5: 9}"]},"execution_count":3026,"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":3027,"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":3028,"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":3029,"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":3030,"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":3031,"id":"f7349218","metadata":{},"outputs":[{"data":{"text/plain":["[25, 21, 21, 13, 4, 6, 18, 8, 10, 3, 21, 25, 9, 3, 16, 3]"]},"execution_count":3031,"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":3032,"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^2, x, x^3, y, yx^2, yx, yx^3 \\quad (N = 8) \\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^2, x, x^3 \\tag{2}\n","$$\n","\n","$$\n","y, yx^2, yx, yx^3 \\tag{3}\n","$$\n","\n","Notice that basis (2) is bit reversed order compared to normal basis, and basis (3) is like basis (2) multiplied by y.\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":3033,"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":3034,"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":3035,"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":3036,"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]=5, f_prime[10*i + 5]=5\n","f[2*i + 20]=8, f_prime[2*i + 20]=8\n","f[5*i + 21]=8, f_prime[5*i + 21]=8\n","f[20*i + 29]=8, f_prime[20*i + 29]=8\n","f[21*i + 26]=19, f_prime[21*i + 26]=19\n","f[29*i + 11]=6, f_prime[29*i + 11]=6\n","f[26*i + 10]=27, f_prime[26*i + 10]=27\n","f[20*i + 2]=1, f_prime[20*i + 2]=1\n","f[5*i + 10]=0, f_prime[5*i + 10]=0\n","f[2*i + 11]=11, f_prime[2*i + 11]=11\n","f[10*i + 26]=10, f_prime[10*i + 26]=10\n","f[11*i + 29]=30, f_prime[11*i + 29]=30\n","f[26*i + 21]=8, f_prime[26*i + 21]=8\n","f[29*i + 20]=10, f_prime[29*i + 20]=10\n","f[21*i + 5]=9, f_prime[21*i + 5]=9\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])"]},{"cell_type":"markdown","id":"f56384cd","metadata":{},"source":["## Circle FRI"]},{"cell_type":"code","execution_count":3037,"id":"d73add2d","metadata":{},"outputs":[],"source":["# Inputs:\n","# f is the polynomial to be folded\n","# D is the domain of f\n","# r is the random number for random linear combination\n","# Outputs:\n","# The first return value is the folded polynomial\n","# The second return value is the new domain\n","def fold(f, D, r, debug=False):\n"," assert len(f) == len(D), \"len(f) != len(D), {} != {}, f={}, D={}\".format(len(f), len(D), f, D)\n","\n"," # divide\n"," N = len(f)\n"," # left is the first half of f, of x from 1 to g^(N/2)\n"," left = f[:N//2]\n"," # right is the second half of f, of x from g^(N-1) to g^(N/2), which corresponds to minus x in left\n"," right = f[:N//2-1:-1]\n"," assert len(left) == len(right), \"len(left) != len(right), {} != {}, left={}, right={}\".format(len(left), len(right), left, right)\n","\n"," for i, x in enumerate(D[:N//2]):\n"," # f == f0 + x * f1\n"," f0 = (left[i] + right[i]) / C31(2)\n"," f1 = (left[i] - right[i]) / (C31(2) * x)\n"," # f[:N//2] stores the folded polynomial\n"," if debug: print(f\"f[{i}] = {f[i]} = ({left[i]} + {right[i]})/2 + {r} * ({left[i]} - {right[i]})/(2 * {x})\")\n"," f[i] = f0 + r * f1\n"," # if debug: print(f\"{f[i]} = ({left[i]} + {right[i]})/2 + {r} * ({left[i]} - {right[i]})/(2 * {x})\")\n"," # reuse f[N//2:] to store new domain\n"," f[N//2 + i] = 2 * x^2 - 1\n","\n"," # return the folded polynomial and the new domain\n"," return f[:N//2], f[N//2:]"]},{"cell_type":"code","execution_count":3038,"id":"1a3ade08","metadata":{},"outputs":[],"source":["from utils import log_2\n","from merlin.merlin_transcript import MerlinTranscript\n","from merkle import MerkleTree\n","\n","def fri(f, D, degree, query_num, transcript, debug=False):\n","\n"," assert len(f) == len(D), \"len(f) != len(D), {} != {}, f={}, D={}\".format(len(f), len(D), f, D)\n"," assert isinstance(transcript, MerlinTranscript)\n","\n"," if debug: print(f\"f={f}, D={D}\")\n","\n"," # first tree\n"," first_tree = MerkleTree(f)\n"," first_oracle = f\n","\n"," # first random number\n"," transcript.append_message(b'first_tree', bytes(first_tree.root, 'ascii'))\n"," r = int.from_bytes(transcript.challenge_bytes(b'r', int(4)), 'big')\n","\n"," # first step (J mapping)\n"," # for f in natural order, we just divide f into 2 parts from the middle\n"," N = len(f)\n"," assert N % 2 == 0, \"N must be even, N={}\".format(N)\n","\n"," left = f[:N//2]\n"," right = f[:N//2-1:-1]\n"," assert len(left) == len(right), \"len(left) != len(right), {} != {}, left={}, right={}\".format(len(left), len(right), left, right)\n"," f = [None for _ in range(N//2)]\n"," for i, (_, y) in enumerate(D[:N//2]):\n"," f0 = (left[i] + right[i]) / C31(2)\n"," f1 = (left[i] - right[i]) / (C31(2) * y)\n"," f[i] = f0 + r * f1\n"," if debug: print(f\"f[{i}] = {f[i]} = ({left[i]} + {right[i]})/2 + {r} * ({left[i]} - {right[i]})/(2 * {y})\")\n","\n"," print(f\"f={f}\")\n","\n"," D = [x for x, _ in D[:N//2]]\n","\n"," # fold\n"," trees = [MerkleTree(f)]\n"," oracles = [f[:]]\n"," # log_2(degree) - 1 because we have already done the first step\n"," for i in range(log_2(degree) - 1):\n"," # random number\n"," transcript.append_message(b'tree', bytes(trees[-1].root, 'ascii'))\n","\n"," r = int.from_bytes(transcript.challenge_bytes(b'r', int(4)), 'big')\n"," f, D = fold(f, D, r, debug)\n","\n"," if debug: print(f\"f={f}, D={D}\")\n","\n"," # merkle tree\n"," trees.append(MerkleTree(f))\n"," oracles.append(f[:])\n","\n"," # query paths\n"," max_size = N\n"," get_q = lambda transcript: int.from_bytes(transcript.challenge_bytes(b'query', int(4)), 'big') % max_size\n"," queries = [get_q(transcript) for _ in range(query_num)]\n","\n"," query_paths = []\n"," for q in queries:\n"," max_size = N\n"," cur_path = []\n"," indices = []\n"," x0 = q\n"," x1 = max_size - q - 1\n"," if x1 < x0:\n"," x0, x1 = x1, x0\n","\n"," if debug: print(f\"max_size={max_size}\")\n"," if debug: print(f\"q={q}\")\n"," if debug: print(f\"x0={x0}, x1={x1}\")\n"," \n"," cur_path.append((first_oracle[x0], first_oracle[x1]))\n"," indices.append(x0)\n"," q = x0\n"," max_size >>= 1\n","\n"," for oracle in oracles:\n"," x0 = q\n"," x1 = max_size - q - 1\n"," if x1 < x0:\n"," x0, x1 = x1, x0\n","\n"," if debug: print(f\"q={q}\")\n"," if debug: print(f\"x0={x0}, x1={x1}\")\n"," \n"," cur_path.append((oracle[x0], oracle[x1]))\n"," indices.append(x0)\n"," q = x0\n"," max_size >>= 1\n","\n"," query_paths.append((cur_path, indices))\n","\n"," # merkle paths\n"," merkle_paths = []\n"," for _, indices in query_paths:\n"," cur_query_paths = []\n"," for i, idx in enumerate(indices):\n"," if i == 0:\n"," if debug: print(f\"i = {i}, idx = {idx}, first_oracle[idx] = {first_oracle[idx]}\")\n"," cur_query_paths.append(first_tree.get_authentication_path(idx))\n"," else:\n"," cur_tree = trees[i - 1]\n"," if debug: print(f\"i = {i}, idx = {idx}, oracles[i-1][idx] = {oracles[i-1][idx]}\")\n"," cur_query_paths.append(cur_tree.get_authentication_path(idx))\n"," merkle_paths.append(cur_query_paths)\n","\n"," return {\n"," 'query_paths': query_paths,\n"," 'merkle_paths': merkle_paths,\n"," 'first_tree': first_tree.root,\n"," 'intermediate_trees': [tree.root for tree in trees],\n"," 'degree_bound': degree,\n"," 'final_value': f[0],\n"," }"]},{"cell_type":"code","execution_count":3039,"id":"98890e94","metadata":{},"outputs":[],"source":["from merkle import verify_decommitment\n","\n","def defri(fri_proof, degree_bound, T, query_num, transcript, debug=False):\n","\n"," assert degree_bound >= fri_proof['degree_bound']\n"," degree_bound = fri_proof['degree_bound']\n","\n"," assert isinstance(transcript, MerlinTranscript)\n","\n"," first_tree = fri_proof['first_tree']\n"," intermediate_trees = fri_proof['intermediate_trees']\n","\n"," transcript.append_message(b'first_tree', bytes(first_tree, 'ascii'))\n"," r = int.from_bytes(transcript.challenge_bytes(b'r', int(4)), 'big')\n","\n"," fold_challenges = [r]\n"," for i in range(log_2(int(degree_bound)) - 1):\n"," transcript.append_message(b'tree', bytes(intermediate_trees[i], 'ascii'))\n"," r = int.from_bytes(transcript.challenge_bytes(b'r', int(4)), 'big')\n"," fold_challenges.append(r)\n","\n"," get_q = lambda transcript: int.from_bytes(transcript.challenge_bytes(b'query', int(4)), 'big') % degree_bound\n"," queries = [get_q(transcript) for _ in range(query_num)]\n","\n"," for q, (cur_path, _), mps in zip(queries, fri_proof['query_paths'], fri_proof['merkle_paths']):\n"," num_vars = degree_bound\n"," for i, mp in enumerate(mps):\n","\n"," assert num_vars > 0, \"num_vars must be positive, num_vars={}\".format(num_vars)\n","\n"," x0 = q\n"," x1 = num_vars - q - 1\n"," if x1 < x0:\n"," x0, x1 = x1, x0\n"," \n"," if debug: print(f\"num_vars={num_vars}\")\n"," if debug: print(f\"q={q}\")\n"," if debug: print(f\"x0={x0}, x1={x1}\")\n","\n"," code_left, code_right = F31(cur_path[i][0]), F31(cur_path[i][1])\n","\n"," table = T[i]\n"," if debug: print(f\"i = {i}, table={table}\")\n"," if i != len(mps) - 1:\n"," alpha = fold_challenges[i]\n"," f_code_folded = cur_path[i + 1][0 if x0 < num_vars / 4 else 1]\n"," assert f_code_folded == (code_left + code_right)/2 + alpha * (code_left - code_right)/(2*table[x0]), f\"{f_code_folded} != ({code_left} + {code_right})/2 + {alpha} * ({code_left} - {code_right})/(2*{table[x0]})\"\n"," else:\n"," assert fri_proof['final_value'] == code_left, f\"{fri_proof['final_value']} != {code_left}\"\n","\n"," if i == 0:\n"," assert verify_decommitment(x0, code_left, mp, fri_proof['first_tree']), f\"verify_decommitment(x0={x0}, code_left={code_left}, mp={mp}, fri_proof['first_tree']={fri_proof['first_tree']})\"\n"," else:\n"," assert verify_decommitment(x0, code_left, mp, fri_proof['intermediate_trees'][i - 1]), f\"verify_decommitment(x0={x0}, code_left={code_left}, mp={mp}, fri_proof['intermediate_trees'][i - 1]={fri_proof['intermediate_trees'][i - 1]})\"\n","\n"," num_vars >>= 1\n"," q = x0"]},{"cell_type":"code","execution_count":3040,"id":"2aa20f33","metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["f=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15], D=[5*i + 10, 11*i + 2, 2*i + 11, 10*i + 5, 10*i + 26, 2*i + 20, 11*i + 29, 5*i + 21, 26*i + 21, 20*i + 29, 29*i + 20, 21*i + 26, 21*i + 5, 29*i + 11, 20*i + 2, 26*i + 10]\n","f[0] = 10 = (0 + 15)/2 + 2502028588 * (0 - 15)/(2 * 5)\n","f[1] = 16 = (1 + 14)/2 + 2502028588 * (1 - 14)/(2 * 11)\n","f[2] = 25 = (2 + 13)/2 + 2502028588 * (2 - 13)/(2 * 2)\n","f[3] = 16 = (3 + 12)/2 + 2502028588 * (3 - 12)/(2 * 10)\n","f[4] = 21 = (4 + 11)/2 + 2502028588 * (4 - 11)/(2 * 10)\n","f[5] = 7 = (5 + 10)/2 + 2502028588 * (5 - 10)/(2 * 2)\n","f[6] = 19 = (6 + 9)/2 + 2502028588 * (6 - 9)/(2 * 11)\n","f[7] = 18 = (7 + 8)/2 + 2502028588 * (7 - 8)/(2 * 5)\n","f=[10, 16, 25, 16, 21, 7, 19, 18]\n","f[0] = 10 = (10 + 18)/2 + 2239804683 * (10 - 18)/(2 * 10)\n","f[1] = 16 = (16 + 19)/2 + 2239804683 * (16 - 19)/(2 * 2)\n","f[2] = 25 = (25 + 7)/2 + 2239804683 * (25 - 7)/(2 * 11)\n","f[3] = 16 = (16 + 21)/2 + 2239804683 * (16 - 21)/(2 * 5)\n","f=[2, 26, 18, 19], D=[13, 7, 24, 18]\n","f[0] = 2 = (2 + 19)/2 + 4177269117 * (2 - 19)/(2 * 13)\n","f[1] = 26 = (26 + 18)/2 + 4177269117 * (26 - 18)/(2 * 7)\n","f=[23, 3], D=[27, 4]\n","f[0] = 23 = (23 + 3)/2 + 2350112906 * (23 - 3)/(2 * 27)\n","f=[2], D=[0]\n","max_size=16\n","q=5\n","x0=5, x1=10\n","q=5\n","x0=2, x1=5\n","q=2\n","x0=1, x1=2\n","q=1\n","x0=0, x1=1\n","q=0\n","x0=0, x1=0\n","max_size=16\n","q=5\n","x0=5, x1=10\n","q=5\n","x0=2, x1=5\n","q=2\n","x0=1, x1=2\n","q=1\n","x0=0, x1=1\n","q=0\n","x0=0, x1=0\n","max_size=16\n","q=14\n","x0=1, x1=14\n","q=1\n","x0=1, x1=6\n","q=1\n","x0=1, x1=2\n","q=1\n","x0=0, x1=1\n","q=0\n","x0=0, x1=0\n","max_size=16\n","q=1\n","x0=1, x1=14\n","q=1\n","x0=1, x1=6\n","q=1\n","x0=1, x1=2\n","q=1\n","x0=0, x1=1\n","q=0\n","x0=0, x1=0\n","i = 0, idx = 5, first_oracle[idx] = 5\n","i = 1, idx = 2, oracles[i-1][idx] = 25\n","i = 2, idx = 1, oracles[i-1][idx] = 26\n","i = 3, idx = 0, oracles[i-1][idx] = 23\n","i = 4, idx = 0, oracles[i-1][idx] = 2\n","i = 0, idx = 5, first_oracle[idx] = 5\n","i = 1, idx = 2, oracles[i-1][idx] = 25\n","i = 2, idx = 1, oracles[i-1][idx] = 26\n","i = 3, idx = 0, oracles[i-1][idx] = 23\n","i = 4, idx = 0, oracles[i-1][idx] = 2\n","i = 0, idx = 1, first_oracle[idx] = 1\n","i = 1, idx = 1, oracles[i-1][idx] = 16\n","i = 2, idx = 1, oracles[i-1][idx] = 26\n","i = 3, idx = 0, oracles[i-1][idx] = 23\n","i = 4, idx = 0, oracles[i-1][idx] = 2\n","i = 0, idx = 1, first_oracle[idx] = 1\n","i = 1, idx = 1, oracles[i-1][idx] = 16\n","i = 2, idx = 1, oracles[i-1][idx] = 26\n","i = 3, idx = 0, oracles[i-1][idx] = 23\n","i = 4, idx = 0, oracles[i-1][idx] = 2\n","------- fri -------\n","num_vars=16\n","q=5\n","x0=5, x1=10\n","i = 0, table=[5, 11, 2, 10, 10, 2, 11, 5, 26, 20, 29, 21, 21, 29, 20, 26]\n","num_vars=8\n","q=5\n","x0=2, x1=5\n","i = 1, table=[10, 2, 11, 5, 26, 20, 29, 21]\n","num_vars=4\n","q=2\n","x0=1, x1=2\n","i = 2, table=[13, 7, 24, 18]\n","num_vars=2\n","q=1\n","x0=0, x1=1\n","i = 3, table=[27, 4]\n","num_vars=1\n","q=0\n","x0=0, x1=0\n","i = 4, table=[0]\n","num_vars=16\n","q=5\n","x0=5, x1=10\n","i = 0, table=[5, 11, 2, 10, 10, 2, 11, 5, 26, 20, 29, 21, 21, 29, 20, 26]\n","num_vars=8\n","q=5\n","x0=2, x1=5\n","i = 1, table=[10, 2, 11, 5, 26, 20, 29, 21]\n","num_vars=4\n","q=2\n","x0=1, x1=2\n","i = 2, table=[13, 7, 24, 18]\n","num_vars=2\n","q=1\n","x0=0, x1=1\n","i = 3, table=[27, 4]\n","num_vars=1\n","q=0\n","x0=0, x1=0\n","i = 4, table=[0]\n","num_vars=16\n","q=14\n","x0=1, x1=14\n","i = 0, table=[5, 11, 2, 10, 10, 2, 11, 5, 26, 20, 29, 21, 21, 29, 20, 26]\n","num_vars=8\n","q=1\n","x0=1, x1=6\n","i = 1, table=[10, 2, 11, 5, 26, 20, 29, 21]\n","num_vars=4\n","q=1\n","x0=1, x1=2\n","i = 2, table=[13, 7, 24, 18]\n","num_vars=2\n","q=1\n","x0=0, x1=1\n","i = 3, table=[27, 4]\n","num_vars=1\n","q=0\n","x0=0, x1=0\n","i = 4, table=[0]\n","num_vars=16\n","q=1\n","x0=1, x1=14\n","i = 0, table=[5, 11, 2, 10, 10, 2, 11, 5, 26, 20, 29, 21, 21, 29, 20, 26]\n","num_vars=8\n","q=1\n","x0=1, x1=6\n","i = 1, table=[10, 2, 11, 5, 26, 20, 29, 21]\n","num_vars=4\n","q=1\n","x0=1, x1=2\n","i = 2, table=[13, 7, 24, 18]\n","num_vars=2\n","q=1\n","x0=0, x1=1\n","i = 3, table=[27, 4]\n","num_vars=1\n","q=0\n","x0=0, x1=0\n","i = 4, table=[0]\n"]}],"source":["f = [i for i in range(len(D_standard_coset_4))]\n","query_num = 4\n","fri_proof = fri(f, D_standard_coset_4, len(D_standard_coset_4), query_num, MerlinTranscript(b'FRI'), debug=True)\n","print(\"------- fri -------\")\n","T = []\n","domain = D_standard_coset_4[:]\n","for i in range(log_2(len(domain)) + 1):\n"," if i != 0:\n"," T.append([x for x, _ in domain])\n"," domain = sq(domain)[:len(domain)//2]\n"," else:\n"," T.append([y for _, y in domain])\n"," domain = domain[:len(domain)//2]\n","\n","defri(fri_proof, len(D_standard_coset_4), T, query_num, MerlinTranscript(b'FRI'), debug=True)"]}],"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} diff --git a/src/circle_fft.ipynb b/src/circle_fft.ipynb deleted file mode 100644 index 689c2bf..0000000 --- a/src/circle_fft.ipynb +++ /dev/null @@ -1 +0,0 @@ -{"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. = 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. = 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^2, x, x^3, y, yx^2, yx, yx^3 \\quad (N = 8) \\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^2, x, x^3 \\tag{2}\n","$$\n","\n","$$\n","y, yx^2, yx, yx^3 \\tag{3}\n","$$\n","\n","Notice that basis (2) is bit reversed order compared to normal basis, and basis (3) is like basis (2) multiplied by y.\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}