多項式細節梳理&模板（多項式）

時間  2019-11-09
標籤 多項式 細節 梳理 模板 简体版
原文   原文鏈接

基礎

好久之前的多項式總結html

如今的碼風又變了。。。c++

FFT和NTT的板子數組

typedef complex<double> C;
const double PI=acos(-1);
void FFT(C*a,R op){
    for(R i=0;i<N;++i)
        if(i<r[i])swap(a[i],a[r[i]]);
    for(R i=1;i<N;i<<=1){
        C wn=C(cos(PI/i),sin(PI/i)*op),w=1,t;
        for(R j=0;j<N;j+=i<<1,w=1)
            for(R k=j;k<j+i;++k,w*=wn)
                t=a[k+i]*w,a[k+i]=a[k]-t,a[k]+=t;
    }
}
const int YL=998244353;
LL Pow(LL b,R k=YL-2,LL a=1){
    for(;k;k>>=1,b=b*b%YL)
        if(k&1)a=a*b%YL;
    return a;
}
void NTT(R*a,R n,R op){
    for(R i=0;i<n;++i)
        if(i<r[i])swap(a[i],a[r[i]]);
    for(R i=1;i<n;i<<=1){
        LL wn=Pow(3,(YL-1)/(i<<1)),w=1,x;
        if(op==-1)wn=Pow(wn);
        for(R j=0;j<n;j+=i<<1,w=1)
            for(R k=j;k<j+i;++k,w=w*wn%YL)
                x=w*a[k+i]%YL,a[k+i]=Mod(a[k]-x+YL),a[k]=Mod(a[k]+x);
    }
    if(op==-1){
        LL x=Pow(n);
        for(R i=0;i<n;++i)a[i]=a[i]*x%YL;
    }
}

任意模數NTT

怎麼看都以爲MTT給人感受很差，用double巨佬說怕掉精度，用longdouble那常數想象一下～
而後就去學了三模NTT
參考Memory of Winter巨佬的題解
細節很少吧，首先要把\(998244353,1004535809,469762049\)背下來，它們都有一個原根爲\(3\)。
而後手動兩兩合併同餘方程，不要直接用CRT公式。
蒟蒻爲了減小碼量直接開二維數組了，實際上像Memory of Winter巨佬同樣封裝結構體會跑得更快。
upd:重寫了一遍，二維數組又長又醜又慢qwq
洛谷P4245 【模板】任意模數NTT函數

#include<bits/stdc++.h>
#define LL long long
#define I inline
#define R register int
#define G if(++ip==ie)if(fread(ip=buf,1,SZ,stdin))
#define Wn(A) Pow(3,Mod((A-1)/(i<<1)*op,A-1),A)
using namespace std;
const int SZ=1<<19,N=1<<18,YL=1e9+7,A=998244353,B=1004535809,C=469762049;
char buf[SZ],*ie=buf+SZ,*ip=ie-1;
inline int in(){
    G;while(*ip<'-')G;
    R x=*ip&15;G;
    while(*ip>'-'){x*=10;x+=*ip&15;G;}
    return x;
}
int L,r[N];
I int Mod(R x,R YL){
    return x+(x>>31&YL);
}
I int Pow(LL b,R k,LL YL,LL a=1){
    for(;k;k>>=1,b=b*b%YL)
        if(k&1)a=a*b%YL;
    return a;
}
I int Inv(LL b,LL YL){
    return Pow(b%YL,YL-2,YL);
}
struct Z{
    int a,b,c;
    Z(){}
    Z(LL a):a(a%A),b(a%B),c(a%C){}
    Z(R a,R b,R c):a(a),b(b),c(c){}
    I Z operator!(){a+=a>>31&A,b+=b>>31&B,c+=c>>31&C;return*this;}
    I Z operator+(const Z&x){return!Z(a+x.a-A,b+x.b-B,c+x.c-C);}
    I Z operator-(const Z&x){return!Z(a-x.a  ,b-x.b  ,c-x.c  );}
    I Z operator*(const Z&x){return Z((LL)a*x.a%A,(LL)b*x.b%B,(LL)c*x.c%C);}
    I int CRT(LL YL){
        static LL AB=(LL)A*B%YL,I0=Inv(A,B),I1=Inv((LL)A*B,C),x;
        x=Mod(b-a,B)*I0%B*A+a;
        return(Mod(c-x%C,C)*I1%C*AB+x)%YL;
    }
}f[N],g[N];
void NTT(Z*a,R op){
    for(R i=0;i<L;++i)
        if(i<r[i])swap(a[i],a[r[i]]);
    for(R i=1;i<L;i<<=1){
        Z wn(Wn(A),Wn(B),Wn(C)),w(1),x;
        for(R j=0;j<L;j+=i<<1,w=1)
            for(R k=j;k<j+i;++k,w=w*wn)
                x=w*a[k+i],a[k+i]=a[k]-x,a[k]=a[k]+x;
    }
    if(op==-1){
        Z w(Inv(L,A),Inv(L,B),Inv(L,C));
        for(R i=0;i<L;++i)a[i]=a[i]*w;
    }
}
int main(){
    R n=in(),m=in(),p=in();
    for(R i=0;i<=n;++i)f[i]=in();
    for(R i=0;i<=m;++i)g[i]=in();
    for(L=1;L<=n+m;L<<=1);
    for(R i=1;i<L;++i)r[i]=(r[i>>1]|L*(1&i))>>1;
    NTT(f,1);NTT(g,1);
    for(R i=0;i<L;++i)f[i]=f[i]*g[i];
    NTT(f,-1);
    for(R i=0;i<=n+m;++i)printf("%d ",f[i].CRT(p));
    return 0;
}

另外，三模NTT不能將三個及以上的多項式一次都乘起來，由於實際值域太大，不能保證模意義下的惟一性。
正確的作法是兩個兩個乘。
好比，任意模數多項式求逆的核心代碼
洛谷P4239 【模板】多項式求逆（增強版）ui

void PolyInv(Z*a,Z*b,Z*a1,R m){
    if(m==1){b[0]=Inv(a[0].b,YL);return;}
    PolyInv(a,b,a1,(m+1)>>1);memcpy(a1,a,12*m);
    Init(2*m);NTT(b,1);NTT(a1,1);
    for(R i=0;i<L;++i)a1[i]=b[i]*a1[i];
    NTT(a1,-1);
    for(R i=0;i<L;++i)a1[i]=Mod(-a1[i].CRT()+2*(i==0),YL);
    NTT(a1,1);
    for(R i=0;i<L;++i)b[i]=b[i]*a1[i];
    NTT(b,-1);memset(a1,0,12*L);memset(b+m,0,12*(L-m));
    for(R i=0;i<m;++i)b[i]=b[i].CRT();
}

多項式運算

分治乘法

若干個總項數不超過\(n\)的多項式的卷積，根據合併果子的方法，能夠在不超過\(O(n\log^2n)\)的時間內完成。
有時候是若干個二項式相乘，還能夠寫成循環的形式，好處是減小Rader排序r數組的預處理次數。
第一次嘗試這種寫法：洛谷CF981H K Pathsthis

void Solve(R x){
    R n=4*e[x].size();
    for(R i=0;i<n;i+=4)
        a[i]=1,a[i+1]=s[e[x][i>>2]];
    for(R i=4;i<n;i<<=1){
        for(R j=1;j<i;++j)
            r[j]=(r[j>>1]>>1)|(1&j?i>>1:0);
        for(R j=0;j+i<n;j+=i<<1){
            R*p=a+j,*q=p+i;
            NTT(p,i,1);NTT(q,i,1);
            for(R k=0;k<i;++k)p[k]=(LL)p[k]*q[k]%YL;
            NTT(p,i,-1);
            memset(q,0,4*i);
        }
    }
}

分治FFT

某些卷積式有一些特色，直接算計算量龐大，卻能夠經過CDQ分治思想，考慮左邊對右邊的貢獻。
洛谷P4721 【模板】分治 FFTspa

void Solve(R l,R r){
    if(l+1==r)return;
    R m=(l+r)>>1;
    Solve(l,m);R n=Init(m-l+r-l);
    memcpy(a,f+l,4*(m-l));memset(a+m-l,0,4*(n-m+l));
    memcpy(b,g  ,4*(r-l));memset(b+r-l,0,4*(n-r+l));
    NTT(a,n,1);NTT(b,n,1);
    for(R i=0;i<n;++i)a[i]=(LL)a[i]*b[i]%YL;
    NTT(a,n,-1);
    for(R i=m;i<r;++i)f[i]=Mod(f[i]+a[i-l]);
    Solve(m,r);
}

泰勒展開

將函數\(f(x)\)在\(x_0\)處展開得
\(f(x)=f(x_0)+\frac{f'(x_0)}{1!}(x-x_0)+\frac{f''(x_0)}{2!}(x-x_0)^2+...+\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n+\xi\)code

牛頓迭代推導倍增求解式

已知多項式函數\(F(x)\)，用倍增法求解\(B\)使得\(F(B)\equiv0(\mod x^{2^k})\)。
假設已經求出\(B,F(B)\equiv0(\mod x^{2^k})\)，如今求\(B_1,F(B_1)\equiv0(\mod x^{2^{k+1}})\)。
\(F(B_1)\)在\(x=B\)處只展開一項獲得
\(F(B_1)=F(B)+F'(B)(B_1-B)\equiv0(\mod x^{2^{k+1}})\)
\(B_1=B-\frac{F(B)}{F'(B)}\)，帶入便可。
接下來的\(A\)爲已知多項式，B爲待求多項式。htm

求逆

\(F(B)=AB-1=0\)
\(B_1=B-\frac{AB-1}{A}=B-B(AB-1)=2B-AB^2\)blog

開方

\(F(B)=B^2-A=0\)
\(B_1=B-\frac{B^2-A}{2B}=\frac{B^2+A}{2B}\)

對數

不必迭代。
\(B=\ln A\)
\(B=\int\frac{A'}{A}\)

指數

\(B=e^A\)
\(\ln B-A=0\)
\(B_1=B-\frac{\ln B-A}{\frac1B}=B(1+A-\ln B)\)

快速冪

\(B=A^k=e^{k\ln A}\)
這常數看着就嚇人。。。

除法&取餘

https://www.luogu.org/problemnew/solution/P4512

代碼模板

寫的時候保證了必定的穩定性（數組清空問題），所以犧牲了一丁點效率。
參數中，a爲已知，b爲待求，a一、b1爲輔助數組。
傳入前需本身確保b、a一、b1爲空。
返回後保證a維持原狀，b存好答案，a一、b1爲空。
還有，線性遞推，多點求值，插值待填坑

const int YL=998244353;
int Inv[N],r[N];//在外面初始化逆元
inline int Mod(const int x){
    return x>=YL?x-YL:x;
}
LL Pow(LL b,R k=YL-2,LL a=1){
    for(;k;k>>=1,b=b*b%YL)
        if(k&1)a=a*b%YL;
    return a;
}
int Init(R m){
    R n=1;while(n<m<<1)n<<=1;
    for(R i=0;i<n;++i)r[i]=(r[i>>1]|(1&i)*n)>>1;
    return n;
}
void NTT(R*a,R n,R op){
    for(R i=0;i<n;++i)
        if(i<r[i])swap(a[i],a[r[i]]);
    for(R i=1;i<n;i<<=1){
        LL wn=Pow(3,(YL-1)/(i<<1)),w=1,x;
        if(op==-1)wn=Pow(wn);
        for(R j=0;j<n;j+=i<<1,w=1)
            for(R k=j;k<j+i;++k,w=w*wn%YL)
                x=w*a[k+i]%YL,a[k+i]=Mod(a[k]-x+YL),a[k]=Mod(a[k]+x);
    }
    if(op==-1){
        LL x=Pow(n);
        for(R i=0;i<n;++i)a[i]=a[i]*x%YL;
    }
}
void PolyRev(R*a,R*b,R n){
    for(R i=0;i<n;++i)b[i]=a[n-i-1];
}
void PolyDer(R*a,R*b,R n){
    for(R i=1;i<n;++i)b[i-1]=(LL)a[i]*i%YL;
    if(a==b)a[n-1]=0;
}
void PolyInt(R*a,R*b,R n){
    for(R i=n;i;--i)b[i]=(LL)a[i-1]*Inv[i]%YL;
    if(a==b)a[0]=0;
}
void PolyInv(R*a,R*b,R*a1,R m){
    if(m==1){b[0]=Pow(a[0]);return;}
    PolyInv(a,b,a1,(m+1)>>1);memcpy(a1,a,4*m);
    R n=Init(m);NTT(b,n,1);NTT(a1,n,1);
    for(R i=0;i<n;++i)b[i]=(YL+2-(LL)b[i]*a1[i]%YL)*b[i]%YL;
    NTT(b,n,-1);memset(a1,0,4*n);memset(b+m,0,4*(n-m));
}
void PolySqrt(R*a,R*b,R*a1,R*b1,R m){
    if(m==1){b[0]=sqrt(a[0]);return;}
    PolySqrt(a,b,a1,b1,(m+1)>>1);PolyInv(b,b1,a1,m);memcpy(a1,a,4*m);
    R n=Init(m);NTT(a1,n,1);NTT(b1,n,1);
    for(R i=0;i<n;++i)a1[i]=(LL)a1[i]*b1[i]%YL;
    NTT(a1,n,-1);
    for(R i=0;i<m;++i)b[i]=(LL)(b[i]+a1[i])*((YL+1)>>1)%YL;
    memset(a1,0,4*n);memset(b1,0,4*n);memset(b+m,0,4*(n-m));
}
void PolyLn(R*a,R*b,R*a1,R m){
    PolyInv(a,b,a1,m);PolyDer(a,a1,m);
    R n=Init(m);NTT(b,n,1);NTT(a1,n,1);
    for(R i=0;i<n;++i)b[i]=(LL)b[i]*a1[i]%YL;
    NTT(b,n,-1);memset(a1,0,4*n);PolyInt(b,b,m);
}
void PolyExp(R*a,R*b,R*a1,R*b1,R m){
    if(m==1){b[0]=1;return;}
    PolyExp(a,b,a1,b1,(m+1)>>1);PolyLn(b,b1,a1,m);memcpy(a1,a,4*m);
    R n=Init(m);NTT(b,n,1);NTT(a1,n,1);NTT(b1,n,1);
    for(R i=0;i<n;++i)b[i]=(LL)b[i]*(YL+1+a1[i]-b1[i])%YL;
    NTT(b,n,-1);memset(a1,0,4*n);memset(b1,0,4*n);memset(b+m,0,4*(n-m));
}
void PolyDiv(R*A,R*a,R*b,R*A1,R*a1,R M,R m){
    PolyRev(a,a1,m);PolyInv(a1,b,A1,M-m+1);PolyRev(A,A1,M);
    R n=Init(M);NTT(b,n,1);NTT(A1,n,1);
    for(R i=0;i<n;++i)b[i]=(LL)b[i]*A1[i]%YL;
    NTT(b,n,-1);memset(A1,0,4*n);memset(a1,0,4*n);
    reverse(b,b+M-m+1);memset(b+M-m+1,0,4*(n-M+m-1));
}
void PolyMod(R*A,R*a,R*b,R*A1,R*a1,R M,R m){
    PolyDiv(A,a,b,A1,a1,M,m);memcpy(A1,A,4*M);memcpy(a1,a,4*m);
    R n=Init(M);NTT(b,n,1);NTT(A1,n,1);NTT(a1,n,1);
    for(R i=0;i<n;++i)b[i]=Mod(A1[i]-(LL)b[i]*a1[i]%YL+YL);
    NTT(b,n,-1);memset(A1,0,4*n);memset(a1,0,4*n);
}
相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程