zoukankan      html  css  js  c++  java
  • 多项式全家桶

    快速傅里叶变换

    多项式的系数表示法

    [f(x)=sum_{i=0}^{n-1}a_ix^i ]

    [f(x) ightarrow {a_0,a_1,a_2dots ,a_{n-1}} ]

    多项式的点值表示法

    注意到一个最高次数为(n-1)的多项式可以由(n)个点值唯一确定
    即将互不相同的({x_0,x_1,dots ,x_{n-1}})带入(f(x))得到(n)个取值({f(x_0),f(x_1),dots ,f(x_{n-1})})
    因此

    [f(x) ightarrow {f(x_0),f(x_1),dots ,f(x_{n-1})} ]

    多项式乘法

    注意到如果用系数表示法,直接做的复杂度为(mathcal O(n^2))(类似于高精度乘法)
    但如果是点值表示法则可以做到(mathcal O(n))
    即,对于两个多项式

    [f(x) ightarrow {f(x_0),f(x_1),dots ,f(x_{n-1})}\ g(x) ightarrow {g(x_0),g(x_1),dots ,g(x_{n-1})} ]

    那么

    [h(x)=f(x)g(x) ightarrow {f(x_0)g(x_0),f(x_1)g(x_1),dots ,f(x_{n-1})g(x_{n-1})} ]

    然而虽然理想美好,但是现实却很残酷
    一般情况下都是用的系数表示法
    于是自然而然地想到将系数表示法转化为点值表示法,乘出来之后再变回来
    显然,暴力得到点值仍然是(mathcal O(n^2))
    于是我们考虑一些特殊的值带入

    复数中的单位根

    注意以下的(n)都是(2)的整数次幂
    对于方程

    [x^n=1 ]

    由代数基本定理,这个方程应该存在(n)个根
    不妨设这(n)个根为(w_n^0,w_n^1,dots ,w_n^{n-1})
    容易发现

    [w_n^k=cos(frac{2pi}nk)+icdot sin(frac{2pi}nk) ]

    可以结合复平面上的单位圆进行理解
    然后单位根有如下神奇性质

    [w_n^a=w_{2n}^{2a}\ w_n^a=-w_n^{a+frac n2}\ w_n^a=w_n^{a-n} ]

    然后就可以开始搞事情了

    进入正题

    注意以下的(n)都是(2)的整数次幂
    不妨设

    [f(x)=a_0+a_1x+a_2x^2dots +a_{n-1}x^{n-1}\ A(x)=a_0+a_2x+a_4x^2dots +a_{n-2}x^{frac{n-2}2}\ B(x)=a_1+a_3x+a_5x^2dots +a_{n-1}x^{frac{n-2}2} ]

    即按照奇偶性分类
    那么我们有

    [f(x)=A(x^2)+xcdot B(x^2) ]

    然后可以注意到

    [f(w_n^k)=A(w_n^{2k})+w_n^kcdot B(w_n^{2k})\ =A(w_{frac n2}^k)+w_n^kcdot B(w_{frac n2}^k) ]

    [f(w_n^{k+frac n2})=A(w_n^{2k+n})+w_n^{k+frac n2}cdot B(w_n^{2k+n})\ =A(w_n^{2k})+w_n^{k+frac n2}cdot B(w_n^{2k})\ =A(w_{frac n2}^k)-w_n^kcdot B(w_{frac n2}^k) ]

    于是我们在计算(f(w_n^k))时就可以顺便计算出(f(w_n^{k+frac n2}))的值

    这样就可以把(mathcal O(n^2))优化到(mathcal O(nlog_2n))

    现在的问题就在于如何把点值再变回系数

    考虑我们现在已经求出(h(x)=f(x)*g(x))(w_n^0,w_n^1,dots w_n^{n-1})处的点值,需要求出(h(x))的各项系数(b_i)

    考虑如此构造函数

    [C(x)=sum_{i=0}^{n-1}h(w_n^i)x^i ]

    那么

    [C(x)=sum_{i=0}^{n-1}sum_{j=0}^{n-1}b_j(w_n^i)^jx^i ]

    考虑将(x=w_n^{-0},w_n^{-1},dots w_n^{-(n-1)})带入(C(x))

    [C(w_n^{-k})=sum_{i=0}^{n-1}sum_{j=0}^{n-1}b_j(w_n^i)^j(w_n^{-k})^i\ =sum_{i=0}^{n-1}sum_{j=0}^{n-1}b_j(w_n^{j-k})^i\ =sum_{j=0}^{n-1}b_jsum_{i=0}^{n-1}(w_n^{j-k})^i ]

    考虑单位根反演定理

    [[n|k]=frac 1nsum_{i=0}^{n-1}(w_n^k)^i ]

    考虑证明此结论

    (nmid k)时,显然(w_n^k=1),因此原式(=frac 1nsum_{i=0}^{n-1}1^i=1)

    (n mid k)时,原式(=frac 1ncdot frac{(w_n^k)^n-1}{w_n^k-1}=frac 1ncdot frac{w_n^{kn}-1}{w_n^k-1}=frac 1ncdot frac{1-1}{w_n^k-1}=0)

    因此

    [C(w_n^{-k})=sum_{j=0}^{n-1}b_jn[nmid j-k]=nb_k ]

    所以我们只用求出(C(x))(x=w_n^{-0},w_n^{-1},dots w_n^{-(n-1)})的所有点值,最后再(/n),就可以计算出(h(x))的各项系数了

    于是我们就在(mathcal O(nlog_2n))的时间内完成了多项式乘法

    但是倘若直接按照这样递归模拟,常数极大

    我们考虑进行优化,直接把每个数交换到对应的位置上

    img

    可以发现:交换后序列的位置为 交换前序列的二进制下 翻转过来

    可以通过如下代码进行实现

    for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
    

    在函数中这样实现

    for(int i=0;i<lim;++i)if(i<rev[i])swap(f[i],f[rev[i]]);
    

    这样就可以化递归为迭代,大大减小常数

    还有其他优化:比如预处理单位根。比较trivial就不再赘述了

    完整代码

    inline void fft(pt*p,int lim,int tp)
    {
    	for(int i=0;i<lim;++i)if(i<rev[i])swap(p[i],p[rev[i]]);
    	for(int mid=1;mid<lim;mid<<=1)
    	{
    		pt w=pt(1.0,0.0),wn=pt(cos(pi/mid),1.0*tp*sin(pi/mid));
    		for(int len=mid<<1,j=0;j<lim;j+=len,w=pt(1.0,0.0))
    			for(int k=0;k<mid;++k,w=w*wn)
    			{
    				pt x=p[j+k],y=w*p[j+mid+k];
    				p[j+k]=x+y,p[j+mid+k]=x-y;
    			}
    	}
    }
    

    其中(tp=1)表示为快速傅里叶变换,(tp=-1)表示快速傅里叶逆变换

    快速数论变换

    存在意义

    普通的(FFT)显然无法支持取模操作

    同时由于运用浮点数进行运算,包括大量使用如(sin,cos)这些函数,精度难免有误差

    于是快速数论变换应运而生

    (a,p)互素,且 (p>1)

    对于 (a^n≡1(mod p))最小(n),我们称之为 (a)(p)的阶,记做 (δ_p(a))

    例如: (δ_7(2)=3)

    原根

    (p)是正整数,(a)是整数,若 (δ_p(a))等于(varphi (p)),则称 (a)为模 (p)的一个原根,记为 (g)

    (δ_7(3)=6=varphi (7)),因此(3)(7)的一个原根

    注意原根的个数是不唯一的

    原根有一个极其重要的性质

    (P)为素数,假设一个数(g)(P)的原根,那么(g^i pmod P (0<i<P))的结果两两不同

    进入正题

    下证

    [w_n^kequiv(g^{frac {p-1}n})^kpmod p ]

    考虑单位根的所有性质

    • (w_n^k=w_{2n}^{2k})

      (Leftrightarrow (g^{frac {p-1}n})^kequiv (g^{frac {p-1}{2n}})^{2k}pmod p)显然成立

    • (w_n^k=-w_n^{k+frac n2})

      (Leftrightarrow (g^{frac {p-1}n})^kequiv -(g^{frac {p-1}n})^{k+frac n2})

      (Leftrightarrow g^{frac {p-1}{2}}equiv -1)

    • (w_n^k=w_n^{k-n})

      (Leftrightarrow (g^{frac {p-1}n})^kequiv (g^{frac {p-1}n})^{k-n})

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

    因此我们就可以用原根来替换单位根

    常用的模数有(998244353,1004535809,469762049),它们的原根都是(3)

    其他地方都和(FFT)一模一样

    代码如下

    inline void ntt(int lim,poly&f,bool tp)
    {
    	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 sw=qpow(tp?ivG:G,(mod-1)/(mid<<1));
    		wn[0]=1;for(int i=1;i<mid;++i)wn[i]=mll(wn[i-1],sw);
    		for(int len=mid<<1,p=0;p+len-1<lim;p+=len)
    			for(int k=0;k<mid;++k)
    			{
    				int x=f[p+k],y=1ll*wn[k]*f[p+mid+k]%mod;
    				f[p+k]=add(x,y),f[p+mid+k]=dec(x,y);
    			}
    	}
    }
    

    多项式求逆

    给出(n-1)次函数(f(x)),求函数(g(x))满足

    [g(x)equivfrac 1{f(x)}pmod {x^n} ]

    系数均对(998244353)取模

    考虑倍增法

    不妨设我们已经求得函数(g_0)满足

    [g_0cdot fequiv 1pmod {x^{frac n2}} ]

    又由于

    [gcdot fequiv 1pmod {x^{frac n2}} ]

    所以

    [fcdot (g-g_0)equiv 0pmod {x^{frac n2}}\ g-g_0equiv 0pmod {x^{frac n2}}\ (g-g_0)^2equiv 0pmod {x^{n}}\ g^2-2gcdot g_0+g_0^2equiv 0pmod {x^{n}}\ fcdot (g^2-2gcdot g_0+g_0^2)equiv 0pmod {x^{n}}\ g-2g_0+fcdot g_0^2equiv 0pmod {x^{n}}\ gequiv 2g_0-fcdot g_0^2pmod {x^{n}}\ ]

    特别地,当(f(x)=c)时,(g(x)=c^{p-1})

    然后递归求解即可

    注意此处常数的优化

    poly ginv(int n,poly f)
    {
    	if(n==1){poly g;g[0]=qpow(f[0],mod-2);return g;}
    	poly g=ginv(n+1>>1,f),h,p=g;
    	int lim=1,l=0;while(lim<n+n)lim<<=1,++l;
    	for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
    	f.pre(n,lim),p.pre(n,lim);
    	ntt(lim,f,0),ntt(lim,p,0);
    	for(int i=0;i<lim;++i)f[i]=mll(f[i],mll(p[i],p[i]));
    	ntt(lim,f,1);int iv=qpow(lim,mod-2);
    	for(int i=0;i<n;++i)h[i]=dec(add(g[i],g[i]),mll(f[i],iv));
    	return h;
    }
    

    多项式( ext{ln})

    给出(n-1)次的多项式函数(f(x)),求函数(g(x))满足

    [g(x)equivln f(x)pmod {x^n} ]

    系数均对(998244353)取模,保证(f(0)=1)

    两边同时求导可以得到

    [g'equiv(ln f)'pmod {x^n}\ g'equivfrac {f'}fpmod {x^n}\ gequiv int frac {f'}fpmod {x^n} ]

    于是多项式求逆+多项式求导+多项式积分即可

    来康康具体实现

    inline void igl(int n,poly&f)
    {
    	for(int i=n-1;i;--i)f[i]=mll(f[i-1],inv[i]);f[0]=0;
    }//积分
    inline void diff(int n,poly&f)
    {
    	for(int i=0,t=f[i+1];i<n-1;++i,t=f[i+1])f[i]=mll(t,i+1);f[n-1]=0;
    }//微分
    inline poly ln(int n,poly f)
    {
    	poly g=ginv(n,f),h;diff(n,f);
    	h=mul(n,n,g,f);igl(n,h);
    	return h;
    }
    

    多项式牛顿迭代

    对于函数方程

    [F(G(x))equiv 0pmod {x^n} ]

    已知(F(x)),求(G(x))

    类似于倍增法

    假设我们已经求出(G_0(x))满足

    [F(G_0(x))equiv 0pmod {x^{frac n2}} ]

    考虑在(G_0(x))处使用泰勒展开

    [F(G(x))=sum_{i=0}^{+infty}frac{F^{(i)}(G_0(x))}{i!}(G(x)-G_0(x))^i ]

    其中(f^{(i)})表示对(f)(i)阶导

    注意这个地方是把(G(x))当做自变量,(G_0(x))当作常数,因此不用考虑链式法则等因素

    容易知道

    [G(x)-G_0(x)equiv 0pmod {x^{frac n2}} ]

    因此对于(forall ige2)

    [(G(x)-G_0(x))^iequiv 0pmod {x^n} ]

    所以

    [F(G(x))equiv F(G_0(x))+F'(G_0(x))(G(x)-G_0(x))pmod {x^n}\ Leftrightarrow F(G_0(x))+F'(G_0(x))G(x)-F'(G_0(x))G_0(x)equiv 0pmod {x^n}\ Leftrightarrow G(x)equiv G_0(x)-frac{F(G_0(x))}{F'(G_0(x))}pmod {x^n} ]

    于是这样不断迭代下去即可求出答案

    但要注意当(n=1)时需要特殊处理(即递归边界)

    那这个东西有什么用呢?以下举两个例子

    多项式(exp)

    给出(n-1)次多项式(f(x)),求(g(x))满足

    [g(x)equiv e^{f(x)}pmod {x^n} ]

    系数均对(998244353)取模,保证(f(0)=0)

    化一下柿子

    [ln g(x)equiv f(x)pmod {x^n}\ Leftrightarrow ln g(x)-f(x)equiv 0pmod {x^n} ]

    带入先前的牛顿迭代表达式可得

    [g(x)equiv g_0(x)-frac{ln g_0(x)-f(x)}{frac 1{g_0(x)}}pmod {x^n}\ g(x)equiv g_0(x)(1-ln g_0(x)+f(x))pmod {x^n}\ ]

    迭代即可

    递归边界(g(x)equiv 1pmod x),因为题目保证函数(f(x))的常数项为(0)

    代码如下

    poly exp(int n,poly f)
    {
    	if(n==1){poly g;g[0]=1;return g;}
    	poly g=exp(n+1>>1,f),h=ln(n,g);
    	for(int i=0;i<n;++i)h[i]=add(dec(0,h[i]),f[i]);
        h[0]=add(h[0],1);
    	return mul(n,n,g,h);
    }
    

    多项式开根

    给出(n-1)次多项式(f(x)),求(g(x))满足

    [g(x)equiv sqrt{f(x)}pmod {x^n} ]

    系数均对(998244353)取模

    仍然推柿子

    [g^2(x)equiv f(x)pmod {x^n}\ g^2(x)-f(x)equiv 0pmod {x^n} ]

    还是带入牛顿迭代表达式可得

    [g(x)equiv g_0(x)-frac{g_0^2(x)-f(x)}{2g_0(x)}pmod {x^n} ]

    多项式乘法+多项式求逆即可

    递归边界(g(x)equivsqrt {f(0)} pmod {998244353})

    于是求出(f(x))常数项的二次剩余即可

    代码如下

    inline poly sqr(int n,poly f)
    {
    	if(n==1){poly g;g[0]=solve(f[0]);return g;}
    	poly g=sqr(n+1>>1,f),h,ivg;
    	ivg=ginv(n,g),h=mul(n,n,g,g);
    	int iv=qpow(2,mod-2);
    	for(int i=0;i<n;++i)ivg[i]=mll(iv,ivg[i]);
    	for(int i=0;i<n;++i)h[i]=dec(h[i],f[i]);
    	h=mul(n,n,h,ivg);
    	for(int i=0;i<n;++i)h[i]=dec(g[i],h[i]);
    	return h;
    }
    

    多项式除法+取模

    这个稍微trivial一点

    给出(n-1)次多项式(f(x))(m-1)次多项式(g(x)),求(q(x),r(x))满足

    [f(x)=q(x)g(x)+r(x) ]

    要求(r(x))的次数小于(m-1)
    系数均对(998244353)取模
    对于多项式

    [F(x)=sum_{i=0}^{n-1}a_ix^i ]

    我们定义

    [F_R(x)=sum_{i=0}^{n-1}a_{n-i-1}x^i ]

    形象地说,将(F(x))的系数倒过来
    回到正题

    [ecause f(x)=q(x)g(x)+r(x)\ herefore f(frac 1x)=q(frac 1x)g(frac 1x)+r(frac 1x)\ herefore x^{n-1}f(frac 1x)=x^{n-1}q(frac 1x)g(frac 1x)+x^{n-1}r(frac 1x)\ herefore x^{n-1}f(frac 1x)=x^{n-m}q(frac 1x)x^{m-1}g(frac 1x)+x^{n-m+1}cdot x^{m-2}r(frac 1x)\ herefore f_R(x)=q_R(x)g_R(x)+x^{n-m+1}r_R(x)\ herefore f_R(x)equiv q_R(x)g_R(x)pmod {x^{n-m+1}}\ herefore q_R(x)equiv frac{f_R(x)}{g_R(x)}pmod {x^{n-m+1}} ]

    于是多项式求逆+乘法就可以求出(q(x))

    那么

    [r(x)=f(x)-q(x)g(x) ]

    问题得以解决

    (q(x))代码如下

    inline poly div(int n,int m,poly f,poly g)
    {
    	rv(n,f),rv(m,g);
    	g=ginv(n-m+1,g);
    	poly q=mul(n-m+1,n-m+1,f,g);
    	rv(n-m+1,q);return q;
    }
    

    多项式快速幂(普通版)

    给定(n-1)次多项式(f(x))和正整数(k),求(g(x))满足

    [g(x)equiv f^k(x)pmod {x^n} ]

    系数均对(998244353)取模且保证(f(0)=1)

    首先对于多项式函数(f(x))满足(f(0)=1)证明如下命题:

    [f^p(x)equiv 1pmod {x^n} ]

    其中(p)是质数并且(n<p)

    原命题等价于

    [(1+sum_{i=1}^{n-1}a_ix^i)^pequiv 1pmod {x^n}\ Leftrightarrow 1^p+sum_{i=1}^{n-1}a_i^px^{ip}equiv 1pmod {x^n}\ Leftrightarrow 1equiv 1pmod {x^n}\ ]

    其中第一步到第二步的变换可以由多项式定理来证明

    于是

    [g(x)equiv f^k(x)equiv f^{k\%p}(x)pmod {x^n} ]

    不妨令(k'=k\%p)

    先两边同时求(ln)

    [ln g(x)equiv ln f^{k'}(x)equiv k'ln f(x)pmod {x^n} ]

    再两边同时求(exp)

    [g(x)equiv e^{k'ln f(x)}pmod {x^n} ]

    于是就做完了

    多项式快速幂(加强版)

    同上,唯一区别在于不保证(f(0)=1)

    这样会有什么影响呢?

    回头一看发现如果不保证这一点那么(ln,exp)都没有办法做

    于是考虑转化

    我们需要找到(f)中最低的系数非零的项,记为(a_tx^t),然后整个多项式除以(a_tx^t),于是我们成功的让(a_0)变成了(1),计算之后我们需要将答案右移(t^k)位,再乘(a_t^k)

    但是还是有很多坑

    • (t^k)如果大于等于(n)需要特判
    • (t^k,a_t^k)中的(k)都应该对(mod-1)取模而不是(mod)

    代码

    inline poly fsp(int n,int k1,int k2,poly f)
    {
    	int u,k;poly ans;
    	for(u=0;u<n&&!f[u];++u);k=qpow(f[u],k2);
    	if(1ll*u*k2>=n)return ans;
    	int iv=qpow(f[u],mod-2);
    	for(int i=u;i<n;++i)f[i]=mll(f[i],iv);
    	for(int i=0,j=u;j<n;++i,++j)ans[i]=f[j];
    	ans=ln(n,ans);
    	for(int i=0;i<n;++i)ans[i]=mll(ans[i],k1);
    	ans=exp(n,ans);
    	for(int i=0;i<u*k2;++i)f[i]=0;
    	for(int i=u*k2,j=0;i<n;++i,++j)f[i]=mll(ans[j],k);
    	return f;
    }
    int main()
    {
    	int n,k1=0,k2=0,k3=0;scanf("%d",&n);pre(n);
    	scanf("%s",ch+1);int len=strlen(ch+1);
    	poly f;read(n,f);
    	for(int i=1;i<=len;++i)
    	{
    		k1=add(mll(10,k1),ch[i]^48),
    		k2=(10ll*k2%(mod-1)+(ch[i]^48))%(mod-1);
    		if(k1>n&&!f[0])//这里要特判
    		{
    			for(int i=0;i<n;++i)printf("0 ");
    			return 0;
    		}
    	}
    	print(n,fsp(n,k1,k2,f));
    	return 0;
    }
    

    封装后所有的完整代码

    #include<bits/stdc++.h>
    using namespace std;
    const int mod=998244353,G=3,ivG=332748118,N=4e6+5;
    inline int in()
    {
        int x=0,f=1;char ch=getchar();
        while(!isdigit(ch)){if(ch=='-') f=-1;ch=getchar();}
        while(isdigit(ch)){x=(x<<3)+(x<<1)+(ch^48);ch=getchar();}
        return x*f;
    } 
    struct poly
    {
    	vector<int>v;
    	inline int&operator[](int x)
    	{
    		while(x>=v.size())v.push_back(0);
    		return v[x];
    	}
    	inline void pre(int x,int lim)
    	{
    		int t=min(lim,(int)v.size());
    		for(int i=x;i<t;++i)v[i]=0;
    	}
    };
    namespace Math
    {
    	int inv[N],pif;
    	inline int add(int x,int y){x+=y;return x>=mod?x-mod:x;}
    	inline int dec(int x,int y){x-=y;return x<0?x+mod:x;}
    	inline int mll(int x,int y){return 1ll*x*y%mod;}
    	struct pt
    	{
    		int x,y;
    		pt(int _x=0,int _y=0){x=_x,y=_y;}
    		inline pt operator*(const pt&rhs)
    		{
    			return pt(add(mll(x,rhs.x),mll(mll(y,rhs.y),pif)),add(mll(x,rhs.y),mll(y,rhs.x)));
    		}
    	};
    	inline int qpow(int x,int y)
    	{
    		int ans=1;
    		for(;y;y>>=1,x=mll(x,x))(y&1)&&(ans=mll(ans,x));
    		return ans;
    	}
    	inline pt qpow(pt x,int y)
    	{
    		pt ans=pt(1,0);
    		for(;y;y>>=1,x=x*x)if(y&1)ans=ans*x;
    		return ans;
    	}
    	inline void pre(int n)
    	{
    		inv[1]=1;
    		for(int i=2;i<=n;++i)inv[i]=mod-mll(mod/i,inv[mod%i]);
    	}
    	inline bool ck(int x){return qpow(x,(mod-1)>>1)==mod-1;}
    	inline int solve(int x)
    	{
    		int ans=0;
    		if(mod==2)return x;
    		if(!x)return 0;
    		else if(qpow(x,(mod-1)>>1)==mod-1)return -1;
    		int a=rand()%mod;
    		while(!a||!ck(dec(mll(a,a),x)))a=rand()%mod;
    		pif=dec(mll(a,a),x);
    		ans=qpow(pt(a,1),(mod+1)>>1).x;
    		return min(ans,mod-ans);
    	}//二次剩余求解
    }
    using namespace Math;
    namespace Bas
    {
    	int rev[N],wn[N];
    	inline void ntt(int lim,poly&f,bool tp)
    	{
    		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 sw=qpow(tp?ivG:G,(mod-1)/(mid<<1));
    			wn[0]=1;for(int i=1;i<mid;++i)wn[i]=mll(wn[i-1],sw);
    			for(int len=mid<<1,p=0;p+len-1<lim;p+=len)
    				for(int k=0;k<mid;++k)
    				{
    					int x=f[p+k],y=1ll*wn[k]*f[p+mid+k]%mod;
    					f[p+k]=add(x,y),f[p+mid+k]=dec(x,y);
    				}
    		}
    	}
    	inline poly mul(int n,int m,poly f,poly g)
    	{
    		int lim=1,l=0;while(lim<n+m)lim<<=1,++l;
    		for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
    		f.pre(n,lim),g.pre(m,lim);
    		ntt(lim,f,0),ntt(lim,g,0);
    		for(int i=0;i<lim;++i)f[i]=mll(f[i],g[i]);
    		ntt(lim,f,1);int iv=qpow(lim,mod-2);
    		for(int i=0;i<lim;++i)f[i]=mll(f[i],iv);
    		return f;
    	}//乘法
    	inline void igl(int n,poly&f)
    	{
    		for(int i=n-1;i;--i)f[i]=mll(f[i-1],inv[i]);f[0]=0;
    	}//积分
    	inline void diff(int n,poly&f)
    	{
    		for(int i=0,t=f[i+1];i<n-1;++i,t=f[i+1])f[i]=mll(t,i+1);f[n-1]=0;
    	}//微分(注意这里写得不当很容易RE)
    	inline void rv(int n,poly&f){reverse(f.v.begin(),f.v.begin()+n);}
    	inline void read(int n,poly&f){for(int i=0;i<n;++i)f[i]=in();}
    	inline void print(int n,poly f){for(int i=0;i<n;++i)printf("%d ",f[i]);puts("");}
    }
    using namespace Bas;
    namespace Poly
    {
    	poly ginv(int n,poly f)
    	{
    		if(n==1){poly g;g[0]=qpow(f[0],mod-2);return g;}
    		poly g=ginv(n+1>>1,f),h,p=g;
    		int lim=1,l=0;while(lim<n+n)lim<<=1,++l;
    		for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
    		f.pre(n,lim),p.pre(n,lim);
    		ntt(lim,f,0),ntt(lim,p,0);
    		for(int i=0;i<lim;++i)f[i]=mll(f[i],mll(p[i],p[i]));
    		ntt(lim,f,1);int iv=qpow(lim,mod-2);
    		for(int i=0;i<n;++i)h[i]=dec(add(g[i],g[i]),mll(f[i],iv));
    		return h;
    	}//多项式求逆
    	inline poly ln(int n,poly f)
    	{
    		poly g=ginv(n,f),h;diff(n,f);
    		h=mul(n,n,f,g);igl(n,h);
    		return h;
    	}//多项式求ln
    	poly exp(int n,poly f)
    	{
    		if(n==1){poly g;g[0]=1;return g;}
    		poly g=exp(n+1>>1,f);
    		poly h=ln(n,g);
    		for(int i=0;i<n;++i)h[i]=add(dec(0,h[i]),f[i]);
            h[0]=add(h[0],1);
    		return mul(n,n,g,h);
    	}//多项式求exp
    	inline poly div(int n,int m,poly f,poly g)
    	{
    		rv(n,f),rv(m,g);
    		g=ginv(n-m+1,g);
    		poly q=mul(n,n-m+1,f,g);
    		rv(n-m+1,q);return q;
    	}//多项式除法求商式
    	inline poly sqr(int n,poly f)
    	{
    		if(n==1){poly g;g[0]=solve(f[0]);return g;}
    		poly g=sqr(n+1>>1,f),h,ivg;
    		ivg=ginv(n,g),h=mul(n,n,g,g);
    		int iv=qpow(2,mod-2);
    		for(int i=0;i<n;++i)ivg[i]=mll(iv,ivg[i]);
    		for(int i=0;i<n;++i)h[i]=dec(h[i],f[i]);
    		h=mul(n,n,h,ivg);
    		for(int i=0;i<n;++i)h[i]=dec(g[i],h[i]);
    		return h;
    	}//多项式开根
    	inline poly fsp(int n,int k1,int k2,poly f)
    	{
    		int u,k;poly ans;
    		for(u=0;u<n&&!f[u];++u);k=qpow(f[u],k2);
    		if(1ll*u*k2>=n)return ans;
    		int iv=qpow(f[u],mod-2);
    		for(int i=u;i<n;++i)f[i]=mll(f[i],iv);
    		for(int i=0,j=u;j<n;++i,++j)ans[i]=f[j];
    		ans=ln(n,ans);
    		for(int i=0;i<n;++i)ans[i]=mll(ans[i],k1);
    		ans=exp(n,ans);
    		for(int i=0;i<u*k2;++i)f[i]=0;
    		for(int i=u*k2,j=0;i<n;++i,++j)f[i]=mll(ans[j],k);
    		return f;
    	}//多项式快速幂
    }
    using namespace Poly;
    char ch[N];
    int main(){return 0;}
    //上面所有的n都表示n-1次多项式 
    
  • 相关阅读:
    bzoj2049 [Sdoi2008]Cave 洞穴勘测——LCT
    洛谷P2679 子串——DP
    bzoj3669 [Noi2014]魔法森林——LCT
    洛谷P3778 [APIO2017]商旅——01分数规划
    bzoj4196 [Noi2015]软件包管理器——树链剖分
    bzoj4881 线段游戏——上升序列方案数
    bzoj1426 (洛谷P4550) 收集邮票——期望
    bzoj1858 [Scoi2010]序列操作——线段树
    bzoj3626 [LNOI2014]LCA——树链剖分
    L The Digits String(没有写完,有空补)
  • 原文地址:https://www.cnblogs.com/zmyzmy/p/14159779.html
Copyright © 2011-2022 走看看