【FFT优化】[UOJ#34]多项式乘法

题目描述

这是一道模板题。

给你两个多项式，请输出乘起来后的多项式。

输入格式

第一行两个整数 $n$ 和 $m$，分别表示两个多项式的次数。

第二行 $n + 1$ 个整数，分别表示第一个多项式的 $0$ 到 $n$ 次项前的系数。

第三行 $m + 1$ 个整数，分别表示第一个多项式的 $0$ 到 $m$ 次项前的系数。

输出格式

一行 $n + m + 1$ 个整数，分别表示乘起来后的多项式的 $0$ 到 $n + m$ 次项前的系数。

样例一

input

1 2
1 2
1 2 1

output

1 4 5 2

explanation

$(1 + 2x) \cdot (1 + 2x + x^2) = 1 + 4x + 5x^2 + 2x^3$。

限制与约定

$0 \leq n, m \leq 10^5$，保证输入中的系数大于等于 $0$ 且小于等于 $9$。

时间限制：$1\texttt{s}$

空间限制：$256\texttt{MB}$

分析

本博客主要用于存$FFT$优化模板.
对于两个实系数的多项式$A(x)$,$B(x)$，长度均为$N$($N$为2的整次数幂)，令

$$P(x)=A(x)+iB(x)\\Q(x)=A(x)-iB(x)$$
$F_p[k],F_q[k]$为 $P(x)$, $Q(x)$进行 $DFT$后的值.
令 $X=\frac{2\pi jk}{N}$
$$\begin{align} F_p[k]&=A(\omega^k_{N})+iB(\omega^k_{N})\\ &=\sum_{j=0}^{N-1}A_j\omega^{jk}_{N}+i\sum_{j=0}^{N-1}B_j\omega^{jk}_{N}\\ &=\sum_{j=0}^{N-1}(A_j+iB_j)*(\cos(X)+i\sin(X)) \end{align}$$
$$\begin{align} F_q[k]&=A(\omega^k_{N})-iB(\omega^k_{N})\\ &=\sum_{j=0}^{N-1}A_j\omega^{jk}_{N}-i\sum_{j=0}^{N-1}B_j\omega^{jk}_{N}\\ &=\sum_{j=0}^{N-1}(A_j-iB_j)*(\cos(X)+i\sin(X))\\ &=\sum_{j=0}^{N-1}A_j*\cos(X)+iA_j\sin(X)-iB_j\cos(X)+B_j\sin(X)\\ &=\sum_{j=0}^{N-1}(A_j\cos(X)+B_j\sin(X))+i(A_j\sin(X)-B_j\cos(X))\\ &=\sum_{j=0}^{N-1}conj((A_j\cos(X)+B_j\sin(X))-i(A_j\sin(X)-B_j\cos(X)))\\ &=\sum_{j=0}^{N-1}conj((A_j\cos(-X)-B_j\sin(-X))+i(A_j\sin(-X)+B_j\cos(-X)))\\ &=\sum_{j=0}^{N-1}conj(A_j(\cos(-X)+i\sin(-X))+iB_j(\cos(-X)+i\sin(-X)))\\ &=\sum_{j=0}^{N-1}conj((A_j+iB_j)*(\cos(-X)+i\sin(-X))\\ &=\sum_{j=0}^{N-1}conj(A_j\omega^{-jk}_{N}+iB_j(\omega^{-jk}_{N}))\\ &=conj(A(\omega^{-k}_{N})+iB(\omega^{-k}_{N}))\\ &=conj(A(\omega^{N-k}_{N})+iB((\omega^{N-k}_{N}))\\ &=conj(F_p[N-k]) \end{align}$$
这样我们就可以通过一次DFT求出 $F_p,F_q$
$$DFT(A[k])=\frac{F_p[k]+F_q[k]}{2}\\ DFT(B[k])=\frac{F_p[k]-F_q[k]}{2i}=-i\frac{F_p[k]-F_q[k]}{2}$$
如果只是做多项式乘法 $C(x)=A(x)*B(x)$，直接使 $N$变为$A$, $B$长度的和，然后 $$DFT(C[k])=(F_p[k]^2-F_q[k]^2)*\frac{-i}{4}=(F_p[k]^2-F_p[N-k]^2)*\frac{-i}{4}$$

代码

#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
#define MAXN 300000
using namespace std;
const double pi=acos(-1);
int n,m,N,ans[MAXN+10];
char s[MAXN+10];
struct cpx{
    double r,i;
    inline cpx(){
    }
    inline cpx(double r,double i):r(r),i(i){
    }
    inline cpx operator+(const cpx &b)const{
        return cpx(r+b.r,i+b.i);
    }
    inline cpx operator-(const cpx &b)const{
        return cpx(r-b.r,i-b.i);
    }
    inline cpx operator*(const cpx &b)const{
        return cpx(r*b.r-i*b.i,r*b.i+i*b.r);
    }
    inline cpx &operator*=(const cpx &b){
        return *this=*this*b;
    }
    inline cpx conj()const{
        return cpx(r,-i);
    }
}a[MAXN+10],b[MAXN+10],t[MAXN+10];
void fft(cpx *a,int N,int f){
    int i,j=0,k,t;
    for(i=1;i<N-1;i++){
        for(t=N;j^=t>>=1,~j&t;);
        if(i<j)
            swap(a[i],a[j]);
    }
    for(i=1;i<N;i<<=1){
        cpx wn(cos(pi/i),f*sin(pi/i));
        t=i<<1;
        for(j=0;j<N;j+=t){
            cpx w(1,0);
            for(k=0;k<i;k++,w*=wn){
                cpx x(a[j+k]),y(w*a[j+k+i]);
                a[j+k]=x+y,a[j+k+i]=x-y;
            }
        }
    }
    if(f==-1)
        for(i=0;i<N;i++)
            a[i].r/=N;
}
void mul(cpx *a,cpx *b,cpx *c,int N){
    int i,j;
    for(i=0;i<N;i++)
        t[i]=cpx(a[i].r,b[i].r);
    fft(t,N,1);
    for(i=0;i<N;i++){
        j=(N-i)&(N-1);
        c[i]=(t[i]*t[i]-(t[j]*t[j]).conj())*cpx(0,-0.25);
    }
    fft(c,N,-1);
    // c refer t
    //a[i]=(c[i]+c[k-i])/2
    //b[i]=-i(c[i]-c[k-i])/2
    //a[i]*b[i]=(c[i]*c[i]-c[k-i]*c[k-i])*((-i)/4）
}
void read(){
    int lenb,lena,i;
    scanf("%s",s);
    lena=strlen(s);
    for(i=0;i<lena;i++)
        a[i].r=s[lena-i-1]-'0';
    scanf("%s",s);
    lenb=strlen(s);
    for(i=0;i<lenb;i++)
        b[i].r=s[lenb-i-1]-'0';
    for(N=1;N<=lena+lenb-2;N<<=1);
}
void print(){
    int i;
    for(i=0;i<N;i++){
        ans[i]+=a[i].r+0.2;
        ans[i+1]+=ans[i]/10;
        ans[i]%=10;
    }
    for(i=N-1;i;i--)
        if(ans[i])
            break;
    for(;i>=0;i--)
        printf("%d",ans[i]);
    puts("");
}
int main()
{
    read();
    mul(a,b,a,N);
    print();
}