zoukankan      html  css  js  c++  java
  • 多项式算法合集

    头函数

    #define poly vector<int>
    #define bg begin
    #define pb push_back
    const int mod=998244353,g=3;
    inline int add(int a,int b){return a+b>=mod?a+b-mod:a+b;}
    inline void Add(int &a,int b){a=add(a,b);}
    inline int dec(int a,int b){return a>=b?a-b:a-b+mod;}
    inline void Dec(int &a,int b){a=dec(a,b);}
    inline int mul(int a,int b){return 1ll*a*b>=mod?1ll*a*b%mod:a*b;}
    inline void Mul(int &a,int b){a=mul(a,b);}
    inline int ksm(int a,int b,int res=1){for(;b;b>>=1,a=mul(a,a))(b&1)&&(res=mul(res,a));return res;}
    

    多项式加减点乘点除

    幼儿园小朋友应该都会了吧

    inline poly operator +(poly a,poly b){
    	poly c;int lim=max(a.size(),b.size());c.resize(lim);
    	a.resize(lim),b.resize(lim);
    	for(int i=0;i<lim;i++)c[i]=add(a[i],b[i]);return c;
    }
    inline poly operator -(poly a,poly b){
    	poly c;int lim=max(a.size(),b.size());c.resize(lim);
    	a.resize(lim),b.resize(lim);
    	for(int i=0;i<lim;i++)c[i]=dec(a[i],b[i]);return c;
    }
    inline poly operator *(poly a,int b){
    	for(int i=0;i<a.size();i++)Mul(a[i],b);return a;
    }
    inline poly operator /(poly a,int b){
    	for(int i=0,inv=ksm(b,mod-2);i<a.size();i++)Mul(a[i],inv);
    	return a;
    }
    

    多项式乘法

    FFT:

    前置

    多项式的点值和系数表示法:

    对于一个nn次多项式f(x)f(x)
    f(x)=a0+a1x+a2x2+a3x3anxnf(x)=a_0+a_1x+a_2x^2+a_3x^3……a_nx^n被称作该多项式的系数表示
    我们可以通过带任意一个xx都可以的到唯一的一个f(x)f(x)
    oioi中一般xx一般都只是一个不定元,不会带入特定值计算,比如用作表示下标之类的

    而如果我们把不同的n+1n+1带入进去得到的n+1n+1个点值叫做点值表示法
    n+1n+1个点也可以还原出唯一一个nn次多项式


    虚数

    i1i,sqrt{-1}
    考虑在平面直角坐标系内
    yy轴用ii来表示

    复数更准确的定义是

    eik=cos(k)+isin(k)e^{ik}=cos(k)+isin(k)

    这样平面上一个点就是(a,bi)(a,bi)的形式
    而2个虚数相乘,就对应平面直角坐标系上2个向量模长相乘,极角相加

    由于C++C++自带的复数很慢
    所以我们一般手写复数结构体

    const double pi=acos(-1);
    struct complex{
        double x,y;
        complex(double _x=0,double _y=0):x(_x),y(_y){}
        friend inline complex operator +(const complex &a,const complex &b){
            return complex(a.x+b.x,a.y+b.y);
        }
        friend inline complex operator -(const complex &a,const complex &b){
            return complex(a.x-b.x,a.y-b.y);
        }
        friend inline complex operator /(const complex &a,const double &b){
            return complex(a.x/b,a.y/b);
        }
        friend inline complex operator *(const complex &a,const complex &b){
            return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
        }
    };
    

    单位根

    fftfft的根本
    nn次单位根指的是满足wn=1w^n=1的复数
    nn次单位根有nn个,分别表示为wnkk=0,1,2,n1w_{n}^{k},k=0,1,2,……n-1
    实际上对应的将平面单位圆周上的nn个点
    这些点将单位圆周均分成nn块,且构成一个正nn边形

    更准确的表示为wnk=e2πik/nw_n^k=e^{2pi ik/n}
    wnk=cos(2πk/n)+isin(2πk/n)w_n^k=cos(2pi k/n)+isin(2pi k/n)表示
    即单位圆上的nn

    单位根几个重要的性质:

    1、消去引理

    对于任何整数n0,k0,d0ngeq0,kgeq 0,dgeq 0
    wdndk=wnkw_{dn}^{dk}=w_{n}^{k}

    证明:wdndk=e2πidk/dn=e2πik/n=wnkw_{dn}^{dk}=e^{2pi idk/dn}=e^{2pi ik/n}=w_{n}^{k}

    2、折半引理

    如果n0ngeq 0nn为偶数,那么
    对所有nn次单位根平方,得到的集合是n/2n/2n/2n/2次单位根的集合
    说白了就是wnk+n2=wnkw_{n}^{k+frac n 2}=-w_n^{k}或者就是wnn2=1w_n^{frac n 2}=-1
    证明:画个单位圆,n2frac n 2就是旋转了180°180°,自然取反

    3、求和引理

    对于任意整数n1ngeq 1n∤kn ot | k
    j=0n1(wnk)j=0sum_{j=0}^{n-1}(w_{n}^k)^j=0

    考虑这其实是一个等比数列,wnk1w_{n}^{k} ot= 1
    原式=(wnk)n1wnk1=(wnn)k1wnk1=0=frac{(w_{n}^{k})^{n}-1}{w_{n}^{k}-1}=frac{(w_{n}^n)^k-1}{w_{n}^{k}-1}=0

    而当nkn|k时,原式显然为nn

    这个性质在后面很重要

    算法

    考虑对于两个多项式
    A(x)=a0+a1x+a2x2+a3x3anxnA(x)=a_0+a_1x+a_2x^2+a_3x^3……a_nx^n
    B(x)=b0+b1x+b2x2+b3x3bmxmB(x)=b_0+b_1x+b_2x^2+b_3x^3……b_mx^m
    我们要求一个n+mn+m次多项式C(x)=A(x)B(x)C(x)=A(x)*B(x)
    更具体的C(x)C(x)满足ci=j=0iajbijc_i=sum_{j=0}^{i}a_j*b_{i-j}
    也就是2个数列倒着乘的和,所谓的卷积

    考虑如果我们直接拆开暴力做是O(n2)O(n^2)
    当然也有一种分治乘法可以做到nlog23n^{log_23}(大概就是大常数nnnsqrt n)

    FFTFFT可以做到O(nlogn)O(nlogn)求出CC


    下面为了方便假设n=mn=mnmn ot =m也没有区别

    考虑如果我们有n+1n+1个值x1,xn+1x_1,……,x_{n+1}
    如果已经求出A(x1)A(xn+1)A(x_1),……A(x_{n+1})
    B(x1),B(xn+1)B(x_1),……B(x_{n+1})
    也就是分别求出x1,xn+1x_1,……,x_n+1的点值

    那么我们可以直接在O(n)O(n)的时间内求出
    C(x1)=A(x1)B(x1)C(xn+1)=A(xn+1)B(xn+1)C(x_1)=A(x_1)*B(x_1),……,C(x_{n+1})=A(x_{n+1})*B(x_{n+1})

    现在我们考虑这样

    先对A(x)A(x),B(x)B(x)求出n+1n+1个点的点值,乘起来得到C(x)C(x)的点值
    又对于一个nn次多项式,如果我们知道其n+1n+1个不同点值
    就一定可以还原出一个唯一的多项式
    于是最后再由点值表示还原出原来的C(x)C(x)

    注意由于实际乘出来C(x)C(x)最高系数是2n12n-1
    所以我们需要带入2n2n个点求值
    于是我们会将A(x),B(x)A(x),B(x)补0到2n2n次项
    实际上由于fftfft的特殊性
    我们会将项数补充到2的整数次幂次

    当然非二的整数次幂项的多项式乘法也是可以做的(见下面补充的混合基fftfftBluesteinBluestein

    第一步将系数转成点值是正变换,称作DFTDFT,单次复杂度O(nlogn)O(nlogn)
    由点值还原系数为逆变换,称作IDFTIDFT,单次复杂度O(nlogn)O(nlogn)
    于是总复杂度就是O(nlogn)O(nlogn)
    整个操作被称为快速傅里叶变换(FFT)(FFT)


    DFT:

    考虑我们带入nn次单位根wnw_n
    A(wnk)=i=0nai(wnk)iA(w_n^k)=sum_{i=0}^na_i(w_n^k)^i
    考虑将下标按照奇偶分类
    A(wnk)=i=0n21a2i(wnk)2i+i=0n21a2i+1(wnk)2i+1A(w_n^k)=sum_{i=0}^{frac n 2-1}a_{2i}(w_n^{k})^{2i}+sum_{i=0}^{frac n2 -1}a_{2i+1}(w_{n}^k)^{2i+1}

    =i=0n21a2i(wn2ki)+wnki=0n21a2i+1(wn2ki)=sum_{i=0}^{frac n 2-1}a_{2i}(w_n^{2ki})+w_{n}^ksum_{i=0}^{frac n2 -1}a_{2i+1}(w_{n}^{2ki})

    =i=0n21a2i(wn2k)+wnki=0n21a2i+1(wn2k)=sum_{i=0}^{frac n 2-1}a_{2i}(w_{frac n 2}^{k})+w_n^ksum_{i=0}^{frac n2 -1}a_{2i+1}(w_{frac n 2}^k)

    =Ao(wn2k)+wnkA1(wn2k)=A_o(w_{frac n 2}^k)+{w_{n}^k}A_1(w_{frac n2 }^k)

    这样A0A_0A1A_1都只有n2frac n2项了
    可以继续递归去做
    尽管现在复杂度并没有变化

    考虑对于A(wnk+n2)A(w_n^{k+frac n 2})
    A(wnk+n2)=i=0n21a2i(wnk+n2)2i+wnk+n2i=0n21a2i+1(wnk+n2)2iA(w_n^{k+frac n 2})=sum_{i=0}^{frac n 2-1}a_{2i}(w_n^{k+frac n 2})^{2i}+w_{n}^{k+frac n 2}sum_{i=0}^{frac n2 -1}a_{2i+1}(w_{n}^{k+frac n2 })^{2i}

    考虑单位根的消去引理

    wnk+n2=wnkw_{n}^{k+frac n 2}=-w_n^k

    A(wnk+n2)=i=0n21a2i(wnk)2i+(wnk)i=0n21a2i+1(wnk)2iA(w_n^{k+frac n 2})=sum_{i=0}^{frac n 2-1}a_{2i}(-w_n^k)^{2i}+(-w_n^k)sum_{i=0}^{frac n2 -1}a_{2i+1}(-w_n^k)^{2i}

    =i=0n21a2i(wnk)2iwnki=0n21a2i+1(wnk)2i=sum_{i=0}^{frac n 2-1}a_{2i}(w_n^k)^{2i}-w_n^ksum_{i=0}^{frac n2 -1}a_{2i+1}(w_n^k)^{2i}

    =i=0n21a2i(wn2k)wnki=0n21a2i+1(wn2k)=sum_{i=0}^{frac n 2-1}a_{2i}(w_{frac n 2}^{k})-w_n^ksum_{i=0}^{frac n2 -1}a_{2i+1}(w_{frac n 2}^k)

    =Ao(wn2k)wnkA1(wn2k)=A_o(w_{frac n 2}^k)-{w_{n}^k}A_1(w_{frac n2 }^k)

    我们发现唯一不同的就是第二项的符号
    也就是说如果我们知道Ao(wn2k)A_o(w_{frac n 2}^k)A1(wn2k)A_1(w_{frac n2 }^k)
    我们就可以同时知道A(wnk+n2)A(w_n^{k+frac n 2})A(wnk)A(w_n^k)

    考虑对于k[0,n1],kin[0,n-1],我们求A(wnk)A(w_n^k)
    就只需要知道k[0,n21],A0(wn2k)kin[0,frac n2-1 ],A_0(w_{frac n2 }^k)A1(wn2k)A_1(w_{frac n 2}^k)
    就可以在O(n)O(n)的时间内得到AA

    A0A_0,A1A_1系数都只有n2frac n 2个,所以规模只有原来的一般
    递归求解即可

    时间复杂度T(n)=2T(n2)+O(n)=O(nlogn)T(n)=2*T(frac n2 )+O(n)=O(nlogn)

    这里也是之所以要把项数补充到2k2^k
    因为每一次都把nn项分成n2frac n 2
    如果nn是奇数,那就没法分开了

    由于递归常数比较大,而一般fftfft有关的题时间瓶颈就在DFTDFT上面
    DFTDFT不知道为什么 常数也很大
    写的差的fftfft甚至可以跑1e61e6的数据TT

    所以我们考虑迭代做
    由于每个数最终在的位置和原来不一样
    所以我们要预处理出最终的位置上

    据说找规律得到了O(n)O(n)预处理的方法
    代码如下:
    没看懂,选择全文背诵
    当然也可以nlognnlogn模拟最终位置

    int rev[N<<2];
    inline void init(int lim){
    	for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)*(lim>>1));
    }
    

    我们先把每个数放到最终应该在的地方然后一步步迭代回去就是了

    inline void DFT(complex f[],int lim){
        for(int i=0;i<lim;i++)if(rev[i]>i)swap(f[i],f[rev[i]]);
        for(int mid=1;mid<lim;mid<<=1){
            complex now=plx(cos(pi/mid),sin(pi/mid));
            for(int i=0;i<lim;i+=(mid<<1)){
                complex w=plx(1,0);
                for(int j=0;j<mid;j++,w=w*now){
                    plx p0=f[i+j],p1=w*f[i+j+mid];
                    f[i+j]=p0+p1,f[i+j+mid]=p0-p1;
                }
            }
        }
    }
    

    IDFT:

    以下摘抄自miskcoomiskcoo神犇markdownmarkdown太难写了QAQ)

    在这里插入图片描述
    I是对角矩阵,即对角线上都是1,其他都是0


    代码实现

    const double pi=acos(-1);
    struct plx{
        double x,y;
        plx(double _x=0,double _y=0):x(_x),y(_y){}
        friend inline plx operator +(const plx &a,const plx &b){
            return plx(a.x+b.x,a.y+b.y);
        }
        friend inline plx operator -(const plx &a,const plx &b){
            return plx(a.x-b.x,a.y-b.y);
        }
        friend inline plx operator /(const plx &a,const double &b){
            return plx(a.x/b,a.y/b);
        }
        friend inline plx operator *(const plx &a,const plx &b){
            return plx(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
        }
    };
    inline void fft(plx f[],int lim,int kd){//kd表示在做正变换还是逆变换
        for(int i=0;i<lim;i++)if(rev[i]>i)swap(f[i],f[rev[i]]);
        for(int mid=1;mid<lim;mid<<=1){
            plx now=plx(cos(pi/mid),kd*sin(pi/mid));
            for(int i=0;i<lim;i+=(mid<<1)){
                plx w=plx(1,0);
                for(int j=0;j<mid;j++,w=w*now){
                    plx p0=f[i+j],p1=w*f[i+j+mid];
                    f[i+j]=p0+p1,f[i+j+mid]=p0-p1;
                }
            }
        }
        if(kd==-1)for(int i=0;i<lim;i++)f[i]=f[i]/lim;
    }
    

    NTT:

    由于fftfft涉及复数和实数运算,实际会出现精度误差
    在整数运算时难免会出锅

    所以我们考虑一个能在模意义下的变换
    这就是nttntt

    首先引入原根的概念

    对于一个素数pp,其原根gg定义为满足g0,g1,gp2g^0,g^1,……,g^{p-2}互不相同的数
    又由于费马定理,对一个素数pp,有ap11(mod p)a^{p-1}equiv 1(mod p)
    这个和单位根很相似

    我们考虑单位根之所以能够做fftfft,是因为ww的三个性质
    考虑如果我们对于一个素数p=2nk+1p=2^n*k+1,我们令gn=gkg_n=g^k
    这样我们就可以满足gn0,gn1,gn2gnn1g_n^0,g_n^1,g_n^2……g_n^{n-1}互不相同且满足这些性质(不想证了)

    但这也就限制了nttntt的模数必须是k2n+1k*2^n+1的形式
    否则就要做任意模数nttntt
    做法就是选出几个模数分别做之后CrtCrt合并答案(并不会)
    代码和fftfft类似
    由于wnk=wnnkw_n^{-k}=w_n^{n-k}
    所以我们可以先做的时候不管正逆变换
    最后把11~n1n-1反转一下就可以了

    代码实现

    inline void ntt(poly &f,int lim,int kd){
        for(int i=0;i<lim;i++)if(i<rev[i])swap(f[i],f[rev[i]]);
        for(int mid=1;mid<lim;mid<<=1){
            int now=ksm(g,(mod-1)/(mid<<1));
            for(int i=0;i<lim;i+=(mid<<1)){
                int w=1;
                for(int j=0;j<mid;j++,w=mul(w,now)){
                    int a0=f[i+j],a1=mul(w,f[i+j+mid]);
                    f[i+j]=add(a0,a1),f[i+j+mid]=dec(a0,a1);
                }
            }
        }
        if(kd==-1&&(reverse(f.bg()+1,f.bg()+lim),1))for(int i=0,inv=ksm(lim,mod-2);i<lim;i++)f[i]=mul(f[i],inv);
    }
    

    由于nttntt中每次乘起来很耗费常数
    我也不知道为什么,但就是很耗时间

    于是我们可以预处理出原根优化常数

    预处理原根

    const int N=1000005,C=21;
    int *w[22];
    inline void init_w(){
    	for(int i=1;i<=C;i++)
    	w[i]=new int[1<<(i-1)];
    	int wn=ksm(g,(mod-1)/(1<<C));
    	w[C][0]=1;
    	for(int i=1;i<(1<<(C-1));i++)w[C][i]=mul(w[C][i-1],wn);
    	for(int i=C-1;i;i--)
    	for(int j=0;j<(1<<(i-1));j++)
    	w[i][j]=w[i+1][j<<1];
    }
    

    速度比不预处理快了差不多25%25\%45%45\%不等

    其实nttntt本身常数不算很大,运算常数大概也只有5、6左右
    不过下标不连续可能会导致慢一些

    inline void ntt(poly &f,int lim,int kd){
    	for(int i=0;i<lim;i++)if(i>rev[i])swap(f[i],f[rev[i]]);
    	for(int mid=1,l=1;mid<lim;mid<<=1,l++)
    	for(int i=0;i<lim;i+=(mid<<1))
    		for(int j=0,a0,a1;j<mid;j++)
    			a0=f[i+j],a1=mul(f[i+j+mid],w[l][j]),
    			f[i+j]=add(a0,a1),f[i+j+mid]=dec(a0,a1);
    	if(kd==-1&&(reverse(f.begin()+1,f.begin()+lim),1))
    		for(int inv=ksm(lim,mod-2),i=0;i<lim;i++)Mul(f[i],inv);
    }
    

    乘法

    可以在比较小的时候暴力加循环展开 优化常数

    inline poly operator *(poly a,poly b){
    	int deg=a.size()+b.size()-1,lim=1;
    	if(deg<=128){
    		poly c(deg,0);
    		for(int i=0;i<a.size();i++)
    		for(int j=0;j<b.size();j++)
    			Add(c[i+j],mul(a[i],b[j]));
    		return c;
    	}
    	while(lim<deg)lim<<=1;init(lim);
    	a.resize(lim),ntt(a,lim,1);
    	b.resize(lim),ntt(b,lim,1);
    	for(int i=0;i<lim;i++)Mul(a[i],b[i]);
    	ntt(a,lim,-1),a.resize(deg);
    	return a;
    }
    

    多项式求逆:

    已知一个n1n-1次多项式A(x)A(x),求多项式B(x)B(x)满足:

    A(x)B(x)1(mod xn)A(x)B(x)equiv 1(mod x^n)

    求解过程

    倍增:

    若已知A(x)B(x)1(mod xn2)A(x)B(x)'equiv 1(mod x^{lceil frac n 2 ceil})

    A(x)B(x)A(x)B(x)0(mod xn2)A(x)B(x)'-A(x)B(x)equiv 0(mod x^{lceil frac n 2 ceil})

    B(x)B(x)0(mod xn2)B(x)'-B(x)equiv 0(mod x^{lceil frac n 2 ceil})

    平方:

    B(x)2+B(x)22B(x)B(x)0(mod xn)B(x)'^2+B(x)^2-2B(x)'B(x)equiv 0(mod x^{ n})
    A(x)(B(x)2+B(x)22B(x)B(x))0(mod xn)A(x)(B(x)'^2+B(x)^2-2B(x)'B(x))equiv 0(mod x^{ n})
    B(x)2B(x)A(x)B(x)2(mod xn)B(x)equiv 2B(x)'-A(x)B(x)'^2(mod x^n)

    复杂度T(n)=T(n2)+O(nlogn)=O(nlogn)T(n)=T(frac n 2)+O(nlogn)=O(nlogn)

    注意次数,最高到3倍,开的4倍

    代码实现

    inline poly Inv(poly a,int deg){
    	poly c,b(1,ksm(a[0],mod-2));
    	for(int lim=4;lim<(deg<<2);lim<<=1){
    		init(lim);
    		c=a,c.resize(lim>>1);
    		c.resize(lim),ntt(c,lim,1);
    		b.resize(lim),ntt(b,lim,1);
    		for(int i=0;i<lim;i++)Mul(b[i],dec(2,mul(b[i],c[i])));
    		ntt(b,lim,-1),b.resize(lim>>1);
    	}b.resize(deg);return b;
    }
    

    多项式开方:

    已知一个n1n-1次多项式A(x)A(x),求一个mod xnmod x^n意义下的多项式B(x)B(x)满足

    B(x)2A(x)(mod xn)B(x)^2equiv A(x)(mod x^n)

    满足A[0]=1A[0]=1

    求解过程

    倍增

    首先当n=1n=1时要满足A[0]=1A[0]=1(否则要二次剩余解,老子不会)

    假设已知B(x)2A(x)(mod xn2)B(x)'^2equiv A(x)(mod x^{lceil frac n 2 ceil})
    由于n2lceilfrac n 2 ceil以上的项是不会有影响的
    所以B(x)B(x)(mod xn2)B(x)'equiv B(x) (mod x^{lceil frac n 2 ceil})
    移项平方得:
    B(x)2+B(x)22B(x)B(x)0(mod xn)B(x)^2+B(x)'^2-2B(x)B(x)'equiv 0(mod x^n)

    A(x)+B(x)22B(x)B(x)(mod xn)A(x)+B(x)'^2equiv 2B(x)B(x)'(mod x^n)
    B(x)A(x)2B(x)+B(x)2(mod xn)B(x)equiv frac{A(x)}{2B(x)'}+frac{B(x)}{2}(mod x^n)

    复杂度T(n)=T(n2)+O(nlogn)=O(nlogn)T(n)=T(frac n 2)+O(nlogn)=O(nlogn)

    代码实现

    inline poly Sqrt(poly a,int deg){
    	poly b(1,1),c,d;int inv2=ksm(2,mod-2);
    	for(int lim=4;lim<(deg<<2);lim<<=1){
    		c=a,c.resize(lim>>1);
    		init_rev(lim),d=Inv(b,lim>>1),
    		c=c*d,b.resize(lim);
    		for(int i=0;i<lim;i++)b[i]=mul(inv2,add(b[i],c[i]));
    	}b.resize(deg);return b;
    }
    

    多项式除法和取模:

    给定一个nn次多项式A(x)A(x)和一个mm次多项式B(x)B(x)
    求一个nmn-m次多项式C(x)C(x)m1m-1次多项式D(x)D(x)满足:
    A(x)=B(x)C(x)+D(x)A(x)=B(x)C(x)+D(x)

    求解过程

    考虑对一个最高次数为nn的多项式B(x)B(x)操作BR(x)=xnB(1x)B^R(x)=x^nB(frac 1 x)
    会发现BR(x)B^R(x)只是B(x)B(x)的系数反转的柿子

    考虑A(x)=B(x)C(x)+D(x)A(x)= B(x)C(x)+D(x)
    A(x)A(x)最高项为nn,B(x)B(x)最高项为mm,则C(x)C(x)最高项为nmn-m,D(x)D(x)m1m-1
    两边同时乘一个xnx^n,并带入1xfrac 1 x
    xnA(x)=xmB(x)xnmC(x)+xnm+1xm1D(x)x^nA(x)= x^mB(x)x^{n-m}C(x)+x^{n-m+1}*x^{m-1}D(x)
    AR(x)=BR(x)CR(x)+xnm+1DR(x)A^R(x)=B^R(x)C^R(x)+x^{n-m+1}D^R(x)
    考虑在mod xnm+1mod x^{n-m+1}意义下,A,BA,B已知,CC最高为nmn-m不影响,而DD被消去
    AR(x)BR(x)CR(x)(mod xnm+1)A^R(x)equiv B^R(x)C^R(x)(mod x^{n-m+1})

    多项式求逆就可以了

    复杂度O(nlogn)O(nlogn)

    代码实现

    inline poly operator /(poly a,poly b){
    	int lim=1,deg=a.size()-b.size()+1;
    	reverse(a.bg(),a.end());
    	reverse(b.bg(),b.end());
    	while(lim<=deg)lim<<=1;
    	b=Inv(b,lim),b.resize(deg);
    	a=a*b,a.resize(deg);
    	reverse(a.bg(),a.end());
    	return a;
    }
    inline poly operator %(poly a,poly b){
    	poly c=a-(a/b)*b;
    	c.resize(b.size()-1);
    	return c;
    }
    

    多项式求导与积分

    我怕是自学的是一个假的微积分
    假装自己会的差不多了

    复杂度O(n)O(n)

    代码实现

    inline poly deriv(poly a){
    	for(int i=0;i<a.size()-1;i++)a[i]=mul(a[i+1],i+1);
    	a.pob();return a;
    }
    inline poly integ(poly a){
    	a.pb(0);
    	for(int i=a.size()-1;i;i--)a[i]=mul(a[i-1],inv[i]);
    	a[0]=0;
    	return a;
    }
    

    多项式Ln

    已知一个n1n-1次多项式g(x)g(x),求一个mod xnmod x^n意义下的多项式f(x)f(x),满足:
    f(x)Ln(g(x))(mod xn)f(x)equiv Ln(g(x))(mod x^n)

    求解过程

    由于有公式
    f(x)=ln(g(x))f(x)=ln(g(x))
    g(x)=f(x)g(x)g'(x)=f'(x)g(x)
    f(x)=g(x)g(x)f'(x)=frac{g'(x)}{g(x)}
    且要满足g[0]=1g[0]=1否则老子不会
    求导求逆最后再积分就可以了

    复杂度O(nlogn)O(nlogn)

    代码实现

    /*
    if f(x)=Ln(g(x))
    then g'(x)=f'(x)g(x)
    */
    inline poly Ln(poly a,int lim){
    	a=integ(deriv(a)*Inv(a,lim)),a.resize(lim);
    	return a;
    }
    

    多项式Exp

    前置知识

    泰勒展开
    牛顿迭代

    泰勒展开

    f(x)=f(x0)+f1(x0)1!(xx0)+f2(x0)2!(xx0)2++fn(x0)n!(xx0)n+f(x)=f(x_0)+frac{f^1{(x_0)}}{1!}(x-x_0)+frac{f^2(x_0)}{2!}(x-x_0)^2+……+frac{f^n(x_0)}{n!}(x-x_0)^n+……

    考虑我们要构造一个函数g(x)g(x)使g(x)g(x)完全拟合f(x)f(x)
    那首先初始点x0x_0的值要和f(x0)f(x_0)一样
    在此基础上只需要保证一阶导数,二阶导数……都完全相同即可
    i,fi(x)=gi(x)forall i,f^i(x)=g^i(x)
    由于求ggii阶导数时为i!ai=fi(x)i!a_i=f^i(x)
    ai=fi(x)i!a_i=frac{f^i(x)}{i!}
    所以得证

    牛顿迭代

    在多项式中一般用来解这类问题:

    假设有函数ff和一个多项式g(x)mod xng(x)mod x^n
    满足f(g(x))0(mod xn)f(g(x))equiv 0 ( mod x^n)
    已知ff,求gg

    说白了就是用来解f(g(x))0(mod xn)f(g(x))equiv 0(mod x^n)之类的方程

    首先在n=1n=1的时候即常数时单独求

    假设已经知道f(g(x))0(mod xn2)f(g'(x))equiv 0(mod x^{lceil frac n2 ceil})
    要求f(g(x))0(mod xn)f(g(x))equiv0(mod x^n)

    考虑将f(g(x))f(g(x))g(x)g'(x)处泰勒展开
    f(g(x))=f(g(x))+f1(g(x))1!(g(x)g(x))+f2(g(x))2!(g(x)g(x))2+f(g(x))=f(g'(x))+frac{f^1(g'(x))}{1!}(g(x)-g'(x))+frac{f^2(g'(x))}{2!}(g(x)-g(x)')^2+……

    首先显然有g(x)g(x)(mod xn2)g(x)equiv g'(x)(mod x^{lceil frac n 2 ceil})
    所以g(x)g(x)g(x)-g'(x)最低项次数一定大于n2lceil frac n 2 ceil

    则在mod xnmod x^n意义下,整个式子从(g(x)g(x))2(g(x)-g(x)')^2开始都为00
    f(g(x))f(g(x))+f1(g(x))(g(x)g(x))(mod xn)f(g(x))equiv f(g'(x))+{f^1(g'(x))}(g(x)-g'(x))(mod x^n)

    f(g(x))0(mod xn)f(g(x))equiv 0 (mod x^n)

    g(x)g(x)f(g(x))f1(g(x))(mod xn)g(x)equiv g'(x)-frac{f(g'(x))}{f^1(g'(x))}(mod x^n)

    这就大功告成了

    例:

    比如多项式开根
    就是解B(x)2A(x)0(mod xn)B(x)^2-A(x)equiv 0(mod x^n)
    假设已知B(x)2A(x)0(mod xn2)B'(x)^2-A(x)equiv 0(mod x^{lceil frac n 2 ceil})
    这时候f(B(x))=B(x)2A(x),f(B(x))=B(x)^2-A(x),f1(B(x))=2B(x)f^1(B(x))=2B(x)
    带入B(x)B(x)B(x)2A(x)2B(x)B(x)2+A(x)2B(x)B(x)equiv B'(x)-frac{B'(x)^2-A(x)}{2B'(x)}equiv frac{B'(x)^2+A(x)}{2B'(x)}
    就是我们推出来的式子


    已知一个n1n-1次多项式g(x)g(x),求一个mod xnmod x^n意义下的多项式f(x)f(x),满足:
    f(x)eg(x)(mod xn)f(x)equiv e^{g(x)}(mod x^n)

    也就是Ln(f(x))g(x)(mod xn)Ln(f(x))equiv g(x)(mod x^n)

    求解过程

    倍增:

    考虑原问题,即求解方程Ln(f(x))g(x)0(mod xn)Ln(f(x))-g(x)equiv 0(mod x^n)

    同样假设已经知道Ln(f(x))g(x)0(mod xn2)Ln(f'(x))-g(x)equiv 0(mod x^{lceil frac n 2 ceil})

    P(f(x))=Ln(f(x))g(x)P(f(x))=Ln(f(x))-g(x),则P1(f(x))=1f(x)P^1(f(x))=frac 1 {f(x)}

    f(x)f(x)Ln(f(x))g(x)1f(x)f(x)f(x)Ln(f(x))+f(x)g(x)f(x)equiv f'(x)-frac{Ln(f'(x))-g(x)}{frac{1}{f'(x)}}equiv f'(x)-f'(x)Ln(f'(x))+f'(x)g(x)

    f(x)(1Ln(f(x))+g(x))(mod xn)equiv f'(x)(1-Ln(f'(x))+g(x)) (mod x^n)

    复杂度T(n)=T(n2)+O(nlogn)=O(nlogn)T(n)=T(frac n 2)+O(nlogn)=O(nlogn)

    代码实现

    inline poly exp(poly a,int deg){
    	poly b(1,1),c;int n=a.size();
    	for(int lim=2;lim<(deg<<1);lim<<=1){
    		c=Ln(b,lim);
    		for(int i=0;i<lim;i++)c[i]=dec(i<n?a[i]:0,c[i]);
    		Add(c[0],1),b=b*c;
    		b.resize(lim);
    	}b.resize(deg);return b;
    }
    

    多项式多点求值

    已知一个nn次多项式f(x)f(x),求f(a1)f(am)f(a_1)……f(a_m)

    求解过程

    考虑构造函数P(x)=i=1m(xai)P(x)=prod_{i=1}^{m}(x-a_i)
    显然i[1,m],P(ai)=0forall iin[1,m],P(a_i)=0

    假设f(x)=P(x)g(x)+A(x)f(x)=P(x)g(x)+A(x)
    A(x)=f(x)%P(x)A(x)=f(x)\%P(x)

    那显然f(ai)=A(ai)f(a_i)=A(a_i)
    但由于P(x)P(x)mm次的,没有起到优化的作用

    而考虑对于kP(x)=(xak)k,P(x)=(x-a_k)
    f(x)%P(x)f(x)\% P(x)后就只剩下一个常数项,即f(ai)f(a_i)的值了
    但是这样一次就nlognnlogn

    考虑分治优化
    P1(x)=i=lmid(xai),P2(x)=i=mid+1r(xai)P_1(x)=prod_{i=l}^{mid}(x-a_i),P_2(x)=prod_{i=mid+1}^{r}(x-a_i)
    i[l,mid]f(x)%P1(x)=f(ai)forall iin[l,mid] , f(x)\%P_1(x)=f(a_i)
    取模之后f(x)f(x)的次数减少了一半,继续递归求解即可

    P(x)P(x)可以先分治nttntt预处理出来
    复杂度T(n)=2T(n2)+O(nlogn)=O(nlog2n)T(n)=2*T(frac n 2)+O(nlogn)=O(nlog^2n)

    代码实现

    第一次写TT掉了,预处理了一波单位根就跑过去了
    也可以在比较小的时候暴力秦九韶展开,然并没写

    poly f[N<<2];
    int a[N];
    int n,m;
    #define lc (u<<1)
    #define rc ((u<<1)|1)
    #define mid ((l+r)>>1)
    void build(int u,int l,int r,int *v){
    	if(l==r){f[u].pb(dec(0,v[l])),f[u].pb(1);return;}
    	build(lc,l,mid,v),build(rc,mid+1,r,v);
    	f[u]=f[lc]*f[rc];
    }
    void calc(int u,int l,int r,poly now,int *v){
    	if(l==r){v[l]=now[0];return;}
    	calc(lc,l,mid,now%f[lc],v),calc(rc,mid+1,r,now%f[rc],v);
    }
    #undef lc
    #undef rc
    #undef mid
    

    多项式快速插值

    写在这里了


    下降幂多项式乘法

    写在这里了


    其他技巧

    多项式快速幂

    直接快速幂要多个loglog而且常数大(虽然lnlnexpexp常数一样大死个仙人)
    A[0]=1A[0]=1时(LnLn要求保证这个),A(x)k=Exp(ln(A(x))k)A(x)^k=Exp(ln(A(x))*k)
    洛谷板子读入时取模

    inline poly ksm(poly a,int deg,int k){
    	a=exp(Ln(a,deg)*k,deg),a.resize(deg);
    	return a;
    }
    

    对于A[0]1A[0] ot=1也可以很简单的做
    写在这里了

    MTTMTT

    写在这里了

    模板合集

    const int mod=998244353,g=3;
    inline int add(int a,int b){return a+b>=mod?a+b-mod:a+b;}
    inline void Add(int &a,int b){a=add(a,b);}
    inline int dec(int a,int b){return a>=b?a-b:a-b+mod;}
    inline void Dec(int &a,int b){a=dec(a,b);}
    inline int mul(int a,int b){return 1ll*a*b>=mod?1ll*a*b%mod:a*b;}
    inline void Mul(int &a,int b){a=mul(a,b);}
    inline int ksm(int a,int b,int res=1){for(;b;b>>=1,a=mul(a,a))(b&1)&&(res=mul(res,a));return res;}
    const int N=100005,C=17;
    int *w[18];
    int rev[N<<2];
    #define poly vector<int> 
    #define pb push_back
    #define bg begin
    inline void init_rev(int lim){
    	for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)*(lim>>1));
    }
    inline void init_w(){
    	for(int i=1;i<=C;i++)
    	w[i]=new int[1<<(i-1)];
    	int wn=ksm(g,(mod-1)/(1<<C));
    	w[C][0]=1;
    	for(int i=1;i<(1<<(C-1));i++)w[C][i]=mul(w[C][i-1],wn);
    	for(int i=C-1;i;i--)
    	for(int j=0;j<(1<<(i-1));j++)
    	w[i][j]=w[i+1][j<<1];
    }
    inline void ntt(poly &f,int lim,int kd){
    	for(int i=0;i<lim;i++)if(i>rev[i])swap(f[i],f[rev[i]]);
    	for(int mid=1,l=1;mid<lim;mid<<=1,l++)
    	for(int i=0;i<lim;i+=(mid<<1))
    		for(int j=0,a0,a1;j<mid;j++)
    			a0=f[i+j],a1=mul(f[i+j+mid],w[l][j]),
    			f[i+j]=add(a0,a1),f[i+j+mid]=dec(a0,a1);
    	if(kd==-1&&(reverse(f.bg()+1,f.bg()+lim),1))
    		for(int inv=ksm(lim,mod-2),i=0;i<lim;i++)Mul(f[i],inv);
    }
    inline poly operator +(poly a,poly b){
    	poly c;int lim=max(a.size(),b.size());c.resize(lim);
    	a.resize(lim),b.resize(lim);
    	for(int i=0;i<lim;i++)c[i]=add(a[i],b[i]);return c;
    }
    inline poly operator -(poly a,poly b){
    	poly c;int lim=max(a.size(),b.size());c.resize(lim);
    	a.resize(lim),b.resize(lim);
    	for(int i=0;i<lim;i++)c[i]=dec(a[i],b[i]);return c;
    }
    inline poly operator *(poly a,int b){
    	for(int i=0;i<a.size();i++)Mul(a[i],b);return a;
    }
    inline poly operator /(poly a,int b){
    	for(int i=0,inv=ksm(b,mod-2);i<a.size();i++)Mul(a[i],inv);
    	return a;
    }
    inline poly operator *(poly a,poly b){
    	int deg=a.size()+b.size()-1,lim=1;
    	if(deg<=128){
    		poly c(deg,0);
    		for(int i=0;i<a.size();i++)
    		for(int j=0;j<b.size();j++)
    			Add(c[i+j],mul(a[i],b[j]));
    		return c;
    	}
    	while(lim<deg)lim<<=1;init(lim);
    	a.resize(lim),ntt(a,lim,1);
    	b.resize(lim),ntt(b,lim,1);
    	for(int i=0;i<lim;i++)Mul(a[i],b[i]);
    	ntt(a,lim,-1),a.resize(deg);
    	return a;
    }
    inline poly Inv(poly a,int deg){
    	poly c,b(1,ksm(a[0],mod-2));
    	for(int lim=4;lim<(deg<<2);lim<<=1){
    		init(lim);
    		c=a,c.resize(lim>>1);
    		c.resize(lim),ntt(c,lim,1);
    		b.resize(lim),ntt(b,lim,1);
    		for(int i=0;i<lim;i++)Mul(b[i],dec(2,mul(b[i],c[i])));
    		ntt(b,lim,-1),b.resize(lim>>1);
    	}b.resize(deg);return b;
    }
    inline poly Sqrt(poly a,int deg){
    	poly b(1,1),c,d;int inv2=ksm(2,mod-2);
    	for(int lim=4;lim<(deg<<2);lim<<=1){
    		c=a,c.resize(lim>>1);
    		init_rev(lim),d=Inv(b,lim>>1),
    		c=c*d,b.resize(lim);
    		for(int i=0;i<lim;i++)b[i]=mul(inv2,add(b[i],c[i]));
    	}b.resize(deg);return b;
    }
    inline poly operator /(poly a,poly b){
    	int lim=1,deg=a.size()-b.size()+1;
    	reverse(a.bg(),a.end());
    	reverse(b.bg(),b.end());
    	while(lim<=deg)lim<<=1;
    	b=Inv(b,lim),b.resize(deg);
    	a=a*b,a.resize(deg);
    	reverse(a.bg(),a.end());
    	return a;
    }
    inline poly operator %(poly a,poly b){
    	poly c=a-(a/b)*b;
    	c.resize(b.size()-1);
    	return c;
    }
    inline poly deriv(poly a){
    	for(int i=0;i<a.size()-1;i++)a[i]=mul(a[i+1],i+1);
    	a.pob();return a;
    }
    inline poly integ(poly a){
    	a.pb(0);
    	for(int i=a.size()-1;i;i--)a[i]=mul(a[i-1],inv[i]);
    	a[0]=0;
    	return a;
    }
    inline poly Ln(poly a,int lim){
    	a=integ(deriv(a)*Inv(a,lim)),a.resize(lim);
    	return a;
    }
    inline poly exp(poly a,int deg){
    	poly b(1,1),c;int n=a.size();
    	for(int lim=2;lim<(deg<<1);lim<<=1){
    		c=Ln(b,lim);
    		for(int i=0;i<lim;i++)c[i]=dec(i<n?a[i]:0,c[i]);
    		Add(c[0],1),b=b*c;
    		b.resize(lim);
    	}b.resize(deg);return b;
    }
    inline poly ksm(poly a,int deg,int k){
    	a=exp(Ln(a,deg)*k,deg),a.resize(deg);
    	return a;
    }
    poly f[N<<2];
    #define mid ((l+r)>>1)
    #define lc (u<<1)
    #define rc ((u<<1)|1)
    inline void build(int u,int l,int r,int *a){
    	if(l==r){
    		f[u].clear();
    		f[u].pb(dec(0,a[l]));
    		f[u].pb(1);return;
    	}build(lc,l,mid,a),build(rc,mid+1,r,a);
    	f[u]=f[lc]*f[rc];
    }
    inline void calc(int u,int l,int r,poly g,int *a){
    	if(l==r){a[l]=g[0];return;}
    	calc(lc,l,mid,g%f[lc],a);
    	calc(rc,mid+1,r,g%f[rc],a);
    }
    inline void getans(poly a,int *b,int num){
    	build(1,1,num,b),calc(1,1,num,a,b);
    }
    #undef mid
    #undef lc
    #undef rc
    
  • 相关阅读:
    AJAX
    前端上传文件 后端PHP获取文件
    PHP基础语法
    JS错误记录
    JS学习笔记
    python利用xlrd读取excel文件始终报错原因
    安装xlwt和xlrd
    编程菜鸟的日记-Linux无处不在
    编程菜鸟的日记-《软件测试》Ron Patton著-读书笔记
    编程菜鸟的日记-初学尝试编程-C++ Primer Plus 第6章编程练习9
  • 原文地址:https://www.cnblogs.com/stargazer-cyk/p/12328281.html
Copyright © 2011-2022 走看看