zoukankan      html  css  js  c++  java
  • 再探快速傅里叶变换(FFT)(其四) 多项式操作

    再探快速傅里叶变换(其四) 多项式操作

    省选前夕爆补一波多项式全家桶,存一些板子

    多项式乘法

    不解释。注意如果式子太复杂,直接全部DFT再按原式子算,最后IDFT。这里封装好是为了简化模板代码,而且可以方便的换成任意模数NTT

    void poly_mul(ll *a,ll *b,ll *c,int n,int m){
        static ll ta[maxn+5],tb[maxn+5];
        int N=1,L=0;
        while(N<n+m-1){
            N*=2;
            L++;
        }
        for(int i=0;i<n;i++) ta[i]=a[i];
        for(int i=n;i<N;i++) ta[i]=0;
        for(int i=0;i<m;i++) tb[i]=b[i];
        for(int i=n;i<N;i++) tb[i]=0;
        for(int i=0;i<N;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
        NTT(ta,N,1);
        NTT(tb,N,1);
        for(int i=0;i<N;i++) c[i]=ta[i]*tb[i]%mod;
        NTT(c,N,-1);
        for(int i=n+m-1;i<N;i++) c[i]=0;
    } 
    

    多项式求逆

    定义:LuoguP4238给出(n-1)次多项式(f(x)),我们要求一个多项式(g(x)),满足(f(x)g(x) equiv 1(mod x^n)).
    系数对NTT模数取模.

    求逆的算法基于倍增,是迭代进行的。假设我们已经求出了一个(h(x)),满足(f(x)h(x) equiv 1 (mod x^{frac{n}{2}}))

    由定义有(f(x)g(x) equiv 1(mod x^{frac{n}{2}}))

    两式做差,的(h(x)-g(x) equiv 1 (mod x^{frac{n}{2}}))

    平方,得(h^2(x)-2h(x)g(x)+g^2(x) equiv 1(mod x^n))

    两边同时乘上(f(x)),将(f(x)g(x)=1)代入,移项可得:

    [g(x) equiv h(x)(2-f(x)h(x)) (mod x^n) ]

    因此递归求出(h(x))后,就可以推出(g(x)).边界条件(n=1)(g_0=f_0^{-1}). 时间复杂度满足递推式(T(n)=T(frac{n}{2})+n log n),根据主定理为(Theta(nlog n))

    void poly_inv(ll *f,ll *g,int n){
        static ll tmp[maxn+5];
        if(n==1){
            g[0]=inv(f[0]);
            return;
        }
        poly_inv(f,g,(n+1)/2);
        poly_mul(f,g,tmp,n,n);
        poly_mul(tmp,g,tmp,n,n);
        for(int i=0;i<n;i++) g[i]=(2*g[i]-tmp[i]+mod)%mod; 
    }
    

    多项式取余

    定义:[LuoguP4512]给定一个(n)次多项式(f(x))和一个(m)次多项式(g(x)),求两个多项式(q(x),r(x)).满足(f(x)=q(x)g(x)+r(x)),且(q(x))次数为(n-m),(r(x))次数小于(m).模数为NTT模数。

    我们先定义一种操作(r),若(f(x)=sum_{i=0}^{n}a_ix^i),则(f_r(x)=sum_{i=0}^{n-1}a_ix^{n-i}).实际上就是把系数顺序倒过来。容易发现(f_r(x)=x^nf(frac{1}{x})).接下来开始推式子

    [egin{aligned} f(x)=q(x) cdot g(x)+r(x) \ Leftrightarrow fleft(frac{1}{x} ight)=qleft(frac{1}{x} ight) cdot gleft(frac{1}{x} ight)+rleft(frac{1}{x} ight) \ Leftrightarrow x^{n} fleft(frac{1}{x} ight)=x^{n-m} qleft(frac{1}{x} ight) cdot x^{m} gleft(frac{1}{x} ight)+x^{n-m+1} cdot x^{m-1} rleft(frac{1}{x} ight) \ Leftrightarrow f_{r}(x)=q_{r}(x) cdot g_{r}(x)+x^{n-m+1} cdot r_{r}(x) \ Leftrightarrow f_{r}(x) equiv q_{r}(x) cdot g_{r}(x) quadleft(mod x^{n-m+1} ight) \ Leftrightarrow q_{r}(x) equiv f_{r}(x) cdot g_{r}^{-1}(x) quadleft(mod x^{n-m+1} ight) end{aligned}]

    我们只需要在(mod x^{n-m+1})下对(g_r)求逆,就能求出(q_r)(q),然后用(r(x)=f(x)-g(x)q(x))求出(r).

    //a就是f,b就是g
    void poly_mod(ll *a,ll *b,ll *q,ll *r,int n,int m){
        static ll ra[maxn+5],rb[maxn+5],invrb[maxn+5];
        for(int i=0;i<n;i++) ra[i]=a[i];
        for(int i=0;i<m;i++) rb[i]=b[i];
        reverse(ra,ra+n);
        reverse(rb,rb+m);
        for(int i=n-m+1;i<n*2;i++) rb[i]=0;
        poly_inv(rb,invrb,n-m+1); //在mod x^(n-m+1)下做 
        poly_mul(ra,invrb,q,n,n-m+1);
        for(int i=n-m+1;i<n*2;i++) q[i]=0;
        reverse(q,q+n-m+1);
        poly_mul(q,b,r,n-m+1,m);
        for(int i=0;i<n;i++) r[i]=(a[i]-r[i]+mod)%mod;
    }
    

    多项式ln

    定义:LuoguP4725给出(n-1)次多项式(f(x)),我们要求一个多项式(g(x)),满足(g(x) equiv ln f(x)(mod x^n)).
    系数对NTT模数取模.保证常数项为1.

    对上式求导,(g'(x)=[ln f(x)]'=ln'(f(x))cdot f'(x)=frac{f'(x)}{f(x)})

    这样就可以用多项式求逆得到(g'(x)),再积分回去就好了。
    保证常数项为1,是保证原多项式为整系数的情况下,(ln)多项式的常数用整数算得出来.并且我们发现求导的过程中实际上把常数项忽略掉了。

    void poly_deriv(ll *f,ll *g,int n){
        for(int i=1;i<n;i++) g[i-1]=f[i]*i%mod;
        g[n-1]=0;
    } 
    void poly_inter(ll *f,ll *g,int n){
        for(int i=n-1;i>=1;i--) g[i]=f[i-1]*inv(i)%mod;
        g[0]=0;
    }
    void poly_ln(ll *f,ll *g,int n){
        static ll invf[maxn+5];
        poly_deriv(f,g,n);
        poly_inv(f,invf,n);
        poly_mul(g,invf,g,n,n);
        poly_inter(g,g,n*2);
        for(int i=n;i<n*2;i++) g[i]=0;
    }
    

    多项式exp

    泰勒展开

    泰勒展开是用多项式来逼近指数函数,三角函数等函数。因为转化成了多项式函数,所以很多时候能简化分析(比如在生成函数问题中)。

    为了逼近(感性理解就是让图像形状更接近)目标函数(f(x)),我们先选点((x_0,f(x_0)))来切入,让多项式函数(g(x))经过这一点。为了让图像更接近,我们让此处的增减性相同(一阶导数相同),凹凸性相同(二阶导数相同),以及更高阶的导数相同...

    实际操作中我们不用无限阶求导,只需算到(n)阶.

    (g(x)=a_0+a_1(x-x_0)+a_2(x-x_0)^2+dots a_n(x-x_0)^n).(每一项里面的(x-x_0)是为了保证过定点)

    由初始点相等,有(a_0=f(0))
    (x_0)(n)阶导数相等,有(f^{(n)}(x_0)=g^{(n)}(x_0)=n!a_n),即(a_n=frac{f^{(n)}}{n!}).(求(n)阶导之后只有第(n)项留下来,(a_nx^n Arr na_nx^{n-1} Arr n(n-1)a_nx^{n-2} Arr dots n!a_n))

    代入,得

    [g(x)=f(0)+frac{f'(x_0)}{1!}(x-x_0)+frac{f''(x_0)}{2!}(x-x_0)^2+dots frac{f^{(n)}(x_0)}{n!}(x-x_0)^n+R_n(x) ]

    这就是著名的泰勒展开式.其中(R_n(x))是佩亚诺(Peano)余项,感性理解就是和真实值的误差。特别地,令(x_0=0)就可以得到麦克劳林展开:

    [g(x)=f(0)+frac{f'(0)}{1!}x+frac{f''(0)}{2!}x^2+dots frac{f^{(n)}(0)}{n!}x^n+R_n(x) ]

    现在给出一些OI中常用函数的泰勒展开。在生成函数问题中经常用到.

    [e^x=1+frac{x}{1!}+frac{x^2}{2!}+dots frac{x^n}{n!}+dots=sum_{n=0}^{infin}frac{x^n}{n!} ag{1} ]

    证明:取(x_0=0),(e^x)求导等于本身,所以每一项都是(e^0=1)

    [frac{1}{1-x}=1+x+x^2+x^3+dots=sum_{n=0}^{infin} x^n ag{2} ]

    证明:取(x_0=0),((frac{1}{1-x})'=frac{1}{(1-x)^2},frac{1}{(1-x)^2}'=frac{1}{(1-x)^3}dots)

    [-ln (1-x) =frac{x}{1}+frac{x^2}{2}+dots frac{x^n}{n}+dots=sum_{n=1}^{infin} frac{x^n}{n} ag{3} ]

    证明: 对上式两边求导,即可得到((2))

    牛顿迭代

    牛顿迭代是一种求方程近似解的方法。

    image.png

    牛顿迭代的几何解释如图,先选取一个初始值(x_0),向上做垂线与(f(x))的图像相交。作(f(x))(x_0)处的切线交(x)轴于(x_1),从(x_1)再向上做垂线,如此迭代,(x_i)会逐渐趋近于方程的解。由导数的几何意义,(x_i)处切线方程为(y=f'(x_i)x+f(x_i)-f'(x_i)x_i).令(y=0),得交点坐标(x_{i+1}=x_{i}-frac{f(x_i)}{f'(x_i)}). 按照这个递推式不断迭代就可以得到近似解。

    其中迭代的实数(x_i)也可以换成一个多项式,并且在模意义下进行,此时就可以求出精确值。

    牛顿迭代也可以用泰勒展开解释,假设我们已经求出(h(x))满足(f(h(x)) equiv 0(mod x^{frac{n}{2}})),要求(f(g(x)) equiv 0(mod x^n))

    (f(g(x)))(h(x))处泰勒展开,有:

    [f(g(x))=f(h(x))+f'(h(x))(g(x)-h(x)+frac{f''(h(x))}{2!}(g(x)-h(x))^2+dots ]

    因为我们是迭代进行的,显然有(g(x)equiv h(x) (mod x^{frac{n}{2}})).那么对于(forall i geq 2)((g(x)-h(x))^i equiv 0(mod x^n)),也就是说
    (f(g(x))=f(h(x))+f'(h(x))(g(x)-h(x))),移项得

    [f(g(x))=h(x)-frac{f(h(x))}{f'(h(x))}(mod x^n) ]

    exp的求解

    定义:LuoguP4726给出(n-1)次多项式(f(x)),我们要求一个多项式(g(x)),满足(g(x) equiv mathrm{e}^{f(x)} (mod x^n)).
    系数对NTT模数取模

    对上式两边求(ln),得(ln(g(x))-f(x) equiv 0(mod x^n))

    (B(t(x))=ln(t(x))-f(x)),把(t(x))看成自变量,用牛顿迭代求使得(B(t(x))=0)(t(x))

    我们设第(i)次迭代后的根为(g_i(x)),根据牛顿迭代的公式有:

    [g_{i+1}(x)=g_{i}(x)-frac{H(g_{i}(x))}{H'(g_{i}(x))}=g_{i}(x)-frac{ln(g_{i}(x))-f(x)}{frac{1}{g_i(x)}}=g_{i}(x)(1+f(x)+ln g_{i}(x)) ]

    注意(H'(g_{i}(x))=frac{1}{g_{i}(x)}),因为我们把(g_{i}(x))看作自变量,(f(x))看作常数。

    每次迭代都使得精度翻倍(指的是最高次翻倍).复杂度(T(n)=T(frac{n}{2})+nlog_2n),即(Theta(nlog n))

    void poly_exp(ll *f,ll *g,int n){
        static ll lng[maxn+5];
        if(n==1){
            g[0]=1;
            return;
        }
        poly_exp(f,g,(n+1)/2);
        poly_ln(g,lng,n);
        for(int i=0;i<n;i++) lng[i]=(f[i]-lng[i]+mod)%mod;
        lng[0]++;
        poly_mul(g,lng,g,n,n);
        for(int i=n;i<n*2;i++) g[i]=0;
    }
    

    多项式开根

    开根求解

    定义: LuoguP5277给定一个 (n-1) 次多项式 (f(x),) 求一个在 (mod x^{n}) 意义下的多项式 (g(x),) 使得 (g^{2}(x) equiv f(x)) (left(mod x^{n} ight),) 答案取正的 (sqrt{f(x)})(正的指模意义下常数项较小的).NTT模数

    多项式开根一般用倍增算法。虽然可以类似求逆直接用定义推式子,但是牛顿迭代更快。另外类似多项式快速幂(ln+exp)也可以,但常数较大.

    定义(B(h(x))=h^2(x)-f(x)),我们要求使得(B(g(x))=0)(g(x)).根据牛顿迭代的式子,假设我们已经在模(x^{frac{n}{2}})的意义下求出了解(h(x)),那么有:

    [g(x)=h(x)-frac{B(h(x))}{B'(h(x))}=h(x)-frac{h^2(x)-f(x)}{2h(x)}=frac{h^2(x)+f(x)}{2h(x)} ]

    套求逆的板子即可。初始值(g(0))需要求(f(0))的二次剩余解,即(g(0)^2 equiv f(0) (mod p)),记得取系数小的.因为二次剩余的算法并不是学习其他多项式算法的前置知识,这里不再赘述,详见此博客

    void poly_sqrt(ll *f,ll *g,int n) {
        static ll invg[maxn+5];
        if(n==1) {
            ll root=(Cipolla::cipolla(f[0])%mod+mod)%mod;
            g[0]=min(root,mod-root);
            return;
        }
        poly_sqrt(f,g,(n+1)/2);
        poly_inv(g,invg,n);
        poly_mul(f,invg,invg,n,n);
        ll inv2=inv(2);
        for(int i=0; i<n; i++) g[i]=(g[i]+invg[i])%mod*inv2%mod;
    }
    

    多项式快速幂

    保证常数项为1

    定义:LuoguP5245给出(n-1)次多项式(f(x))和正整数(k),我们要求一个多项式(g(x)),满足(g(x) equiv f^k(x) (mod x^n)).保证常数项为1
    系数对NTT模数取模,(n leq 10^5,k leq 10^{10^5})

    显然有$f^k(x)=exp(k ln f(x)) ((常数项为1是为了求)ln()。套前面的板子,复杂度)O(nlog n)( 注意)k(很大需要取模,因为我们一直在对系数进行操作所以取模的是)p(而不是)p-1$.(这里和费马小定理没有关系!)

    void fast_pow(ll *a,ll *b,ll k,int n){
        static ll lna[maxn+5];
        poly_ln(a,lna,n);
        for(int i=0;i<n;i++) lna[i]=lna[i]*k%mod; 
        poly_exp(lna,b,n);//不能直接带入a,否则会破坏初始值
    }
    

    不保证常数项为1

    定义:LuoguP5273给出(n-1)次多项式(f(x))和正整数(k),我们要求一个多项式(g(x)),满足(g(x) equiv f^k(x) (mod x^n)).不保证常数项为1
    系数对NTT模数取模,(n leq 10^5,k leq 10^{10^5})

    既然不保证常数项为1,我们可以把常数项除掉。但是常数项也可能为0.那么我们就要找到第一个非0的项,设(f(x))的最低次项为(ax^d).那么有:

    [f^{k}(x)=a^{k} x^{k d}left(frac{f(x)}{a x^{d}} ight)^{k} ]

    对应到系数表达式上,就是把系数先左移(d)位,除以(a),按普通版的方法求(k)次幂,然后乘上(a^k),最后再右移(kd)位(如果(kd>n)就直接输出0).注意做多项式(k)次幂和右移(kd)位的时候(k)取模的是(p),但是求(a^k)的时候(k)取模的是(p-1).

    void fast_pow(ll *a,ll *b,ll k1,ll k2,int n){
        static ll lna[maxn+5];
        ll pos=n;
        for(int i=0;i<n;i++){
            if(a[i]!=0){//找到第一个系数不为0的位置.然后整体左移pos,相当于除法
                pos=i;
                break;
            }
        }  
        if(pos*k1>=n){
            for(int i=0;i<n;i++) printf("0 ");
            return;
        }
    
        ll d=a[pos];//模板的多项式ln要求常数项a[0]=1,所以要整体除掉常数
        ll invd=inv(d),powd=fast_pow(d,k2);
        for(int i=pos;i<n;i++) a[i-pos]=a[i]*invd%mod; 
        poly_ln(a,lna,n-pos); 
        for(int i=0;i<n-pos;i++) lna[i]=lna[i]*k1%mod; 
        poly_exp(lna,b,n-pos);//不能直接带入a,否则会破坏初始值
        for(ll i=0;i<min(pos*k1,n*1ll);i++) printf("0 ");
        for(ll i=pos*k1;i<n;i++) printf("%lld ",b[i-pos*k1]*powd%mod); //k2在指数上,模mod-1
    }
    

    此题细节较多,放上完整代码

    多项式多点求值

    多项式快速插值

  • 相关阅读:
    一行代码搞定Dubbo接口调用
    测试周期内测试进度报告规范
    jq 一个强悍的json格式化查看工具
    浅析Docker容器的应用场景
    HDU 4432 Sum of divisors (水题,进制转换)
    HDU 4431 Mahjong (DFS,暴力枚举,剪枝)
    CodeForces 589B Layer Cake (暴力)
    CodeForces 589J Cleaner Robot (DFS,或BFS)
    CodeForces 589I Lottery (暴力,水题)
    CodeForces 589D Boulevard (数学,相遇)
  • 原文地址:https://www.cnblogs.com/birchtree/p/13263252.html
Copyright © 2011-2022 走看看