在Python中计算n次单位根

Question

我正在尝试写一个叫做croot(k, n)的算法，它可以返回单位根的第k个根，条件是n等于n。我得到的答案大部分是对的，但有些数字的表示方式看起来很奇怪，似乎不太对。这里有一个例子。

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)

输出是：

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

哇哇哇。当k等于2，n等于8时，结果是错的，应该是i，表示成1j、j或者1.00000j等等。有人能帮我一下吗？我这样做是因为我想写一个快速傅里叶变换（FFT）的算法。我对复数和Python不是很熟悉，所以可能犯了简单的错误。

谢谢，

如果你们需要更多信息，随时问我。

Answer 1

我在使用Python 3.5.2的时候，尝试计算1200个单位根时，numpy.roots耗尽了内存，导致我的Chromebook崩溃。崩溃发生在它构建多项式的伴随矩阵时，所以我觉得它没有使用稀疏表示。从文档中可以看到：

这个算法依赖于计算伴随矩阵的特征值。

如果你只是想计算单位根，使用三角函数的方法在时间和空间复杂度上都更好：

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

如果你不介意放弃并行计算，可以使用生成器，这样每增加一个根所需的空间和时间都是固定的，而第一种方法必须先计算出所有的n个根才能返回结果：

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

在我的机器上，不使用并行计算，这些方法计算第10,000,000个根的时间，和numpy.roots计算第1,000个根的时间差不多：

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

Answer 2

这里介绍了单位立方根和四次根的使用示例。输入的数组应该被理解为多项式的系数。

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

---

编辑： 多项式的系数在输入数组 p 中以以下顺序给出，传递给 np.roots(p)：

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

举个例子，如果你想得到单位根的 n 次根，也就是方程 1 * x**n - 1 == 0 的解，你可以使用这样的输入：p = [1] + [0] * (n - 1) + [-1]。

Answer 3

看看这个数字

(6.12303176911e-17+1j)

6.12303176911e-17 等于 0.0000000000000000612303176911，这个数字非常小（接近于零）。你看到的这个情况是因为浮点数的表示方式有限，导致了四舍五入的误差。

这个误差就像是你测量到太阳的距离，误差大约在10微米左右。如果你在处理来自现实世界的数据时进行快速傅里叶变换（FFT），那么测量误差通常会比这个大得多。