Агуулгыг алгасах

Хурдан Фурье хувиргалт

Энэ өгүүлэлд бид $n$ урттай хоёр олон гишүүнтийг $O(n \log n)$ хугацаанд үржүүлэх боломж олгодог алгоритмыг авч үзэх бөгөөд энэ нь $O(n^2)$ хугацаа шаарддаг тривиал үржүүлэлтээс сайн юм. Мэдээж урт хоёр тоог үржүүлэхийг мөн олон гишүүнт үржүүлэх бодлого болгон шилжүүлж болох тул урт хоёр тоог мөн $O(n \log n)$ хугацаанд үржүүлж болно (энд $n$ нь тоонуудын цифрийн тоо).

Хурдан Фурье хувиргалт (FFT)-ыг нээсэн гавьяаг 1965 онд алгоритмаа нийтэлсэн Кули, Тьюки нарт хамааруулдаг. Гэвч үнэндээ FFT-г өмнө нь олон удаа нээсэн боловч орчин үеийн компьютер зохион бүтээгдэхээс өмнө түүний ач холбогдлыг ойлгоогүй байсан. Зарим судлаачид FFT-г нээсэн гавьяаг 1924 онд Рунге, Кёниг нарт хамааруулдаг. Гэвч үнэндээ Гаусс ийм аргыг аль хэдийн 1805 онд боловсруулсан боловч хэзээ ч нийтлээгүй.

Энд танилцуулж буй FFT алгоритм $O(n \log n)$ хугацаанд ажилладаг боловч дурын том коэффициенттэй дурын том олон гишүүнтийг үржүүлэх, эсвэл дурын том бүхэл тоог үржүүлэхэд ажиллахгүй гэдгийг анхаараарай. Энэ нь бага коэффициенттэй $10^5$ хэмжээтэй олон гишүүнт, эсвэл $10^6$ хэмжээтэй хоёр тоог үржүүлэхийг амархан боловсруулж чадах бөгөөд энэ нь программчлалын тэмцээний бодлого бодоход ихэвчлэн хангалттай байдаг. $10^6$ биттэй тоог үржүүлэх хэмжээнээс хэтэрвэл тооцоолол хийх явцад ашиглагдах хөвөгч цэгтэй тооны хүрээ ба нарийвчлал нь эцсийн үр дүнг үнэн зөв гаргахад хангалтгүй болно, гэхдээ дурын том олон гишүүнт/бүхэл тооны үржүүлэлт хийж чадах илүү нийлмэл хувилбарууд байдаг. Жишээ нь 1971 онд Шёнхаге, Штрассен нар цагираг бүтцэд FFT-г рекурсивээр хэрэглэдэг, $O(n \log n \log \log n)$-д ажилладаг дурын том тоо үржүүлэх хувилбарыг боловсруулсан. Мөн саяхан (2019 онд) Харви, ван дер Ховен нар жинхэнэ $O(n \log n)$-д ажилладаг алгоритм нийтэлсэн.

Дискрет Фурье хувиргалт

$n - 1$ зэрэгтэй олон гишүүнт байг:

$$A(x) = a_0 x^0 + a_1 x^1 + \dots + a_{n-1} x^{n-1}$$

Ерөнхий чанарыг алдагдуулалгүйгээр бид $n$ буюу коэффициентийн тоог $2$-ын зэрэг гэж үзнэ. Хэрэв $n$ нь $2$-ын зэрэг биш бол бид зүгээр л дутуу гишүүд $a_i x^i$-г нэмж, $a_i$ коэффициентүүдийг $0$ болгоно.

Комплекс тооны онол нь $x^n = 1$ тэгшитгэл $n$ комплекс шийдтэй ($n$ дугаар нэгжийн язгуур гэж нэрлэдэг) бөгөөд шийдүүд нь $k = 0 \dots n-1$-тэйгээр $w_{n, k} = e^{\frac{2 k \pi i}{n}}$ хэлбэртэй болохыг хэлдэг. Түүнчлэн эдгээр комплекс тоо маш сонирхолтой шинж чанартай: жишээ нь үндсэн $n$ дугаар язгуур $w_n = w_{n, 1} = e^{\frac{2 \pi i}{n}}$-г бусад бүх $n$ дугаар язгуурыг илэрхийлэхэд ашиглаж болно: $w_{n, k} = (w_n)^k$.

$A(x)$ олон гишүүнтийн (эсвэл эквивалентаар коэффициентийн вектор $(a_0, a_1, \dots, a_{n-1})$-ийн) дискрет Фурье хувиргалт (DFT)-ыг $x = w_{n, k}$ цэгүүд дэх олон гишүүнтийн утгууд гэж тодорхойлно, өөрөөр хэлбэл энэ нь дараах вектор:

$$\begin{align} \text{DFT}(a_0, a_1, \dots, a_{n-1}) &= (y_0, y_1, \dots, y_{n-1}) \\ &= (A(w_{n, 0}), A(w_{n, 1}), \dots, A(w_{n, n-1})) \\ &= (A(w_n^0), A(w_n^1), \dots, A(w_n^{n-1})) \end{align}$$

Үүнтэй адилаар урвуу дискрет Фурье хувиргалт-ыг тодорхойлно: Олон гишүүнтийн утгууд $(y_0, y_1, \dots, y_{n-1})$-ийн урвуу DFT нь олон гишүүнтийн коэффициентүүд $(a_0, a_1, \dots, a_{n-1})$ юм.

$$\text{InverseDFT}(y_0, y_1, \dots, y_{n-1}) = (a_0, a_1, \dots, a_{n-1})$$

Тиймээс шууд DFT нь $n$ дугаар язгуур дахь цэгүүдэд олон гишүүнтийн утгыг тооцоолдог бол урвуу DFT нь тэдгээр утгыг ашиглан олон гишүүнтийн коэффициентүүдийг сэргээж чадна.

DFT-ийн хэрэглээ: олон гишүүнтийг хурдан үржүүлэх

$A$ ба $B$ гэсэн хоёр олон гишүүнт байг. Бид тус бүрийн DFT-г тооцоолно: $\text{DFT}(A)$ ба $\text{DFT}(B)$.

Хэрэв бид эдгээр олон гишүүнтийг үржүүлбэл юу болох вэ? Мэдээж цэг бүрт утгууд нь зүгээр л үржигдэнэ, өөрөөр хэлбэл

$$(A \cdot B)(x) = A(x) \cdot B(x).$$

Энэ нь хэрэв бид $\text{DFT}(A)$ ба $\text{DFT}(B)$ векторуудыг үржүүлбэл — нэг векторын элемент бүрийг нөгөө векторын харгалзах элементээр үржүүлэх замаар — бид $\text{DFT}(A \cdot B)$ олон гишүүнтийн DFT-ээс өөр юу ч авахгүй гэсэн үг:

$$\text{DFT}(A \cdot B) = \text{DFT}(A) \cdot \text{DFT}(B)$$

Эцэст нь урвуу DFT-г хэрэглэснээр бид дараахыг авна:

$$A \cdot B = \text{InverseDFT}(\text{DFT}(A) \cdot \text{DFT}(B))$$

Баруун талд хоёр DFT-ийн үржвэр гэдгээр бид векторын элементүүдийн хос хосоор нь үржүүлэхийг хэлж байна. Үүнийг $O(n)$ хугацаанд тооцоолж болно. Хэрэв бид DFT ба урвуу DFT-г $O(n \log n)$-д тооцоолж чадвал хоёр олон гишүүнтийн (улмаар хоёр урт тооны) үржвэрийг ижил time complexity-тэйгээр тооцоолж чадна.

Хоёр олон гишүүнт ижил зэрэгтэй байх ёстойг тэмдэглэх нь зүйтэй. Эс бөгөөс DFT-ийн хоёр үр дүнгийн вектор өөр өөр урттай болно. Үүнийг бид $0$ утгатай коэффициент нэмэх замаар шийдэж болно.

Мөн хоёр олон гишүүнтийн үржвэрийн үр дүн нь $2 (n - 1)$ зэрэгтэй олон гишүүнт байдаг тул бид олон гишүүнт бүрийн зэргийг хоёр дахин нэмэгдүүлэх ёстой (дахин $0$-ээр дүүргэх замаар). $n$ утгатай вектороос бид $2n - 1$ коэффициенттэй хүссэн олон гишүүнтээ сэргээж чадахгүй.

Хурдан Фурье хувиргалт

Хурдан Фурье хувиргалт нь DFT-г $O(n \log n)$ хугацаанд тооцоолох боломж олгодог арга юм. FFT-ийн үндсэн санаа нь хуваа ба ялагтун аргыг хэрэглэх явдал юм. Бид олон гишүүнтийн коэффициентийн векторыг хоёр вектор болгон хувааж, тус бүрийн DFT-г рекурсивээр тооцоолж, бүтэн олон гишүүнтийн DFT-г тооцоолохын тулд үр дүнг нэгтгэнэ.

Тиймээс $n - 1$ зэрэгтэй $A(x)$ олон гишүүнт байг, энд $n$ нь $2$-ын зэрэг бөгөөд $n > 1$:

$$A(x) = a_0 x^0 + a_1 x^1 + \dots + a_{n-1} x^{n-1}$$

Бид үүнийг хоёр жижиг олон гишүүнт болгон хуваана: нэг нь зөвхөн тэгш байрлалын коэффициентүүдийг, нөгөө нь сондгой байрлалын коэффициентүүдийг агуулна:

$$\begin{align} A_0(x) &= a_0 x^0 + a_2 x^1 + \dots + a_{n-2} x^{\frac{n}{2}-1} \\ A_1(x) &= a_1 x^0 + a_3 x^1 + \dots + a_{n-1} x^{\frac{n}{2}-1} \end{align}$$

Дараахыг харахад амархан

$$A(x) = A_0(x^2) + x A_1(x^2).$$

$A_0$ ба $A_1$ олон гишүүнтүүд нь $A$ олон гишүүнтийн зөвхөн хагас хэмжээний коэффициенттэй. Хэрэв бид $\text{DFT}(A_0)$ ба $\text{DFT}(A_1)$-г ашиглан $\text{DFT}(A)$-г шугаман хугацаанд тооцоолж чадвал time complexity-ийн хувьд $T_{\text{DFT}}(n) = 2 T_{\text{DFT}}\left(\frac{n}{2}\right) + O(n)$ рекуррент хамаарал гарах ба мастер теорем-оор $T_{\text{DFT}}(n) = O(n \log n)$ болно.

Үүнийг хэрхэн хийхийг сурцгаая.

Бид $\left(y_k^0\right)_{k=0}^{n/2-1} = \text{DFT}(A_0)$ ба $\left(y_k^1\right)_{k=0}^{n/2-1} = \text{DFT}(A_1)$ векторуудыг тооцоолсон гэж үзье. $\left(y_k\right)_{k=0}^{n-1} = \text{DFT}(A)$-ийн илэрхийллийг олъё.

Эхний $\frac{n}{2}$ утгын хувьд бид өмнө тэмдэглэсэн $A(x) = A_0(x^2) + x A_1(x^2)$ тэгшитгэлийг зүгээр ашиглаж болно:

$$y_k = y_k^0 + w_n^k y_k^1, \quad k = 0 \dots \frac{n}{2} - 1.$$

Гэвч хоёр дахь $\frac{n}{2}$ утгын хувьд бид арай өөр илэрхийлэл олох хэрэгтэй:

$$\begin{align} y_{k+n/2} &= A\left(w_n^{k+n/2}\right) \\ &= A_0\left(w_n^{2k+n}\right) + w_n^{k + n/2} A_1\left(w_n^{2k+n}\right) \\ &= A_0\left(w_n^{2k} w_n^n\right) + w_n^k w_n^{n/2} A_1\left(w_n^{2k} w_n^n\right) \\ &= A_0\left(w_n^{2k}\right) - w_n^k A_1\left(w_n^{2k}\right) \\ &= y_k^0 - w_n^k y_k^1 \end{align}$$

Энд бид дахин $A(x) = A_0(x^2) + x A_1(x^2)$ ба $w_n^n = 1$, $w_n^{n/2} = -1$ гэсэн хоёр адилтгалыг ашигласан.

Тиймээс бид бүтэн вектор $(y_k)$-г тооцоолох хүссэн томьёогоо авна:

$$\begin{align} y_k &= y_k^0 + w_n^k y_k^1, &\quad k = 0 \dots \frac{n}{2} - 1, \\ y_{k+n/2} &= y_k^0 - w_n^k y_k^1, &\quad k = 0 \dots \frac{n}{2} - 1. \end{align}$$

(Энэ $a + b$ ба $a - b$ хэв маягийг заримдаа эрвээхэй гэж нэрлэдэг.)

Ингээд бид DFT-г $O(n \log n)$ хугацаанд хэрхэн тооцоолохыг сурлаа.

Урвуу FFT

$(y_0, y_1, \dots y_{n-1})$ вектор буюу $x = w_n^k$ цэгүүд дэх $n - 1$ зэрэгтэй $A$ олон гишүүнтийн утгууд өгөгдсөн байг. Бид олон гишүүнтийн коэффициентүүд $(a_0, a_1, \dots, a_{n-1})$-г сэргээхийг хүсэж байна. Энэ мэдэгдэж буй бодлогыг интерполяц гэж нэрлэдэг бөгөөд үүнийг бодох ерөнхий алгоритмууд байдаг. Гэвч энэ тусгай тохиолдолд (учир нь бид нэгжийн язгуур дахь цэгүүдийн утгыг мэддэг) бид хамаагүй энгийн алгоритм (практикт шууд FFT-тэй ижил) авч болно.

Бид DFT-г тодорхойлолтынх нь дагуу матриц хэлбэрээр бичиж болно:

$$ \begin{pmatrix} w_n^0 & w_n^0 & w_n^0 & w_n^0 & \cdots & w_n^0 \\ w_n^0 & w_n^1 & w_n^2 & w_n^3 & \cdots & w_n^{n-1} \\ w_n^0 & w_n^2 & w_n^4 & w_n^6 & \cdots & w_n^{2(n-1)} \\ w_n^0 & w_n^3 & w_n^6 & w_n^9 & \cdots & w_n^{3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ w_n^0 & w_n^{n-1} & w_n^{2(n-1)} & w_n^{3(n-1)} & \cdots & w_n^{(n-1)(n-1)} \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{pmatrix} = \begin{pmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1} \end{pmatrix} $$

Энэ матрицыг Вандермондын матриц гэж нэрлэдэг.

Тиймээс бид $(y_0, y_1, \dots y_{n-1})$ векторыг зүүнээс нь матрицын урвуугаар үржүүлэх замаар $(a_0, a_1, \dots, a_{n-1})$ векторыг тооцоолж болно:

$$ \begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{pmatrix} = \begin{pmatrix} w_n^0 & w_n^0 & w_n^0 & w_n^0 & \cdots & w_n^0 \\ w_n^0 & w_n^1 & w_n^2 & w_n^3 & \cdots & w_n^{n-1} \\ w_n^0 & w_n^2 & w_n^4 & w_n^6 & \cdots & w_n^{2(n-1)} \\ w_n^0 & w_n^3 & w_n^6 & w_n^9 & \cdots & w_n^{3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ w_n^0 & w_n^{n-1} & w_n^{2(n-1)} & w_n^{3(n-1)} & \cdots & w_n^{(n-1)(n-1)} \end{pmatrix}^{-1} \begin{pmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1} \end{pmatrix} $$

Хурдан шалгахад матрицын урвуу нь дараах хэлбэртэй болохыг харж болно:

$$ \frac{1}{n} \begin{pmatrix} w_n^0 & w_n^0 & w_n^0 & w_n^0 & \cdots & w_n^0 \\ w_n^0 & w_n^{-1} & w_n^{-2} & w_n^{-3} & \cdots & w_n^{-(n-1)} \\ w_n^0 & w_n^{-2} & w_n^{-4} & w_n^{-6} & \cdots & w_n^{-2(n-1)} \\ w_n^0 & w_n^{-3} & w_n^{-6} & w_n^{-9} & \cdots & w_n^{-3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ w_n^0 & w_n^{-(n-1)} & w_n^{-2(n-1)} & w_n^{-3(n-1)} & \cdots & w_n^{-(n-1)(n-1)} \end{pmatrix} $$

Тиймээс бид дараах томьёог авна:

$$a_k = \frac{1}{n} \sum_{j=0}^{n-1} y_j w_n^{-k j}$$

Үүнийг $y_k$-ийн томьёотой харьцуулбал

$$y_k = \sum_{j=0}^{n-1} a_j w_n^{k j},$$

эдгээр бодлого бараг ижил болохыг анзаарах бөгөөд тиймээс $a_k$ коэффициентүүдийг шууд FFT-тэй ижил хуваа ба ялагтун алгоритмаар олж болно, зөвхөн $w_n^k$-ийн оронд $w_n^{-k}$ ашиглаж, эцэст нь гарсан коэффициентүүдийг $n$-д хуваах хэрэгтэй.

Тиймээс урвуу DFT-г тооцоолох нь шууд DFT-г тооцоолохтой бараг ижил бөгөөд үүнийг мөн $O(n \log n)$ хугацаанд гүйцэтгэж болно.

Implementation

Here we present a simple recursive implementation of the FFT and the inverse FFT, both in one function, since the difference between the forward and the inverse FFT are so minimal. To store the complex numbers we use the complex type in the C++ STL.

using cd = complex<double>;
const double PI = acos(-1);

void fft(vector<cd> & a, bool invert) {
    int n = a.size();
    if (n == 1)
        return;

    vector<cd> a0(n / 2), a1(n / 2);
    for (int i = 0; 2 * i < n; i++) {
        a0[i] = a[2*i];
        a1[i] = a[2*i+1];
    }
    fft(a0, invert);
    fft(a1, invert);

    double ang = 2 * PI / n * (invert ? -1 : 1);
    cd w(1), wn(cos(ang), sin(ang));
    for (int i = 0; 2 * i < n; 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;
    }
}

The function gets passed a vector of coefficients, and the function will compute the DFT or inverse DFT and store the result again in this vector. The argument $\text{invert}$ shows whether the direct or the inverse DFT should be computed. Inside the function we first check if the length of the vector is equal to one, if this is the case then we don't have to do anything. Otherwise we divide the vector $a$ into two vectors $a0$ and $a1$ and compute the DFT for both recursively. Then we initialize the value $wn$ and a variable $w$, which will contain the current power of $wn$. Then the values of the resulting DFT are computed using the above formulas.

If the flag $\text{invert}$ is set, then we replace $wn$ with $wn^{-1}$, and each of the values of the result is divided by $2$ (since this will be done in each level of the recursion, this will end up dividing the final values by $n$).

Using this function we can create a function for multiplying two polynomials:

vector<int> multiply(vector<int> const& a, vector<int> const& b) {
    vector<cd> fa(a.begin(), a.end()), fb(b.begin(), b.end());
    int n = 1;
    while (n < a.size() + b.size()) 
        n <<= 1;
    fa.resize(n);
    fb.resize(n);

    fft(fa, false);
    fft(fb, false);
    for (int i = 0; i < n; i++)
        fa[i] *= fb[i];
    fft(fa, true);

    vector<int> result(n);
    for (int i = 0; i < n; i++)
        result[i] = round(fa[i].real());
    return result;
}

This function works with polynomials with integer coefficients, however you can also adjust it to work with other types. Since there is some error when working with complex numbers, we need round the resulting coefficients at the end.

Finally the function for multiplying two long numbers practically doesn't differ from the function for multiplying polynomials. The only thing we have to do afterwards, is to normalize the number:

    int carry = 0;
    for (int i = 0; i < n; i++)
        result[i] += carry;
        carry = result[i] / 10;
        result[i] %= 10;
    }

Since the length of the product of two numbers never exceed the total length of both numbers, the size of the vector is enough to perform all carry operations.

Improved implementation: in-place computation

To increase the efficiency we will switch from the recursive implementation to an iterative one. In the above recursive implementation we explicitly separated the vector $a$ into two vectors - the element on the even positions got assigned to one temporary vector, and the elements on odd positions to another. However if we reorder the elements in a certain way, we don't need to create these temporary vectors (i.e. all the calculations can be done "in-place", right in the vector $A$ itself).

Note that at the first recursion level, the elements whose lowest bit of the position was zero got assigned to the vector $a_0$, and the ones with a one as the lowest bit of the position got assigned to $a_1$. In the second recursion level the same thing happens, but with the second lowest bit instead, etc. Therefore if we reverse the bits of the position of each coefficient, and sort them by these reversed values, we get the desired order (it is called the bit-reversal permutation).

For example the desired order for $n = 8$ has the form:

$$a = \bigg\{ \Big[ (a_0, a_4), (a_2, a_6) \Big], \Big[ (a_1, a_5), (a_3, a_7) \Big] \bigg\}$$

Indeed in the first recursion level (surrounded by curly braces), the vector gets divided into two parts $[a_0, a_2, a_4, a_6]$ and $[a_1, a_3, a_5, a_7]$. As we see, in the bit-reversal permutation this corresponds to simply dividing the vector into two halves: the first $\frac{n}{2}$ elements and the last $\frac{n}{2}$ elements. Then there is a recursive call for each halve. Let the resulting DFT for each of them be returned in place of the elements themselves (i.e. the first half and the second half of the vector $a$ respectively.

$$a = \bigg\{ \Big[y_0^0, y_1^0, y_2^0, y_3^0\Big], \Big[y_0^1, y_1^1, y_2^1, y_3^1 \Big] \bigg\}$$

Now we want to combine the two DFTs into one for the complete vector. The order of the elements is ideal, and we can also perform the union directly in this vector. We can take the elements $y_0^0$ and $y_0^1$ and perform the butterfly transform. The place of the resulting two values is the same as the place of the two initial values, so we get:

$$a = \bigg\{ \Big[y_0^0 + w_n^0 y_0^1, y_1^0, y_2^0, y_3^0\Big], \Big[y_0^0 - w_n^0 y_0^1, y_1^1, y_2^1, y_3^1\Big] \bigg\}$$

Similarly we can compute the butterfly transform of $y_1^0$ and $y_1^1$ and put the results in their place, and so on. As a result we get:

$$a = \bigg\{ \Big[y_0^0 + w_n^0 y_0^1, y_1^0 + w_n^1 y_1^1, y_2^0 + w_n^2 y_2^1, y_3^0 + w_n^3 y_3^1\Big], \Big[y_0^0 - w_n^0 y_0^1, y_1^0 - w_n^1 y_1^1, y_2^0 - w_n^2 y_2^1, y_3^0 - w_n^3 y_3^1\Big] \bigg\}$$

Thus we computed the required DFT from the vector $a$.

Here we described the process of computing the DFT only at the first recursion level, but the same works obviously also for all other levels. Thus, after applying the bit-reversal permutation, we can compute the DFT in-place, without any additional memory.

This additionally allows us to get rid of the recursion. We just start at the lowest level, i.e. we divide the vector into pairs and apply the butterfly transform to them. This results with the vector $a$ with the work of the last level applied. In the next step we divide the vector into vectors of size $4$, and again apply the butterfly transform, which gives us the DFT for each block of size $4$. And so on. Finally in the last step we obtained the result of the DFTs of both halves of $a$, and by applying the butterfly transform we obtain the DFT for the complete vector $a$.

using cd = complex<double>;
const double PI = acos(-1);

int reverse(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<cd> & a, bool invert) {
    int n = a.size();
    int lg_n = 0;
    while ((1 << lg_n) < n)
        lg_n++;

    for (int i = 0; i < n; i++) {
        if (i < reverse(i, lg_n))
            swap(a[i], a[reverse(i, lg_n)]);
    }

    for (int len = 2; len <= n; len <<= 1) {
        double ang = 2 * PI / len * (invert ? -1 : 1);
        cd wlen(cos(ang), sin(ang));
        for (int i = 0; i < n; i += len) {
            cd w(1);
            for (int j = 0; j < len / 2; j++) {
                cd 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 (cd & x : a)
            x /= n;
    }
}

At first we apply the bit-reversal permutation by swapping the each element with the element of the reversed position. Then the $\log n - 1$ states of the algorithm we compute the DFT for each block of the corresponding size $\text{len}$. For all those blocks we have the same root of unity $\text{wlen}$. We iterate all blocks and perform the butterfly transform on each of them.

We can further optimize the reversal of the bits. In the previous implementation we iterated all bits of the index and created the bitwise reversed index. However we can reverse the bits in a different way.

Suppose that $j$ already contains the reverse of $i$. Then by to go to $i + 1$, we have to increment $i$, and we also have to increment $j$, but in a "reversed" number system. Adding one in the conventional binary system is equivalent to flip all tailing ones into zeros and flipping the zero right before them into a one. Equivalently in the "reversed" number system, we flip all leading ones, and the also the next zero.

Thus we get the following implementation:

using cd = complex<double>;
const double PI = acos(-1);

void fft(vector<cd> & a, bool invert) {
    int n = 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);
        cd wlen(cos(ang), sin(ang));
        for (int i = 0; i < n; i += len) {
            cd w(1);
            for (int j = 0; j < len / 2; j++) {
                cd 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 (cd & x : a)
            x /= n;
    }
}

Additionally we can precompute the bit-reversal permutation beforehand. This is especially useful when the size $n$ is the same for all calls. But even when we only have three calls (which are necessary for multiplying two polynomials), the effect is noticeable. Also we can precompute all roots of unity and their powers.

Тооны онолын хувиргалт

Одоо бид зорилгоо бага зэрэг өөрчилье. Бид хоёр олон гишүүнтийг $O(n \log n)$ хугацаанд үржүүлэхийг хүссэн хэвээр байгаа ч энэ удаад коэффициентүүдийг ямар нэг анхны тоо $p$ модулиар тооцоолохыг хүсэж байна. Мэдээж энэ ажилд бид ердийн DFT ашиглаж, үр дүнд нь модулийн операторыг хэрэглэж болно. Гэвч ингэх нь дугуйруулах алдаанд хүргэж болзошгүй, ялангуяа том тоотой ажиллах үед. Тооны онолын хувиргалт (NTT) нь зөвхөн бүхэл тоотой ажилладаг тул үр дүн нь зөв байх баталгаатай гэсэн давуу талтай.

Дискрет Фурье хувиргалт нь комплекс тоо ба $n$ дугаар нэгжийн язгуурт тулгуурладаг. Үүнийг үр ашигтай тооцоолохын тулд бид язгуурын шинж чанарыг өргөнөөр ашигладаг (жишээ нь зэрэгт дэвшүүлэх замаар бусад бүх язгуурыг үүсгэдэг нэг язгуур байдаг).

Гэвч модулийн арифметик дэх $n$ дугаар нэгжийн язгуурын хувьд мөн ижил шинж чанар биелнэ. Анхдагч талбар дахь $n$ дугаар нэгжийн язгуур гэдэг нь дараах нөхцөлийг хангадаг $w_n$ тоо юм:

$$\begin{align} (w_n)^n &= 1 \pmod{p}, \\ (w_n)^k &\ne 1 \pmod{p}, \quad 1 \le k < n. \end{align}$$

Бусад $n-1$ язгуурыг $w_n$ язгуурын зэрэг хэлбэрээр авч болно.

Үүнийг хурдан Фурье хувиргалтын алгоритмд хэрэглэхийн тулд $2$-ын зэрэг болох ямар нэг $n$-ийн хувьд, мөн бүх бага зэргүүдийн хувьд язгуур оршин байх шаардлагатай. Бид дараах сонирхолтой шинж чанарыг анзаарч болно:

$$\begin{align} (w_n^2)^m = w_n^n &= 1 \pmod{p}, \quad \text{энд } m = \frac{n}{2}\\ (w_n^2)^k = w_n^{2k} &\ne 1 \pmod{p}, \quad 1 \le k < m. \end{align}$$

Тиймээс хэрэв $w_n$ нь $n$ дугаар нэгжийн язгуур бол $w_n^2$ нь $\frac{n}{2}$ дугаар нэгжийн язгуур болно. Улмаар хоёрын бүх бага зэргийн хувьд шаардлагатай зэрэгтэй язгуур оршин байх ба тэдгээрийг $w_n$ ашиглан тооцоолж болно.

Урвуу DFT-г тооцоолохын тулд бидэнд $w_n$-ийн урвуу $w_n^{-1}$ хэрэгтэй. Гэвч анхны модулийн хувьд урвуу үргэлж оршдог.

Тиймээс комплекс язгуураас бидэнд хэрэгтэй бүх шинж чанар нь $n$ дугаар нэгжийн язгуур оршин байх хангалттай том модуль $p$ байгаа тохиолдолд модулийн арифметикт мөн бэлэн байна.

Жишээ нь бид дараах утгуудыг авч болно: модуль $p = 7340033$, $w_{2^{20}} = 5$. Хэрэв энэ модуль хангалтгүй бол бид өөр хос олох хэрэгтэй. Бид $p = c 2^k + 1$ хэлбэрийн модулийн хувьд ($p$ анхны байхад) $2^k$ дугаар нэгжийн язгуур үргэлж оршдог гэсэн баримтыг ашиглаж болно. $g$ нь $p$-ийн анхдагч язгуур байхад $g^c$ нь ийм $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 = 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)(1LL * wlen * 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)(1LL * a[i+j+len/2] * 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)(1LL * w * wlen % mod);
            }
        }
    }

    if (invert) {
        int n_1 = inverse(n, mod);
        for (int & x : a)
            x = (int)(1LL * x * n_1 % mod);
    }
}

Here the function inverse computes the modular inverse (see Modular Multiplicative Inverse). The constants mod, root, root_pw determine the module and the root, and root_1 is the inverse of root modulo mod.

In practice this implementation is slower than the implementation using complex numbers (due to the huge number of modulo operations), but it has some advantages such as less memory usage and no rounding errors.

Дурын модультай үржүүлэлт

Энд бид өмнөх хэсэгтэй ижил зорилгод хүрэхийг хүсэж байна. $A(x)$ ба $B(x)$ хоёр олон гишүүнтийг үржүүлж, коэффициентүүдийг ямар нэг тоо $M$ модулиар тооцоолох. Тооны онолын хувиргалт нь зөвхөн тодорхой анхны тоонуудын хувьд ажилладаг. Модуль нь хүссэн хэлбэртэй биш байх тохиолдолд яах вэ?

Нэг сонголт бол $c 2^k + 1$ хэлбэрийн өөр өөр анхны тоогоор олон удаа тооны онолын хувиргалт хийж, дараа нь эцсийн коэффициентүүдийг тооцоолохын тулд Хятадын үлдэгдлийн теорем-ыг хэрэглэх явдал юм.

Өөр нэг сонголт бол $A(x)$ ба $B(x)$ олон гишүүнтүүдийг тус бүрийг нь хоёр жижиг олон гишүүнт болгон хуваарилах явдал юм

$$\begin{align} A(x) &= A_1(x) + A_2(x) \cdot C \\ B(x) &= B_1(x) + B_2(x) \cdot C \end{align}$$

энд $C \approx \sqrt{M}$.

Тэгвэл $A(x)$ ба $B(x)$-ийн үржвэрийг дараах байдлаар илэрхийлж болно:

$$A(x) \cdot B(x) = A_1(x) \cdot B_1(x) + \left(A_1(x) \cdot B_2(x) + A_2(x) \cdot B_1(x)\right)\cdot C + \left(A_2(x) \cdot B_2(x)\right)\cdot C^2$$

$A_1(x)$, $A_2(x)$, $B_1(x)$ ба $B_2(x)$ олон гишүүнтүүд нь зөвхөн $\sqrt{M}$-ээс бага коэффициент агуулах тул гарч ирэх бүх үржвэрийн коэффициентүүд $M \cdot n$-ээс бага байх ба энэ нь ердийн хөвөгч цэгтэй төрлөөр боловсруулахад ихэвчлэн хангалттай бага байдаг.

Тиймээс энэ арга нь бага коэффициенттэй олон гишүүнтүүдийн үржвэрийг тооцоолохыг (ердийн FFT ба урвуу FFT ашиглан) шаарддаг бөгөөд дараа нь анхны үржвэрийг модулийн нэмэх, үржүүлэх үйлдлээр $O(n)$ хугацаанд сэргээж болно.

Хэрэглээ

DFT-г эхлээд харахад олон гишүүнт үржүүлэхтэй ямар ч холбоогүй мэт олон төрлийн бусад бодлогод ашиглаж болно.

Бүх боломжит нийлбэр

Бидэнд $a[]$ ба $b[]$ хоёр массив өгөгдсөн. Бид бүх боломжит нийлбэр $a[i] + b[j]$-г олж, нийлбэр бүр хэдэн удаа гарч ирэхийг тоолох ёстой.

Жишээ нь $a = [1,~ 2,~ 3]$ ба $b = [2,~ 4]$-ийн хувьд бид дараахыг авна: $3$ нийлбэрийг $1$ аргаар, $4$ нийлбэрийг мөн $1$ аргаар, $5$-г $2$, $6$-г $1$, $7$-г $1$ аргаар авч болно.

Бид $a$ ба $b$ массивуудад $A$ ба $B$ хоёр олон гишүүнт байгуулна. Массивын тоонууд нь олон гишүүнт дэх зэрэг илтгэгчийн үүрэг гүйцэтгэнэ ($a[i] \Rightarrow x^{a[i]}$); мөн энэ гишүүний коэффициент нь тухайн тоо массивт хэдэн удаа гарч ирэхийг илэрхийлнэ.

Дараа нь эдгээр хоёр олон гишүүнтийг $O(n \log n)$ хугацаанд үржүүлснээр бид $C$ олон гишүүнт авах бөгөөд зэрэг илтгэгч нь ямар нийлбэр авч болохыг, коэффициент нь хэдэн удаа болохыг хэлж өгнө. Үүнийг жишээн дээр үзүүлбэл:

$$(1 x^1 + 1 x^2 + 1 x^3) (1 x^2 + 1 x^4) = 1 x^3 + 1 x^4 + 2 x^5 + 1 x^6 + 1 x^7$$

Бүх боломжит скаляр үржвэр

Бидэнд $n$ урттай $a[]$ ба $b[]$ хоёр массив өгөгдсөн. Бид $a$-г $b$-ийн циклик шилжилт бүртэй үржүүлэхийг тооцоолох ёстой.

Бид $2n$ хэмжээтэй хоёр шинэ массив үүсгэнэ: $a$-г урвуулж, түүн рүү $n$ тэг нэмнэ. Мөн $b$-г өөр дээр нь нэмнэ. Эдгээр хоёр массивыг олон гишүүнт мэт үржүүлээд, $c$ үржвэрийн $c[n-1],~ c[n],~ \dots,~ c[2n-2]$ коэффициентүүдийг харвал бид дараахыг авна:

$$c[k] = \sum_{i+j=k} a[i] b[j]$$

Мөн $i \ge n$-ийн хувьд бүх элемент $a[i] = 0$ тул:

$$c[k] = \sum_{i=0}^{n-1} a[i] b[k-i]$$

Энэ нийлбэр нь $a$ векторын $b$-ийн $(k - (n - 1))$ дугаар циклик зүүн шилжилттэй скаляр үржвэр болохыг харахад амархан. Тиймээс эдгээр коэффициент нь бодлогын хариу бөгөөд бид үүнийг $O(n \log n)$ хугацаанд авч чадсан. Энд $c[2n-1]$ нь мөн $n$ дугаар циклик шилжилтийг өгдөг боловч энэ нь $0$ дугаар циклик шилжилттэй ижил тул бид үүнийг хариудаа тусад нь авч үзэх шаардлагагүй гэдгийг анхаараарай.

Хоёр судал

Бидэнд $a$ ба $b$ хоёр Булийн судал ($0$ ба $1$ утгуудын циклик массив) өгөгдсөн. Бид эхний судлын $1$ нь хоёр дахь судлын $1$-ийн хажууд байхгүй байхаар эхний судлыг хоёр дахьд залгах бүх аргыг олохыг хүсэж байна.

Энэ бодлого өмнөх бодлогоос үнэндээ төдийлөн ялгардаггүй. Хоёр судлыг залгана гэдэг нь зүгээр л хоёр дахь массив дээр циклик шилжилт хийхийг хэлэх ба хоёр массивын скаляр үржвэр $0$ бол бид хоёр судлыг залгаж болно.

Тэмдэгт мөр тааруулах

Бидэнд жижиг үсгээс бүрдсэн текст $T$ ба хэв маяг $P$ гэсэн хоёр тэмдэгт мөр өгөгдсөн. Бид текст дэх хэв маягийн бүх тохиолдлыг тооцоолох ёстой.

Бид тэмдэгт мөр бүрд олон гишүүнт үүсгэнэ ($T[i]$ ба $P[I]$ нь цагаан толгойн $26$ үсэгт харгалзах $0$-ээс $25$ хоорондох тоо):

$$A(x) = a_0 x^0 + a_1 x^1 + \dots + a_{n-1} x^{n-1}, \quad n = |T|$$

энд

$$a_i = \cos(\alpha_i) + i \sin(\alpha_i), \quad \alpha_i = \frac{2 \pi T[i]}{26}.$$

Мөн

$$B(x) = b_0 x^0 + b_1 x^1 + \dots + b_{m-1} x^{m-1}, \quad m = |P|$$

энд

$$b_i = \cos(\beta_i) - i \sin(\beta_i), \quad \beta_i = \frac{2 \pi P[m-i-1]}{26}.$$

$P[m-i-1]$ илэрхийлэл нь хэв маягийг тодорхой урвуулж байгааг анхаараарай.

$C(x) = A(x) \cdot B(x)$ хоёр олон гишүүнтийн үржвэрийн $(m-1+i)$ дугаар коэффициент нь текстийн $i$ байрлалд хэв маяг гарч ирж байгаа эсэхийг хэлж өгнө.

$$c_{m-1+i} = \sum_{j = 0}^{m-1} a_{i+j} \cdot b_{m-1-j} = \sum_{j=0}^{m-1} \left(\cos(\alpha_{i+j}) + i \sin(\alpha_{i+j})\right) \cdot \left(\cos(\beta_j) - i \sin(\beta_j)\right)$$

энд $\alpha_{i+j} = \frac{2 \pi T[i+j]}{26}$ ба $\beta_j = \frac{2 \pi P[j]}{26}$

Хэрэв тохирол байвал $T[i+j] = P[j]$ болох ба тиймээс $\alpha_{i+j} = \beta_j$. Эндээс (Пифагорын тригонометрийн адилтгалыг ашиглан):

$$\begin{align} c_{m-1+i} &= \sum_{j = 0}^{m-1} \left(\cos(\alpha_{i+j}) + i \sin(\alpha_{i+j})\right) \cdot \left(\cos(\alpha_{i+j}) - i \sin(\alpha_{i+j})\right) \\ &= \sum_{j = 0}^{m-1} \cos(\alpha_{i+j})^2 + \sin(\alpha_{i+j})^2 = \sum_{j = 0}^{m-1} 1 = m \end{align}$$

Хэрэв тохирол байхгүй бол дор хаяж нэг тэмдэгт өөр байх ба энэ нь $a_{i+1} \cdot b_{m-1-j}$ үржвэрүүдийн нэг нь $1$-тэй тэнцүү биш байхад хүргэж, улмаар $c_{m-1+i} \ne m$ коэффициент гарна.

Орлуулагч тэмдэгттэй тэмдэгт мөр тааруулах

Энэ бол өмнөх бодлогын өргөтгөл юм. Энэ удаад бид хэв маяг нь боломжит үсэг бүртэй тохирч чадах орлуулагч тэмдэгт $\*$ агуулахыг зөвшөөрнө. Жишээ нь $a*c$ хэв маяг нь $abccaacc$ текстэд яг гурван байрлалд, индекс $0$, индекс $4$ ба индекс $5$-д гарч ирнэ.

Бид яг ижил олон гишүүнт үүсгэх ба зөвхөн $P[m-i-1] = *$ бол $b_i = 0$ гэж тавина. Хэрэв $x$ нь $P$ дахь орлуулагч тэмдэгтийн тоо бол $c_{m-1+i} = m - x$ үед $T$-д $i$ индекст $P$-ийн тохирол байна.

Дасгал бодлогууд