feat: add SageMath

README.md

@@ -15,6 +15,11 @@
   <!-- [![Documentation](https://docs.rs/ronkathon/badge.svg)](https://docs.rs/ronkathon) -->
   </div>
 
+## Math
+To see computations used in the background, go to the `math/` directory.
+From there, you can run the `.sage` files in a SageMath environment.
+In particular, the `math/field.sage` computes roots of unity in the `PlutoField` which is of size 101.
+
 ## License
 Licensed under your option of either:
 - Apache License, Version 2.0 ([LICENSE-APACHE](LICENSE-APACHE) or http://www.apache.org/licenses/LICENSE-2.0)

math/field.sage

@@ -0,0 +1,44 @@
+# Define a finite field, our `PlutoField` which is 
+# a field of 101 elements (101 being prime).
+F = GF(101)
+print(F)
+
+# Generate primitive element
+primitive_element = F.primitive_element()
+print("The primitive element is: ", primitive_element)
+
+######################################################################
+# Let's find a mth root of unity (for m = 5)
+# First, check that m divides 101 - 1 = 100
+m = 5
+assert (101 - 1) % m == 0
+quotient = (101 - 1) // m
+print("The quotient is: ", quotient)
+
+# Find a primitive root of unity using the formula:
+# omega = primitive_element^quotient
+omega_m = primitive_element^quotient
+print("The ", m, "th root of unity is: ", omega_m)
+
+# Check that this is actually a root of unity:
+assert omega_m^m == 1
+print(omega_m, "^", m, " = ", omega_m^m)
+######################################################################
+
+######################################################################
+# Let's find a mth root of unity (for n = 25)
+# First, check that m divides 101 - 1 = 100
+n = 25
+assert (101 - 1) % n == 0
+quotient = (101 - 1) // n
+print("The quotient is: ", quotient)
+
+# Find a primitive root of unity using the formula:
+# omega = primitive_element^quotient
+omega_n = primitive_element^quotient
+print("The ", n, "th root of unity is: ", omega_n)
+
+# Check that this is actually a root of unity:
+assert omega_n^n == 1
+print(omega_n, "^", n, " = ", omega_n^n)
+######################################################################