The nth roots of unity are the solutions to the polynomial equation x^n = 1. Is there a good known algorithm for n = 2^k, for some k (i.e. where n is a power of 2)?
There are many algorithms for calculating nth roots of unity. For example, this is a Python implementation using Numpy's "roots" function to return an array of roots.
import numpy as np;
def omegas(n):
#make the n-degree polynomial and solve for its roots
poly = [0]*(n+1)
poly[0] = -1; #constant
poly[-1] = 1; #highest degree
return np.roots(poly)
You can also use trig functions:
import numpy as np
import cmath
def trig_omegas(n):
return np.array([cmath.rect(1,x*np.pi) for x in range(n)])
But the accuracy leaves me wanting. This is what the answer SHOULD be, for n=4:
array([-1.+0.j, 0.+1.j, -0.-1.j, 1.+0.j])
#or, in counterclockwise order
array([ 1.+0.j, 0.+1.j, -1.+0.j, -0.-1.j])
And this is the result of the above functions.
>>> omegas(4)
array([ -1.00000000e+00+0.j, 8.32667268e-17+1.j, 8.32667268e-17-1.j,
1.00000000e+00+0.j])
>>> trig_omegas(4)
array([ 1. +0.00000000e+00j, -1. +1.22464680e-16j, 1. -2.44929360e-16j,
-1. +3.67394040e-16j])
The trailing zeros indicate that there's a small error at the end. The first entry of omegas(4) is actually a little smaller than -1.
omegas(4)[0] (-1.0000000000000004+0j)
Is there a better way to get roots of unity of powers of 2?
You might like the library Sympy. Try solving for the fifth roots of unity
solve(x**5-1)
in the console at http://live.sympy.org/
To convert the exact solutions to floating point approximations, do [x.evalf() for x in solve(x**5-1)]