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

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

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

目录

• 再探快速傅里叶变换(其四) 多项式操作
  • 多项式乘法
  • 多项式求逆
  • 多项式取余
  • 多项式ln
  • 多项式exp
    • 泰勒展开
    • 牛顿迭代
    • exp的求解
  • 多项式开根
    • 开根求解
  • 多项式快速幂
    • 保证常数项为1
    • 不保证常数项为1
  • 多项式多点求值
  • 多项式快速插值

多项式乘法

不解释。注意如果式子太复杂，直接全部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))

牛顿迭代

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

牛顿迭代的几何解释如图,先选取一个初始值(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
}

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

多项式多点求值

多项式快速插值