3

因此，我正在尝试编写一个算法 croot(k, n)，它返回 n == n 的第 k 个统一根。我得到的答案大多是正确的，但它给了我非常奇怪的表示，对于某些数字来说似乎是错误的。这是一个例子。

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 不是很有经验，所以我很可能犯了一个简单的错误。

谢谢，

如果你们需要任何其他信息，请询问。

4

3 回答 3

7

这是一个单位的立方根和第四个根的用法示例。输入数组应被解释为多项式系数。

>>> 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方程 的解的 th 单位根1 * x**n - 1 == 0，您可以使用类似 的输入p = [1] + [0] * (n - 1) + [-1]。

于 2013-03-15T04:07:38.567 回答
6

看看这个数字

(6.12303176911e-17+1j)

6.12303176911e-17=0.0000000000000000612303176911非常小（接近于零）。您看到的是由于浮点数的有限表示而导致的舍入错误

误差相当于测量到太阳的距离在10微米左右。如果您对来自现实世界的数据运行 FFT，则测量误差通常远大于此。

于 2013-03-15T04:01:39.637 回答
2

使用 Python 3.5.2，numpy.roots当我尝试确定第 1,200 个统一根时，内存不足并导致我的 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))

在我的机器上并且不使用并行化，这些方法在计算第 1,000 个根所需的时间内计算第 10,000,000numpy.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
于 2017-01-09T23:49:08.890 回答