Дискретное преобразование Фурье

Давайте научимся быстро умножать многочлены.

Пусть у нас есть два многочлена $A(x)$ и $B(x)$ степени $n$ .

\align{ A(x) &= \sum\limits_{k=0}^n a_k \cdot x^k = a_{n} \cdot x^n + a_{n-1} \cdot x^{n-1} + \dotsb + a_2 \cdot x^2 + a_1 \cdot x + a_0 \\ B(x) &= \sum\limits_{k=0}^n b_k \cdot x^k = b_{n} \cdot x^n + b_{n-1} \cdot x^{n-1} + \dotsb + b_2 \cdot x^2 + b_1 \cdot x + b_0 }

В школе изучался неэффективный метод умножения «в столбик»

A(x) \cdot B(x) = \sum\limits_{k=0}^n x^k \sum\limits_{i+j=k} a_i b_j

Этот метод требует $(n+1)(n+2)/2 = \Theta(n^2)$ умножений.

Каждый многочлен $P(x) = \sum\limits_{k=0}^n a_k \cdot x^k$ степени $n$ может быть представлен списком своих $n+1$ коэффициентов $a_0, a_1, \dotsc, a_n$ , а может быть представлен списком значений в $n+1$ точках $\bigl( x_0, P(x_0) \bigr), \bigl( x_1, P(x_1) \bigr), \dotsc, \bigl( x_n, P(x_n) \bigr)$ .

Если два многочлена представлены значениями, а не коэффициентами, то их произведение считается быстро. Пусть многочлены $A(x)$ и $B(x)$ имеют в точках $x_0, x_1, \dotsc, x_{2n+1}$ значения $A_0, A_1, \dotsc, A_{2n+1}$ и $B_0, B_1, \dotsc, B_{2n+1}$ . Тогда многочлен $A(x) \cdot B(x)$ имеет в этих точках значения $A_0 B_0, A_1 B_1, \dotsc, A_{2n+1} B_{2n+1}$ .

Такой способ вычисления требует $2n+1 = O(n)$ умножений. Получаем следующий код

input const int n
input array[int] a[n+1]
input array[int] b[n+1]

array[int] x[2*n+1] = [???, ???, ..., ???]

array[int] A[2*n+1] = evaluate(a, x)
array[int] B[2*n+1] = evaluate(b, x)

array[int] C[2*n+1]
for int i = 0; i < 2*n+1; i++:
    C[i] = A[i] * B[i]

array[int] c[2*n+1] = interpolate(x, C)

Нам осталось придумать эффективную реализацию операций evaluate(p, x) — вычисление значений многочлена, заданного коэффициентами p, в точках x и interpolate(x, P) — восстановление коэффициентов многочлена по его значениям P в точках x.

Стандартными методами можно реализовать операцию evaluate(p, x) за $n+1$ умножений с помощью схемы Горнера и операцию interpolate(x, P) за $O(n^3)$ методами обычной интерполяции. Можно, конечно, попробовать оптимизировать интерполяцию до $O(n^2)$ , но это тоже не подходит, потому что школьный метод всё еще эффективнее.

Мы вольны выбирать узлы $x_j$ , и эта способность является ключом к построению эффективного алгоритма.

Дискретное преобразование Фурье

Выберем число $N \ge n+1$ , чтобы иметь запас на будущее (мы из двух многочленов степени $n$ хотим сделать многочлен степени $2n$ ) и чтобы было проще считать.

Для того, чтобы все операции делать эффективно, можно выбрать узлы $x_j$ следующим образом

x_j = e^{2 \pi i j / N} = w^j

где $w = e^{2 \pi / N}$ — первый комплексный корень из единицы.

Дискретное преобразование Фурье многочлена $P(x)$ , заданного коэффициентами $P(x) = \sum\limits_{j=0}^{N-1} p_j \cdot x^j$ , называется вычисление его значений в $N$ точках $w^0, w^1, w^2, \dotsc, w^{N-1}$ . При этом список коэффициентов $a_0, a_1, a_2, \dotsc, a_n$ дополнен нулями до длины $n$ .

Можно сформулировать задачу поиска дискретного преобразования Фурье в матричных терминах

\pmatrix{ ~ 1 & 1 & 1 & 1 & \cdots & 1 & 1\\[0.4em]~ 1 & w & w^2 & w^3 & \cdots & w^{N-2} & w^{N-1}\\[0.4em]~ 1 & w^2 & w^4 & w^6 & \cdots & w^{2 (N-2)} & w^{2 (N-1)}\\[0.4em]~ 1 & w^3 & w^6 & w^9 & \cdots & w^{3 (N-2)} & w^{3 (N-1)} \\ ~ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\[0.4em]~ 1 & w^{(N-2)} & w^{(N-2) \cdot 3} & w^{(N-2) \cdot 4} & \cdots & w^{(N-2) \cdot (N-2)} & w^{(N-2) \cdot (N-1)}\\[0.4em]~ 1 & w^{(N-1)} & w^{(N-1) \cdot 3} & w^{(N-1) \cdot 4} & \cdots & w^{(N-1) \cdot (N-2)} & w^{(N-1) \cdot (N-1)} } \cdot \pmatrix{p_0 \\ p_1 \\ p_2 \\ p_3 \\ \vdots \\ p_{N-2} \\ p_{N-1}} = \pmatrix{P_0 \\ P_1 \\ P_2 \\ P_3 \\ \vdots \\ P_{N-2} \\ P_{N-1}}

Алгоритм Кули-Тьюки

Алгоритм Кули-Тьюки находит дискретного преобразования Фурье многочлена $P(x) = \sum\limits_{j=0}^{n} p_j \cdot x^j$ .

Перед всеми манипуляциями нам надо дополнить массив коэффициентов $p_j$ нулями так, чтобы его длина равнялась степени двойки. То есть, если исходная степень многочлена $P(x)$ равнялась $n$ , то после дополнения его степень будет $N = 2^{\lceil \log_2 n \rceil}$ , а коэффициенты $p_j = 0$ при $j > n$ . Исходный многочлен теперь выглядит так

P(x) = \sum\limits_{j=0}^{N-1} p_j \cdot x^j

Представим многочлен в виде $P(x) = S(x^2) + x \cdot T(x^2)$ , где $S(x)$ состоит из коэффициентов при чётных степенях $P(x)$ , а $T(x)$ состоит из коэффициентов при нечётных степенях $P(x)$ .

S(x) = \sum\limits_{j=0}^{N/2-1} p_{2j} \cdot x^j \qquad T(x) = \sum\limits_{j=0}^{N/2-1} p_{2j+1} \cdot x^j

При $N = 2k$ можно использовать равенство

w^{2j} = w^{2j \bmod 2k} = w^{2 \cdot (j \bmod k)}

С помощью этого можно разделить задачу вычисления дискретного преобразования Фурье $P(x)$ к двум вдвое меньшим задачам

P(w^j) = S(w^{2j}) + w^j \cdot T(w^{2j}) = S \big( (w^2)^{j \bmod k} \big) + w^j \cdot T \big( (w^2)^{j \bmod k} \big)

function fft(array[int] p[N], complex w) -> array[N]:
    if N == 1: return P

    const int k = N / 2

    array[int] s[k], t[k]

    for int i = 0; i < N; i++:
        if i % 2 == 0:
            s[i/2] = p[i]
        else:
            t[i/2] = p[i]

    S = fft(s, w * w)
    T = fft(t, w * w)

    array[int] result[N]
    complex wt = 1

    for int i = 0; i < N; i++:
        result[i] = S[i % k] + wt * T[i % k]
        wt *= w

    return result

Другие основания

Что нам нужно было от $w^j$ ? Только то, что все они являются корнями из единицы, то есть образуют циклическую группу порядка $n$ . Мы выбрали первый попавшийся объект — дискретную подгруппу группы $\CC / \RR$ , имеющую порядок $n$ .