android 快速傅里叶变换库快速傅里叶变换计算量

转载

flyingsmiling 2024-07-31 17:05:14

文章标签 android 快速傅里叶变换库多项式 bc 2d 文章分类 Android 移动开发

快速傅里叶变换($ \rm Fast\ Fourier\ Transformation $), 简称 $\rm FFT$, 用于在 $ \Theta(n\log n) $ 时间内求两个多项式的乘积. 快速数论变换($ \rm Fast\ Number\ Theoretic\ Transforms$), 简称 $\rm NTT$, 用于在 $ \Theta(n\log n) $ 时间内求两个多项式的乘积, 系数对 $p$ 取模.

快速傅里叶变换(FFT)学习笔记

简介

快速傅里叶变换($ \rm Fast\ Fourier\ Transformation $), 简称 $\rm FFT$, 用于在 $ \Theta(n\log n) $ 时间内求两个多项式的乘积.

快速数论变换($ \rm Fast\ Number\ Theoretic\ Transforms$), 简称 $\rm NTT$, 用于在 $ \Theta(n\log n) $ 时间内求两个多项式的乘积, 系数对 $p$

前置技能

卷积

(注: 以下所有多项式均只包含 $ x $ 一个变量)

一个 $ n - 1 $ 次 $ n $ 项式 $ f(x) $ 可以表示为 $ f(x) = \sum_{i=0}^{n - 1} a_ix^i $.

设有多项式 $ f(x) = \sum_{i=0}^{n - 1} a_i x^i $ 和 $ g(x) = \sum_{i = 0}^{n - 1} b_ix^i $, 则有:

\[f(x)\cdot g(x) = (f \cdot g)(x) = h(x) = \sum_{i = 0}^{n - 1} a_ix^ig(x) = \sum_{i = 0}^{n - 1} a_ix^i\left(\sum_{i = 0}^{n - 1} b_ix^i\right) = \sum_{k = 0}^{2n - 2}\sum_{i + j = k} a_i b_j x^{i+j} \]

我们称 $ (f \cdot g)(x) $ 为 $ f(x) $ 和 $ g(x) $ 的卷积.

复数

定义

定义 $ i = \sqrt{-1} $, 则出现了形如 $ a+bi(a,b\in \mathbb{R}) $ 的数.

$ a + bi $ 可以被表示为复平面上 $(a, b)$ 的点, 也可以表示为向量 $(a, b)$.

复数的模(长):

$|a + bi| = \sqrt{a^2 + b^2} $.

为该复数所表示的点在复平面上到原点的距离.

复数的幅角:

一个复数的幅角为该数在复平面上与实轴正半轴的夹角(逆时针)记作 $ arg(a + bi) $, 显然复数的幅角有无穷多个, 每个幅角都相差 $2\pi$, 若 $ a + bi $ 的幅角 $ \theta$ 满足 $ -\pi \le \theta \le \pi $, 则称 $ \theta $ 是 $ a + bi $ 的幅角主值.

$0$

运算

加(减)法:

$ (a + bi) + (c + di) = (a + c) + (b + d)i h(x) = f(x)\cdot g(x) $.

乘法:

$ (a + bi)(c + di) = ac + adi + bci - bd = (ac - bd) + (bc + ad)i $.

android 快速傅里叶变换库快速傅里叶变换计算量_多项式

如图, $ A \cdot B = C $, 复数的乘法有以下性质: 模长相乘, 幅角相加.

模长相乘可以以如下方法证明:

设 $ z_1 = a + bi $, $ z_2 = c + di $, $ z_3 = z_1z_2 = (ac - bd) + (ad + bc) i $;

有:

$|z_1|\cdot |z_2| = \sqrt{a^2 + b^2}\cdot \sqrt{c^2 + d^2} = \sqrt{(a^2 + b^2)(c2 + d^2)} = \sqrt{a^2c2 + a^2d2 + b^2c2 + b^2d2} $

$|z_3| = \sqrt{(ac - bd) ^ 2 + (ad + bc) ^ 2} = \sqrt{a^2c^2 + a^2d^2 + b^2c^2 + b^2d^2} = |z_1|\cdot |z_2|$

共轭:

$ \overline{a + bi} = a - bi $, 易证任何复数乘上它的共轭都为实数, $ (a + bi)(a - bi) = a ^ 2 + b ^2 $.

除法:

\[\frac{a + bi}{c + di} = \frac{(a + bi)(c - di)}{c^2 + d^2} = \frac{ac + bd}{c^2 + d^2} + \frac{ad + bc}{c^2 + d^2}i \]

欧拉公式

求证 $e^{i\theta} = \cos \theta + i \sin \theta$;

设

\[f(\theta) = \frac{e^{i\theta}}{ \cos \theta + i \sin \theta} \]

可得

\[f'(\theta) = \frac{(e^{i\theta})'(\cos \theta + i \sin \theta) - e^{i\theta}(\cos \theta + i \sin \theta)'}{(\cos \theta + i \sin \theta)^2} \]

\[= \frac{ie^{i\theta}(\cos \theta + i \sin \theta) - e^{i\theta}(-\sin \theta + i \cos \theta)}{(\cos \theta + i \sin \theta)^2} \]

\[= \frac{e^{i\theta} i \cos \theta - e^{i\theta} \sin \theta + e^{i\theta}\sin \theta - e^{i\theta}i \cos \theta}{(\cos \theta + i \sin \theta)^2} = 0 \]

所以, $ f(\theta) $ 是常数.

求 $ f(0) $, 得 $ f(\theta) = f(0) = \frac{1}{1 + 0} = 1 $.

所以, $ e^{i\theta} = \cos \theta + i \sin \theta $;

证毕.

单位根

定义

定义 $ n(n\ge 1) $ 次单位根的 $ n $ 次幂等于 $ 1 $.

显然单位根的模长等于 $ 1 $, 否则怎么乘都不会是 $ 1 $.

所以, 单位根都在单位圆上.

考虑设单位根幅角为 $ \alpha $, 则有 $ n\alpha = 2k\pi, (k\in \mathbb N) $.

解得:

\[\alpha = \frac{2k\pi}{n} = \frac{2\pi}{n}k \]

所以, 单位根 $ n $ 等分单位圆.

android 快速傅里叶变换库快速傅里叶变换计算量_bc_02

注意到, 当 $ k = n $ 时, 进入了循环, 所以只有 $n$

对于 $ n $ 次的幅角为 $ \frac{2\pi}{n}k $ 的单位根我们记作 $ \omega_n^k $.

显然, $ \omega_n^k = \omega_n^{k\bmod n} $.

得到幅角 $ \alpha $ 后, 我们可以很轻易的得到该单位根所对应的数.

引理

android 快速傅里叶变换库快速傅里叶变换计算量_2d_03

如上图, $ A'B = OA \cdot \sin \alpha = \sin \alpha $;

$ OB = OA \cdot \cos \alpha = \cos \alpha $.

所以点 $ A' $ 所代表的数为 $ \cos \alpha + i\sin \alpha $.

所以

\[\begin{array}{} \omega_n^k & = & \cos \frac{2\pi}{n}k + i\sin \frac{2\pi}{n}k \\ & = & e^{i\frac{2\pi}{n}k} = \left(e^{i\frac{2\pi}{n}}\right)^k \\ & = & (\omega_n^1)^k \end{array} \]

由上式, 有:

引理 1:

$ \omega_n^j\omega_nk = (\omega_n¹⁾j(\omega_n¹⁾k = (\omega_n¹⁾ =\omega_n^{j + k} $;

引理 2:

同上, 显然有 $ (\omega_n^j)k = \omega_n^{kj} $;

引理 3 (消去引理):

$ \omega_{dn}^{dk} = e^{i\frac{2\pi}{dn}dk} = e^{i\frac{2\pi}{n}k} = \omega_n^k $;

引理 4 (折半引理):

由消去引理可得 $ \omega_{2n}^n = \omega_2^1 = -1 $;

$ \omega_{2n}^{k + n} = \omega_{2n}^k\omega_{2n}n = -\omega_{2n}^k $;

引理 5:

当 $ n $ 为偶数时.

根据消去引理, 有:

$ (\omega_n^k)2 = \omega_n^{2k} = \omega^k_{\frac{n}{2}} $

快速傅里叶变换

离散傅立叶变换(DFT)

多项式的点值表示

定义

显然, $ n + 1 $ 个点确定一个 $ n $ 项式.

所以我们可以用形如 $ (x_0, y_0),(x_1,y_1),(x_2,y_2),\dots,(x_n,y_n) $ 来表示一个多项式.

点值表示下的卷积

显然, 设有 $ (x_0, y_0),(x_1,y_1),(x_2,y_2),\dots,(x_n,y_n) $ 表示$ f(x) $.

有 $ (x_0, y_0'),(x_1,y_1'),(x_2,y_2'),\dots,(x_n,y_n') $ 表示$ g(x) $.

则由于 $ h(x) = (f\cdot g)(x) = f(x)\cdot g(x) $, 很自然的就有:

$ (x_0, y_0\cdot y_0'),(x_1, y_1\cdot y_1'),(x_2, y_2\cdot y_2'),\cdots,(x_n, y_n\cdot y_n') $ 表示 $ h(x) $.

我们就有很快的 $ \Theta(n) $ 算法来解决已经求出点值表示的两个多项式相乘.

需要注意的是, 由于 $ h $ 理论上有 $ 2n + 1 $ 项, 所以需要在先前多造几组点.

求值

(根据系数表示求点值表示叫做求值)

(以下所有 $ n $ 都保证 $ n = 2^t, t\in \mathbb N $, 实际情况可以补一些系数为 $0$

设:

\[f(x) = \sum_{i = 0}^{n-1} a_i x^i \]

我们要求出其点值表示.

于是我们想(luan)到(gao), 给 f 数组代入单位根, 并把奇数项和偶数项分开.

设 $ f_0(x) = \sum_{i = 0} ^ {\frac{n}{2} - 1} a_{2i} x^i $

$ f_1(x) = \sum_{i = 0} ^ {\frac{n}{2} - 1} a_{2i + 1} x^i $

显然有: $f(x) = f_0(x^2) + xf_1(x^2)$.

我们代入 $\omega_n^k$, 得:

容易得到:

\[\begin{array}{} f(\omega_n^k) & = & f_0((\omega_n^k)^2) + \omega_n^kf_1((\omega_n^k)^2) \\ & = & f_0\left(\omega_{\frac{n}{2}}^k\right) + \omega_n^k f_1\left(\omega_{\frac{n}{2}}^k\right) \end{array} \]

(其中第二行用了上面的引理 5)

然而我们还有 $\frac{n}{2}$

又设 $1 \le k < \frac{n}{2}$.

有:

\[\begin{array}{} f(\omega_n^{\frac{n}{2} + k}) & = & f_0((\omega_n^{\frac{n}{2} + k})^2) + \omega_n^{\frac{n}{2} + k}f_1((\omega_n^{\frac{n}{2} + k})^2) \\ & = & f_0(\omega_n^{n + 2k}) + \omega_n^{\frac{n}{2} + k}f_1(\omega_n^{n + 2k})\\ & = & f_0(\omega_n^{2k}) + \omega_n^{\frac{n}{2} + k}f_1(\omega_n^{2k})\\ & = & f_0\left(\omega_{\frac{n}{2}}^k\right) + \omega_n^{\frac{n}{2} + k}f_1\left(\omega_{\frac{n}{2}}^k\right) \\ & = & f_0\left(\omega_{\frac{n}{2}}^k\right) - \omega_n^ kf_1\left(\omega_{\frac{n}{2}}^k\right) \end{array} \]

(其中第二行用到了引理 2, 第三行用到了 $\omega_n^k = \omega_n^{k\bmod n}$, 第四行用到了消去引理, 第五行用到了折半引理)

我们发现, 两个式子都包含了 $f_0\left(\omega_{\frac{n}{2}}^k\right), \omega_n^k f_1\left(\omega_{\frac{n}{2}}^k\right)$ , 原问题被分解成两个更小的子问题, 采用分治即可, 时间复杂度为 $\Theta(n \log n)$.

代码实现

void DFT(complex f[], int len) {
    if(!len) return;
    complex fl[len + 5], fr[len + 5];
    for(int i = 0; i < len; i++) {
        fl[i] = f[i << 1];
        fr[i] = f[i << 1 | 1];
    }
    DFT(fl, len >> 1, flag);
    DFT(fr, len >> 1, flag);
    complex ur(cos(pi / len), sin(pi / len)), tmp(1, 0);
    for(int i = 0; i < len; i++) {
        f[i] = fl[i] + tmp * fr[i];
        f[i + len] = fl[i] - tmp * fr[i];
        tmp *= ur;
    }
}

离散傅里叶逆变换(IDFT)

将插值转化成矩阵的形式

(根据点值表示求系数表示叫做插值)

根据暴力求值的公式, 我们可以推出:

\[\left[\begin{matrix}{} (\omega_n^0)^0 & (\omega_n^0)^1 & (\omega_n^0)^2 & \cdots & {(\omega_n^0)^{n-1}} \\ (\omega_n^1)^0 & (\omega_n^1)^1 & (\omega_n^1)^2 & \cdots & {(\omega_n^1)^{n-1}} \\ (\omega_n^2)^0 & (\omega_n^2)^1 & (\omega_n^2)^2 & \cdots & {(\omega_n^2)^{n-1}} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ (\omega_n^{n - 1})^0 & (\omega_n^{n - 1})^1 & (\omega_n^{n - 1})^2 & \cdots & {(\omega_n^{n - 1})^{n-1}} \\ \end{matrix}\right] \left[\begin{matrix}{} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n-1} \\ \end{matrix}\right] = \left[\begin{matrix}{} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{n-1} \\ \end{matrix}\right] \]

求该矩阵的逆

可以发现这是个范德蒙矩阵

设该矩阵为 $D$.

前人已经告诉我们了逆矩阵:

\[D^{-1} = \frac{1}{n} \times \left[\begin{matrix}{} (\omega_n^0)^0 & (\omega_n^0)^1 & (\omega_n^0)^2 & \cdots & {(\omega_n^0)^{n-1}} \\ (\omega_n^{-1})^0 & (\omega_n^{-1})^1 & (\omega_n^{-1})^2 & \cdots & {(\omega_n^{-1})^{n-1}} \\ (\omega_n^{-2})^0 & (\omega_n^{-2})^1 & (\omega_n^{-2})^2 & \cdots & {(\omega_n^{-2})^{n-1}} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ (\omega_n^{-(n - 1)})^0 & (\omega_n^{-(n - 1)})^1 & (\omega_n^{-(n - 1)})^2 & \cdots & {(\omega_n^{-(n - 1)})^{n-1}} \\ \end{matrix}\right] \]

虽然我不会找逆矩阵, 但是我可以证明一下这个逆矩阵: (我太菜了)

于是设 $E = DD^{-1}$

\[E_{jk} = \sum_{p = 0}^{n - 1} (\omega_n^j)^p \frac{(\omega_n^{-p})^k}{n} \\ = \frac{\sum_{p = 0}^{n - 1} \omega_n^{jp} \omega_n^{-kp}}{n} \\ = \frac{\sum_{p = 0}^{n - 1} \omega_n^{(j-k)p}}{n} \\ \]

当 $j \not =k$

\[E_{jk} = \frac{\omega_n^{(j-k)n} - 1}{n (\omega^{j - k}_n - 1)} = \frac{\omega_n^{0} - 1}{n (\omega^{j - k}_n - 1)} = 0 \]

当 $j = k$

\[E_{jk} = \frac{\sum_{p = 0}^{n - 1} \omega_n^{0}}{n} = \frac{\sum_{p = 0}^{n - 1}1}{n} = 1 \]

所以, $E$

代码实现

void IDFT(complex f[], int len, int flag) {
    if(!len) return;
    complex fl[len + 5], fr[len + 5];
    for(int i = 0; i < len; i++) {
        fl[i] = f[i << 1];
        fr[i] = f[i << 1 | 1];
    }
    IDFT(fl, len >> 1, flag);
    IDFT(fr, len >> 1, flag);
    complex ur(cos(pi / len), -sin(pi / len)), tmp(1, 0);
    for(int i = 0; i < len; i++) {
        f[i] = fl[i] + tmp * fr[i];
        f[i + len] = fl[i] - tmp * fr[i];
        tmp *= ur;
    }
}

简化代码

发现 IDFT 的代码其实和 DFT 相差不大.

发现:

我们需要求:

\[a_k = \frac{\sum_{j = 0}^{n - 1} y_j(\omega_n^{-j})^k}{n} = \frac{\sum_{j = 0}^{n - 1} y_j[(\omega_n^{-1})^{-j}]^k}{n} \]

发现 DFT 为:

\[y_k = \sum_{j = 0} ^ {n - 1} a_j(\omega^j)^k \]

然后发现这个式子和 DFT 要求的式子很像, 只不过是把 $y$

可以发现 $\omega_n^{-1}$ 和 $\omega_n^1$

所以我们就只需要写一个 DFT 就可以了.

递归版代码实现

#include <bits/stdc++.h>

using std::cin;
using std::cout;

typedef double f64;

const f64 pi = acos(-1);
const f64 eps = 1e-5;
const int MAXN = 1e5 + 5;

struct complex {
    f64 real, imag;

    complex(f64 rpt = 0, f64 ipt = 0) {real = rpt; imag = ipt;}

    complex operator + (const complex &x) const {
        complex ans(real + x.real, imag + x.imag);
        return ans;
    }
    complex operator - (const complex &x) const {
        complex ans(real - x.real, imag - x.imag);
        return ans;
    }
    complex operator * (const complex &x) const {
        complex ans(real * x.real - imag * x.imag, real * x.imag + imag * x.real);
        return ans;
    }
    complex operator += (const complex &x) {return *this = *this + x;}
    complex operator -= (const complex &x) {return *this = *this - x;}
    complex operator *= (const complex &x) {return *this = *this * x;}
} fa[MAXN], fb[MAXN];

int n, m;

void FFT(complex f[], int len, int flag) {
    if(!len) return;
    complex fl[len + 5], fr[len + 5];
    for(int i = 0; i < len; i++) {
        fl[i] = f[i << 1];
        fr[i] = f[i << 1 | 1];
    }
    FFT(fl, len >> 1, flag);
    FFT(fr, len >> 1, flag);
    complex ur(cos(pi / len), sin(pi / len) * flag), tmp(1, 0);
    for(int i = 0; i < len; i++) {
        f[i] = fl[i] + tmp * fr[i];
        f[i + len] = fl[i] - tmp * fr[i];
        tmp *= ur;
    }
}

int main() {
    scanf("%d%d", &n, &m);
    for(int i = 0; i <= n; i++)
    	scanf("%lf", &fa[i].real);
    for(int i = 0; i <= m; i++)
    	scanf("%lf", &fb[i].real);
    int len = 1;
    for(n += m; len <= n; len <<= 1);
    FFT(fa, len >> 1, 1); FFT(fb, len >> 1, 1);
    for(int i = 0; i < len; i++)
    	fa[i] *= fb[i];
    FFT(fa, len >> 1, -1);
    for(int i = 0; i <= n; i++) printf("%.0lf ", fabs(fa[i].real) / len);
    return 0;
}

蝴蝶算法

发现规律

我们发现, 直接用递归写, 常数非常大.

$ 0, 1, 2, 3, 4, 5, 6, 7 $

$ 0, 2, 4, 6|1, 3, 5, 7 $

$ 0, 4|2, 6|1, 5|3, 7 $

$ 0|4|2|6|1|5|3|7 $

发现最后的数的二进制表示是:
$000, 100, 110, 001, 101, 011, 111$

是原来的序号的二进制反序.

所以先求出最后的数组, 在原数组上迭代就可以了.

迭代实现

void FFT(complex f[], int len, int flag) {
    for(int i = 0; i < len; i++)
        if(i < rev[i]) std::swap(f[i], f[rev[i]]);
    for(int i = 1; i < len; i <<= 1) {
        complex ur(cos(PI / i), flag * sin(PI / i));
        for(int j = 0; j < len; j += (i << 1)) {
            complex tmp(1, 0);
            for(int k = 0; k < i; k++, tmp *= ur) {
                complex fr = f[i + j + k], fl = f[j + k];
                f[j + k] = fl + tmp * fr;
                f[i + j + k] = fl - tmp * fr;
            }
        }
    }
}

代码实现

#include <bits/stdc++.h>

using std::cin;
using std::cout;

typedef double f64;

const f64 PI = acos(-1);
const int MAXN = 4e6 + 5;

struct complex {
    f64 real, imag;
    complex(f64 rpt = 0, f64 ipt = 0) {real = rpt; imag = ipt;}
    complex operator + (const complex &x) const {
        return complex(real + x.real, imag + x.imag);
    }
    complex operator - (const complex &x) const {
        return complex(real - x.real, imag - x.imag);
    }
    complex operator * (const complex &x) const {
        return complex(real * x.real - imag * x.imag, real * x.imag + imag * x.real);
    }
    complex operator += (const complex &x) {return *this = *this + x;}
    complex operator -= (const complex &x) {return *this = *this - x;}
    complex operator *= (const complex &x) {return *this = *this * x;}
} fa[MAXN], fb[MAXN];

int n, m, rev[MAXN];

void FFT(complex f[], int len, int flag) {
    for(int i = 0; i < len; i++)
        if(i < rev[i]) std::swap(f[i], f[rev[i]]);
    for(int i = 1; i < len; i <<= 1) {
        complex ur(cos(PI / i), flag * sin(PI / i));
        for(int j = 0; j < len; j += (i << 1)) {
            complex tmp(1, 0);
            for(int k = 0; k < i; k++, tmp *= ur) {
                complex fr = f[i + j + k], fl = f[j + k];
                f[j + k] = fl + tmp * fr;
                f[i + j + k] = fl - tmp * fr;
            }
        }
    }
}

int read() {
    int x = 0; char c = getchar();
    while(c < '0' || c > '9') c = getchar();
    while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
    return x;
}

void write(int x) {
    if(x / 10) write(x / 10);
    putchar(x % 10 + '0');
}

int main() {
    n = read(); m = read();
    for(int i = 0; i <= n; i++) fa[i].real = (f64)read();
    for(int i = 0; i <= m; i++) fb[i].real = (f64)read();
    int len = 1, maxBit = 0;
    for(n += m; len <= n; len <<= 1, maxBit++);
    for(int i = 0; i < len; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (maxBit - 1));
    FFT(fa, len, 1); FFT(fb, len, 1);
    for(int i = 0; i < len; i++)
        fa[i] *= fb[i];
    FFT(fa, len, -1);
    for(int i = 0; i <= n; i++)
        write((int)(fabs(fa[i].real) / len + 0.5)), putchar(' ');
    return 0;
}

快速数论变换

原根与阶

阶

最小的 $ t $ 使得 $ a^t \equiv 1 \pmod p$, 则 $\delta_p (a) = t$, 称 $ t $ 为 $a$

原根

若 $ \delta_p(a) = \varphi(p) $, 则 $a$ 是 $p$

求原根

原根通常很小, 模数通常不会超过 $10^9$, 暴力枚举即可.

常见模数 $998244353 = 7 \times 17\times 2^{13} + 1$ 的原根为 $3$.

快速数论变换

根据相关定理, 可以证明, 模 $p$ 意义下, $e^{\frac{-2\pi i}{n}} \equiv g^{\frac{p - 1}{n}} \pmod p$.

怎么证明? 我也不会.

代码实现

#include <bits/stdc++.h>

typedef long long ll;

const ll mod = 7 * 17 * (1 << 23) + 1;
const ll g = 3, gInv = 332748118;
const ll MAXN = 4e6 + 5;

ll n, m, rev[MAXN], fa[MAXN], fb[MAXN];

ll qPow(ll a, ll b) {
	ll ans = 1, base = a;
	while(b) {
		if(b & 1) ans = ans * base % mod;
		base = base * base % mod;
		b >>= 1;
	}
	return ans;
}

void FFT(ll f[], int len, int flag) {
    for(int i = 0; i < len; i++)
        if(i < rev[i]) std::swap(f[i], f[rev[i]]);
    for(int i = 1; i < len; i <<= 1) {
        ll ur = qPow(flag ? g : gInv, (mod - 1) / (i << 1));
        for(int j = 0; j < len; j += (i << 1)) {
            ll tmp = 1;
            for(int k = 0; k < i; k++, tmp = tmp * ur % mod) {
                ll fr = f[i + j + k], fl = f[j + k];
                f[j + k] = (fl + tmp * fr % mod) % mod;
                f[i + j + k] = (fl - tmp * fr % mod + mod) % mod;
            }
        }
    }
}

ll read() {
    ll x = 0; char c = getchar();
    while(c < '0' || c > '9') c = getchar();
    while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
    return x;
}

void write(ll x) {
    if(x / 10) write(x / 10);
    putchar(x % 10 + '0');
}

int main() {
    n = read(); m = read();
    for(int i = 0; i <= n; i++) fa[i] = read();
    for(int i = 0; i <= m; i++) fb[i] = read();
    int len = 1, maxBit = 0;
    for(n += m; len <= n; len <<= 1, maxBit++);
    for(int i = 0; i < len; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (maxBit - 1));
    FFT(fa, len, 1); FFT(fb, len, 1);
    for(int i = 0; i < len; i++)
        fa[i] = fa[i] * fb[i] % mod;
    FFT(fa, len, 0);
    ll lenInv = qPow(len, mod - 2);
    for(int i = 0; i <= n; i++)
        write(fa[i] * lenInv % mod), putchar(' ');
    return 0;
}

例题

题目

输入两个数 $a, b$ 求 $a\times b$. $0\le a,b \le 10^{100000}$.

分析

发现普通高精乘的复杂度为 $\Theta((\lg n)^2)$

可以把两个数看成两个多项式相乘, 最后处理一下进位即可.

本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。

上一篇：在Python中LDA的 passes python lda模型

下一篇：java跟linux怎么交互 java和linux

提问和评论都可以，用心的回复会被更多人看到评论

发布评论

相关文章

官方博客	全部文章	热门标签	班级博客
了解我们	网站地图	意见反馈

鸿蒙开发者社区	51CTO学堂
51CTO	软考资讯