Свёртка: что на самом деле умеет это преобразование

Если ты умеешь свести задачу к свёртке, знаешь, чем NTT отличается от комплексного БПФ и почему на маленьких массивах преобразование проигрывает, листай до «Теперь сам».

Как думаешь, почему это называется не «умножением многочленов»

Многочлены в прикладных задачах встречаются нечасто. Зато выражение вида

c[k] = a[0]·b[k] + a[1]·b[k−1] + ... + a[k]·b[0]

встречается везде: сумма двух случайных величин, скользящее среднее, размытие картинки, число способов набрать сумму k двумя наборами монет. Это и есть свёртка, а умножение многочленов - одна из её записей. Вопрос: как часто встречаются задачи, которые сводятся к свёртке, если специально их искать?

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

Логарифм. Умножать многозначные числа вручную тяжело, складывать - легко; логарифм переводит умножение в сложение, там ты работаешь дешёвой операцией, потом возвращаешься обратно. Преобразование делает то же самое для свёртки: в мире значений тяжёлая свёртка превращается в дешёвое поточечное умножение. Всё, что делает БПФ, - это дешёвый переход туда и обратно.

Механизм: как свести задачу к свёртке

СВЁРТКА:   c[k] = Σ a[i] · b[j]  по всем парам с i + j = k

           то же одной переменной:  c[k] = Σ_i a[i] · b[k−i],
           где i пробегает только те значения, при которых ОБА индекса
           существуют: max(0, k − len(b) + 1) <= i <= min(k, len(a) − 1)
           (в Python про это забывают и получают b[−1] - элемент С КОНЦА
            массива, то есть тихо неверный ответ)

в лоб:                  Θ(n·m)
через преобразование:   Θ(n·log n)

    1. дополнить a и b нулями до длины n ≥ len(a)+len(b)−1, степень двойки
    2. A = fft(a),  B = fft(b)                  два прямых преобразования
    3. C[j] = A[j] · B[j]                       поточечно, линия
    4. c = ifft(C)                              одно обратное

ЧТО СВОДИТСЯ К СВЁРТКЕ

умножение длинных чисел     цифры = коэффициенты, потом разнести переносы
все суммы двух множеств     индикаторы множеств, c[k] > 0 ⇔ k представим
распределение суммы         свёртка распределений двух величин
скользящее среднее          свёртка с массивом из w чисел 1/w (w - ширина окна)
размытие картинки           двумерная свёртка с ядром
поиск с джокерами           свёртка текста и перевёрнутого шаблона
число пар с разностью d     свёртка индикаторов со сдвигом

Три строки стоит расшифровать. Ядро - это маленький массив коэффициентов, с которым сворачивают данные: для размытия это, например, девять чисел 1/9 в квадрате 3 × 3, и каждый пиксель заменяется средним по своей окрестности; двумерная свёртка - та же сумма произведений, только индексы двумерные. Скользящее среднее - тот же приём в одном измерении: ядро из w чисел 1/w. Поиск с джокерами (его же называют поиском подстроки с масками) - это когда в шаблоне есть позиции «здесь подойдёт любой символ»; тогда совпадение в каждой позиции текста выражается суммой произведений, то есть свёрткой текста с перевёрнутым шаблоном.

ДВА ВИДА ПРЕОБРАЗОВАНИЯ

комплексное БПФ (FFT):      ω = e^(2πi/n), числа с плавающей точкой
                            + работает при любой длине-степени-двойки
                            − накапливает погрешность, целые надо округлять

теоретико-числовое (NTT):   ω - корень из единицы по простому модулю p
                            + арифметика целая и ТОЧНАЯ
                            − нужен p вида k·2^m + 1, длина ограничена 2^m
                            примеры: 998244353 = 7·17·2^23 + 1, корень 3
                                      65537 = 2^16 + 1,        корень 3

Правило

Любая задача вида c[k] = Σ a[i]·b[k−i] считается за Θ(n·log n): дополнить оба массива нулями до степени двойки не меньше len(a)+len(b)−1, перевести в значения, перемножить поточечно, вернуть обратно. Длину надо брать с запасом на ВСЮ длину результата - иначе старшие коэффициенты «завернутся» на младшие. Для целых ответов берут NTT по простому модулю (точная арифметика), для сигналов и вещественных данных - комплексное БПФ. Порог выгоды - сотни-тысячи элементов: у преобразования большая константа.

Разбор: умножение длинных чисел через свёртку

Число 1234 - это многочлен 4 + 3x + 2x² + 1x³ в точке x = 10. Умножение двух чисел - это умножение многочленов плюс разнесение переносов.

1234 · 5678

коэффициенты (младший первый):   a = [4, 3, 2, 1]   b = [8, 7, 6, 5]

свёртка c[k] = Σ a[i]·b[k−i]:
c[0] = 4·8 = 32
c[1] = 4·7 + 3·8 = 28 + 24 = 52
c[2] = 4·6 + 3·7 + 2·8 = 24 + 21 + 16 = 61
c[3] = 4·5 + 3·6 + 2·7 + 1·8 = 20 + 18 + 14 + 8 = 60
c[4] = 3·5 + 2·6 + 1·7 = 15 + 12 + 7 = 34
c[5] = 2·5 + 1·6 = 16
c[6] = 1·5 = 5

переносы (справа налево по разрядам):
32 → цифра 2, перенос 3
52+3 = 55 → цифра 5, перенос 5
61+5 = 66 → цифра 6, перенос 6
60+6 = 66 → цифра 6, перенос 6
34+6 = 40 → цифра 0, перенос 4
16+4 = 20 → цифра 0, перенос 2
5+2  = 7  → цифра 7

читаем от старшего: 7 006 652

Проверка: 1234 · 5678 = 7 006 652. ✓

Обрати внимание, что сама свёртка не знает про основание системы счисления: она посчитала коэффициенты, а переносы - отдельный линейный проход. Поэтому основание можно брать любое, и в реальных библиотеках длинной арифметики берут не 10, а 10⁴ или 2¹⁶, чтобы в один коэффициент уходило сразу несколько цифр, а произведение всё ещё помещалось в машинное слово.

Разбор: почему для целых берут NTT, а не комплексное БПФ

Комплексное БПФ считает в числах с плавающей точкой. Корень из единицы e^(2πi/n) - иррациональное число, его значение округляется, и ошибка накапливается по log n уровням. Для сигналов это не важно: там данные и так приближённые. Для целых - важно, потому что ответ надо округлить до целого, и если погрешность дошла до 0.5, округление даст неверную цифру.

Числа с плавающей точкой (тип double) хранят значение как мантиссу - 53 двоичные цифры значащей части - и отдельно показатель степени. Поэтому любое число известно с относительной точностью около 2⁻⁵³ ≈ 1.1·10⁻¹⁶; эту величину и называют машинным эпсилоном (обычно её берут как 2.2·10⁻¹⁶). Каждое сложение и умножение добавляет свою такую ошибку, и по log n уровням преобразования они накапливаются.

оценка погрешности комплексного БПФ:
    ошибка ~ n · log n · максимальное_значение · машинный_эпсилон

при double (эпсилон ~ 2.2e−16) и коэффициентах до 10⁶:
    n = 2¹⁰  → ошибка ~ 10⁻⁶      безопасно
    n = 2²⁰  → ошибка ~ 10⁻³      ещё терпимо
    коэффициенты до 10⁹ и n = 2²⁰ → ошибка около единицы   УЖЕ НЕТ

Второй вид преобразования обходится без вещественных чисел вовсе, но для этого надо на минуту заглянуть в арифметику по модулю.

АРИФМЕТИКА ПО МОДУЛЮ p: работаем только с остатками от деления на p
    складываем и умножаем как обычно, а результат каждый раз берём по модулю
    если p ПРОСТОЕ, у каждого остатка кроме нуля есть обратный: a · a⁻¹ ≡ 1 (mod p),
    поэтому делить тоже можно - на этом и держится «поделить на n» в обратном переходе

КОРЕНЬ ИЗ ЕДИНИЦЫ ПО МОДУЛЮ: остаток ω, у которого ω^n ≡ 1, а меньшие степени - нет
    p = 65537, n = 4, ω = 256:   256² = 65536 ≡ −1,   256⁴ ≡ 1      ✓
    свойства те же, что у комплексных: квадраты дают корни порядка n/2,
    а ω^(n/2) ≡ −1 - значит бабочка работает слово в слово

«корень 3» в таблице выше - это ПЕРВООБРАЗНЫЙ корень: степени тройки дают все
    ненулевые остатки по модулю p. Из него получают корень любого нужного
    порядка: ω = 3^((p−1)/n). Для p = 65537 и n = 4 это 65281 - тоже корень
    четвёртой степени (65281² ≡ −1), просто другой из двух

В NTT ничего не округляется: все вычисления - целые по модулю p. Корень из единицы существует ТОЧНО, если p − 1 делится на длину n. Отсюда требование к модулю: p = k·2^m + 1, и тогда любая длина n = 2^t при t ≤ m годится.

Цена - ограничения. Во-первых, ответ получается по модулю p: если настоящие коэффициенты больше, нужно либо брать модуль побольше, либо посчитать всё два-три раза по РАЗНЫМ модулям и собрать настоящий ответ по остаткам: если известны остатки числа по нескольким взаимно простымРазбирается в разделе 18, глава «Модульная арифметика, НОД и обратный элемент»: Алгоритм Евклида: не больше пяти кругов на цифру модулям, то само число восстанавливается однозначно, пока оно меньше их произведения (это и есть китайская теорема об остатках). Во-вторых, m ограничивает максимальную длину: у 65537 это m = 16, то есть до 65536 элементов; у 998244353 - m = 23, до восьми миллионов.

В задаче этой главы модуль 65537 выбран не случайно: движок курса считает в тех же числах, что браузер (double, 53 бита мантиссы), и произведение двух чисел по модулю 998244353 вышло бы порядка 10¹⁸ - за пределом точности. С модулем 65537 произведение не превосходит 4.3·10⁹ и остаётся точным. Это ровно та ловушка, в которую попадают, перенося код NTT из C++ в Python или JavaScript без проверки: код выглядит правильным, а результат тихо неверный.

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

Первое: взять длину без запаса. Если n меньше len(a)+len(b)−1, старшие коэффициенты «завернутся» и прибавятся к младшим - получится циклическая свёртка вместо обычной, ответ неверный без единой ошибки в коде. Второе: применять преобразование к маленьким массивам. Дело не в числе умножений, а в константе: три преобразования подряд, бит-реверсная перестановка, таблица корней, обращения к памяти вразброс. Замер на движке курса (задача этой главы): на 128 коэффициентах свёртка в лоб стоит 51 206 шагов против 63 847 у преобразования - лоб дешевле; на 256 уже 200 710 против 137 088 - дешевле преобразование. Порог у каждой реализации свой, и его ЗАМЕРЯЮТ, а не выводят. Третье: забыть про точность (см. разбор выше). Четвёртое: считать, что БПФ ускоряет любое перемножение массивов: поточечное умножение и так линейно, ускоряется ровно свёртка. Пятое: думать, что комплексное БПФ даёт целые числа: оно даёт 2.9999999997, и округление - обязательный шаг.

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

Свёртка - основа цифровой обработки сигналов и изображений: любой фильтр (размытие, повышение резкости, эквалайзер) - это свёртка с ядром, и для больших ядер её считают через БПФ. Спектральный анализ - это прямое применение преобразования: увидеть, из каких частот состоит звук или вибрация подшипника. В свёрточных нейросетях большие ядра тоже иногда считают через БПФ. В постквантовой криптографии NTT прописан прямо в стандартах решётчатых схем (ML-KEM, он же Kyber, и ML-DSA): там перемножают многочлены, и преобразование - основная операция. А вот в RSA и на эллиптических кривых числа слишком короткие для этого порога: там работают столбиком, Карацубой и умножением по Монтгомери. Зато в библиотеках длинной арифметики вроде GMP преобразование включается на очень больших числах - как последняя ступень той самой лестницы. В прикладном коде готовые библиотеки закрывают всё: numpy.fft, scipy.signal.fftconvolve, FFTW. Знать устройство нужно, чтобы понимать, почему fftconvolve быстрее convolve только на больших данных, откуда берутся комплексные числа в результате и почему массив дополняют нулями.

Теперь сам

Дано два множества целых: A = {1, 3, 4} и B = {2, 3}. Нужно узнать, какие суммы a + b достижимы и сколькими способами. Сведи к свёртке и посчитай.

A = {1, 3, 4}   →   a = [0, 1, 0, 1, 1]        индексы 0..4
B = {2, 3}      →   b = [0, 0, 1, 1]           индексы 0..3

c[k] = Σ a[i]·b[k−i]   - число пар с суммой k

c[3] = a[1]·b[2] = 1                           1+2
c[4] = a[1]·b[3] = 1                           1+3
c[5] = a[3]·b[2] = 1                           3+2
c[6] = a[3]·b[3] + a[4]·b[2] = 1 + 1 = 2       3+3 и 4+2
c[7] = a[4]·b[3] = 1                           4+3
остальные c[k] = 0