Accurate roots of unity for powers of 2?

Question

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?

Answer 1

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)]

Answer 2

Here's a recursive algorithm that generates the n roots by taking the n/2 roots and the points in between. If n is 4 or lower, it hardcodes the result (because you'll never find the two midpoints of -1 and 1 on the complex plane). Otherwise, it finds the n/2 roots, takes every two consecutive roots, finds a point on their angle bisector, and scales that point to the unit circle.

import numpy as np
#power of 2 roots of unity
#computed by recursively taking the averages of the n/2 roots
def omegas_2(n):
    #base cases
    if n == 1:
        return np.array([1])
    if n == 2:
        return np.array([1,-1])
    if n == 4: # need three points to specify the plane
        return np.array([1,1j,-1,-1j])

    halfroots = omegas_2(n//2)
    result = np.empty(n)
    for i in range(n//2):
        result[2*i] = halfroots[i];

        #take the center of the i'th and the i+1'th n/2-root
        result[2*i+1] = (halfroots[i]+halfroots[(i+1)%(n//2)]);
        result[2*i+1] /= abs(result[2*i+1]); #normalize to length 1
    return result

Answer 3

Maxima implementation:

omegan(n):=block( [ hr, hn, q, j ],
    local(ret),
    if n=1 then return([1]),
    if n=2 then return([1, -1]),
    if n=4 then return([1, %i, -1, -%i]),
    if oddp(n) then error ("length must be power of 2; found ", n),
    hn:fix(n/2),
    hr: omegan(hn),
    for i:1 thru hn do (
        ret[2*i]:hr[i],
        j: if i+1>hn then i+1-hn else i+1,
        q:hr[i]+hr[j],          
        ret[2*i+1]:rectform(q/abs(q))
    ),  
    listarray (ret)     
)

Answer 4

Pretty Print using LaTeX:

import numpy as np 
import cmath
from IPython.display import display, Markdown
from qiskit.visualization import array_to_latex
 
# Change as required:
N = 8

root = []
text = ""
A_deg = round(2*cmath.pi/N * 180/cmath.pi, 2)
A_rad = round(2.0/N,4)
text = r"For N={N}, the Complex plane is divided into {N} equal parts by an angle of ${A_rad} \pi$ or ${A_deg} ^{{\circ}}$ .".format(N=N,A_rad=A_rad,A_deg=A_deg)
display(Markdown(text))

for k in range(N):
    root.append(cmath.exp(2*cmath.pi*1j*k/N))
    
_complex = 0
_real =0
_imaginary = 0
    
for i in range(N):
    number_type = ''
    if round(root[i].real,8) == 0: 
        number_type = 'PURE IMAGINARY'
        _imaginary += 1
        
    elif round(root[i].imag,8) == 0: 
        number_type = 'PURE REAL'
        _real += 1
    else: 
        number_type = 'COMPLEX'
        _complex += 1
    
    #print("Root-",i," is \t",np.round(root[i],3)," \twhich is: ",number_type )
    temp_text = format_latex(array_to_latex( np.array( [cmath.exp(2*cmath.pi*1j*i/N)] ) , source=True) ) 
    text = r'Root-{i}: $\omega^{{{i}}} = '.format(i=i) + temp_text +  r'$' + ' which is ' + number_type
    display(Markdown(text))

print("Number of Real Roots: ",_real) 
print("Number of Imaginary Roots: ",_imaginary)  
print("Number of Complex Roots: ",_complex)  

Output (Jupyter Notebook):