Опечатки, баги и конструктивная критика принимаются (Антон Тимофеев).
На лекции забыли ввести умножение в алгебраической форме: (a + ib)(с + id) = (ac - bd, ad + bc)
typedef double dbl;
struct Complex
{
dbl a, b; /* Действительная и мнимая части числа */
explicit Complex( dbl a = 0, dbl b = 0 ) : a(a), b(b) {}
Complex operator+( Complex c ) const { return Complex(a + c.a, b + c.b); }
Complex operator-( Complex c ) const { return Complex(a - c.a, b - c.b); }
Complex operator*( Complex c ) const { return Complex(a * c.a - b * c.b, a * c.b + b * c.a); }
void operator+=( Complex c ) { a += c.a; b += c.b; }
void operator*=( Complex c ) { dbl t = a; a = a * c.a - b * c.b; b = t * c.b + b * c.a; }
void operator*=( dbl m ) { a *= m, b *= m; }
void operator=( Complex const& c ) { a = c.a, b = c.b; }
};
Use cases:
dbl ang = 2 * pi / n;
Complex w1(cos(ang), sin(ang)), z;
z = Complex(30);
...
Complex arr[maxn];
vector<Complex> vec;
fill(arr, arr + maxn, Complex()); // заполнение массива
fill(vec.begin(), vec.end(), Complex()); // заполнение вектора
dbl va[] = {1, 2, 3, 4, 5}; // инициализация
vector<Complex> a(va, va + 5); // vector<Complex> от vector<dbl>
Input: a = {1, 2, 3, 4, 5}
Next pow of 2: n = 8
Complex array: a = {(1, 0); (2, 0); (3, 0); (4, 0); (5, 0); (0, 0); (0, 0); (0, 0)}
fft(a, false):
0. reverse bits reordering:
a = {(1.0, 0.0); (5.0, 0.0); (3.0, 0.0); (0.0, 0.0); (2.0, 0.0); (0.0, 0.0); (4.0, 0.0); (0.0, 0.0)}
1. merged intervals with length 1:
a = {(6.00, 0.00);(-4.00, 0.00);(3.00, 0.00);(3.00, 0.00);(2.00, 0.00);(2.00, 0.00);(4.00, 0.00);(4.00, 0.00)}
2. merged intervals with length 2:
a = {(9.00, 0.00);(-4.00, 3.00);(3.00, 0.00);(-4.00, -3.00);(6.00, 0.00);(2.00, 4.00);(-2.00, 0.00);(2.00, -4.00)}
3. merged intervals with length 4:
a = {(15.00, 0.00);(-5.41, 7.24);(3.00, -2.00);(-2.59, 1.24);(3.00, 0.00);(-2.59, -1.24);(3.00, 2.00);(-5.41, -7.24)}
fft(a, true):
0. reverse bits reordering:
a = {(15.00, 0.00);(3.00, 0.00);(3.00, -2.00);(3.00, 2.00);(-5.41, 7.24);(-2.59, -1.24);(-2.59, 1.24);(-5.41, -7.24)}
1. merged intervals with length 1:
a = {(18.00, 0.00);(12.00, 0.00);(6.00, 0.00);(0.00, -4.00);(-8.00, 6.00);(-2.83, 8.49);(-8.00, -6.00);(2.83, 8.49)}
2. merged intervals with length 2:
a = {(24.00, 0.00);(8.00, 0.00);(12.00, 0.00);(16.00, 0.00);(-16.00, 0.00);(5.66, 5.66);(0.00, 12.00);(-11.31, 11.31)}
3. merged intervals with length 4:
a = {(8.00, 0.00);(16.00, 0.00);(24.00, 0.00);(32.00, 0.00);(40.00, 0.00);(0.00, 0.00);(0.00, 0.00);(0.00, 0.00)}
4. Divide all on n = 8 (because of inverse transform):
a = {(1.00, 0.00);(2.00, 0.00);(3.00, 0.00);(4.00, 0.00);(5.00, 0.00);(0.00, 0.00);(0.00, 0.00);(0.00, 0.00)}
a[j], b[j] — многочлены. На этот раз нужно посчитать свёртку c[k] = ∑ajbk-j взятую по простому модулю p.
Для любого простого числа p существует примитивный (первообразный) корень g, т.е. такое число, что
g0, g1, g2, ..., gp-2 по модулю p — различные числа из множества {1, 2, ..., p-1}.
По малой теореме Ферма gp-1 = 1 = g0 (mod p), значит все остальные степени будут циклически повторять первые p-1.
Недоказанное утверждение состоит в том, что среди множества первообразных корней для конкретного p существуют маленькие первообразные корни, т.е. можно перебрать числа, начиная с 2, проверя свойство примитивности с помощью следующей теоремы:
Th: g — примитивный корень тогда и только тогда, когда для любого несобственного делителя q числа p-1 верно, что g(p-1)/q > 1. (Иначе, цикл получится длины (p-1)/q вместо p-1).
Если число (p-1) делится на n = 2k, то числа вида ω(n, k) = gk(p-1)/n удовлетворяют всем свойствам комплексных корней из 1, требуемых для FFT (свойства 1-3 с лекции). Осталось только находить число обратное числу nпо модулю p (чтобы делить в обратном преобразовании).
Таким образом, для простого p = t2k + 1 можно считать свёртку от векторов с длиной, кратной n = 2k.
n = 220 = 1048576. Тогда p = 7*220 + 1 = 7340033 — простое.
Один из примитивных кореней g = 3. Числа для FFT, в которых будут вычислены многолены по модулю p имеют вид
ω(n, k) = 37k*220 / n (mod p).