Корни из единицы: почему деление пополам работает

Если ты умеешь вывести разбиение на чётные и нечётные, знаешь, почему ω^{k+n/2} = −ω^k, и помнишь, чем отличается обратное преобразование, листай до «Теперь сам».

Как думаешь, какие точки выбрать

Вычислить многочлен степени n−1 в одной точке стоит O(n). В n точках - O(n²), и это ровно то, что мы хотим обойти.

Обойти можно только одним способом: сделать так, чтобы вычисления в разных точках ПЕРЕИСПОЛЬЗОВАЛИ друг друга. Значит точки должны быть связаны между собой. Вопрос: какая связь позволит считать половину работы один раз на две точки?

На что это похоже

Циферблат часов. Двенадцать делений расставлены не случайно: каждое второе - это тоже правильный набор, только из шести делений; каждое четвёртое - из трёх. Корни из единицы устроены так же: возведи их в квадрат - и n точек схлопнутся ровно в n/2 точек того же вида. Именно эта самоподобность и делает деление пополам возможным.

Механизм: чётные, нечётные и бабочка

Сначала о самих числах. Точки, которые нам понадобятся, - КОМПЛЕКСНЫЕ: к обычным числам добавляют одно новое, i, у которого i² = −1, и дальше считают по обычным правилам. Отсюда сразу i³ = i²·i = −i и i⁴ = 1. Комплексное число a + bi удобно представлять точкой на плоскости с координатами (a, b); умножение на i - это поворот этой точки на 90 градусов против часовой стрелки, а e^(2πi/n) - точка на окружности единичного радиуса, повёрнутая на 1/n полного оборота.

КОРНИ ИЗ ЕДИНИЦЫ - n точек, равномерно расставленных по окружности радиуса 1

ω = e^(2πi/n)      поворот на 1/n оборота
ω^k                поворот на k/n оборота
ω^n = 1            n поворотов по 1/n - это полный круг, вернулись в 1

n = 4:  ω = e^(2πi/4) - поворот на четверть круга, то есть ровно i
        точки: ω⁰ = 1,  ω¹ = i,  ω² = i² = −1,  ω³ = i³ = −i
n = 8:  восемь точек через 45 градусов: 1, (1+i)/√2, i, (−1+i)/√2, −1, ... 

Проверять формулы можно и не думая о геометрии: достаточно правил i² = −1 и ω^n = 1.

ДВА СВОЙСТВА, НА КОТОРЫХ ДЕРЖИТСЯ ВСЁ

(ω^k)² = ω^(2k)        квадраты n корней дают n/2 корней порядка n/2
ω^(k + n/2) = −ω^k     половина корней - «минус» другой половины

Второе свойство выводится в одну строку: ω^(n/2) = e^(2πi/n · n/2) = e^(πi), то есть поворот на ПОЛОВИНУ круга - а это −1. Значит ω^(k + n/2) = ω^k · ω^(n/2) = −ω^k. При n = 4 это видно и руками: ω^(1+2) = i³ = −i = −ω¹.

РАЗБИЕНИЕ

A(x) = a₀ + a₁x + a₂x² + ... 
     = (a₀ + a₂x² + a₄x⁴ + ...) + x·(a₁ + a₃x² + a₅x⁴ + ...)
     = E(x²) + x·O(x²)

где E - многочлен из ЧЁТНЫХ коэффициентов, O - из НЕЧЁТНЫХ, оба степени n/2 − 1

БАБОЧКА: считаем сразу пару точек

A(ω^k)         = E(ω^{2k}) + ω^k · O(ω^{2k})
A(ω^{k+n/2})   = E(ω^{2k}) − ω^k · O(ω^{2k})       потому что ω^{k+n/2} = −ω^k

то есть из двух чисел E и O получаются ДВА значения - одно сложением, другое вычитанием

цена:  T(n) = 2·T(n/2) + n   →   Θ(n·log n)
обратное преобразование: то же самое с ω^{−1} вместо ω и делением всего на n

Правило

Точки выбираются как корни из единицы ω^k, где ω^n = 1. Их самоподобность ((ω^k)² пробегает корни порядка n/2) позволяет разрезать многочлен на чётные и нечётные коэффициенты: A(x) = E(x²) + x·O(x²). Тогда значения в n точках получаются из значений двух половин в n/2 точках, причём пара A(ω^k) и A(ω^{k+n/2}) считается из одних и тех же E и O сложением и вычитанием - это БАБОЧКА. Отсюда T(n) = 2·T(n/2) + n и Θ(n·log n). Обратное преобразование устроено так же, но с ω^{−1} и делением результата на n. Длина должна быть степенью двойки - иначе деление пополам не доходит до конца.

Разбор: преобразование четырёх коэффициентов руками

Возьмём A(x) = 1 + 2x + 3x² + 4x³ и n = 4. Корни четвёртой степени из единицы: 1, i, −1, −i.

Сначала честно, в лоб:

A(1)  = 1 + 2 + 3 + 4 = 10
A(i)  = 1 + 2i − 3 − 4i = −2 − 2i
A(−1) = 1 − 2 + 3 − 4 = −2
A(−i) = 1 − 2i − 3 + 4i = −2 + 2i

Теперь через разбиение. Чётные коэффициенты - [1, 3], нечётные - [2, 4]:

E(y) = 1 + 3y        O(y) = 2 + 4y

нужные значения E и O - в квадратах наших точек:
    1² = 1,  i² = −1,  (−1)² = 1,  (−i)² = −1     → всего ДВЕ точки: 1 и −1

E(1) = 4    E(−1) = −2
O(1) = 6    O(−1) = −2

Собираем бабочками:

A(1)  = E(1)  + 1·O(1)  = 4 + 6   = 10            ← пара k = 0
A(−1) = E(1)  − 1·O(1)  = 4 − 6   = −2            ← и k + n/2 = 2

A(i)  = E(−1) + i·O(−1) = −2 + i·(−2) = −2 − 2i   ← пара k = 1
A(−i) = E(−1) − i·O(−1) = −2 + 2i                 ← и k + n/2 = 3

Все четыре значения совпали с честным счётом. Посчитаем экономию по умножениям. В лоб: четыре точки, в каждой многочлен степени 3 по Горнеру - три умножения, итого 12. Через разбиение: два многочлена степени 1 (E и O) в двух точках каждый - это 4 вычисления по одному умножению, плюс четыре бабочки по одному умножению ω^k · O, итого 8. На четырёх точках выигрыш скромный, и это нормально: он берётся не с одного уровня, а с рекурсии - на тысяче точек разница уже между 10⁶ и 10⁴.

Разбор: почему обратное преобразование - почти то же самое

Прямое преобразование - это умножение вектора коэффициентов на матрицу из степеней ω. Матрица здесь - просто таблица чисел, а «умножить на неё вектор» значит посчитать n сумм: j-е значение равно Σ_k W[j][k] · a_k. Проверь, что это ровно наше определение: A(ω^j) = Σ_k a_k · (ω^j)^k.

значения = W · коэффициенты,     где W[j][k] = ω^{jk}

Обратная матрица W⁻¹ - такая таблица, что умножение на неё возвращает исходный вектор: сначала W, потом W⁻¹, и получили то, с чего начали. В общем случае её ищут решением системы, но здесь она известна заранее.

Обратное - это умножение на обратную матрицу, и у неё поразительно простой вид:

W⁻¹[j][k] = ω^{−jk} / n

То есть обратное преобразование - это ТО ЖЕ САМОЕ преобразование, только с ω^{−1} вместо ω, а в конце всё делится на n. Никакого решения системы, никакой интерполяции Лагранжа - тот же алгоритм, запущенный второй раз.

Почему так выходит: сумма Σ_k ω^{(j−m)k} по всем k равна n, если j = m, и нулю иначе - потому что при j ≠ m это сумма геометрической прогрессии с знаменателем, отличным от единицы, а все корни из единицы в сумме дают ноль. Именно эта ортогональность и делает обратный переход бесплатным.

Практический вывод: в коде прямое и обратное преобразование - одна функция с флагом. Ровно так устроена и готовая функция в задаче этой главы: вызов «вперёд» и вызов «назад» отличаются двумя строчками. Работает она не с комплексными числами, а с их целочисленным двойником - корнями из единицы по простому модулю; об этом следующий урокРазбирается в разделе 13, глава «Быстрое преобразование Фурье»: Свёртка: что на самом деле умеет это преобразование.

Что ломается без этого

Первое: брать длину, не равную степени двойки, - деление пополам не дойдёт до конца, и алгоритм либо сломается, либо потребует отдельной обработки остатка. Второе: забыть поделить на n в обратном преобразовании - коэффициенты будут в n раз больше нужного. Третье: перепутать ω и ω^{−1} местами: результат окажется «перевёрнутым» по индексам. Четвёртое: считать корни приближённо и накопить ошибку: у комплексного БПФ значения вычисляются в числах с плавающей точкой, и на больших n требуется аккуратность - именно поэтому для целочисленных задач берут преобразование по модулю, где корень из единицы существует точно. Пятое: думать, что бабочка - это оптимизация: без неё n значений собирались бы независимо, и лишний множитель n вернулся бы.

Зачем это в работе

Прямо писать БПФ почти никогда не приходится - оно есть в любой библиотеке (numpy.fft, scipy.signal, FFTW). Понимать устройство нужно, чтобы правильно им пользоваться: почему длину дополняют до степени двойки, почему результат комплексный, откуда берётся деление на n, что значит «частота». Кстати о ней: если массив - это записанный звук, то j-е значение преобразования показывает, насколько сильно в этом звуке присутствует колебание, укладывающееся в запись ровно j раз. Это и называют частотой; ноль отвечает за постоянную составляющую, а середина массива - за самые быстрые колебания, какие вообще различимы при такой частоте записи. А идея самоподобных точек встречается шире: на ней же стоят быстрые преобразования Уолша-Адамара (свёртка по исключающему или), преобразования в конечных полях и даже некоторые схемы в криптографии. Общий принцип, который стоит унести: если хочешь применить деление пополам к вычислению в n точках, подбирай точки так, чтобы они были самоподобны.

Теперь сам

Дан многочлен A(x) = 5 + 0x + 1x² + 0x³, n = 4, корни 1, i, −1, −i. Посчитай все четыре значения через разбиение на чётные и нечётные.

E(y) = 5 + y       O(y) = 0

квадраты точек: 1² = 1, i² = −1, (−1)² = 1, (−i)² = −1  → точки 1 и −1

E(1) = 6      E(−1) = 4
O(1) = 0      O(−1) = 0

A(1)  = E(1)  + 1·O(1)  = 6 + 0 = 6
A(−1) = E(1)  − 1·O(1)  = 6 − 0 = 6
A(i)  = E(−1) + i·O(−1) = 4 + 0 = 4
A(−i) = E(−1) − i·O(−1) = 4 − 0 = 4