[学习笔记]MTT

\[\color{red}{\text{校长者，真神人也，左马桶，右永神，会执利笔破邪炁，何人当之？}} \\ \begin{array}{|} \hline \color{pink}{\text{The principal is really a god}} \\ \color{pink}{\text{with a closestool on the left and Yongshen on the right}} \\ \color{pink}{\text{holding a sharp pen to pierce the truth}} \\ \color{pink}{\text{Who can resist him? }} \\ \hline \end{array} \\ \begin{array}{|} \hline \color{green}{\text{校長は本当に神であり、左側にトイレ、右側にヨンシェンがあり}} \\ \color{green}{\text{鋭いペンを持って真実を突き刺している。誰が彼に抵抗できるだろうか？ }} \\ \hline \end{array} \\ \begin{array}{|} \hline \color{lightblue}{\text{Le principal est vraiment un dieu}} \\ \color{lightblue}{\text{avec des toilettes à gauche et Yongshen à droite}} \\ \color{lightblue}{\text{tenant un stylo pointu pour percer la vérité}} \\ \color{lightblue}{\text{Qui peut lui résister ? }} \\ \hline \end{array} \\ \begin{array}{|} \hline \color{purple}{\text{Der Direktor ist wirklich ein Gott}} \\ \color{purple}{\text{mit einer Toilette links und Yongshen rechts}} \\ \color{purple}{\text{der einen spitzen Stift hält}} \\ \color{purple}{\text{um die Wahrheit zu durchdringen.}} \\ \color{purple}{\text{Wer kann ihm widerstehen? }} \\ \hline \end{array} \\ \begin{array}{|} \hline \color{cyan}{\text{Principalis deus est, Yongshen a dextris cum latrina}} \\ \color{cyan}{\text{acuto stylo ad perforandum veritatem: quis resistet ei? }} \\ \hline \end{array} \\ \color{red}{\text{对曰：“无人，狗欲当之，还请赐教！”}} \\ \newcommand\brak[1]{\left({#1}\right)} \newcommand\Brak[1]{\left\{{#1}\right\}} \newcommand\d[0]{\text{d}} \newcommand\string[2]{\genfrac{\{}{\}}{0pt}{}{#1}{#2}} \newcommand\down[2]{{#1}^{\underline{#2}}} \newcommand\ddiv[2]{\left\lfloor\frac{#1}{#2}\right\rfloor} \newcommand\udiv[2]{\left\lceil\frac{#1}{#2}\right\rceil} \newcommand\lcm[0]{\operatorname{lcm}} \newcommand\set[1]{\left\{{#1}\right\}} \newcommand\ceil[1]{\left\lceil{#1}\right\rceil} \newcommand\floor[1]{\left\lfloor{#1}\right\rfloor} \]

---

  简单记录一下 \(\rm MTT\) 的推导和使用。

7 次 DFT 的 MTT

  最朴素，并且最容易理解。

  令 \(\text{bas} = \sqrt p\)，那么，若计算 \(\Lambda=F\circ G\)，我们可以考虑将每一项拆成两个部分：

\[F(x)=\sum_{i}f(i)x^i=\sum_{i}(c_1(i)bas+r_1(i))x^i \\ G(x) = \sum_ig(i)x^i=\sum_i(c_2(i)bas+r_2(i))x^i \]

  可得

\[\begin{aligned} (F\cdot G)(x)&=\brak{\sum_{i}(c_1(i)bas+r_1(i))x^i}\brak{\sum_i(c_2(i)bas+r_2(i))x^i} \\ &=\sum_{i}\sum_j\Big(c_1(i)c_2(j)\text{bas}^2+(c_1(i)r_2(j)+r_1(i)c_2(j))\text{bas}+r_1(i)r_2(j)\Big)x^{i+j} \end{aligned} \]

  根据 \(\text{bas}\) 的次数区分，我们需要做 \(4\) 次 \(\rm DFT\)，针对 \(c_1,c_2,r_1,r_2\)，然后做三次 \(\rm IDFT\)，针对 \(c_1\circ c_2,\;c_1\circ r_2+c_2\circ r_1,\;r_1\circ r_2\)，常数比较大，但是易于推导与理解。

5 次 DFT 的 MTT

  事实上和 4 次只有一步之遥了。

  第一种方法为什么常数那么大？因为每一次变换的时候，我们将 \((h_i,0)\)（其中 \(h\) 为任意一个函数） 作为一个复数进行变换，每次的变换，初始的虚部均为 \(0\)，这个无用却不得不计算的部分带给我们很大的开销，考虑能否将这个东西用起来。

  我们还是考虑计算 \(\Lambda=(c_1\text{bas}+r_1)(x)+(c_2\text{bas}+r_2)(x)\)，显然，我们其实想要求 \(c_1\circ c_2,\;c_1\circ r_2,\;r_1\circ c_2\;,r_1\circ r_2\) 这四个卷积，并不需要将他们单独求出来，可以构造一个结构使其出现这四个玩意，考虑到 \((a+bi)(c+di)=ac-bd+(ad+bc)i\)，并且，如果在 \(b\) 之前添一个负号，同理可得 \(ac+bd+(ad-bc)i\)，这两个东西上下相加，实部和虚部均只剩下一个，上下相减也只会剩下一个，于是我们就可以求出这四个玩意。具体地，令

\[P(x)=c_1(x)+ir_1(x) \\ Q(x)=c_2(x)+ir_2(x) \]

  则

\[\begin{align} P(x)\circ Q(x)=(c_1\circ c_2-r_1\circ r_2)(x)+i(c_1\circ r_2+r_1\circ c_2)(x)\tag{1} \\ \bar P(x)\circ Q(x)=(c_1\circ c_2+r_1\circ r_2)(x)+i(c_1\circ r_2-r_1\circ c_2)(x)\tag{2} \end{align} \]

  于是

\[(1)+(2)=2(c_1\circ c_2)(x)+2i(c_1\circ r_2)(x) \\ (1)-(2)=-2(r_1\circ r_2)(x)+2i(r_1\circ c_2)(x) \]

  这四个东西分别单独出现在实部和虚部，显然可以直接求的。

  接下来统计一下我们需要的变换次数：三次正变换，针对 \(P,\bar P,Q\)，两次逆变换，针对 \(P(x)\circ Q(x)\) 和 \(\bar P(x)\circ Q(x)\). 常数大大减小了！

  思路确实简单，不过难点在想到如此构造。

typedef complex<long double> cplx;

const int maxn = 1e5;
const long double Pi = acos(-1.0);
int mod;

namespace __poly {

int rev[maxn * 6 + 5], n;

inline void prepare(int len) {
    for (n = 1; n < len; n <<= 1);
    for (int i = 0; i < n; ++i)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1)? n >> 1: 0);
    return ;
}
inline void fft(vector<cplx>& f, const int opt) {
    f.resize(n);
    for (int i = 0; i < n; ++i) if (i < rev[i])
        swap(f[i], f[rev[i]]);
    for (int p = 2; p <= n; p <<= 1) {
        int len = p >> 1;
        cplx w(cos(Pi / len), opt * sin(Pi / len));
        for (int k = 0; k < n; k += p) {
            cplx buf(1, 0), tmp;
            for (int i = k; i < k + len; ++i, buf *= w) {
                tmp = f[i + len] * buf;
                f[i + len] = f[i] - tmp;
                f[i] = f[i] + tmp;
            }
        }
    }
    if (opt == -1) for (int i = 0; i < n; ++i) f[i] /= n;
    return ;
}

vector<cplx> P, _P, Q, mul1, mul2;
const int bas = 32000;
inline void mtt(vector<int>& f, const vector<int>& g) {
    P.resize(f.size()), _P.resize(f.size()), Q.resize(g.size());
    for (int i = 0; i < f.size(); ++i) {
        P[i] = cplx(f[i] / bas, f[i] % bas);
        _P[i] = conj(P[i]);
    }
    for (int i = 0; i < g.size(); ++i)
        Q[i] = cplx(g[i] / bas, g[i] % bas);
    prepare(f.size() + g.size());
    fft(P, 1), fft(_P, 1), fft(Q, 1);
    transform(P.begin(), P.end(), Q.begin(), P.begin(), multiplies<cplx>());
    transform(_P.begin(), _P.end(), Q.begin(), _P.begin(), multiplies<cplx>());
    fft(P, -1), fft(_P, -1), f.resize(n);
    for (int i = 0; i < n; ++i) {
        int c1c2 = (ll)((P[i].real() + _P[i].real()) * 0.5 + 0.5) % mod;
        int c1r2 = (ll)((P[i].imag() + _P[i].imag()) * 0.5 + 0.5) % mod;
        int r1c2 = (ll)((P[i].imag() - _P[i].imag()) * 0.5 + 0.5) % mod;
        int r1r2 = (ll)((_P[i].real() - P[i].real()) * 0.5 + 0.5) % mod;
        f[i] = 1ll * c1c2 * bas % mod * bas % mod;
        f[i] = (f[i] + 1ll * (c1r2 + r1c2) * bas % mod) % mod;
        f[i] = (f[i] + r1r2) % mod;
        f[i] = (f[i] + mod) % mod;
    }
    return ;
}

} // using namespace __poly;

4 次 DFT 的 MTT

  还可以更快！事实上，我们可以直接从 \(P(x)\) 得到 \(\bar P(x)\)，先把 \(P(x)\) 变换出来，对于其共轭，假设 \(p_k=a+bi,\bar p_k=a-bi\)，不难发现 \(\bar p_k\) 经过变换之后就是 \(p_{n-k}\) 经过变换得到的东西，因为共轭复数与原来的复数关于实轴对称，不过一个特例的是 \(\bar p_0=p_0\)，特殊处理一下就好了。

typedef complex<long double> cplx;

const int maxn = 1e5;
const long double Pi = acos(-1.0);
int mod;

namespace __poly {

int rev[maxn * 6 + 5], n;

inline void prepare(int len) {
    for (n = 1; n < len; n <<= 1);
    for (int i = 0; i < n; ++i)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1)? n >> 1: 0);
    return ;
}
inline void fft(vector<cplx>& f, const int opt) {
    f.resize(n);
    for (int i = 0; i < n; ++i) if (i < rev[i])
        swap(f[i], f[rev[i]]);
    for (int p = 2; p <= n; p <<= 1) {
        int len = p >> 1;
        cplx w(cos(Pi / len), opt * sin(Pi / len));
        for (int k = 0; k < n; k += p) {
            cplx buf(1, 0), tmp;
            for (int i = k; i < k + len; ++i, buf *= w) {
                tmp = f[i + len] * buf;
                f[i + len] = f[i] - tmp;
                f[i] = f[i] + tmp;
            }
        }
    }
    if (opt == -1) for (int i = 0; i < n; ++i) f[i] /= n;
    return ;
}

vector<cplx> P, _P, Q, mul1, mul2;
const int bas = 32000;
inline void mtt(vector<int>& f, const vector<int>& g) {
    P.resize(f.size()), Q.resize(g.size());
    for (int i = 0; i < f.size(); ++i)
        P[i] = cplx(f[i] / bas, f[i] % bas);
    for (int i = 0; i < g.size(); ++i)
        Q[i] = cplx(g[i] / bas, g[i] % bas);
    prepare(f.size() + g.size());
    fft(P, 1), fft(Q, 1);
    _P.resize(n);
    for (int i = 0; i < n; ++i) _P[i] = conj(P[i? n - i: 0]);
    transform(P.begin(), P.end(), Q.begin(), P.begin(), multiplies<cplx>());
    transform(_P.begin(), _P.end(), Q.begin(), _P.begin(), multiplies<cplx>());
    fft(P, -1), fft(_P, -1), f.resize(n);
    for (int i = 0; i < n; ++i) {
        int c1c2 = (ll)((P[i].real() + _P[i].real()) * 0.5 + 0.5) % mod;
        int c1r2 = (ll)((P[i].imag() + _P[i].imag()) * 0.5 + 0.5) % mod;
        int r1c2 = (ll)((P[i].imag() - _P[i].imag()) * 0.5 + 0.5) % mod;
        int r1r2 = (ll)((_P[i].real() - P[i].real()) * 0.5 + 0.5) % mod;
        f[i] = 1ll * c1c2 * bas % mod * bas % mod;
        f[i] = (f[i] + 1ll * (c1r2 + r1c2) * bas % mod) % mod;
        f[i] = (f[i] + r1r2) % mod;
        f[i] = (f[i] + mod) % mod;
    }
    return ;
}

} // using namespace __poly;