3

So, I'm trying to write an algorithm croot(k, n), that returns the kth root of unity with n == n. I'm getting mostly the right answer, however it's giving me really weird representations that seem wrong for certain numbers. Here is an example.

import cmath

def croot(k, n):
    if n<=0:
        return None
    return cmath.exp((2 * cmath.pi * 1j * k) / n)


for k in range(8):
    print croot(k, 8)

Output is:

(1+0j)
(0.70710...+0.70710...j)
(6.12323399574e-17+1j)

Whoa whoa whoa. So the root when k = 2 and n = 8 is wrong, as it should just be i, which would be represented like 1j, or j, or 1.00000j, etc. Could somebody help me here? I'm doing this because I'm trying to write an FFT algorithm. I'm not very experienced with complex numbers and Python so I could very well be making a simple mistake.

Thanks,

If you guys need any additional information just ask.

robins35
  • 653
  • 12
  • 37
  • 1
    Just so you know, you already have this functionality in numpy.. – wim Mar 15 '13 at 04:04
  • Wait, did you mean functionality for roots of unity, or for fft, because this is for an assignment and I need to write the fft part myself. – robins35 Mar 15 '13 at 04:11

4 Answers4

7

Look at this number

(6.12303176911e-17+1j)

6.12303176911e-17 = 0.0000000000000000612303176911 which is really small (close to zero). What you are seeing is rounding errors due to the limited representation of floating point numbers

The error is equivalent to measuring the distance to the sun to within 10 microns or so. If you're running FFTs on data from the real world the measurement errors are usually far larger than this.

John La Rooy
  • 295,403
  • 53
  • 369
  • 502
  • What? 6.12303176911e = 16.644 according to my calculator – robins35 Mar 15 '13 at 04:06
  • 1
    @Scriptonaut: The "e" here means "times ten to the power of", not the base of the natural logarithms. It's ~6.123 * 10^(-17). – DSM Mar 15 '13 at 04:07
  • 2
    Try typing `1e-3` and `1e-4` into interactive interpreter and you'll get the idea ;) – wim Mar 15 '13 at 04:11
7

Here's cube roots of unity and 4th roots for a usage example. The input array should be interpreted as polynomial coefficients.

>>> import numpy as np
>>> np.roots([1, 0, 0, -1])
array([-0.5+0.8660254j, -0.5-0.8660254j,  1.0+0.j       ])
>>> np.roots([1, 0, 0, 0, -1])
array([ -1.00000000e+00+0.j,   5.55111512e-17+1.j,   5.55111512e-17-1.j,
         1.00000000e+00+0.j])

edit: The polynomial coefficients are given in the input array p to np.roots(p) in the following order:

p[0] * x**n + p[1] * x**(n-1) + ... + p[n-1]*x + p[n]

So, for example, to return the nth roots of unity, which are the solutions of the equation 1 * x**n - 1 == 0, you would use an input like p = [1] + [0] * (n - 1) + [-1].

wim
  • 338,267
  • 99
  • 616
  • 750
3

Using Python 3.5.2, numpy.roots ran out of memory and crashed my Chromebook when I tried to determine the 1,200th roots of unity. The crash occurred when they construct the companion matrix to the polynomial, so I don't think they're using a sparse representation. From the docs:

The algorithm relies on computing the eigenvalues of the companion matrix

If you just need to compute the roots of unity, a trigonometric approach is asymptotically better in both time and space complexity:

def nthRootsOfUnity1(n): # linear space, parallelizable
    from numpy import arange, pi, exp
    return exp(2j * pi / n * arange(n))

If you're willing to give up on parallelization, you can use generators to achieve constant space and constant time for each additional root, whereas the first method must compute all n roots before returning:

def nthRootsOfUnity2(n): # constant space, serial
    from cmath import exp, pi
    c = 2j * pi / n
    return (exp(k * c) for k in range(n))

On my machine and without using parallelization, these methods compute the 10,000,000th roots in the time it takes numpy.roots to compute the 1,000th roots:

In [3]: n = 1000

In [4]: %time _ = sum(numpy.roots([1] + [0]*(n - 1) + [-1]))
CPU times: user 40.7 s, sys: 811 ms, total: 41.5 s
Wall time: 10.8 s

In [5]: n = 10000000

In [6]: %time _ = sum(nthRootsOfUnity1(n))
CPU times: user 4.98 s, sys: 256 ms, total: 5.23 s
Wall time: 5.23 s

In [7]: %time _ = sum(nthRootsOfUnity2(n))
CPU times: user 11.6 s, sys: 2 ms, total: 11.6 s
Wall time: 11.6 s
GeoffChurch
  • 395
  • 1
  • 13
1

If you happen to use python on Jupyter Notebook you can have some sleek LaTeX rendering using Qiskit.visualization > array_to_latex() like this.

Sample:

enter image description here

Rajesh Swarnkar
  • 601
  • 1
  • 6
  • 18