MAXimal

добавлено: 10 Jun 2008 19:04
редактировано: 9 Nov 2012 12:38

Быстрое преобразование Фурье за O (N log N). Применение к умножению двух полиномов или длинных чисел

Здесь мы рассмотрим алгоритм, который позволяет перемножить два полинома длиной $n$ за время $O(n \log n)$ , что значительно лучше времени $O(n^2)$ , достигаемого тривиальным алгоритмом умножения. Очевидно, что умножение двух длинных чисел можно свести к умножению полиномов, поэтому два длинных числа также можно перемножить за время $O(n \log n)$ .

Изобретение Быстрого преобразования Фурье приписывается Кули (Coolet) и Таки (Tukey) — 1965 г. На самом деле БПФ неоднократно изобреталось до этого, но важность его в полной мере не осознавалась до появления современных компьютеров. Некоторые исследователи приписывают открытие БПФ Рунге (Runge) и Кёнигу (Konig) в 1924 г. Наконец, открытие этого метода приписывается ещё Гауссу (Gauss) в 1805 г.

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

Пусть имеется многочлен $n$ -ой степени:

$A(x) = a_0 x^0 + a_1 x^1 + \ldots + a_{n-1} x^{n-[...]$

Не теряя общности, можно считать, что $n$ является степенью 2. Если в действительности $n$ не является степенью 2, то мы просто добавим недостающие коэффициенты, положив их равными нулю.

Из теории функций комплексного переменного известно, что комплексных корней $n$ -ой степени из единицы существует ровно $n$ . Обозначим эти корни через $w_{n,k}, k = 0 \ldots {n-1}$ , тогда известно, что $w_{n,k} = e^{ i \frac{ 2 \pi k }{ n } }$ . Кроме того, один из этих корней $w_n = w_{n,1} = e^{ i \frac{ 2 \pi }{ n } }$ (называемый главным значением корня $n$ -ой степени из единицы) таков, что все остальные корни являются его степенями: $w_{n,k} = (w_n)^k$ .

Тогда дискретным преобразованием Фурье (ДПФ) (discrete Fourier transform, DFT) многочлена $A(x)$ (или, что то же самое, ДПФ вектора его коэффициентов $(a_0, a_1, \dots, a_{n-1})$ ) называются значения этого многочлена в точках $x = w_{n,k}$ , т.е. это вектор:

${\rm DFT}(a_0, a_1, \ldots, a_{n-1}) = (y_0, y_1,[...]$
$= (A(w_n^0), A(w_n^1), \ldots, A(w_n^{n-1})).$

Аналогично определяется и обратное дискретное преобразование Фурье (InverseDFT). Обратное ДПФ для вектора значений многочлена $(y_0, y_1, \ldots y_{n-1})$ — это вектор коэффициентов многочлена $(a_0, a_1, \ldots, a_{n-1})$ :

${\rm InverseDFT}(y_0, y_1, \ldots, y_{n-1}) = (a_[...]$

Таким образом, если прямое ДПФ переходит от коэффициентов многочлена к его значениям в комплексных корнях $n$ -ой степени из единицы, то обратное ДПФ — наоборот, по значениям многочлена восстанавливает коэффициенты многочлена.

Применение ДПФ для быстрого умножения полиномов

Пусть даны два многочлена $A$ и $B$ . Посчитаем ДПФ для каждого из них: ${\rm DFT}(A)$ и ${\rm DFT}(B)$ — это два вектора-значения многочленов.

Теперь, что происходит при умножении многочленов? Очевидно, в каждой точке их значения просто перемножаются, т.е.

$(A \times B)(x) = A(x) \times B(x).$

Но это означает, что если мы перемножим вектора ${\rm DFT}(A)$ и ${\rm DFT}(B)$ , просто умножив каждый элемент одного вектора на соответствующий ему элемент другого вектора, то мы получим не что иное, как ДПФ от многочлена $A \times B$ :

${\rm DFT} (A \times B) = {\rm DFT} (A) \times {\r[...]$

Наконец, применяя обратное ДПФ, получаем:

$A \times B = {\rm InverseDFT}( {\rm DFT} (A) \tim[...]$

где, повторимся, справа под произведением двух ДПФ понимается попарные произведения элементов векторов. Такое произведение, очевидно, требует для вычисления только $O(n)$ операций. Таким образом, если мы научимся вычислять ДПФ и обратное ДПФ за время $O(n \log n)$ , то и произведение двух полиномов (а, следовательно, и двух длинных чисел) мы сможем найти за ту же асимптотику.

Следует заметить, что, во-первых, два многочлена следует привести к одной степени (просто дополнив коэффициенты одного из них нулями). Во-вторых, в результате произведения двух многочленов степени $n$ получается многочлен степени $2n-1$ , поэтому, чтобы результат получился корректным, предварительно нужно удвоить степени каждого многочлена (опять же, дополнив их нулевыми коэффициентами).

Быстрое преобразование Фурье

Быстрое преобразование Фурье (fast Fourier transform) — это метод, позволяющий вычислять ДПФ за время $O(n \log n)$ . Этот метод основывается на свойствах комплексных корней из единицы (а именно, на том, что степени одних корней дают другие корни).

Основная идея БПФ заключается в разделении вектора коэффициентов на два вектора, рекурсивном вычислении ДПФ для них, и объединении результатов в одно БПФ.

Итак, пусть имеется многочлен $A(x)$ степени $n$ , где $n$ — степень двойки, и $n>1$ :

$A(x) = a_0 x^0 + a_1 x^1 + \ldots + a_{n-1} x^{n-[...]$

Разделим его на два многочлена, один — с чётными, а другой — с нечётными коэффициентами:

$A_0(x) = a_0 x^0 + a_2 x^1 + \ldots + a_{n-2} x^{[...]$
$A_1(x) = a_1 x^0 + a_3 x^1 + \ldots + a_{n-1} x^{[...]$

Нетрудно убедиться, что:

$A(x) = A_0(x^2) + x A_1(x^2). ~~~~~~~(1)$

Многочлены $A_0$ и $A_1$ имеют вдвое меньшую степень, чем многочлен $A$ . Если мы сможем за линейное время по вычисленным ${\rm DFT}(A_0)$ и ${\rm DFT}(A_1)$ вычислить ${\rm DFT}(A)$ , то мы и получим искомый алгоритм быстрого преобразования Фурье (т.к. это стандартная схема алгоритма "разделяй и властвуй", и для неё известна асимптотическая оценка $O(n \log n)$ ).

Итак, пусть мы имеем вычисленные вектора $\{ y_k^0 \}_{k=0}^{n/2-1} = {\rm DFT}(A_0)$ и $\{ y_k^1 \}_{k=0}^{n/2-1} = {\rm DFT}(A_1)$ . Найдём выражения для $\{ y_k \}_{k=0}^{n-1} = {\rm DFT}(A)$ .

Во-первых, вспоминая (1), мы сразу получаем значения для первой половины коэффициентов:

$y_k = y_k^0 + w_n^k y_k^1, ~~~~k = 0 \ldots n/2-1[...]$

Для второй половины коэффициентов после преобразований также получаем простую формулу:

$y_{k+n/2} = A(w_n^{k+n/2}) = A_0(w_n^{2k+n}) + w_[...]$
$= A_0(w_n^{2k}) - w_n^k A_1(w_n^{2k}) = y_k^0 - w[...]$

(Здесь мы воспользовались (1), а также тождествами $w_n^n = 1$ , $w_n^{n/2} = -1$ .)

Итак, в результате мы получили формулы для вычисления всего вектора $\{ y_k \}$ :

$y_k = y_k^0 + w_n^k y_k^1, \ \ \ \ k = 0 \ldots n[...]$
$y_{k+n/2} = y_k^0 - w_n^k y_k^1, \ \ \ \ k = 0 \l[...]$

(эти формулы, т.е. две формулы вида $a+bc$ и $a-bc$ , иногда называют "преобразование бабочки" ("butterfly operation"))

Тем самым, мы окончательно построили алгоритм БПФ.

Обратное БПФ

Итак, пусть дан вектор $(y_0, y_1, \ldots, y_{n-1})$ — значения многочлена $A$ степени $n$ в точках $x = w_n^k$ . Требуется восстановить коэффициенты $(a_0, a_1, \ldots, a_{n-1})$ многочлена. Эта известная задача называется интерполяцией, для этой задачи есть и общие алгоритмы решения, однако в данном случае будет получен очень простой алгоритм (простой тем, что он практически не отличается от прямого БПФ).

ДПФ мы можем записать, согласно его определению, в матричном виде:

$\pmatrix{ w_n^0 & w_n^0 & w_n^0 & w_n^0 & \cdots [...]$

Тогда вектор $(a_0, a_1, \ldots, a_{n-1})$ можно найти, умножив вектор $(y_0, y_1, \ldots, y_{n-1})$ на обратную матрицу к матрице, стоящей слева (которая, кстати, называется матрицей Вандермонда):

$\pmatrix{ a_0 \cr a_1 \cr a_2 \cr a_3 \cr \vdots [...]$

Непосредственной проверкой можно убедиться в том, что эта обратная матрица такова:

$\frac{1}{n} \pmatrix{ w_n^0 & w_n^0 & w_n^0 & w_n[...]$

Таким образом, получаем формулу:

$a_k = \frac{1}{n} \sum_{j=0}^{n-1} y_j w_n^{-kj}.[...]$

Сравнивая её с формулой для $y_k$ :

$y_k = \sum_{j=0}^{n-1} a_j w_n^{kj},$

мы замечаем, что эти две задачи почти ничем не отличаются, поэтому коэффициенты $a_k$ можно находить таким же алгоритмом "разделяй и властвуй", как и прямое БПФ, только вместо $w_n^k$ везде надо использовать $w_n^{-k}$ , а каждый элемент результата надо разделить на $n$ .

Таким образом, вычисление обратного ДПФ почти не отличается от вычисления прямого ДПФ, и его также можно выполнять за время $O(n \log n)$ .

Реализация

Рассмотрим простую рекурсивную реализацию БПФ и обратного БПФ, реализуем их в виде одной функции, поскольку различия между прямым и обратным БПФ минимальны. Для хранения комплексных чисел воспользуемся стандартным в C++ STL типом complex (определённым в заголовочном файле <complex>).

typedef complex<double> base;
 
void fft (vector<base> & a, bool invert) {
	int n = (int) a.size();
	if (n == 1)  return;
 
	vector<base> a0 (n/2),  a1 (n/2);
	for (int i=0, j=0; i<n; i+=2, ++j) {
		a0[j] = a[i];
		a1[j] = a[i+1];
	}
	fft (a0, invert);
	fft (a1, invert);
 
	double ang = 2*PI/n * (invert ? -1 : 1);
	base w (1),  wn (cos(ang), sin(ang));
	for (int i=0; i<n/2; ++i) {
		a[i] = a0[i] + w * a1[i];
		a[i+n/2] = a0[i] - w * a1[i];
		if (invert)
			a[i] /= 2,  a[i+n/2] /= 2;
		w *= wn;
	}
}

В аргумент $\rm a$ функции передаётся входной вектор коэффициентов, в нём же и будет содержаться результат. Аргумент $\rm invert$ показывает, прямое или обратное ДПФ следует вычислить. Внутри функции сначала проверяется, что если длина вектора $\rm a$ равна единице, то ничего делать не надо - он сам и является ответом. Иначе вектор $\rm a$ разделяется на два вектора $\rm a0$ и $\rm a1$ , для которых рекурсивно вычисляется ДПФ. Затем вычисляется величина $w_n$ , и заводится переменная $w$ , содержащая текущую степень $w_n$ . Затем вычисляются элементы результирующего ДПФ по вышеописанным формулам.

Если указан флаг $\rm invert = true$ , то $w_n$ заменяется на $w_n^{-1}$ , а каждый элемент результата делится на 2 (учитывая, что эти деления на 2 произойдут в каждом уровне рекурсии, то в итоге как раз получится, что все элементы поделятся на $n$ ).

Тогда функция для перемножения двух многочленов будет выглядеть следующим образом:

void multiply (const vector<int> & a, const vector<int> & b, vector<int> & res) {
	vector<base> fa (a.begin(), a.end()),  fb (b.begin(), b.end());
	size_t n = 1;
	while (n < max (a.size(), b.size()))  n <<= 1;
	n <<= 1;
	fa.resize (n),  fb.resize (n);
 
	fft (fa, false),  fft (fb, false);
	for (size_t i=0; i<n; ++i)
		fa[i] *= fb[i];
	fft (fa, true);
 
	res.resize (n);
	for (size_t i=0; i<n; ++i)
		res[i] = int (fa[i].real() + 0.5);
}

Эта функция работает с многочленами с целочисленными коэффициентами (хотя, понятно, теоретически ничто не мешает ей работать и с дробными коэффициентами). Однако здесь проявляется проблема большой погрешности при вычислении ДПФ: погрешность может оказаться значительной, поэтому округлять числа лучше самым надёжным способом — прибавлением 0.5 и последующим округлением вниз (внимание: это будет работать неправильно для отрицательных чисел, если таковые могут появиться в вашем применении).

Наконец, функция для перемножения двух длинных чисел практически ничем не отличается от функции для перемножения многочленов. Единственная особенность — что после выполнения умножения чисел как многочлены их следует нормализовать, т.е. выполнить все переносы разрядов:

	int carry = 0;
	for (size_t i=0; i<n; ++i) {
		res[i] += carry;
		carry = res[i] / 10;
		res[i] %= 10;
	}

(Поскольку длина произведения двух чисел никогда не превзойдёт суммарной длины чисел, то размера вектора $\rm res$ хватит, чтобы выполнить все переносы.)

Улучшенная реализация: вычисления "на месте" без дополнительной памяти

Для увеличения эффективности откажемся от рекурсии в явном виде. В приведённой выше рекурсивной реализации мы явно разделяли вектор $\rm a$ на два вектора — элементы на чётных позициях отнесли к одному временно созданному вектору, а на нечётных — к другому. Однако, если бы мы переупорядочили элементы определённым образом, то необходимость в создании временных векторов тогда бы отпала (т.е. все вычисления мы могли бы производить "на месте", прямо в самом векторе $a$ ).

Заметим, что на первом уровне рекурсии элементы, младшие (первые) биты позиций которых равны нулю, относятся к вектору $a_0$ , а младшие биты позиций которых равны единице — к вектору $a_1$ . На втором уровне рекурсии выполняется то же самое, но уже для вторых битов, и т.д. Поэтому если мы в позиции $i$ каждого элемента $a[i]$ инвертируем порядок битов, и переупорядочим элементы массива $a$ в соответствии с новыми индексами, то мы и получим искомый порядок (он называется поразрядно обратной перестановкой (bit-reversal permutation)).

Например, для $n=8$ этот порядок имеет вид:

$a = \biggl\{ \Bigl[ (a_0,a_4), (a_2, a_6) \Bigr] [...]$

Действительно, на первом уровне рекурсии (окружено фигурными скобками) обычного рекурсивного алгоритма происходит разделение вектора на две части: $[a_0,a_2,a_4,a_6]$ и $[a_1,a_3,a_5,a_7]$ . Как мы видим, в поразрядно обратной перестановке этому соответствует просто разделение вектора на две половинки: первые $n/2$ элементов, и последние $n/2$ элементов. Затем происходит рекурсивный вызов от каждой половинки; пусть результирующее ДПФ от каждой из них было возвращено на месте самих элементов (т.е. в первой и второй половинах вектора $a$ соответственно):

$a = \biggl\{ \Bigl[ y_0^0,\ y_1^0,\ y_2^0,\ y_3^0[...]$

Теперь нам надо выполнить объединение двух ДПФ в одно для всего вектора. Но элементы встали так удачно, что и объединение можно выполнить прямо в этом массиве. Действительно, возьмём элементы $y_0^0$ и $y_0^1$ , применим к ним преобразование бабочки, и результат поставим на их месте — и это место и окажется тем самым, которое и должно было получиться:

$a = \biggl\{ \Bigl[ y_0^0+w_n^0y_0^1,\ y_1^0,\ y_[...]$

Аналогично, применяем преобразование бабочки к $y_1^0$ и $y_1^1$ и результат ставим на их место, и т.д. В итоге получаем:

$a = \biggl\{ \Bigl[ y_0^0+w_n^0y_0^1,\ y_1^0+w_n^[...]$
$~~~~~~~~ \Bigl[ y_0^0-w_n^0y_0^1,\ y_1^0-w_n^1y_1[...]$

Т.е. мы получили именно искомое ДПФ от вектора $a$ .

Мы описали процесс вычисления ДПФ на первом уровне рекурсии, но понятно, что те же самые рассуждения верны и для всех остальных уровней рекурсии. Таким образом, после применения поразрядно обратной перестановки вычислять ДПФ можно на месте, без привлечения дополнительных массивов.

Но теперь можно избавиться и от рекурсии в явном виде. Итак, мы применили поразрядно обратную перестановку элементов. Теперь выполним всю работу, выполняемую нижним уровнем рекурсии, т.е. вектор $a$ разделим на пары элементов, для каждого применим преобразование бабочки, в результате в векторе $a$ будут находиться результаты работы нижнего уровня рекурсии. На следующем шаге разделим вектор $a$ на четвёрки элементов, к каждой применим преобразование бабочки, в результате получим ДПФ для каждой четвёрки. И так далее, наконец, на последнем шаге мы, получив результаты ДПФ для двух половинок вектора $a$ , применим к ним преобразование бабочки и получим ДПФ для всего вектора $a$ .

Итак, реализация:

typedef complex<double> base;
 
int rev (int num, int lg_n) {
	int res = 0;
	for (int i=0; i<lg_n; ++i)
		if (num & (1<<i))
			res |= 1<<(lg_n-1-i);
	return res;
}
 
void fft (vector<base> & a, bool invert) {
	int n = (int) a.size();
	int lg_n = 0;
	while ((1 << lg_n) < n)  ++lg_n;
 
	for (int i=0; i<n; ++i)
		if (i < rev(i,lg_n))
			swap (a[i], a[rev(i,lg_n)]);
 
	for (int len=2; len<=n; len<<=1) {
		double ang = 2*PI/len * (invert ? -1 : 1);
		base wlen (cos(ang), sin(ang));
		for (int i=0; i<n; i+=len) {
			base w (1);
			for (int j=0; j<len/2; ++j) {
				base u = a[i+j],  v = a[i+j+len/2] * w;
				a[i+j] = u + v;
				a[i+j+len/2] = u - v;
				w *= wlen;
			}
		}
	}
	if (invert)
		for (int i=0; i<n; ++i)
			a[i] /= n;
}

Вначале к вектору $a$ применяется поразрядно обратная перестановка, для чего вычисляется количество значащих бит ( $\rm lg\_n$ ) в числе $n$ , и для каждой позиции $i$ находится соответствующая ей позиция, битовая запись которой есть битовая запись числа $i$ , записанная в обратном порядке. Если получившаяся в результате позиция оказалась больше $i$ , то элементы в этих двух позициях надо обменять (если не это условие, то каждая пара обменяется дважды, и в итоге ничего не произойдёт).

Затем выполняется $\lg n - 1$ стадий алгоритма, на $k$ -ой из которых ( $k=2 \ldots \lg n$ ) вычисляются ДПФ для блоков длины $2^k$ . Для всех этих блоков будет одно и то же значение первообразного корня $w_{2^k}$ , которое и запоминается в переменной $\rm wlen$ . Цикл по $i$ итерируется по блокам, а вложенный в него цикл по $j$ применяет преобразование бабочки ко всем элементам блока.

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

Например, пусть $j$ — уже подсчитанное число, равное обратной перестановке битов числа $i$ . Тогда, при переходе к следующему числу $i+1$ мы должны и к числу $j$ прибавить единицу, но прибавить её в такой "инвертированной" системе счисления. В обычной двоичной системе счисления прибавить единицу — значит удалить все единицы, стоящие на конце числа (т.е. группу младших единиц), а перед ними поставить единицу. Соответственно, в "инвертированной" системе мы должны идти по битам числа, начиная со старших, и пока там стоят единицы, удалять их и переходить к следующему биту; когда же встретится первый нулевой бит, поставить в него единицу и остановиться.

Итак, получаем такую реализацию:

typedef complex<double> base;
 
void fft (vector<base> & a, bool invert) {
	int n = (int) a.size();
 
	for (int i=1, j=0; i<n; ++i) {
		int bit = n >> 1;
		for (; j>=bit; bit>>=1)
			j -= bit;
		j += bit;
		if (i < j)
			swap (a[i], a[j]);
	}
 
	for (int len=2; len<=n; len<<=1) {
		double ang = 2*PI/len * (invert ? -1 : 1);
		base wlen (cos(ang), sin(ang));
		for (int i=0; i<n; i+=len) {
			base w (1);
			for (int j=0; j<len/2; ++j) {
				base u = a[i+j],  v = a[i+j+len/2] * w;
				a[i+j] = u + v;
				a[i+j+len/2] = u - v;
				w *= wlen;
			}
		}
	}
	if (invert)
		for (int i=0; i<n; ++i)
			a[i] /= n;
}

Дополнительные оптимизации

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

Предпосчитать реверс битов для всех чисел в некоторой глобальной таблице. Особенно легко это, когда размер $n$ при всех вызовах одинаков.
Эта оптимизация становится заметной при большом количестве вызовов $fft()$ . Впрочем, эффект от неё можно заметить даже при трёх вызовах (три вызова — наиболее распространённая ситуация, т.е. когда требуется один раз перемножить два многочлена).
Отказаться от использования $\rm vector$ (перейти на обычные массивы).
Эффект от этого зависит от конкретного компилятора, однако обычно он присутствует и составляет примерно 10%-20%.
Предпосчитать все степени числа . В самом деле, в этом цикле алгоритма раз за разом производится проход по всем степеням числа от до :
```
		for (int i=0; i<n; i+=len) {
			base w (1);
			for (int j=0; j<len/2; ++j) {
				[...]
				w *= wlen;
			}
		}
```
Соответственно, перед этим циклом мы можем предпосчитать в некотором массиве все требуемые степени, и избавиться тем самым от лишних умножений во вложенном цикле.
Ориентировочное ускорение — 5-10%.
Избавиться от обращений к массивам по индексам, использовать вместо этого указатели на текущие элементы массивов, продвигая их на 1 вправо на каждой итерации.
На первый взгляд, оптимизирующие компиляторы должны быть способны самостоятельно справиться с этим, однако на практике оказывается, что замена обращений к массивам $a[i+j]$ и $a[i+j+len/2]$ на указатели ускоряет программу в распространённых компиляторах. Выигрыш составляет 5-10%.
Отказаться от стандартного типа комплексных чисел $\rm complex$ , переписав его на собственную реализацию.
Опять же, это может показаться удивительным, но даже в современных компиляторах выигрыш от такого переписывания может составлять до нескольких десятков процентов! Это косвенно подтверждает расхожее утверждение, что компиляторы хуже справляются с шаблонными типами данных, оптимизируя работу с ними гораздо хуже, чем с не-шаблонными типами.
Другой полезной оптимизацией является отсечение по длине: когда длина рабочего блока становится маленькой (скажем, 4), вычислять ДПФ для него "вручную". Если расписать эти случаи в виде явных формул при длине, равной $4/2$ , то значения синусов-косинусов примут целочисленные значения, за счёт чего можно получить прирост скорости ещё на несколько десятков процентов.

Приведём здесь реализацию с описанными улучшениями (за исключением двух последних пунктов, которые приводят к чрезмерному разрастанию кода):

int rev[MAXN];
base wlen_pw[MAXN];
 
void fft (base a[], int n, bool invert) {
	for (int i=0; i<n; ++i)
		if (i < rev[i])
			swap (a[i], a[rev[i]]);
 
	for (int len=2; len<=n; len<<=1) {
		double ang = 2*PI/len * (invert?-1:+1);
		int len2 = len>>1;
 
		base wlen (cos(ang), sin(ang));
		wlen_pw[0] = base (1, 0);
		for (int i=1; i<len2; ++i)
			wlen_pw[i] = wlen_pw[i-1] * wlen;
 
		for (int i=0; i<n; i+=len) {
			base t,
				*pu = a+i,
				*pv = a+i+len2, 
				*pu_end = a+i+len2,
				*pw = wlen_pw;
			for (; pu!=pu_end; ++pu, ++pv, ++pw) {
				t = *pv * *pw;
				*pv = *pu - t;
				*pu += t;
			}
		}
	}
 
	if (invert)
		for (int i=0; i<n; ++i)
			a[i] /= n;
}
 
void calc_rev (int n, int log_n) {
	for (int i=0; i<n; ++i) {
		rev[i] = 0;
		for (int j=0; j<log_n; ++j)
			if (i & (1<<j))
				rev[i] |= 1<<(log_n-1-j);
	}
}

На распространённых компиляторах данная реализация быстрее предыдущего "улучшенного" варианта в 2-3 раза.

Дискретное преобразование Фурье в модульной арифметике

В основе дискретного преобразования Фурье лежат комплексные числа, корни $n$ -ой степени из единицы. Для эффективного его вычисления использовались такие особенности корней, как существование $n$ различных корней, образующих группу (т.е. степень одного корня — всегда другой корень; среди них есть один элемент — генератор группы, называемый примитивным корнем).

Но то же самое верно и в отношении корней $n$ -ой степени из единицы в модульной арифметике. Точнее, не для любого модуля $p$ найдётся $n$ различных корней из единицы, однако такие модули всё же существуют. По-прежнему нам важно найти среди них примитивный корень, т.е.:

$(w_n)^n = 1 \pmod p,$
$(w_n)^k \ne 1 {\pmod p}, ~~~~~1 \le k < n.$

Все остальные $n-1$ корней $n$ -ой степени из единицы по модулю $p$ можно получить как степени примитивного корня $w_n$ (как и в комплексном случае).

Для применения в алгоритме Быстрого преобразования Фурье нам было нужно, чтобы примивный корень существовал для некоторого $n$ , являвшегося степенью двойки, а также всех меньших степеней. И если в комплексном случае примитивный корень существовал для любого $n$ , то в случае модульной арифметики это, вообще говоря, не так. Однако, заметим, что если $n = 2^k$ , т.е. $k$ -ая степень двойки, то по модулю $m = 2^{k-1}$ имеем:

$(w_n^2)^m = (w_n)^n = 1 \pmod p,$
$(w_n^2)^k = w_n^{2k} \ne 1 {\pmod p}, ~~~~~1 \le [...]$

Таким образом, если $w_n$ — примитивный корень $n=2^k$ -ой степени из единицы, то $w_n^2$ — примитивный корень $2^{k-1}$ -ой степени из единицы. Следовательно, для всех степеней двойки, меньших $n$ , примитивные корни нужной степени также существуют, и могут быть вычислены как соответствующие степени $w_n$ .

Последний штрих — для обратного ДПФ мы использовали вместо $w_n$ обратный ему элемент: $w_n^{-1}$ . Но по простому модулю $p$ обратный элемент также всегда найдётся.

Таким образом, все нужные нам свойства соблюдаются и в случае модульной арифметики, при условии, что мы выбрали некоторый достаточно большой модуль $p$ и нашли в нём примитивный корень $n$ -ой степени из единицы.

Например, можно взять такие значения: модуль $p = 7340033$ , $w_{2^{20}} = 5$ . Если этого модуля будет недостаточно, для нахождения другой пары можно воспользоваться фактом, что для модулей вида $c 2^k + 1$ (но по-прежнему обязательно простых) всегда найдётся примитивный корень степени $2^k$ из единицы.

const int mod = 7340033;
const int root = 5;
const int root_1 = 4404020;
const int root_pw = 1<<20;
 
void fft (vector<int> & a, bool invert) {
	int n = (int) a.size();
 
	for (int i=1, j=0; i<n; ++i) {
		int bit = n >> 1;
		for (; j>=bit; bit>>=1)
			j -= bit;
		j += bit;
		if (i < j)
			swap (a[i], a[j]);
	}
 
	for (int len=2; len<=n; len<<=1) {
		int wlen = invert ? root_1 : root;
		for (int i=len; i<root_pw; i<<=1)
			wlen = int (wlen * 1ll * wlen % mod);
		for (int i=0; i<n; i+=len) {
			int w = 1;
			for (int j=0; j<len/2; ++j) {
				int u = a[i+j],  v = int (a[i+j+len/2] * 1ll * w % mod);
				a[i+j] = u+v < mod ? u+v : u+v-mod;
				a[i+j+len/2] = u-v >= 0 ? u-v : u-v+mod;
				w = int (w * 1ll * wlen % mod);
			}
		}
	}
	if (invert) {
		int nrev = reverse (n, mod);
		for (int i=0; i<n; ++i)
			a[i] = int (a[i] * 1ll * nrev % mod);
	}
}

Здесь функция $\rm reverse$ находит обратный к $n$ элемент по модулю $\rm mod$ (см. Обратный элемент в поле по модулю). Константы $\rm mod$ , $\rm root$ $\rm root\_pw$ определяют модуль и примитивный корень, а $\rm root\_1$ — обратный к $\rm root$ элемент по модулю $\rm mod$ .

Как показывает практика, реализация целочисленного ДПФ работает даже медленней реализации с комплексными числами (из-за огромного количества операций взятия по модулю), однако она имеет такие преимущества, как меньшее использование памяти и отсутствие погрешностей округления.

Некоторые применения

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

Всевозможные суммы

Задача: даны два массива $a[]$ и $b[]$ . Требуется найти всевозможные числа вида $a[i]+b[j]$ , и для каждого такого числа вывести количество способов получить его.

Например, для $a = (1,2,3)$ и $b = (2,4)$ получаем: число 3 можно получить 1 способом, 4 — также одним, 5 — 2, 6 — 1, 7 — 1.

Построим по массивам $a$ и $b$ два многочлена $A$ и $B$ . В качестве степеней в многочлене будут выступать сами числа, т.е. значения $a[i]$ ( $b[i]$ ), а в качестве коэффициентов при них — сколько раз это число встречается в массиве $a$ ( $b$ ).

Тогда, перемножив эти два многочлена за $O(n \log n)$ , мы получим многочлен $C$ , где в качестве степеней будут всевозможные числа вида $a[i]+b[j]$ , а коэффициенты при них будут как раз искомыми количествами

Всевозможные скалярные произведения

Даны два массива $a[]$ и $b[]$ одной длины $n$ . Требуется вывести значения каждого скалярного произведения вектора $a$ на очередной циклический сдвиг вектора $b$ .

Инвертируем массив $a$ и припишем к нему в конец $n$ нулей, а к массиву $b$ — просто припишем самого себя. Затем перемножим их как многочлены. Теперь рассмотрим коэффициенты произведения $c[n \ldots 2n-1]$ (как всегда, все индексы в 0-индексации). Имеем:

$c[k] = \sum_{i+j=k} a[i] b[j].$

Поскольку все элементы $a[i]=0,\ i=n \ldots 2n-1$ , то мы получаем:

$c[k] = \sum_{i=0}^{n-1} a[i] b[k-i].$

Нетрудно увидеть в этой сумме, что это именно скалярное произведение вектора $a$ на $k-n-1$ -ый циклический сдвиг. Таким образом, эти коэффициенты (начиная с $n-1$ -го и закачивая $2n-2$ -ым) — и есть ответ на задачу.

Решение получилось с асимптотикой $O (n \log n)$ .

Две полоски

Даны две полоски, заданные как два булевских (т.е. числовых со значениями 0 или 1) массива $a[]$ и $b[]$ . Требуется найти все такие позиции на первой полоске, что если приложить, начиная с этой позиции, вторую полоску, ни в каком месте не получится $\rm true$ сразу на обеих полосках. Эту задачу можно переформулировать таким образом: дана карта полоски, в виде 0/1 — можно вставать в эту клетку или нет, и дана некоторая фигурка в виде шаблона (в виде массива, в котором 0 — нет клетки, 1 — есть), требуется найти все позиции в полоске, к которым можно приложить фигурку.

Эта задача фактически ничем не отличается от предыдущей задачи — задачи о скалярном произведении. Действительно, скалярное произведение двух 0/1 массивов — это количество элементов, в которых одновременно оказались единицы. Наша задача в том, чтобы найти все циклические сдвиги второй полоски так, чтобы не нашлось ни одного элемента, в котором бы в обеих полосках оказались единицы. Т.е. мы должны найти все циклические сдвиги второго массива, при которых скалярное произведение равно нулю.

Таким образом, эту задачу мы решили за $O(n \log n)$ .