добавлено: 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%.

  • Предпосчитать все степени числа wlen. В самом деле, в этом цикле алгоритма раз за разом производится проход по всем степеням числа wlen от 0 до len/2-1:

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

 

 


powered by e-maxx.ru