多项式Ⅱ

前言

因为之前的 Ⅰ 写的太屑了,所以又写了篇更屑的Ⅱ,其实是不同时间的个人垃圾理解了..

牛顿迭代

  • 要求解(f(x)=0)

    [frac{f(x_0)}{x_0-x}=f'(x_0)\ x=x_0-frac{f(x_0)}{f'(x_0)} ]

    不断往下迭代可以找到近似解。

多项式牛顿迭代

和普通的牛顿迭代其实一样

[F=F_0-frac{G(F_0)}{G'(F_0)} ]

如何理解近似解?

(mod x^i)开始向上(mod x^{i+1}...)

  • 多项式求逆

    [A*Bequiv 1pmod {x^n} ]

    构造

    [G(B)=A*B-1 ]

    牛顿迭代找根

[ egin{aligned} B&=B_0-frac{AB_0-1}{A}\ &=B_0-B_0(AB_0-1)\ &=2B_0-AB_0^2 end{aligned} ]

向上倍增即可。

  • 多项式开根

    [B^2equiv Apmod {x^n} ]

    构造

    [G(B)=A-B^2 ]

    牛顿迭代找根

    [egin{aligned} B&=B_0-frac{A-B_0^2}{-2B_0}\ &=frac{A}{2B_0}+frac{B_0}{2} end{aligned} ]

    向上倍增即可

  • 多项式exp

    [Bequiv e^A pmod {x^n} ]

    并没有弄清楚这个咋定义的,因此还不会直接用。

    [ln B=A ]

    构造

    [G(B)=ln B-A ]

    [egin{aligned} B&=B_0-frac{ln B_0-A}{frac{1}{B_0}}\ &=B_0(1-ln B_0+A) end{aligned} ]

多项式取ln

  • 多项式积分和求导

    [egin{aligned} F(x)&=sum_{i=0}^na_ix^i\ frac{mathrm dF}{mathrm dx}&=sum_{i=0}^{n-1}a_{i+1}(i+1)x^i\ int F(x)&=sum_{i=1}^nfrac{a_{i-1}}{i}x^i end{aligned} ]

  • 多项式取ln

    [egin{aligned} B&=ln Apmod {x^n}\ B'&=ln'(A)A'pmod {x^n}\ &=frac{A'}{A}pmod {x^n} end{aligned} ]

  • 多项式快速幂

    [F^k=e^{kln n}pmod {x^n} ]

    常数无敌的(nlog n)

拉格朗日反演和快速插值什么的,我还没学...


放一个exp的代码

提几个细节

  1. 不知道为什么递归版效率优于非递归版,而递归版又好写一点,因此以后就写递归版了
  2. 每个需要倍增的,每次要外面开两个自己的数组,不许改传进来的值,注意这两个数组一定需要清空
  3. 要注意一下每次NTT的长度什么的
  4. 不要封装卷积乘法,因为慢,可以把多个多项式乘起来全部先转换成点表达,然后一起乘,再转回系数,这样需要注意长度。
#include <cstdio>
#include <algorithm>
const int N=(1<<18)+10;
const int mod=998244353,Gi=332748118;
#define mul(a,b) (1ll*(a)*(b)%mod)
#define add(a,b) ((a+b)%mod)
int invA[N],invB[N],lnA[N],lnB[N],exA[N],exB[N],ans[N],f[N],turn[N],n;
int qp(int d,int k){int f=1;while(k){if(k&1)f=mul(f,d);d=mul(d,d),k>>=1;}return f;}
void NTT(int *a,int typ,int len)
{
    int L=-1;for(int i=1;i<len;i<<=1) ++L;
    for(int i=1;i<len;i++)
    {
        turn[i]=turn[i>>1]>>1|(i&1)<<L;
        if(i<turn[i]) std::swap(a[i],a[turn[i]]);
    }
    for(int le=1;le<len;le<<=1)
    {
        int wn=qp(typ?3:Gi,(mod-1)/(le<<1));
        for(int p=0;p<len;p+=le<<1)
        {
            int w=1;
            for(int i=p;i<p+le;i++,w=mul(w,wn))
            {
                int tx=a[i],ty=mul(w,a[i+le]);
                a[i]=add(tx,ty);
                a[i+le]=add(tx,mod-ty);
            }
        }
    }
    if(!typ)
    {
        int inv=qp(len,mod-2);
        for(int i=0;i<len;i++) a[i]=mul(a[i],inv);
    }
}
void polyinv(int *a,int *b,int len)
{
    if(len==1) {b[0]=qp(a[0],mod-2);return;}
    polyinv(a,b,len>>1);
    for(int i=0;i<len<<1;i++) invA[i]=invB[i]=0;//
    for(int i=0;i<len;i++) invA[i]=a[i],invB[i]=b[i];
    NTT(invA,1,len<<1),NTT(invB,1,len<<1);
    for(int i=0;i<len<<1;i++) invA[i]=mul(invB[i],add(2,mod-mul(invA[i],invB[i])));
    NTT(invA,0,len<<1);
    for(int i=0;i<len;i++) b[i]=invA[i];
}
void direv(int *a,int n){for(int i=0;i<n;i++)a[i]=mul(a[i+1],i+1);a[n]=0;}
void inter(int *a,int n){for(int i=n;i;i--)a[i]=mul(a[i-1],qp(i,mod-2));a[0]=0;}
void polyln(int *a,int *b,int len)
{
    for(int i=0;i<len;i++) lnA[i]=a[i],lnB[i]=b[i];
    polyinv(lnA,lnB,len);
    direv(lnA,len-1);
    NTT(lnA,1,len<<1),NTT(lnB,1,len<<1);
    for(int i=0;i<len<<1;i++) lnA[i]=mul(lnA[i],lnB[i]);
    NTT(lnA,0,len<<1);
    inter(lnA,len-1);
    for(int i=0;i<len;i++) b[i]=lnA[i];
}
void polyexp(int *a,int *b,int len)
{
    if(len==1){b[0]=1;return;}//默认a[0]=0
    polyexp(a,b,len>>1);
    for(int i=0;i<len<<1;i++) exA[i]=exB[i]=0;//
    polyln(b,exA,len);
    for(int i=0;i<len;i++) exA[i]=add(a[i]+(i==0),mod-exA[i]);
    for(int i=0;i<len;i++) exB[i]=b[i];
    NTT(exA,1,len<<1),NTT(exB,1,len<<1);
    for(int i=0;i<len<<1;i++) exA[i]=mul(exA[i],exB[i]);
    NTT(exA,0,len<<1);
    for(int i=0;i<len;i++) b[i]=exA[i];
}
int main()
{
    scanf("%d",&n);
    for(int i=0;i<n;i++) scanf("%d",f+i);
    int len=1;while(len<n) len<<=1;
    polyexp(f,ans,len);
    for(int i=0;i<n;i++) printf("%d ",ans[i]);
    return 0;
}

2018.12.29

原文地址:https://www.cnblogs.com/butterflydew/p/10196095.html