Немного теории. Дискретное преобразование Фурье - это линейное преобразование, которое задается матрицей Вандермонда. Матрица состоит из степеней примитивного корня из единицы степени n
, где n
- длина векторов и размерность квадратной матрицы. Проще показать, как она строится, чем описывать словами.
Матрица обратного преобразования порождается степенями элемента, обратного к выбранному примитивному корню из единицы.
Преобразование Фурье над конечным числовым полем почему-то непопулярно у кодеров в Роиссе, но оно есть и используется в криптографии. Разница в том, что корень из единицы извлекается не в комплексном поле, а в поле классов вычетов. Остальное все то же самое.
На свёртку можно смотреть как на операцию умножения многочленов с коэффициентами из поля или конечного кольца классов вычетов. Если применить нормализацию (перекидывать перенос при переполнении разрядов), то получается умножение длинных чисел. Как я понял, свёртка и преобразование Фурье достаточно умны, чтобы при умножении чисел расширять конечное поле, если в нем не хватает цифр. Тут надо сделать важное пояснение примером.
Хотим найти преобразование Фурье от вектора (4, 3, 2, 1)
над числовым полем Галуа GF(5)
характеристики 5. Матрица Вандермонда порождается примитивным корнем этого поля, т.е. двойкой.
Замечание: не всякий элемент, взаимно простой с модулем, по которому строится поле, является примитивным корнем. Русская педивикия врет. Например, степени двойки не покроют все элементы поля GF(7)
.
В результате действия этой матрицей на вектор (4, 3, 2, 1)
получается вектор (0, 1, 2, 3)
. Но что, если мы захотим найти преобразование Фурье от вектора (7, 1, 7, 1)
? В поле GF(5)
нет элемента 7, но есть класс эквивалентности [2]
, содержащий семерку. Как в этом случае отработает преобразование?
Аналогичная ситуация возникнет, если будем умножать числа при помощи этих преобразований. Например, умножить число 7777 (вектор (7, 7, 7, 7)
) на себя. Здесь снова используется поле GF(5)
, в котором опять нет семерки.
Свёртка работает так: дополняем два данных вектора нулями (их длина при этом увеличивается в два раза), применяем к ним прямое преобразование Фурье, два полученных вектора умножаем почленно и к результату применяем обратное преобразование Фурье. Эта операция работает как умножение двух многочленов с заданными коэффициентами из поля.
Пример. Найдем свертку (2, 1)
и (1, 1)
. Дополним векторы нулями, получим a' = (2, 1, 0, 0)
и b' = (1, 1, 0, 0)
. Применим к ним прямое преобразование: F(a') = (3, 4, 1, 0)
и F(b') = (2, 3, 0, 4)
. Умножим последние векторы почленно: (1, 2, 0, 0)
. Применим к полученному вектору обратное преобразование: (2, 3, 1, 0)
.
Моя реализация преобразований Фурье, свёртки и умножения длинных чисел:
#include <iostream>
#include <vector>
// Возведение a в степень b по модулю p
int powmod (int a, int b, int p)
{
int res = 1;
while (b)
if (b & 1)
res = int ((long long) res * a % p), --b;
else
a = int ((long long) a * a % p), b >>= 1;
return res;
}
// Вычисление обратного к a элемента по модулю n
int inverse(int a, int n)
{
int b0 = n, t, q;
int x0 = 0, x1 = 1;
if (n == 1) return 1;
while (a > 1) {
q = a / n;
t = n, n = a % n, a = t;
t = x0, x0 = x1 - q * x0, x1 = t;
}
if (x1 < 0) x1 += b0;
return x1;
}
// Вычисление примитивного корня числового поля
// (Стандартный DFT использует комплексные числа, поэтому изменим алгоритм для работы в Zp)
// Ограничение: p - простое
int generator (int p)
{
std::vector<int> fact;
int phi = p - 1, n = phi;
for (int i = 2; i * i <= n; ++i)
if (n % i == 0)
{
fact.push_back (i);
while (n % i == 0)
n /= i;
}
if (n > 1)
fact.push_back (n);
for (int res = 2; res <= p; ++res) {
bool ok = true;
for (size_t i = 0; i < fact.size() && ok; ++i)
ok &= powmod (res, phi / fact[i], p) != 1;
if (ok) return res;
}
return -1;
}
// Реализация прямого DFT над Zp
// Принимает вектор длины n, эта длина используется при вычислении
// примитивного корня поля
// Вычисления приводятся по mod (n + 1), т.к. длина вектора на 1 меньше p,
// где p: Zp
std::vector<int> forward_dft (std::vector<int>& a)
{
std::vector<int> result;
int prim_root = generator (a.size()); // Корень степени n из единицы
for (int i = 0; i < a.size(); i++)
{
int sum = 0;
int power = 0;
for (int j = 0; j < a.size(); j++)
{
int power_of_root = powmod (prim_root, power, a.size() + 1);
sum += power_of_root * a[j]; // Накапливаем сумму произведений элементов строки матрицы на вектор
power += i; // Увеличиваем степень, в которую возводится примитивный корень из единицы
}
result.push_back(sum % (a.size() + 1)); // Набрали сумму, приводим по mod p, p = длина вектора + 1
}
return result;
}
std::vector<int> inversed_dft (std::vector<int>& a)
{
std::vector<int> result;
int prim_root = generator (a.size()); // Корень степени n из единицы
int inv_prim_root = inverse (prim_root, a.size() + 1); // Обратный к найденному примитивному корню
int inv_n = inverse (a.size(), a.size() + 1); // Обратный к n по mod (n + 1)
for (int i = 0; i < a.size(); i++)
{
int sum = 0;
int power = 0;
for (int j = 0; j < a.size(); j++)
{
int power_of_inv_root = powmod (inv_prim_root, power, a.size() + 1);
sum += inv_n * power_of_inv_root * a[j]; // Накапливаем сумму произведений элементов строки матрицы на вектор
power += i; // Увеличиваем степень, в которую возводится примитивный корень из единицы
}
result.push_back(sum % (a.size() + 1)); // Набрали сумму, приводим по mod p, p = длина вектора + 1
}
return result;
}
// Свертка двух векторов при помощи дискретного преобразования Фурье
std::vector<int> convolution (std::vector<int> a, std::vector<int> b)
{
a.insert(a.end(), a.size(), 0);
b.insert(b.end(), b.size(), 0);
a = forward_dft(a);
b = forward_dft(b);
std::vector<int> result(a.size(), 0);
for(int i = 0; i < result.size(); i++)
result[i] = a[i] * b[i];
result = inversed_dft(result);
return result;
}
// Умножение двух больших чисел
std::vector<int> multiply (std::vector<int> a, std::vector<int> b)
{
std::vector<int> result = convolution(a, b);
// Нормализация, выполнение переносов
int carry = 0;
for (int i = 0; i < result.size(); i++)
{
result[i] += carry;
carry = result[i] / 10;
result[i] %= 10;
}
return result;
}
int main()
{
int a[] = {1, 2, 1, 2, 1, 2, 1, 2};
int b[] = {7, 5, 7, 5, 7, 5, 7, 5};
std::vector<int> u(std::begin(a), std::end(a));
std::vector<int> v(std::begin(b), std::end(b));
std::vector<int> result = multiply(u, v);
for(int i = 0; i < result.size(); i++)
std::cout << result[i] << " ";
/*std::vector<int> u(std::begin(a), std::end(a));
std::vector<int> v(std::begin(b), std::end(b));
std::vector<int> result = convolution(u, v);
for(int i = 0; i < result.size(); i++)
std::cout << result[i] << " ";*/
/*int a[] = {4, 3, 2, 1};
std::vector<int> v(std::begin(a), std::end(a));
std::vector<int> result = forward_dft(v);
std::cout << "[DEBUG]: Forward DFT" << std::endl;
for(int i = 0; i < result.size(); i++)
std::cout << result[i] << " ";
std::cout << std::endl;
result = inversed_dft(result);
std::cout << "[DEBUG]: Inversed DFT" << std::endl;
for(int i = 0; i < result.size(); i++)
std::cout << result[i] << " ";*/
}
Умножение длинных чисел в функции main
работает неправильно. Сначала я думал, что это связано с тем, что алгоритм работает над конечным полем, в котором нет некоторых "цифр" большого числа. Но видно, что при умножении чисел 12121212 на 75757575 используется поле достаточного размера, чтобы содержать там соответствие для цифры 7. При возведении числа 12 в квадрат алгоритм работает почти правильно, т.е. получается вектор (1, 4, 4, 0)
, но возникает лишний ноль из-за размерности векторов и дополнения нулями. Где я ошибся?