任意模数多项式乘法MTT（可拆系数FFT）详解

更好的阅读体验

任意模数多项式乘法

前言：
在教练讲的时候脑子并不清醒，所以没听懂。后来自己看博客学会了，但目前只学了一种方法：可拆系数FFT。为了方便日后复习，决定先写下这个的笔记，关于三模数NTT下次再补。

建议：准备好演算纸和笔，本篇含有大量推算部分。

为什么不直接使用NTT/FFT

此处的模数是任意的，所以我们使用NTT时，有局限性。只有当模数满足下列情况时才可使用NTT。

模数为$P$，$P=r\times 2^k+1$，其中$k$足够大，满足$2^k\ge N$，其中$N$为多项式扩充后的项数，多项式乘法都需要将项数扩充到$2$的幂。

如果不满足上述条件，就不能直接使用NTT。

那为什么不使用FFT带模运算？
首先不考虑模数的情况下，多项式结果的系数不能超过$double$能表示的精度（一般在$10^{16}$）。超过后，$double$所表示的结果将不再精确。

那怎么办？目前可以给出两种解决方法

1. 可拆系数FFT
2. 三模数NTT<s>（可是我还没学懂QAQ）</s>

所以，向可拆系数FFT进发！

可拆系数FFT

怎么用可拆系数FFT

首先我们还是先假设一个情景：

求：$c_1$与$c_2$的卷积

我们可以把一个数拆成$a\times 2^{15}+b$的形式（不一定是$2^{15}$，大概在$\sqrt{值域}$内）。

$$\begin{gathered} c_1=a_1\times 2^{15}+a_2 \\ {c}_2=b_1\times 2^{15}+b_2 \end{gathered}$$

那这俩的积为

$$\begin{array}{rl}{c}_1\times c_2& =(a_1\times 2^{15}+a_2)\times (b_1\times 2^{15}+b_2)\\ & =a_1{b}_1\times 2^{30}+(a_1{b}_2+a_2{b}_1)\times 2^{15}+a_2{b}_2\end{array}$$

好耶！那我们可以直接做4次FFT（$a_1{b}_1$、$a_1{b}_2$、$a_2{b}_1$、$a_2{b}_2$）！
<s>然后你发现正逆变换总共做了8次常数爆炸然后炸了</s>

所以，我们需要优化！

优化可拆系数FFT

注意：我们这里的优化会用到复数，你可能会害怕得逃走，但是你无需害怕，因为<s>（我也是这样过来的）</s>我会简单地讲一讲。

小资料（可能不太学术规范）：
啥是复数？

其实复数，是一种含实部和虚部的一种数。我们知道$i=\sqrt{?1}$，$i$就是虚数。那我们以实部建立$x$轴、虚部建立$y$轴。

那我们假设在这个平面直角坐标系上有一个点$A(2,3)$，那这个点的复数表示为$2+3i$。

那你就基本知道什么是复数了，让我们学一下基本运算吧。（这个自己记一记好了）

我们可以利用一下复数的乘法运算：

$$\begin{gathered} (a+bi)(c+di)=(ac?bd)+(ad+bc)i \\ (a?bi)(c+di)=(ac+bd)+(ad?bc)i \end{gathered}$$

现在令$P_1=a_1+a_{2}i$、$P_2=a_1?a_{2}i$、$Q=b_1+b_{2}i$。
计算：

$$\begin{array}{rl}{P}_{1}Q+P_{2}Q& =(a_1{b}_1?a_2{b}_2)+(a_1{b}_2+a_2{b}_1)i+(a_1{b}_1+a_2{b}_2)+(a_1{b}_2?a_2{b}_1)i\\ & =2(a_1{b}_1+a_1{b}_{2}i)\end{array}$$

如果我们将上式除以$2$，那我们可以得到$a_1{b}_1,a_1{b}_2$（分别通过“实部相加/2”、“虚部相加/2”可得）。

我们把得到的$a_1{b}_1$代入$P_{2}Q$的实部，可得$a_2{b}_2$。
类似地，将$a_1{b}_2$代入$P_{1}Q$的虚部，可得$a_2{b}_1$。

（注：这里的运算请自己用演算纸推一下）

我们就可以带回$c_1\times c_2$了。

那我们的任务就完成啦！！

代码实现

注：此处的FFT部分与FFT版题的部分是一模一样的，可参照本人之前所写的FFT笔记。

板子：P4245 【模板】任意模数多项式乘法

code:

#include<bits/stdc++.h>
using namespace std;

#define ll long long
#define rp(i,o,p) for(ll i=o;i<=p;++i)
#define pr(i,o,p) for(ll i=o;i>=p;--i)
#define cp complex<long double>

const ll MAXN=1e5+5,P=1e9+7;
const ll lim=32000; // lim = sqrt(值域) <- 1e9
const long double PI=acos(-1.0);

cp p1[MAXN<<2],p2[MAXN<<2],q[MAXN<<2];
ll n,m,limit;
cp p[MAXN<<2],omg[MAXN<<2],inv[MAXN<<2];

void init() {
    for (ll i = 0; i < limit; ++i) {
        omg[i] = cp(cos(2 * PI * i / limit), sin(2 * PI * i / limit));
        inv[i] = conj(omg[i]);
    }
}

void fft(cp *a, cp *omg) {
    for (ll i = 0, j = 0; i < limit; ++i) {
        if (i > j)
            swap(a[i], a[j]);
        for (ll l = limit >> 1; (j ^= l) < l; l >>= 1)
            ;
    }
    for (ll l = 2; l <= limit; l <<= 1) {
        ll m = l >> 1;
        for (cp *p = a; p != a + limit; p += l) {
            rp(i, 0, m - 1) {
                cp t = omg[limit / l * i] * p[i + m];
                p[i + m] = p[i] - t;
                p[i] += t;
            }
        }
    }
}

int main()
{
    scanf("%lld%lld",&n,&m);
    rp(i,0,n)
    {
        ll x;
        scanf("%lld",&x);
        p1[i]=cp(x/lim,x%lim);
        p2[i]=cp(x/lim,-x%lim);
    }
    rp(i,0,m)
    {
        ll x;
        scanf("%lld",&x);
        q[i]=cp(x/lim,x%lim);
    }

    limit=1;
    while(limit<=n+m) limit<<=1;

    init();
    fft(p1,omg),fft(p2,omg),fft(q,omg);
  
    rp(i,0,limit-1)
    {
    	long double xx,xy;
  
        xx=p1[i].real()*q[i].real(),yy=p1[i].imag()*q[i].imag();
        xy=p1[i].real()*q[i].imag(),yx=p1[i].imag()*q[i].real();
        p1[i]=cp(xx/limit-yy/limit,xy/limit+yx/limit);
      
        xx=p2[i].real()*q[i].real(),yy=p2[i].imag()*q[i].imag();
        xy=p2[i].real()*q[i].imag(),yx=p2[i].imag()*q[i].real();
        p2[i]=cp(xx/limit-yy/limit,xy/limit+yx/limit);
    }
    /*
    	上面的循环等价于这两行循环，但是因为c++中complex模板为const类型，
    	单独的real()或imag()值不可以直接修改，只可以两个同时赋值来修改，
    	所以上面的循环还用到了复数的除法运算，自己看看吧
      	rp(i,0,limit-1) q[i].real()/=limit,q[i].imag()/=limit;
      	rp(i,0,limit-1) p1[i]*=q[i],p2[i]*=q[i];
    */
    fft(p1,inv),fft(p2,inv);

    rp(i,0,n+m)
    {
        ll a1b1=(ll)((p1[i].real()+p2[i].real())/2+0.5)%P;
        ll a1b2=(ll)((p1[i].imag()+p2[i].imag())/2+0.5)%P;
        ll a2b1=(ll)((p1[i].imag()+0.5)-a1b2)%P;
        ll a2b2=(ll)((p2[i].real()+0.5)-a1b1)%P;
        ll ans=(a1b1*lim*lim+(a1b2+a2b1)*lim+a2b2)%P;

        ans=(ans%P+P)%P;
        printf("%lld ",ans);
    }

    return 0;
}

后记

三模数NTT的内容会尽快补上的。