پرش به محتویات

تبدیل فوریه سریع

در این مقاله، الگوریتمی را مورد بحث قرار می‌دهیم که به ما اجازه می‌دهد دو چندجمله‌ای به طول $n$ را در زمان $O(n \log n)$ ضرب کنیم، که بهتر از ضرب معمولی با زمان $O(n^2)$ است. بدیهی است که ضرب دو عدد طولانی نیز می‌تواند به ضرب چندجمله‌ای‌ها کاهش یابد، بنابراین دو عدد طولانی نیز می‌توانند در زمان $O(n \log n)$ ضرب شوند (که در آن $n$ تعداد ارقام اعداد است).

کشف تبدیل فوریه سریع (FFT) به کولی و توکی نسبت داده می‌شود که در سال ۱۹۶۵ الگوریتمی را منتشر کردند. اما در واقع FFT بارها پیش از آن کشف شده بود، اما اهمیت آن قبل از اختراع کامپیوترهای مدرن درک نشده بود. برخی محققان کشف FFT را به رونگه و کونیگ در سال ۱۹۲۴ نسبت می‌دهند. اما در واقع، گاوس چنین روشی را در سال ۱۸۰۵ توسعه داده بود، اما هرگز آن را منتشر نکرد.

توجه داشته باشید که الگوریتم FFT ارائه‌شده در اینجا در زمان $O(n \log n)$ اجرا می‌شود، اما برای ضرب چندجمله‌ای‌های دلخواه بزرگ با ضرایب دلخواه بزرگ یا برای ضرب اعداد صحیح دلخواه بزرگ کار نمی‌کند. این الگوریتم به راحتی می‌تواند چندجمله‌ای‌هایی با اندازه $10^5$ و ضرایب کوچک را مدیریت کند یا دو عدد با اندازه $10^6$ را ضرب کند، که معمولاً برای حل مسائل برنامه‌نویسی رقابتی کافی است. فراتر از مقیاس ضرب اعداد با $10^6$ بیت، محدوده و دقت اعداد ممیز شناور مورد استفاده در طول محاسبات برای ارائه نتایج نهایی دقیق کافی نخواهد بود، اگرچه نسخه‌های پیچیده‌تری وجود دارند که می‌توانند ضرب‌های چندجمله‌ای/اعداد صحیح دلخواه بزرگ را انجام دهند. به عنوان مثال، در سال ۱۹۷۱، شونهاگه و اشتراسر نسخه‌ای برای ضرب اعداد دلخواه بزرگ توسعه دادند که FFT را به صورت بازگشتی در ساختارهای حلقه‌ای با زمان اجرای $O(n \log n \log \log n)$ به کار می‌برد. و اخیراً (در سال ۲۰۱۹) هاروی و فان در هوون الگوریتمی را منتشر کردند که دقیقاً در زمان $O(n \log n)$ اجرا می‌شود.

تبدیل فوریه گسسته

یک چندجمله‌ای از درجه $n-1$ را در نظر بگیرید:

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

بدون از دست دادن کلیت، فرض می‌کنیم که $n$ - تعداد ضرایب - توانی از ۲ است. اگر $n$ توانی از ۲ نباشد، به سادگی جملات недостающие $a_i x^i$ را اضافه کرده و ضرایب $a_i$ را برابر با ۰ قرار می‌دهیم.

نظریه اعداد مختلط به ما می‌گوید که معادله $x^n = 1$ دارای $n$ جواب مختلط است (که به آنها ریشه‌های $n$-ام واحد گفته می‌شود) و جواب‌ها به شکل $w_{n, k} = e^{\frac{2 k \pi i}{n}}$ برای $k = 0 \dots n-1$ هستند. علاوه بر این، این اعداد مختلط خواص بسیار جالبی دارند: برای مثال، ریشه اصلی $n$-ام $w_n = w_{n, 1} = e^{\frac{2 \pi i}{n}}$ می‌تواند برای توصیف تمام ریشه‌های $n$-ام دیگر استفاده شود: $w_{n, k} = (w_n)^k$.

تبدیل فوریه گسسته (DFT) چندجمله‌ای $A(x)$ (یا به طور معادل بردار ضرایب $(a_0, a_1, \dots, a_{n-1})$) به عنوان مقادیر چندجمله‌ای در نقاط $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})$، ضرایب چندجمله‌ای $(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)$ را - با ضرب هر عنصر از یک بردار در عنصر متناظر از بردار دیگر - ضرب کنیم، چیزی جز DFT چندجمله‌ای $\text{DFT}(A \cdot B)$ به دست نمی‌آوریم:

$$\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)$ محاسبه کنیم، آنگاه می‌توانیم حاصل‌ضرب دو چندجمله‌ای (و در نتیجه دو عدد طولانی) را با همان پیچیدگی زمانی محاسبه کنیم.

باید توجه داشت که دو چندجمله‌ای باید درجه یکسانی داشته باشند. در غیر این صورت، دو بردار حاصل از DFT طول‌های متفاوتی خواهند داشت. ما می‌توانیم این کار را با افزودن ضرایب با مقدار ۰ انجام دهیم.

و همچنین، از آنجایی که حاصل‌ضرب دو چندجمله‌ای، چندجمله‌ای با درجه $2(n-1)$ است، باید درجه هر چندجمله‌ای را دو برابر کنیم (دوباره با پر کردن با ۰). از یک بردار با $n$ مقدار نمی‌توانیم چندجمله‌ای مورد نظر با $2n-1$ ضریب را بازسازی کنیم.

تبدیل فوریه سریع

تبدیل فوریه سریع روشی است که امکان محاسبه DFT را در زمان $O(n \log n)$ فراهم می‌کند. ایده اصلی FFT استفاده از روش تقسیم و حل است. ما بردار ضرایب چندجمله‌ای را به دو بردار تقسیم می‌کنیم، به صورت بازگشتی DFT را برای هر کدام محاسبه می‌کنیم، و نتایج را برای محاسبه DFT چندجمله‌ای کامل ترکیب می‌کنیم.

پس یک چندجمله‌ای $A(x)$ با درجه $n-1$ را در نظر بگیرید، که در آن $n$ توانی از ۲ است و $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)$ را در زمان خطی با استفاده از $\text{DFT}(A_0)$ و $\text{DFT}(A_1)$ محاسبه کنیم، آنگاه رابطه بازگشتی $T_{\text{DFT}}(n) = 2 T_{\text{DFT}}\left(\frac{n}{2}\right) + O(n)$ را برای پیچیدگی زمانی به دست می‌آوریم، که طبق قضیه اصلی (Master Theorem) به $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$ گاهی اوقات پروانه (butterfly) نامیده می‌شود.)

بنابراین یاد گرفتیم که چگونه DFT را در زمان $O(n \log n)$ محاسبه کنیم.

تبدیل فوریه معکوس

فرض کنید بردار $(y_0, y_1, \dots y_{n-1})$ - مقادیر چندجمله‌ای $A$ از درجه $n-1$ در نقاط $x = w_n^k$ - داده شده است. ما می‌خواهیم ضرایب $(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} $$

این ماتریس، ماتریس وندرموند نامیده می‌شود.

بنابراین می‌توانیم بردار $(a_0, a_1, \dots, a_{n-1})$ را با ضرب بردار $(y_0, y_1, \dots y_{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)$ انجام شود.

پیاده‌سازی

در اینجا یک پیاده‌سازی بازگشتی ساده از FFT و FFT معکوس را ارائه می‌دهیم، هر دو در یک تابع، زیرا تفاوت بین FFT مستقیم و معکوس بسیار ناچیز است. برای ذخیره اعداد مختلط از نوع complex در کتابخانه استاندارد ++C استفاده می‌کنیم.

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;
    }
}

تابع یک بردار از ضرایب را دریافت می‌کند و DFT یا DFT معکوس را محاسبه کرده و نتیجه را دوباره در همین بردار ذخیره می‌کند. آرگومان invert نشان می‌دهد که آیا DFT مستقیم یا معکوس باید محاسبه شود. درون تابع ابتدا بررسی می‌کنیم که آیا طول بردار برابر با یک است، اگر این‌طور باشد، نیازی به انجام کاری نداریم. در غیر این صورت، بردار a را به دو بردار a0 و a1 تقسیم می‌کنیم و DFT را برای هر دو به صورت بازگشتی محاسبه می‌کنیم. سپس مقدار wn و متغیر w را که توان فعلی wn را نگه می‌دارد، مقداردهی اولیه می‌کنیم. سپس مقادیر DFT حاصل با استفاده از فرمول‌های بالا محاسبه می‌شوند.

اگر پرچم invert تنظیم شده باشد، wn را با wn^{-1} جایگزین می‌کنیم و هر یک از مقادیر نتیجه بر ۲ تقسیم می‌شود (از آنجایی که این کار در هر سطح از بازگشت انجام می‌شود، در نهایت مقادیر نهایی بر $n$ تقسیم خواهند شد).

با استفاده از این تابع می‌توانیم تابعی برای ضرب دو چندجمله‌ای ایجاد کنیم:

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;
}

این تابع با چندجمله‌ای‌هایی با ضرایب صحیح کار می‌کند، اما می‌توانید آن را برای کار با انواع دیگر داده‌ها نیز تنظیم کنید. از آنجایی که هنگام کار با اعداد مختلط مقداری خطا وجود دارد، باید در انتها ضرایب حاصل را گرد کنیم.

در نهایت، تابع ضرب دو عدد طولانی عملاً با تابع ضرب چندجمله‌ای‌ها تفاوتی ندارد. تنها کاری که باید بعد از آن انجام دهیم، نرمال‌سازی عدد است:

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

از آنجایی که طول حاصل‌ضرب دو عدد هرگز از مجموع طول هر دو عدد بیشتر نمی‌شود، اندازه بردار برای انجام تمام عملیات نقلی کافی است.

پیاده‌سازی بهبودیافته: محاسبه درجا

برای افزایش کارایی، از پیاده‌سازی بازگشتی به یک پیاده‌سازی تکراری تغییر می‌دهیم. در پیاده‌سازی بازگشتی بالا، ما به صراحت بردار a را به دو بردار - عناصری که در موقعیت‌های زوج قرار داشتند به یک بردار موقت و عناصری که در موقعیت‌های فرد قرار داشتند به بردار موقت دیگر - تقسیم کردیم. اما اگر عناصر را به روش خاصی مرتب کنیم، دیگر نیازی به ایجاد این بردارهای موقت نداریم (یعنی تمام محاسبات می‌تواند "درجا" و مستقیماً در خود بردار A انجام شود).

توجه داشته باشید که در اولین سطح بازگشت، عناصری که کم‌ارزش‌ترین بیت موقعیتشان صفر بود به بردار a_0 و آنهایی که کم‌ارزش‌ترین بیت موقعیتشان یک بود به a_1 اختصاص داده شدند. در سطح دوم بازگشت، همین اتفاق می‌افتد، اما این بار با دومین بیت کم‌ارزش و غیره. بنابراین، اگر بیت‌های موقعیت هر ضریب را معکوس کنیم و آنها را بر اساس این مقادیر معکوس شده مرتب کنیم، ترتیب مورد نظر را به دست می‌آوریم (این ترتیب، جایگشت بیت-معکوس نامیده می‌شود).

برای مثال، ترتیب مورد نظر برای $n=8$ به این شکل است:

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

در واقع، در اولین سطح بازگشت (که با آکولادها احاطه شده است)، بردار به دو بخش $[a_0, a_2, a_4, a_6]$ و $[a_1, a_3, a_5, a_7]$ تقسیم می‌شود. همانطور که می‌بینیم، در جایگشت بیت-معکوس این معادل تقسیم ساده بردار به دو نیمه است: $\frac{n}{2}$ عنصر اول و $\frac{n}{2}$ عنصر آخر. سپس برای هر نیمه یک فراخوانی بازگشتی وجود دارد. فرض کنید DFT حاصل برای هر یک از آنها در جای خود عناصر بازگردانده شود (یعنی به ترتیب نیمه اول و دوم بردار a).

$$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\}$$

حالا می‌خواهیم دو DFT را در یک DFT برای کل بردار ترکیب کنیم. ترتیب عناصر ایده‌آل است و ما می‌توانیم این ادغام را نیز مستقیماً در این بردار انجام دهیم. می‌توانیم عناصر $y_0^0$ و $y_0^1$ را برداشته و تبدیل پروانه‌ای را روی آنها اعمال کنیم. محل دو مقدار حاصل همان محل دو مقدار اولیه است، بنابراین به دست می‌آوریم:

$$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\}$$

به طور مشابه می‌توانیم تبدیل پروانه‌ای $y_1^0$ و $y_1^1$ را محاسبه کرده و نتایج را در جای خودشان قرار دهیم، و به همین ترتیب ادامه دهیم. در نتیجه به دست می‌آوریم:

$$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\}$$

بنابراین DFT مورد نیاز را از بردار a محاسبه کردیم.

در اینجا ما فرآیند محاسبه DFT را فقط در اولین سطح بازگشت توصیف کردیم، اما همین فرآیند به وضوح برای تمام سطوح دیگر نیز کار می‌کند. بنابراین، پس از اعمال جایگشت بیت-معکوس، می‌توانیم DFT را به صورت درجا و بدون هیچ حافظه اضافی محاسبه کنیم.

این کار به ما اجازه می‌دهد تا از بازگشت خلاص شویم. ما فقط از پایین‌ترین سطح شروع می‌کنیم، یعنی بردار را به جفت‌ها تقسیم کرده و تبدیل پروانه‌ای را روی آنها اعمال می‌کنیم. این کار منجر به بردار a با اعمال کار آخرین سطح می‌شود. در مرحله بعد، بردار را به بردارهایی با اندازه ۴ تقسیم می‌کنیم و دوباره تبدیل پروانه‌ای را اعمال می‌کنیم، که به ما DFT را برای هر بلوک با اندازه ۴ می‌دهد. و به همین ترتیب ادامه می‌دهیم. در نهایت، در آخرین مرحله، نتیجه DFTهای هر دو نیمه a را به دست می‌آوریم و با اعمال تبدیل پروانه‌ای، DFT را برای کل بردار 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;
    }
}

ابتدا با جابجایی هر عنصر با عنصر موقعیت معکوس شده، جایگشت بیت-معکوس را اعمال می‌کنیم. سپس در $\log n - 1$ مرحله الگوریتم، DFT را برای هر بلوک با اندازه متناظر len محاسبه می‌کنیم. برای تمام این بلوک‌ها، ریشه واحد wlen یکسانی داریم. ما روی تمام بلوک‌ها تکرار می‌کنیم و تبدیل پروانه‌ای را روی هر یک از آنها انجام می‌دهیم.

می‌توانیم معکوس‌سازی بیت‌ها را بیشتر بهینه کنیم. در پیاده‌سازی قبلی، ما روی تمام بیت‌های اندیس تکرار کرده و اندیس معکوس بیتی را ایجاد کردیم. اما می‌توانیم بیت‌ها را به روش دیگری معکوس کنیم.

فرض کنید $j$ از قبل حاوی معکوس $i$ است. سپس برای رفتن به $i+1$، باید $i$ را افزایش دهیم، و همچنین باید $j$ را نیز افزایش دهیم، اما در یک سیستم عددی "معکوس". اضافه کردن یک در سیستم باینری معمولی معادل تغییر تمام یک‌های انتهایی به صفر و تغییر صفر درست قبل از آنها به یک است. به طور معادل، در سیستم عددی "معکوس"، ما تمام یک‌های ابتدایی و همچنین صفر بعدی را تغییر می‌دهیم.

بنابراین پیاده‌سازی زیر را به دست می‌آوریم:

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;
    }
}

علاوه بر این، می‌توانیم جایگشت بیت-معکوس را از قبل محاسبه کنیم. این کار به ویژه زمانی مفید است که اندازه $n$ برای تمام فراخوانی‌ها یکسان باشد. اما حتی زمانی که فقط سه فراخوانی داریم (که برای ضرب دو چندجمله‌ای لازم است)، تأثیر آن قابل توجه است. همچنین می‌توانیم تمام ریشه‌های واحد و توان‌های آنها را از قبل محاسبه کنیم.

تبدیل نظریه اعداد

حالا هدف را کمی تغییر می‌دهیم. ما هنوز می‌خواهیم دو چندجمله‌ای را در زمان $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$ به دست آورد.

برای اعمال آن در الگوریتم تبدیل فوریه سریع، نیاز داریم که یک ریشه برای یک $n$ که توانی از ۲ است، و همچنین برای تمام توان‌های کوچکتر وجود داشته باشد. می‌توانیم به خاصیت جالب زیر توجه کنیم:

$$\begin{align} (w_n^2)^m = w_n^n &= 1 \pmod{p}, \quad \text{with } 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^{-1}$ از $w_n$ نیاز داریم. اما برای یک پیمانه اول، معکوس همیشه وجود دارد.

بنابراین تمام خواصی که از ریشه‌های مختلط نیاز داریم، در حساب پیمانه‌ای نیز موجود است، به شرطی که یک پیمانه $p$ به اندازه کافی بزرگ داشته باشیم که برای آن یک ریشه $n$-ام واحد وجود داشته باشد.

برای مثال، می‌توانیم مقادیر زیر را در نظر بگیریم: پیمانه $p = 7340033$ و $w_{2^{20}} = 5$. اگر این پیمانه کافی نباشد، باید یک جفت دیگر پیدا کنیم. می‌توانیم از این واقعیت استفاده کنیم که برای پیمانه‌هایی به شکل $p = c 2^k + 1$ (و $p$ اول است)، همیشه ریشه $2^k$-ام واحد وجود دارد. می‌توان نشان داد که $g^c$ چنین ریشه $2^k$-ام واحدی است، که در آن $g$ یک ریشه اولیه از $p$ است.

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);
    }
}

در اینجا تابع inverse معکوس پیمانه‌ای را محاسبه می‌کند (به معکوس ضربی پیمانه‌ای مراجعه کنید). ثابت‌های mod، root و root_pw پیمانه و ریشه را تعیین می‌کنند و root_1 معکوس root به پیمانه mod است.

در عمل، این پیاده‌سازی کندتر از پیاده‌سازی با استفاده از اعداد مختلط است (به دلیل تعداد زیاد عملیات پیمانه)، اما مزایایی مانند استفاده کمتر از حافظه و عدم وجود خطاهای گرد کردن دارد.

ضرب با پیمانه دلخواه

در اینجا می‌خواهیم به همان هدفی برسیم که در بخش قبل داشتیم. ضرب دو چندجمله‌ای $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$ به ۱ روش، مجموع $4$ نیز به ۱ روش، $5$ به ۲ روش، $6$ به ۱ روش و $7$ به ۱ روش قابل دستیابی است.

ما برای آرایه‌های $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$$

تمام ضرب‌های داخلی ممکن

دو آرایه $a[]$ و $b[]$ به طول $n$ به ما داده شده است. باید حاصل‌ضرب $a$ با هر شیفت دایره‌ای از $b$ را محاسبه کنیم.

دو آرایه جدید به اندازه $2n$ ایجاد می‌کنیم: $a$ را معکوس کرده و $n$ صفر به آن اضافه می‌کنیم. و $b$ را به خودش الحاق می‌کنیم. وقتی این دو آرایه را به عنوان چندجمله‌ای ضرب می‌کنیم و به ضرایب $c[n-1],~ c[n],~ \dots,~ c[2n-2]$ حاصل‌ضرب $c$ نگاه می‌کنیم، به دست می‌آوریم:

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

و از آنجایی که تمام عناصر $a[i] = 0$ برای $i \ge n$ هستند:

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

به راحتی می‌توان دید که این مجموع همان ضرب داخلی بردار $a$ با شیفت دایره‌ای چپ $(k - (n - 1))$-ام از $b$ است. بنابراین این ضرایب پاسخ مسئله هستند و ما همچنان توانستیم آن را در زمان $O(n \log n)$ به دست آوریم. توجه داشته باشید که $c[2n-1]$ نیز به ما شیفت دایره‌ای $n$-ام را می‌دهد اما این همان شیفت دایره‌ای صفرم است، بنابراین نیازی به در نظر گرفتن آن به طور جداگانه در پاسخ خود نداریم.

دو نوار

دو نوار بولی (آرایه‌های دایره‌ای از مقادیر $0$ و $1$) $a$ و $b$ به ما داده شده است. می‌خواهیم تمام راه‌های اتصال نوار اول به نوار دوم را پیدا کنیم، به طوری که در هیچ موقعیتی یک $1$ از نوار اول در کنار یک $1$ از نوار دوم قرار نگیرد.

این مسئله در واقع تفاوت چندانی با مسئله قبلی ندارد. اتصال دو نوار به این معنی است که ما یک شیفت دایره‌ای روی آرایه دوم انجام می‌دهیم و می‌توانیم دو نوار را متصل کنیم، اگر ضرب داخلی دو آرایه $0$ باشد.

تطابق رشته

دو رشته به ما داده شده است، یک متن $T$ و یک الگو $P$ که از حروف کوچک تشکیل شده‌اند. باید تمام رخدادهای الگو را در متن محاسبه کنیم.

برای هر رشته یک چندجمله‌ای ایجاد می‌کنیم ($T[i]$ و $P[I]$ اعدادی بین $0$ تا $25$ متناظر با $26$ حرف الفبا هستند):

$$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]$ الگو به صراحت معکوس می‌شود.

ضریب $(m-1+i)$-ام حاصل‌ضرب دو چندجمله‌ای $C(x) = A(x) \cdot B(x)$ به ما خواهد گفت که آیا الگو در متن در موقعیت $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$ منجر می‌شود.

تطابق رشته با کاراکترهای Wildcard

این تعمیمی از مسئله قبلی است. این بار اجازه می‌دهیم که الگو شامل کاراکتر wildcard * باشد که می‌تواند با هر حرف ممکنی تطابق داشته باشد. برای مثال، الگوی $a*c$ در متن $abccaacc$ دقیقاً در سه موقعیت، در اندیس $0$، اندیس $4$ و اندیس $5$ ظاهر می‌شود.

ما دقیقاً همان چندجمله‌ای‌ها را ایجاد می‌کنیم، با این تفاوت که اگر $P[m-i-1] = *$ باشد، $b_i = 0$ قرار می‌دهیم. اگر $x$ تعداد کاراکترهای wildcard در $P$ باشد، آنگاه ما یک تطابق از $P$ در $T$ در اندیس $i$ خواهیم داشت اگر $c_{m-1+i} = m - x$ باشد.

مسائل تمرینی