Nth roots of unity

Justus Perlwitz has a nice interactive post which calculates (and plots) the nth roots of unity: www.justus.pw/garden/nt… and it made me want to do something similar in R which is well known for its statistical capabilities.

The formula Justus shows for the complex roots of z^n = 1 is

z=cis(2*k*pi/n) for 0 <= k < n

and despite R having great support for complex numbers (as I explored previously) and trigonometry functions, there’s no specific cis() so I’ll have to define one (easily enough)

cis <- function(theta) cos(theta) + 1i * sin(theta)

I can then calculate the roots by evaluating this function with the arguments 2*pi*(0:(n-1))/n thanks in part to R’s vectorisation

roots <- function(n) cis(2*(0:(n-1))*pi/n)

Just to make this look a bit nicer, I can also use zapsmall() to round any tiny values towards 0, and sort() the results (first by the real part, then the imaginary part, by default for complex numbers)

roots <- function(n) sort(zapsmall(cis(2*(0:(n-1))*pi/n)))

Of course, there’s ways to actually solve the equation z^n = 1 in R… If I set up a polynomial

p(x) = z1 + z2*x + z3*x^2 + ... + zn*x^{n-1}

then rearrange it to

-1 + ... + zn*x^{n-1} = 0

then that’s equivalent to z^n = 1 and can be solved. I just need the coefficients of the left side, which are -1, n-1 lots of 0, and a 1, i.e.

c(-1, rep(0, n-1), 1)

I can pass that to the base function polyroot() and it will return the complex roots. Again, I’ll use zapsmall() and sort() for the sake of comparison

solver <- function(n) sort(zapsmall(polyroot(c(-1, rep(0, n-1), 1))))

One last option is to leverage the fact that while R doesn’t have a cis(), Euler’s formula says that it’s equivalent to an exponential, exp(ix) = cos(x) + i sin(x) so I could just use that

euler <- function(n) sort(zapsmall(exp((2*pi*1i*(0:(n-1)))/n)))

So, does it work? You bet

roots(2)
#> [1] -1+0i  1+0i
solver(2)
#> [1] -1+0i  1+0i
euler(2)
#> [1] -1+0i  1+0i

roots(3)
#> [1] -0.5-0.8660254i -0.5+0.8660254i  1.0+0.0000000i
solver(3)
#> [1] -0.5-0.8660254i -0.5+0.8660254i  1.0+0.0000000i
euler(3)
#> [1] -0.5-0.8660254i -0.5+0.8660254i  1.0+0.0000000i

roots(4)
#> [1] -1+0i  0-1i  0+1i  1+0i
solver(4)
#> [1] -1+0i  0-1i  0+1i  1+0i
euler(4)
#> [1] -1+0i  0-1i  0+1i  1+0i

roots(5)
#> [1] -0.809017-0.5877853i -0.809017+0.5877853i  0.309017-0.9510565i
#> [4]  0.309017+0.9510565i  1.000000+0.0000000i
solver(5)
#> [1] -0.809017-0.5877853i -0.809017+0.5877853i  0.309017-0.9510565i
#> [4]  0.309017+0.9510565i  1.000000+0.0000000i
euler(5)
#> [1] -0.809017-0.5877853i -0.809017+0.5877853i  0.309017-0.9510565i
#> [4]  0.309017+0.9510565i  1.000000+0.0000000i

I was also tempted to try this out in Julia, which does have a cis() built-in and also deals well with complex numbers

nroots(n) = round.([cis(2π * k / n) for k in 0:(n-1)]; digits=3)
nroots (generic function with 1 method)

nroots(4)
4-element Vector{ComplexF64}:
  1.0 + 0.0im
  0.0 + 1.0im
 -1.0 + 0.0im
 -0.0 - 1.0im

nroots(5)
5-element Vector{ComplexF64}:
    1.0 + 0.0im
  0.309 + 0.951im
 -0.809 + 0.588im
 -0.809 - 0.588im
  0.309 - 0.951im

Simply beautiful, isn’t it?