It will be another blog with nothing but the board

Although there may be a pile of things in front of the board, it's all nonsense..

Convolution and DFT

The formula like this is called convolution:

Where \ (\ circ \) is any operation. When it is \ (\), the convolution is polynomial convolution. It is not difficult to find that \ (\ left\{c\right \} \) is the coefficient of the polynomial obtained by multiplying the polynomials with coefficients of \ (\ left\{a\right \} \) and \ (\ left\{b\right \} \) respectively.

The direct calculation of convolution is \ (O(n^2) \). Consider transforming the three sequences so that the transformed \ (\ hat a \), \ (\ hat b \), \ (\ hat c \) satisfy the following formula

That is, the convolution is converted to parity multiplication, so that the convolution can be completed \ (O(n) \). Finally, the transformation can be mapped back in some way.

A convenient method is discrete Fourier transform DFT. It converts the coefficient vector \ (\ VEC a = (a_0, a_1, \ ldots, a {n-2}, a {n-1}) \) of \ (n-1 \) degree polynomial \ (f \ (\ Omega _n ^ 0), f (\ Omega _n ^ 1), \ ldots, f (\ Omega _n ^ {n-2}), f (\ Omega _n ^ {n-1})) \). Where \ (\ omega_n \) is the sub unit root of \ (n \), that is, the solution of \ (x^n=1 \) in the complex field. Generally, the one with the smallest argument is selected, i.e. \ (\ cos(\frac{2\pi}{n})+\sin(\frac{2\pi}{n})i \).

The unit root is chosen as the abscissa of the point value because of its many excellent properties. No watch here, No. Can refer to Weekly guidance . (in fact, most of the things in this blog are learned here

FFT

An algorithm to accelerate the discrete Fourier transform process by using the property of unit root in complex field.

The specific nature will not. You can learn from Zhou's guidance.

Because I can't understand it, at first, the backplane was really painful in the FFT part. Add a note without logic here.

void fft(complex *a,int n){
    for(int i=0;i<n;i++)
        if(r[i]>i) swap(a[i],a[r[i]]); //The sequence to be iterated is proposed, and r[i] is processed first
    for(int t=n>>1,d=1;d<n;t>>=1,d<<=1) //d is half of the length of the iteration interval, and t is the coefficient for calculating the number of unit roots
        for(int i=0;i<n;i+=(d<<1)) //i is the left end point of the current iteration interval
            for(int j=0;j<d;j++){ //j is the location of the current process
                complex tmp=w[t*j]*a[i+j+d];
                a[i+j+d]=a[i+j]-tmp;
                a[i+j]=a[i+j]+tmp;
            }
}

It can't help understanding, but it can help memory

Rogue P3803 [template] polynomial multiplication

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

namespace IO{
    typedef long long LL;
    typedef double DB;
    int read(){
        int x=0,f=0; char ch=getchar();
        while(ch>'9'||ch<'0'){ f|=(ch=='-'); ch=getchar(); }
        while(ch>='0'&&ch<='9'){ x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); }
        return f?-x:x;
    } char output[50];
    void write(int x,char sp){
        int len=0;
        if(x<0) putchar('-'), x=-x;
        do{ output[len++]=x%10+'0'; x/=10; }while(x);
        for(int i=len-1;~i;i--) putchar(output[i]); putchar(sp);
    }
    void ckmax(int& x,int y){ x=x>y?x:y; }
    void ckmin(int& x,int y){ x=x<y?x:y; }
} using namespace IO;

const int NN=1<<22;
int dega,degb;

namespace Poly_Calculation{
    const DB PI=acos(-1.0);
    struct complex{
        DB r,i;
        complex(){ r=i=0; }
        complex(DB x,DB y){ r=x; i=y; }
        complex operator+(const complex& t)const{ return complex(r+t.r,i+t.i); }
        complex operator-(const complex& t)const{ return complex(r-t.r,i-t.i); }
        complex operator*(const complex& t)const{ return complex(r*t.r-i*t.i,r*t.i+i*t.r); }
    }a[NN],b[NN],w[NN];
    int l,len,r[NN];
    void prework(){
        for(len=1;len<=dega+degb;len<<=1) ++l;
        for(int i=0;i<len;i++){
            r[i]=(r[i>>1]>>1)|((i&1)<<l-1);
            w[i]=complex(cos(2.0*i*PI/len),sin(2.0*i*PI/len));
        }
    }
    void fft(complex *a,int n){
        for(int i=0;i<n;i++)
            if(r[i]>i) swap(a[i],a[r[i]]);
        for(int t=n>>1,d=1;d<n;t>>=1,d<<=1)
            for(int i=0;i<n;i+=(d<<1))
                for(int j=0;j<d;j++){
                    complex tmp=w[t*j]*a[i+j+d];
                    a[i+j+d]=a[i+j]-tmp;
                    a[i+j]=a[i+j]+tmp;
                }
    }
    void poly_mul(complex *a,complex *b){
        fft(a,len); fft(b,len);
        for(int i=0;i<len;i++)
            a[i]=a[i]*b[i], w[i].i=-w[i].i;
        fft(a,len);
        for(int i=0;i<len;i++) a[i].r=a[i].r/len;
    }
} using namespace Poly_Calculation;

signed main(){
    dega=read()+1; degb=read()+1;
    for(int i=0;i<dega;i++) a[i].r=read();
    for(int i=0;i<degb;i++) b[i].r=read();
    prework(); poly_mul(a,b);
    --dega; --degb;
    for(int i=0;i<=dega+degb;i++)
        write(int(a[i].r+0.5),' ');
    return puts(""),0;
}

NTT

Although FFT can complete fast polynomial convolution, it also has some disadvantages. Such as the problem of accuracy, and the difficulty of convenient operation when it comes to module taking, etc.

NTT is required at this time. Generally speaking, the principles of NTT and FFT are completely consistent, except that the original root (integer power of) in the modular sense is used to replace the unit root in the complex field, avoiding floating-point calculation.

After some derivation, if the original root in the modular sense is \ (g \), then the original \ (w_n \) can be replaced by \ (g^{\frac{mod-1}{n}} \).

Because the polynomial length takes the integer power of \ (2 \), the exponent of \ (g \) is required to be an integer, and there are some requirements for modulus.

Two common modules are \ (998244353 \), \ (1004535809 \), and the original root is \ (3 \).

Almost as like as two peas in FFT.

Rogue P3803 [template] polynomial multiplication

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

namespace IO{
    typedef long long LL;
    typedef double DB;
    LL read(){
        LL x=0,f=0; char ch=getchar();
        while(ch>'9'||ch<'0'){ f|=(ch=='-'); ch=getchar(); }
        while(ch>='0'&&ch<='9'){ x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); }
        return f?-x:x;
    } char output[50];
    void write(LL x,char sp){
        int len=0;
        if(x<0) putchar('-'), x=-x;
        do{ output[len++]=x%10+'0'; x/=10; }while(x);
        for(int i=len-1;~i;i--) putchar(output[i]); putchar(sp);
    }
    void ckmax(int& x,int y){ x=x>y?x:y; }
    void ckmin(int& x,int y){ x=x<y?x:y; }
} using namespace IO;

const int NN=1<<22,mod=998244353;
int dega,degb;
LL qpow(LL a,LL b=mod-2,LL res=1){
    for(;b;b>>=1,a=a*a%mod)
        if(b&1) res=res*a%mod;
    return res;
}

namespace Poly_Calculation{
    LL l,len,inv,g[NN],r[NN],a[NN],b[NN];
    void prework(){
        for(len=1;len<=dega+degb;len<<=1) ++l;
        for(int i=0;i<len;i++)
            r[i]=(r[i>>1]>>1)|((i&1)<<l-1);
        g[0]=1; g[1]=qpow(3,(mod-1)/len);
        for(int i=2;i<len;i++) g[i]=g[i-1]*g[1]%mod;
    }
    void ntt(LL *a,int n){
        for(int i=0;i<n;i++)
            if(r[i]>i) swap(a[r[i]],a[i]);
        for(int t=n>>1,d=1;d<n;d<<=1,t>>=1)
            for(int i=0;i<n;i+=(d<<1))
                for(int j=0;j<d;j++){
                    LL tmp=g[t*j]*a[i+j+d]%mod;
                    a[i+j+d]=(a[i+j]-tmp+mod)%mod;
                    a[i+j]=(a[i+j]+tmp)%mod;
                }
    }
    void poly_mul(LL *a,LL *b){
        ntt(a,len); ntt(b,len);
        for(int i=0;i<len;i++)
            a[i]=a[i]*b[i]%mod, g[i]=qpow(g[i]);
        ntt(a,len);
        inv=qpow(len);
        for(int i=0;i<len;i++) a[i]=a[i]*inv%mod;
    }
} using namespace Poly_Calculation;

signed main(){
    dega=read()+1; degb=read()+1;
    for(int i=0;i<dega;i++) a[i]=read();
    for(int i=0;i<degb;i++) b[i]=read();
    prework(); poly_mul(a,b);
    --dega; --degb;
    for(int i=0;i<=dega+degb;i++)
        write(a[i],' ');
    return puts(""),0;
}