投稿
  1. 高性能架构探索首页
  2. C/C++

多项式乘法

admin • • C/C++ • 阅读 138

\(part\, 0\) 关于优化多项式乘法

很自然的，一个n次多项式乘一个m阶多项式，其最终的结果是n+m阶，在空间上来说是线性的，然而为什么我们就会必须要一个\(O(nm)\)的二次级别的时间复杂度来进行运算呢，能不能把这个多项式乘法优化到线性级别或者是接近线性级别的呢

虽然我有心无力，但在这个神仙辈出的时代，定有人能顺利解答这个问题。很自然地，在多项式乘法广泛应用的今天，\(fft\&ntt\)横空出世，堪称将数学应用于信息学的一大经典

\(part\,1\) 点值表示法

命题：对于一个任意的n阶多项式\(f(x)=\sum\limits_{i=0}^nA_ix^i\)，必定有n+1个互不相同的点\((a_i,f(a_i))\)使这个多项式唯一确定

\(p.f.\)

\(\text{显然我们会有：}\)

\[\begin{pmatrix}1&a_0&\cdots&a_0^n\\1&a_1&\cdots&a_1^n\\1&a_2&\cdots&a_2^n\\\vdots&\vdots&\ddots&\vdots\\1&a_n&\cdots&a_n^n\end{pmatrix}\begin{pmatrix}A_0\\A_1\\A_2\\\vdots\\A_n\end{pmatrix}=\begin{pmatrix}f(a_0)\\f(a_1)\\f(a_2)\\f(a_3)\\\vdots\\f(a_n)\end{pmatrix}
\]

\(\text{不妨我们设}\)

\[\Lambda=\begin{pmatrix}1&a_0&\cdots&a_0^n\\1&a_1&\cdots&a_1^n\\1&a_2&\cdots&a_2^n\\\vdots&\vdots&\ddots&\vdots\\1&a_n&\cdots&a_n^n\end{pmatrix}，\Alpha=\begin{pmatrix}A_0\\A_1\\A_2\\\vdots\\A_n\end{pmatrix}，\alpha=\begin{pmatrix}f(a_0)\\f(a_1)\\f(a_2)\\f(a_3)\\\vdots\\f(a_n)\end{pmatrix}
\]

\(\text{发现}|\Lambda|恰好是一个Vandermonde\text{行列式，于是}\)

\[|\Lambda|=\prod\limits_{0\leqslant i<j\leqslant n}(a_j-a_i)\neq0
\]

\(\text{那么}\Lambda^{-1}\text{存在，于是}\)

\[\Alpha=\Lambda^{-1}\alpha
\]

\(\text{系数向量}\Alpha\text{被唯一确定}\)，证毕

关于为什么要介绍点值表示法，是因为大家会发现如果我们用点值表示法来代表一个多项式的话，这样的转移是线性的

\(e.g.\)(A(x)为n次多项式，B(x)为m次多项式)

\[\text{若}C(x)=A(x)B(x)，\text{那么}C(x_i)=A(x_i)B(x_i)
\]

对应的我们遍取\(1\leqslant i\leqslant n+m\)，那么我们就可以在\(O(n+m)\)的时间内由A与B的点值表达式推到C的点值表达式

\(attention:\) 这种思想还可以进一步得到更广泛的应用，一定要重点了解

\(part~2~fft\)

这里我们的想法是引入单位根，而为什么这么想？我也不知道。

先把单位根放在一边，我们考虑朴素的从系数表达式往点值转换，如果我们的多项式自变量是随意决定的n个数的话，那么我们每次计算\(f(a_i)\)的复杂度都是\(O(n)\)，n次之和即为\(O(n^2)\)，显然已经不符合我们的优化目的，因此我们需要一个更加优秀的方法能一次性算出某n个数对应的多项式数值，因此这n个数如果不互相有点特殊关系那可是很难达到我们的目的

在这一点上，\(Fourier\)想到了一个很好的解决办法——单位根（当然我更愿意认为这是\(Euler\)的功劳），而一切思想的基础都建立在我们的\(Euler\)公式上：

\[e^{i\theta}=\cos\theta+i\sin\theta
\]

其中最特殊的情况也就是\(\theta=\pi\)时我们可以得到我们最熟悉的数学中最美丽的公式：

\[e^{i\pi}=-1\Leftrightarrow e^{i\pi}+1=0
\]

我们现在定义单位根\(w_n^k=\cos \dfrac{2k}{n}\pi+i\sin\dfrac{2k}{n}\pi\)，那么\(w_n^k=e^{i\frac{2k}{n}\pi}\)，顺其自然地我们将得到一些有趣的性质，进而我们会发现求值其实是一件很简单的事

首先假设我们有多项式\(f(x)\)，我们现在的目的是求出n个单位根对应的点值\(f(w_n^k)\)，为方便起见，我们假设n恰好为2的正整数次幂，同时我们假设\(f(x)\)的项数（注意不是次数）也是对应的2的正整数次幂，也就是n（简单的道理，不足n的项数我们可以以\(0x^{m}\)的形式补齐），那么我们有：

\[f(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}
\]

考虑将\(f(x)\)奇偶项交错分开，那么

\[f(x)=a_0+a_2x^2+\cdots+a_{n-2}x^{n-2}+a_1x+a_3x^3+\cdots+a_{n-1}x^{n-1}
\]

假设我们有

\(\begin{cases}f_1(x)=a_0+a_2x+\cdots+a_{n-2}x^{\frac{n}{2}-1}\\f_2(x)=a_1+a_3x+\cdots+a_{n-1}x^{\frac{n}{2}-1}\end{cases}\)

显然\(f_1(x)\)与\(f_2(x)\)都是\(\dfrac{n}{2}\)项，于是我们嗅到了一丝分治的气息，事实上我们可以得到

\[f(x)=f_1(x^2)+xf_2(x^2)
\]

接下来我们可以来推式子了

\[\begin{aligned}f(w_n^k)=f(e^{i\frac{2k}{n}\pi})&=f_1(e^{i\frac{4k}{n}\pi})+e^{i\frac{2k}{n}\pi}f_2(e^{i\frac{4k}{n}\pi})\\&=f_1(e^{i\frac{2k}{n/2}\pi})+e^{i\frac{2k}{n}\pi}f_2(e^{i\frac{2k}{n/2}\pi})\\&=f_1(w^k_{\frac{n}{2}})+e^{i\frac{2k}{n}\pi}f_2(w_{\frac{n}{2}}^{k})&\cdots&\cdots① \end{aligned}
\]

可见我们比较轻松的把括号内的单位根的下标给降为了原来的二分之一，那么剩下的就是对\(e^{i\frac{2k}{n}\pi}\)的处理问题以及把\(k\)的取值范围给对应减半的问题，其实这两个问题很好解决

首先考虑n次单位根\(w_n^k\)中k的取值范围，显然对k来说每n是一个周期，因此很明显对应的k任取的范围从简应该是\(k\in [0,n)\bigcap\Z\)，总共n个取值. 那么我们对应的要把k降级求其\(k\in[0,\dfrac{n}{2})\bigcap\Z\)任取的对应单位根.

对此，我们考虑将\(k\in[\dfrac{n}{2},n)\bigcap\Z\)中的任取转化为\(k\in[0,\dfrac{n}{2})\bigcap\Z\)中的情况，我们假设\(k\in[0,\dfrac{n}{2})\bigcap\Z\)，并将\(\dfrac{n}{2}+k\)代入①式，那么就相当于\(k\in[\dfrac{n}{2},n)\bigcap\Z\)的情况：

\[\begin{aligned}f(w_n^{\frac{n}{2}+k})&=f_1(e^{i\frac{n+2k}{n/2}\pi})+e^{i\frac{n+2k}{n}\pi}f_2(e^{i\frac{n+2k}{n/2}\pi})\\&=f_1(e^{i\frac{2k}{n/2}\pi}e^{2i\pi})+e^{i\frac{2k}{n}\pi}e^{i\pi}f_2(e^{i\frac{2k}{n/2}\pi}e^{2i\pi})\\&=f_1(e^{i\frac{2k}{n/2}\pi})-e^{i\frac{2k}{n}\pi}f_2(e^{i\frac{2k}{n/2}\pi})\\&=f_1(w_{\frac{n}{2}}^k)-w_{n}^kf_2(w_{\frac{n}{2}}^k) \end{aligned}
\]

而且我们又有

\[f_1(w_n^k)=f_1(w_{\frac{n}{2}}^k)+w_{n}^kf_2(w_{\frac{n}{2}}^k)
\]

那么整个式子的\(n\)个点值便可以通过求出\(f_1\)与\(f_2\)对应的\(n\)个点值来求出

如此一来我们发现只需要求出对应的\(k\in[0,\dfrac{n}{2})\bigcap\Z\)中任取的对应单位根，即可递归调用来分治解决，显然这个时候我们的时间复杂度已经降到一个可以接受的范围

时间复杂度会满足一个我们十分熟悉的递推式：

\[T(n)=2T(\frac{n}{2})+O(n)
\]

通过数学归纳法我们便可得到我们最终求点值的时间复杂度：

\[T(n)=O(n\log_2n)
\]

那么在系数表达式转化到点值表达式这一步我们就达到目的了，将原本朴素的\(O(n^2)\)带到了\(O(n\log_2n)\)

以上想法可以很快地将我们n个单位根对应的点值给求出，假设我们原来的系数数组为\(f_i\)，转换为点值后替换原来的系数数组，那么递归形式的模板如下：

const int xx=(1<<20)+1;
const double pi=acos(-1.0);
struct Complex
{
    double a,b;
    inline Complex operator +(const Complex x){return (Complex){a+x.a,b+x.b};}
    inline Complex operator -(const Complex x){return (Complex){a-x.a,b-x.b};}
    inline Complex operator *(const Complex x){return (Complex){a*x.a-b*x.b,a*x.b+b*x.a};}
};
inline void fft(Complex *f,int n,const int op)
{
    if(n==1) return;
    n>>=1;
    Complex f1[n],f2[n];
    for(register int i=0;i<n;++i)
    {
        f1[i]=f[i<<1];
        f2[i]=f[(i<<1)+1];
    }
    fft(f1,n,op),fft(f2,n,op);
    const Complex w=(Complex){cos(pi/n),sin(pi/n)*op};
    Complex k=(Complex){1,0};
    for(register int i=0;i<n;++i,k=k*w) 
    {
        f[i]=f1[i]+k*f2[i];
        f[i+n]=f1[i]-k*f2[i];
    }
}

接下来便是一个十分细节的操作了，因为我们得到了点值表达式之后还需要将点值表达式转回系数表达式才能最后输出得到我们常用的形式，那么怎么转换？对此前人想到了一个极为巧妙的办法，将\(w_n^k\)的对应点值作为一个新的系数表达式再将\(w_n^{-k}\)的任取依次代入得到对应的点值，那么推一下式子之后我们可以惊奇地发现以下结论：

\[设原来我们欲求的系数表达式为f(x)=a_0+a_1x+\cdots+a_{n-1}x^{n-1}=\sum\limits_{i=0}^{n-1}a_ix^i\\对应地将w_n^k代入可以得到f(w_n^k)=\sum\limits_{i=0}^{n-1}a_iw_n^{ik}
\]

那么将\(w_n^k\)依次代入得到新的系数表达式：

\[\begin{aligned}g(x)&=f(w_n^0)+f(w_n^1)x+f(w_n^2)x^2+\cdots+f(w_n^{n-1})x^{n-1}\\&=\sum\limits_{j=0}^{n-1}f(w_n^j)x^j=\sum\limits_{j=0}^{n-1}\sum\limits_{i=0}^{n-1}a_iw_n^{ij}x^j\\&=\sum\limits_{i=0}^{n-1}(a_i\sum\limits_{j=0}^{n-1}w_n^{ij}x^j)\\&=\sum\limits_{i=0}^{n-1}(a_i\sum\limits_{j=0}^{n-1}(w_n^{i}x)^j) \end{aligned}
\]

将我们前文提到的的\(w_n^{-k}\)任取依次代入，那么我们可以得到

\[g(w_n^{-k})=\sum\limits_{i=0}^{n-1}(a_i\sum\limits_{j=0}^{n-1}(w_n^{i-k})^j)
\]

认真考虑这个式子\(\sum\limits_{j=0}^{n-1}(w_n^{i-k})^j\)，我们可以从单位根的定义惊讶地发现一些十分有趣的结论

\[\sum\limits_{j=0}^{n-1}(w_n^{i-k})^j=\begin{cases}\sum\limits_{j=0}^{n-1}(w_n^0)^j=n&,i=k\\\sum\limits_{j=0}^{n-1}w_n^j=0&,i\not=k \end{cases}
\]

于是乎

\[g(w_n^{-k})=na_k
\]

之后我们将\(n\)除回去即可得到对应的\(a_k\)：

fft(a,n,1),fft(b,n,1);
for(register int i=0;i<n;++i)
    a[i]=a[i]*b[i];
fft(a,n,-1);
for(register int i=0;i<=m;++i)
    printf("%d ",(int)(a[i].a/n+0.5));

由于精度问题，输出时还是需要考虑四舍五入

当然，现在思路是出来了，但追求速度的我们怎么可能把如此大的常数放在这里不管呢.因为我们的主要工作程序是递归实现的，必定会导致极大的资源浪费与常数消耗，我们考虑怎么直接迭代递推解决

目的是要求出每次系数都进行奇偶分序后，最后原编号会到达哪个位置，考虑我们分序的原理：偶数放左，奇数放右——等同于我们将所有数以二进制排列，以低位往高位开始，越早出现1的在右，越早出现0的在左，那么最后我们简单分析可以知道新的位置编号的二进制相当于原编号二进制倒序，于是我们可以有如下两种递推出新位置的递推式：

法一：根据“越早出现1的在右，越早出现0的在左”从低位至高位依次加上某个数位上的1，复杂度\(O(n\log_2n)\)

inline void init(int n)
{
    for(sup=1;sup<=n;sup<<=1);
    for(register int k=sup>>1;k;k>>=1)
    {
        n=(sup/k)>>1;
        for(register int i=k;i<sup;i+=(k<<1))
            for(register int j=0;j<k;++j)
                rev[i+j]+=n;
    }
}

法二：根据“新的位置编号的二进制相当于原编号二进制倒序”直接进行数位运算，复杂度\(O(n)\)

inline void init(int n)
{
    for(sup=1,all=0;sup<=n;sup<<=1,++all);
	for(rg int i=0;i<sup;i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(all-1));
}

在程序中\(rev[]\)不需要更新时两种方法差别不大，可一旦需要较多更新时第二种方法会出现惊人的优化优势

当我们处理完\(rev[]\)后我们就可以有一个递推的模板了，常数会小约\(\dfrac14\)~\(\dfrac12\)：

inline void fft(Complex *x,int n,const int op)
{
    for(register int i=0;i<n;++i)
        if(i<rev[i])
            swap(x[i],x[rev[i]]);
    for(register int k=1;k<n;k<<=1)
    {
        const Complex w1=(Complex){cos(pi/k),op*sin(pi/k)};
        for(register int i=0;i<n;i+=(k<<1))
        {
            Complex wj=(Complex){1,0};
            for(register int j=0;j<k;++j,wj=wj*w1)
            {
                Complex l=x[i+j],r=x[i+k+j]*wj;
                x[i+j]=l+r;
                x[i+k+j]=l-r;
            }
        }
    }
}

\(part~3~ntt\)

从快速\(Fourier\)变换到快速数论变换，不得不说，这真的是一个神仙想法

首先说一说为什么我们需要快速数论变换

想必大家刚才也都发现了在\(fft\)的代码中我们用到了\(double\)类型的浮点数，既然有浮点数那就意味着我们的计算会伴随误差，最终的结果可能会在一些较大的运算或者故意构造的数据中产生一些不可避免的偏差.对此，前人就创造了一种方法使得我们的计算能够全在模某个质数的意义下进行整数范围内的运算，有效避免了浮点误差

接着是一些必须了解的数论概念：

\(·~\phi(x)为Euler函数(Euler永远滴神)，表示小于x的正整数中与x互质的数的个数（包括1）\\~ ~那么显然，对于一个质数p而言，\phi(p)=p-1\\·~对于某个数a，如果a^k\equiv1\pmod p，且对\forall i\in[1,k)\bigcap\Z，我们有a^k\not\equiv1\pmod p \\~ ~那么我们称k为a模p意义下的阶数，记作ord_p a=k\\·~如果对于某个数g满足ord_p g=\phi(p)，我们称g为p的一个原根\)

接下来我们考虑一点特殊的：当p为质数时我们看看其一个原根g会怎么样

同时我们不妨设\(g_n\equiv g^{\frac{p-1}{n}}\pmod p\)，并且我们的\(p-1\)在题目所要求的范围内一定是2的整数次幂倍数（放心，这样的质数有很多），那么

最首先的我们会得到\(g^i\not\equiv g^j\pmod p，i\not=j\and i,j\in[1,p)\bigcap\Z\)

\(p.f.\)

\(根据Fermat小定理，显然我们会有\)

\[g^{p-1}\equiv g^{\phi(p)}\equiv1\pmod p
\]

\(同时根据阶数的定义我们可以得到\)

\[g^i\not\equiv 1\pmod p，i\in[1,p-1)\bigcap\Z
\]

\(那么考虑反证法，假设 ~ \exist ~ i,j\in[1,p)\bigcap\Z ~ ~ s.t.g^i\equiv g^j\pmod p\)

\(不妨设i< j，那么根据数论意义下乘法的性质我们会有\)

\[g^{j-i}\equiv1\pmod p
\]

\(显然j-i\in[1,p-1)\bigcap\Z，与之前结论矛盾\)

\(故假设不成立，原命题成立\)，证毕

既然如此，那么显然\(g^k\)在模p意义下当k遍取1 ~ p-1时\(g^k\)同样遍取1 ~ p-1

接下来\(ntt\)的想法其实就是直接用我们的\(g_n^k\)来代替我们的\(w_n^k\)，是不是感觉为什么可以这样，确实很不可思议，但事实如此

总结一下我们用单位根来进行点值替换和重新求系数的两个特别的性质：

\(① ~ w_n^{2k}=w_{n/2}^k \\② ~ w_n^{k+\frac{n}2}=-w_n^{k}，k\in[0,\frac{n}2)\)

显然，如果我们的\(g_n\)也满足这样的性质，那么进行与\(part2\)一模一样的运算是显然可以的

那么接下来我们就来证明这两个对应的性质（方便起见，最后的“（mod p）”全部略去）

\[\begin{aligned} g_n^{2k}&\equiv (g^{\frac{p-1}{n}})^{2k}\equiv (g^{\frac{p-1}{n/2}})^k\equiv g_{n/2}^k \end{aligned}
\]

对于第二个，我们先证明一个引理：\(g_n^{n/2}\equiv -1\)

\(p.f.\)

\(显然\)

\[g_n^{n/2}\equiv g^{\frac{p-1}2}
\]

\(而根据Fermat小定理\)

\(g^{p-1}\equiv 1\)

\(因此\)

\[g^{\frac{p-1}{2}}\equiv \pm 1
\]

\(而由之前得出的原根的性质\)

\[g^{\frac{p-1}2}\not\equiv g^{p-1}\equiv 1
\]

\(因此\)

\[g_n^{n/2}\equiv g^{\frac{p-1}2}\equiv -1
\]

证毕

而由以上结论我们可以很轻松地得到性质二：

\[\begin{aligned} g_n^{k+\frac{n}2}\equiv g_n^kg_n^{n/2}\equiv -g_n^k \end{aligned}
\]

因此我们就可以顺利地用\(fft\)的思路得到\(ntt\)的算法啦

当然还有最后一个细节：得到满足我们最初假设条件的质数. 当然对于这个我们的想法不会是自己去硬算暴搞，一个好的办法是写一个线性筛算法去暴力枚举看看是否满足我们的2的整数次幂的条件，大概\(12000ms\)我们就能得到\(10^9\)以内的所有满足条件（\(2^{20}|p-1\)）的质数如下：

23068673
69206017
81788929
104857601
111149057
113246209
132120577
136314881
138412033
155189249
163577857
167772161
169869313
186646529
199229441
211812353
230686721
249561089
257949697
270532609
274726913
377487361
383778817
387973121
415236097
459276289
463470593
469762049
576716801
595591169
597688321
635437057
639631361
645922817
648019969
666894337
683671553
710934529
715128833
740294657
754974721
786432001
799014913
824180737
880803841
897581057
899678209
918552577
924844033
935329793
943718401
950009857
962592769
975175681
985661441
998244353

与之对应的原根和原根的逆（数论中用逆来代替负一次方）如下：

3 7689558
5 27682407
7 46736531
3 34952534
3 37049686
7 16178030
5 52848231
3 45438294
5 83047220
6 25864875
23 28448323
3 55924054
5 101921588
3 62215510
3 66409814
3 70604118
6 38447787
3 83187030
5 103179879
22 110672431
3 91575638
7 53926766
5 153511527
6 64662187
5 166094439
11 41752390
3 154490198
3 156587350
6 96119467
3 198530390
11 54335302
11 288835026
6 106605227
3 215307606
17 266831752
5 266757735
3 227890518
17 460016460
3 238376278
3 246764886
11 617706590
7 449389715
13 491701485
5 329672295
26 711418487
3 299193686
7 257050917
5 367421031
5 554906420
3 311776598
7 269633829
7 135715694
7 550053011
17 803085855
3 328553814
3 332748118

可以观察到原根一般都比较小，这也是我们暴力枚举起来效率还比较高的原因

而对于所有满足条件的质数，我们在日常用得比较多的是998244353，这是因为这个数是在满足相乘后不超过64位级整数的最大满足条件的质数，我们通常取其做模数以便完成一些特殊的运算

\(ntt\)的一种代码如下（\(ps\)：该模板满足输入的系数在0~9范围内，如果不满足需要重改并替换一个大一点的模数）：

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

typedef long long ll;
const int g=3,_g=7689558;
const int mod=23068673;
const int xx=(1<<21)+1;

ll a[xx],b[xx];
int rev[xx];

inline int read()
{
    char ch=getchar();
    for(;!isalnum(ch);ch=getchar());
    return ch-'0';
}
inline ll power(ll x,int y)
{
    ll ans=1;
    for(;y;(x*=x)%=mod,y>>=1)
        if(y&1)
            (ans*=x)%=mod;
    return ans;
}
inline void ntt(ll *x,int n,const bool op)
{
    for(register int i=0;i<n;++i)
        if(i<rev[i])
            swap(x[i],x[rev[i]]);
    for(register int k=1;k<n;k<<=1)
    {
        const ll g1=power(op?g:_g,(mod-1)/(k<<1));
        for(register int i=0;i<n;i+=(k<<1))
        {
            ll gj=1;
            for(register int j=0;j<k;++j,(gj*=g1)%=mod)
            {
                const ll l=x[i+j],r=(x[i+k+j]*gj)%mod;
                x[i+j]=(l+r)%mod;
                x[i+k+j]=(l-r+mod)%mod;
            }
        }
    }
}

int main()
{
    int n,m;
    scanf("%d%d",&n,&m);
    for(register int i=0;i<=n;++i) a[i]=read();
    for(register int i=0;i<=m;++i) b[i]=read();
    m=n+m;
    for(n=1;n<=m;n<<=1);
    for(register int k=n>>1;k;k>>=1)
        for(register int i=k;i<n;i+=(k<<1))
            for(register int j=0;j<k;++j)
                rev[i+j]+=(n/k/2);
    ntt(a,n,1),ntt(b,n,1);
    for(register int i=0;i<n;++i) (a[i]*=b[i])%=mod;
    ntt(a,n,0);
    for(register ll i=0,j=power(n,mod-2);i<=m;++i) printf("%lld ",(a[i]*j)%mod);
    putchar('\n');
    system("pause");
    return 0;
}

原文链接: https://www.cnblogs.com/cjsyh/p/17061885.html

欢迎关注

微信关注下方公众号，第一时间获取干货硬货；公众号内回复【pdf】免费获取数百本计算机经典书籍

    多项式乘法

原创文章受到原创版权保护。转载请注明出处：https://www.ccppcoding.com/archives/313279

非原创文章文中已经注明原地址，如有侵权，联系删除

关注公众号【高性能架构探索】，第一时间获取最新文章

转载文章受原作者版权保护。转载请注明原作者出处！

赞 (0)
adminadmin
0 0
上一篇 2023年2月16日 下午12:45
下一篇 2023年2月16日 下午12:45

相关推荐