Multiply Large Numbers with the Discrete Fourier Transform

A bit of theory. The discrete Fourier transform is a linear transform defined by a Vandermonde matrix. The matrix consists of degrees of a primitive root of a unit of degree n , where n is the length of the vectors and the dimension of the square matrix. It's easier to show how it is built, than to describe in words.

The inverse transformation matrix is ​​generated by the powers of the element inverse to the selected primitive root of one.

Fourier transform over a finite number field is somehow unpopular with Roisse coders, but it is and is used in cryptography. The difference is that the root of the unit is not extracted in the complex field, but in the field of residue classes. The rest is all the same.

Convolution can be viewed as the multiplication of polynomials with coefficients from a field or a finite ring of classes of residues. If we apply the normalization (transfer the transfer when the digits overflow), then we get the multiplication of long numbers. As I understand it, convolutions and Fourier transforms are smart enough to expand the final field when multiplying numbers if there are not enough digits in it. Here it is necessary to make an important explanation of the example.

We want to find the Fourier transform of the vector (4, 3, 2, 1) over the Galois number field GF(5) characteristic 5. The Vandermonde matrix is ​​generated by the primitive root of this field, i.e. two.

Note: not every element that is mutually simple with the module on which the field is built is a primitive root. Russian pedivikiya is lying. For example, powers of two will not cover all elements of the GF(7) field.

As a result of the action of this matrix on the vector (4, 3, 2, 1) , a vector is obtained (0, 1, 2, 3) . But what if we want to find the Fourier transform of the vector (7, 1, 7, 1) ? In the field GF(5) there is no element 7, but there is an equivalence class [2] containing a seven. How does the conversion work in this case?

A similar situation will arise if we multiply the numbers with the help of these transformations. For example, multiply the number 7777 (vector (7, 7, 7, 7) ) by itself. Here again, the GF(5) field is used, in which again there is no seven.

The convolution works like this: we supplement the two vector data with zeros ( their length doubles ), apply the Fourier transform to them, multiply the two vectors obtained by term, and apply the inverse Fourier transform to the result. This operation works as a multiplication of two polynomials with given coefficients from the field.

Example. Find the convolution (2, 1) and (1, 1) . Let's add vectors with zeros, we get a' = (2, 1, 0, 0) and b' = (1, 1, 0, 0) . Apply the direct transformation to them: F(a') = (3, 4, 1, 0) and F(b') = (2, 3, 0, 4) . Multiply the last vectors by term: (1, 2, 0, 0) . Apply to the resulting vector the inverse transform: (2, 3, 1, 0) .

My implementation of Fourier transforms, convolutions and multiplications of long numbers:

 #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] << " ";*/ } 

Multiplying long numbers in the main function does not work correctly. At first I thought that this was due to the fact that the algorithm was working on a finite field, in which there were not some “digits” of a large number. But it can be seen that when multiplying the numbers 12121212 by 75757575, a field of sufficient size is used to contain the correspondence for the digit 7. There, when the number 12 is squared, the algorithm works almost correctly, i.e. a vector is obtained (1, 4, 4, 0) , but an extra zero occurs due to the dimension of the vectors and the padding with zeros. Where am I wrong?