Вычисление nth Roots of Unity в Python

Question

Итак, я пытаюсь написать алгоритм croot (k, n), который возвращает корень kth из единицы с 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 и т. Д. Может ли кто-нибудь помочь мне здесь? Я делаю это, потому что я пытаюсь написать алгоритм БПФ. Я не очень опытен с комплексными числами и Python, поэтому я мог бы очень ошибиться.

Спасибо,

Если вы, ребята, нужна дополнительная информация просто спросить.

Answer 1

Посмотрите на этот номер

(6.12303176911e-17+1j) 

6.12303176911e-17 = 0.0000000000000000612303176911 который действительно мал (близок к нулю). То, что вы видите, - это ошибки округления из-за ограниченного представления чисел с плавающей запятой.

Ошибка эквивалентна измерению расстояния до солнца до 10 мкм или около того. Если вы используете БПФ на данных из реального мира, ошибки измерения обычно намного больше, чем это.

Answer 2

Вот кубические корни единства и 4-й корни для примера использования. Входной массив следует интерпретировать как полиномиальные коэффициенты.

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

Использование Python 3.5.2, numpy.roots закончилось с памятью и разбило мой Chromebook, когда я попытался определить 1200-й корень единства. Сбой произошел, когда они построили матрицу компаньонов для полинома, поэтому я не думаю, что они используют разреженное представление.Из документов:

Алгоритм опирается на вычисления собственных компаньона матрицы

Если вам просто нужно вычислить корни из единицы, тригонометрический подход асимптотически лучше во времени и пространстве сложности :

def nthRootsOfUnity1(n): # linear space, parallelizable 
    from numpy import arange, pi, exp 
    return exp(2j * pi/n * arange(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 вычислить корни: 1000-

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