Немного теории. Дискретное преобразование Фурье - это линейное преобразование, которое
задается матрицей Вандермонда. Матрица состоит из степеней примитивного корня из единицы
степени 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
#include
// Возведение 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 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 forward_dft (std::vector& a)
{
std::vector 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 inversed_dft (std::vector& a)
{
std::vector 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 convolution (std::vector a, std::vector 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 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 multiply (std::vector a, std::vector b)
{
std::vector 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 u(std::begin(a), std::end(a));
std::vector v(std::begin(b), std::end(b));
std::vector result = multiply(u, v);
for(int i = 0; i < result.size(); i++)
std::cout << result[i] << " ";
/*std::vector u(std::begin(a), std::end(a));
std::vector v(std::begin(b), std::end(b));
std::vector result = convolution(u, v);
for(int i = 0; i < result.size(); i++)
std::cout << result[i] << " ";*/
/*int a[] = {4, 3, 2, 1};
std::vector v(std::begin(a), std::end(a));
std::vector 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), но возникает лишний ноль из-за размерности векторов и дополнения нулями. Где
я ошибся?
Ответы
Ответ 1
У вас в программе вообще отсутствует поиск простого числа, в поле которого дальше
идёт дискретное преобразование Фурье - это довольно странно. Я предположил, что вы
считаете, что размер массива плюс единица является простым числом (что в вашем примере
действительно так - 2*8+1=17 - простое).
Далее у вас неконсистентно используется функция generator - объявлена она с параметром-простым
числом p, а вызывается от a.size() (что, по моему предположению, является p-1). При
этом, что забавно, первообразный корень тоже находится - он существует по в том числе
по модулям, являющимися степенями двойки.
Последняя и фатальная проблема заключается в том, что модуля 17 недостаточно для
умножения чисел длины 8 по основанию 10. При умножении чисел мы переходим к умножению
многочленов (например, 4321 переходит в многочлен 4x^3+3*x^2+2*x+1), потом умножаем
многочлены в поле, а потом превращаем многочлен обратно в число. Если в процессе умножения
хотя бы один коэффициент многочлена "переполнился" (т.е. стал больше или равен модуля),
то мы всё получим корректный результат умножения многочленов, но он уже не будет соответствовать
корректному числу.
Например, умножая 4321 на 321 (т.е. (4x^3+3*x^2+2*x+1)*(3*x^2+2*x+1)) мы получим
многочлен 12*x^5+17*x^4+16*x^3+10*x^2+4*x+1, что при подстановке x=10 даёт как раз
число 1387041. Однако если известны коэффициенты многочлена лишь по модулю 17, то коэффициент
при x^4 зануляется и становится невозможно восстановить, чему он был равен на самом
деле. Вам следует выбирать поле достаточно большого размера, чтобы при перемножении
двух многочленов коэффициенты результата не превосходили модуль.
