「总结」多项式生成函数相关

Posted 2021-03-23 lrefrain

实在是太毒瘤了。

大纲。

多项式生成函数相关

默认前置：微积分，各种数和各种反演，FFT，NTT，各种卷积，基本和式变换。

主要内容：
泰勒展开，级数求和，牛顿迭代，主定理。　　　　//例题：在美妙的数学王国中畅游,礼物
多项式全家桶：乘法，求逆，求导，积分，分治，ln，exp，fwt,MTT。 //城市规划，图的价值，染色，遗失的答案，按位或，随机游走。
生成函数：普通型生成函数，指数型生成函数计数原理。 //猎人杀，遗忘的集合，生成树计数
例题。

一、泰勒展开和级数求和

1.泰勒展开

即对于任何函数(f(x)),如果在(x_0)处(n)阶可导，那么满足如下公式：

[f(x)=sumlimits_{i=0}^{n}frac{f^{(i)}(x_0)(x-x_0)^i}{i!}]

这里当(x_0)为0的时候被称作麦克劳林公式。

先推导麦克劳林公式

即：

[f(x)=sumlimits_{i=0}^{n}frac{f^{(i)}(0)x^i}{i!}]

这里只证明多项式函数的正确性（其实是因为任意函数太难证了吧）。

设多项式函数:

[f(x)=sumlimits_{i=0}^{n}a_ix^i]

那么：

[f^{(i)}(x)=sumlimits_{j=0}^{n}a_{j+i}x^{j}prodlimits_{k=j+1}^{j+i+1}k]

而我们的(x_0=0)

也就是说每阶导数都只有常数项。

即：

[f^{(i)}(x)=i!a_i]

所以

[f(x)=sumlimits_{i=0}^{n}frac{f^{(i)}x^i}{i!}]

那如果(x_0!=0)呢？

这时候就是真正的泰勒展开了。

我们照旧设多项式函数：

[f(x)=sumlimits_{i=0}^{n}a_ix^i]

将之更换为关于(w=x-x_0)的多项式。

[g(w)=f(x)=sumlimits_{i=0}^{n}b_i(x-x_0)^i]

用麦克劳林公式展开：

[g(w)=sumlimits_{i=0}^{n}frac{g^{(i)}(0)w^i}{i!}]

也就是说:[b_i=frac{g^{(i)}}{i!}]

那么:

[g(w)=f(w+x_0),g^{(i)}(w)=f^{(i)}(w+x_0)]

[g(0)=f(x_0),g^{(i)}(0)=f^{(i)}(x_0)]

所以:

[b_i=frac{f^{(i)}(x_0)}{i!}]

这样的话就是说：

[f(x)=sumlimits_{i=0}^{n}b_i(x-x_0)^i=sumlimits_{i=0}^{n}frac{f^{(i)}(x_0)(x-x_0)^i}{i!}]

证毕。

2.一些级数求和的公式

其实个人认为泰勒展开是级数求和的逆运算。

[egin{array}{rcl}f(x)&=&e^x&=&limlimits_{n ightarrow +infty}sumlimits_{i=0}^{n}frac{f^{(i)}(0)x^i}{i!}&=&limlimits_{n ightarrow +infty}sumlimits_{i=0}^{n}frac{e^0x^i}{i!}&=&limlimits_{n ightarrow +infty}sumlimits_{i=0}^{n}frac{x^i}{i!}end{array}]

[egin{array}{rcl}f(x)&=&frac{1}{1-x}&=&limlimits_{n ightarrow +infty}sumlimits_{i=0}^{n}frac{f^{(i)}(0)x^i}{i!}&=&limlimits_{n ightarrow +infty}sumlimits_{i=0}^{n}frac{frac{i!}{1-0}x^i}{i!}&=&limlimits_{n ightarrow +infty}sumlimits_{i=0}^{n}x^iend{array}]

[xin(-1,1)]

[egin{array}{rcl}f(x)&=&frac{1}{1+x}&=&limlimits_{n ightarrow +infty}sumlimits_{i=0}^{n}frac{f^{(i)}(0)x^i}{i!}&=&limlimits_{n ightarrow +infty}sumlimits_{i=0}^{n}frac{frac{(-1)^i}{1-0}x^i}{i!}&=&limlimits_{n ightarrow +infty}sumlimits_{i=0}^{n}(-1)^ix^iend{array}]

[xin(-1,1)]

[egin{array}{rcl}f(x)&=&sin(x)&=&frac{sin(0)x^0}{0!}+frac{cos(0)x^1}{1!}+frac{-sin(0)x^2}{2!}+frac{-cos(0)x^3}{3!}…&=&limlimits_{n ightarrow +infty}sumlimits_{i=1}^{n}(-1)^{n-1}frac{x^{2i-1}}{(2i-1)!}end{array}]

[egin{array}{rcl}f(x)&=&cos(x)&=&frac{cos(0)x^0}{0!}+frac{-sin(0)x^1}{1!}+…&=&limlimits_{n ightarrow +infty}sumlimits_{i=0}^{n}(-1)^ifrac{x^{2i}}{(2i)!}end{array}]

[egin{array}{rcl}f(x)&=&ln(1+x)&=&frac{ln(1+0)x^0}{0!}+frac{frac{1}{1+0}x^1}{1!}+frac{-frac{1!}{(1+0)^2}x^2}{2!}+…&=&limlimits_{n ightarrow +infty}sumlimits_{i=1}^{n}(-1)^{i-1}frac{x^i}{i}end{array}]

二、牛顿迭代

1.牛顿迭代

本来这个东西是用来求连续函数的。

先讲一下连续函数是怎么搞得。

对于一个函数(f(x))

我们期望求解(f(x)=0)的解(x_z)

首先需要一个接近(x_z)的值(x_0)，保证在区间([x_z,x_0])范围内，(f^{(2)}(x))的符号不变。

而我们期望得到一个更加接近(x_z)的解(x_1)。

怎么得到呢？

将(f(x_1))在(x_0)处泰勒展开。

得到：

[f(x_1)=sumlimits_{i=0}^{+infty}frac{f^{(i)}(x_0)(x_1-x_0)^i}{i!}]

我们只保留其线性部分。

个人猜测牛顿这样做的原因是因为这样迭代一次之后(x_1)仅有一个解，如果保留二次项部分，虽然会让迭代的次数减少，不过此时有两个解，剩余部分的迭代增加量变成了指数级别的。

那么也就是说我们需要求：

[f(x_1)=f(x_0)+f'(x_0)(x_1-x_0)=0]

此时(x_1=x_0-frac{f(x_0)}{f'(x_0)})

不断迭代下去就可以不断逼近解(x_z)

考虑为什么。

由于二阶导符号不变，那么也就是说函数的凹凸性不变，这样的话一阶导的就是单调的了。

其实就是函数切线的斜率单调，我们求出的(x_1)就是切线的解，这样就可以不断的逼近(x_z)。

2.多项式牛顿迭代

本来这个东西使用来求连续函数的。

现在我们把它搞到离散数学里面来，利用相似的思想。

求解如下方程。

[F(G(x))equiv 0(mod x^n)]

按照套路来。

我们假设已经求出了(G'(x))

满足：

[F(G'(x))equiv 0(mod x^{frac{n}{2}})]

对函数(F(G(x)))在(G'(x))处泰勒展开。

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

可以发现：

[G(x)-G'(x)equiv 0(mod x^{frac{n}{2}})]

所以：

[(G(x)-G'(x))^iequiv 0(mod x^{frac{ni}{2}}),i>=2]

这样也就是说，不准确的部分只有一阶导部分。

也就是：

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

每次迭代精度翻倍了。

可以递归求解。

复杂度是：(T(n)=nlogn+T(frac{n}{2})=nlogn)。

可能要疑问这个是干什么用的。

下面的全家桶就知道了。

讲个题：https://www.cnblogs.com/Lrefrain/p/11921361.html.

三、主定理

主定理用来证明一些分治递归算法的复杂度。

即。

对于如下的复杂度：

[T(n)=aT(frac{n}{b})+f(n)]

设(f(n)=n^d)

有：

[T(n)=egin{cases}n^d& a<n^d ^dlog_bn& a=n^d ^{log_ba}& a>n^d end{cases}]

偷一张算导的图：

是递归示意图了。

那么可以发现总的复杂度可以这样来表示：

[egin{array}{rcl}T(n)&=&sumlimits_{i=0}^{log_bn}left(frac{n}{b^i} ight)^da^{i}&=&n^dsumlimits_{i=0}^{log_bn}left(frac{a}{b^d} ight)^iend{array}]

讨论一下。

1.(a<n^d)后面那个东西是小数，所以说是收敛的，是一个常数，所以(T(n)=n^d)

2.(a=n^d)后面全是1，所以(T(n)=n^dlog_bn)

3.(a>n^d)后面东西，对于最后一项来说，前面的全是常数，所以舍去前面的项。

推一下：

[egin{array}{rcl}T(n)&=&n^dleft(frac{a}{b^d} ight)^{log_bn}&=&n^dfrac{a^{log_bn}}{n^{dlog_bn}}&=&n^dfrac{a^{log_bn}}{n^d}&=&n^{log_ba}end{array}]

四、多项式全家桶

1.导数(deriv)
就和普通导数一样，对于第(i)项，(A[i]=A[i+1]*(i+1))常数项舍去。

2.不定积分(integ)
就和普通不定积分一样，对于第(i)项，(A[i]=frac{A[i-1]}{i})最高次项舍去。

3.求逆(inv)
这个稍微推一下。
对于一个多项式(A(x))，我们需要的(B(x))是满足：
[A(x)B(x)equiv 1(mod x^n)]
假设我们已经求出了一个(B'(x))使得：
[A(x)B_0(x)equiv 1(mod x^{frac{n}{2}})]
那么：
[A(x)(B(x)-B_0(x))equiv 0(mod x^{frac{n}{2}})]
因为他俩前(frac{n}{2})项一样。
那么：
[(B(x)-B_0(x))^2equiv 0(mod x^n)]
也就是：
[B^2(x)+B_0^2(x)-2B(x)B_0(x)equiv 0(mod x^n)]
两边同乘(A(x))
[B(x)+A(x)B_0^2(x)-2B_0(x)equiv 0(mod x^n)]
[B(x)equiv 2B_0(x)-A(x)B_0^2(x)(mod x^n)]
这样我们就得到了一个递归算法，使得每次回溯精度翻倍，而在最底层只有常数项，直接费马小定理即可。
注意这里我们求出的逆可能不仅有(n)项，但是我们只保留(n)项，这样虽然不精确，不过更高的项是没有用的。

4.(ln)
对于一个多项式(A(x))(必须保证常数项为1)，我们要求的(B(x))是满足：
[ln(A(x))equiv B(x)(mod x^n)]
两侧求导。
[B'(x)equiv frac{A'(x)}{A(x)}(mod x^n)]
这样可通过求逆和求导得到(B'(x))
再对(B'(x))不定积分即可得到(B(x))
同样省略大于(n)的高次项。

5.开方(sqr)
必须保证(A(0)=1)
首先，
对于任何一个多项式运算我们均可以设计一种多项式函数,使得运算结果(B(x))为方程：
[F(B(x))equiv 0(mod x^n)]
的解。
开方就是这样的。
我们要设计的函数就是：
[F(B(x))=B^2(x)-A(x)]
那么我们要求的就是：
[F(B(x))equiv 0(mod x^n)]
这样延续套路，假设我们求出了:
[F(B_0(x))equiv 0(mod x^{frac{n}{2}})]
我们对(F(B(x)))在(F(B_0(x)))处泰勒展开。
[F(B(x))equiv sumlimits_{i=0}^{n}frac{F^{(i)}(B_0(x))(B(x)-B_0(x))^i}{i!}(mod x^n)]
由于(B_0(x))和(B(x))的前(x^{frac{n}{2}})项是相同的，那么包括2次及以上的项均为0。
那么只剩下常数项和一阶导了。
[F(B(x))equiv F(B_0(x))+F'(B_0(x))(B(x)-B_0(x))equiv 0(mod x^n)]
消序展开得到：
[B(x)equiv frac{F(B_0(x))}{F'(B_0(x))}]
这种方法可以求解任何多项式函数。
也被称作多项式牛顿迭代。
在开方中就是：
[B(x)equiv B_0(x)-frac{A(x)}{2B_0(x)}]
这样的到了递归算法，回溯一次精度翻倍，递归到最后一层常数项为1。

6.(exp)
对于一个函数(A(x))(常数项为0)。
求：
[B(x)equiv e^{A(x)}(mod x^n)]
构造[F(B(x))=ln(B(x))-A(x)]
那么要求：
[F(B(x))equiv 0(mod x^n)]
利用牛顿迭代公式：
[B(x)equiv B_0-frac{F(B_0(x))}{F'(B_0(x))}equiv B_0(x)(1-ln(B_0(x))+A(x))]
同样是递归算法，递归一次精度翻倍，递归到最后一层常数项为1。

7.快速幂
不是直接倍增的。
要求：
[B(x)equiv A^k(x)(mod x^n)]
[lnB(x)equiv klnA(x)(mod x^n)]
这样我们求个(ln)，然后给他乘上(k)，然后再(exp)回去即可。

8.分治(FFT)
对于这样一个东西：
[f(n)=sumlimits_{i=0}^{n-1}f(i)g(n-i)]
求全部的(f(x))
朴素的算法是(O(n^2))的。
这里我们运用(CDQ)分治的方法求这个东西，就是(O(nlog^2n))的了，具体实现你们自己去思考。

9.(MTT)
我们有时候要求任意模数，这个时候有一种古老的方法是三模数(NTT),学长的话说：”这已经是时代的眼泪了“。
现在介绍(毛TT)，也就是所谓拆系数(FFT)
我们怕他爆(ll)，因为在(1e5)卷积的情况下，极限可以到达(10^{23})的地步，所以(MTT)出现了。
我们一般设一个常数(M=1<<15),接近(sqrt{1e9})的数。
[f(x)=sumlimits_{i=0}^{n}(A1(i)M+B1(i))x^i]
[g(x)=sumlimits_{i=0}^{n}(A2(i)M+B2(i))x^i]
那么((f*g))其实就被分开了:
[f_ig_{n-i}=(A1(i)M+B1(i))(A2(i)M+B2(i))=A1(i)A2(i)M^2+(A1(i)B2(i)+A2(i)B1(i))M+B1(i)B2(i)]
这样我们对于[A1,A2,B1,B2]分别做(4)次(FFT),对于上面三项分别求出,再做(3)次(FFT)。
这样每一项的大小都不会超过(10^{13}),再分别乘上相对应的(M)即可。
虽然乘上(M)还是爆了，不过这个时候是可以直接取模的。

另外：
对于如上的全部递归算法(求逆，exp,开方)复杂度均为(O(nlogn))
证明：
我们设复杂度为：(T(n))
那么：
[T(n)=T(frac{n}{2})+nlogn]
根据主定理得到：
[T(n)=nlogn]

给个板子（乘法，求导，积分，求逆，(ln)，(exp)，开方，另外一种分治(FFT),快速幂）。

int add(int x,int y) {return x+y>=mod?x+y-mod:x+y;}
int mul(int x,int y) {return 1LL*x*y%mod;}
int dic(int x,int y) {return x-y<0?x-y+mod:x-y;}
int qw(int a,int b)
{
    int ans=1;
    for(;b;b>>=1,a=mul(a,a)) if(b&1) ans=mul(ans,a);
    return ans;
}
struct Poly{
    int r[maxn],t1[maxn],t2[maxn],t3[maxn],f[maxn],g[maxn],h[maxn],A[maxn],B[maxn],C[maxn];
    pair<int,int> getlen(int n) {int mx,tim;for(mx=1,tim=0;mx<=n<<1;mx<<=1,++tim);return make_pair(mx,tim);}
    void clear(int *a,int n) {for(int i=0;i<n;++i) a[i]=0;}
    void NTT(int mx,int tim,int *a,int op)
    {
        for(int i=0;i<mx;++i) r[i]=(r[i>>1]>>1)|((i&1)<<tim-1);
        for(int i=0;i<mx;++i) if(i<r[i]) swap(a[i],a[r[i]]);
        for(int mid=1;mid<mx;mid<<=1)
            for(int i=0,len=mid<<1,wn=qw(3,(mod-1)/len);i<mx;i+=len)
                for(int j=0,w=1,x=a[i+j],y=mul(w,a[i+mid+j]);j<mid;++j,w=mul(w,wn),x=a[i+j],y=mul(w,a[i+mid+j]))
                    a[i+j]=add(x,y),a[i+mid+j]=dic(x,y);
        if(op==1) return ;
        reverse(a+1,a+mx);int re=qw(mx,mod-2);
        for(int i=0;i<mx;++i) a[i]=mul(a[i],re);
    }
    void Mul(int n,int *a,int *b,int *c)
    {
         pair<int,int> dr=getlen(n);
        clear(t3,dr.first);clear(h,dr.first);
        for(int i=0;i<n;++i) t3[i]=a[i],h[i]=b[i];
        NTT(dr.first,dr.second,t3,1);
        NTT(dr.first,dr.second,h,1);
        for(int i=0;i<dr.first;++i) c[i]=mul(t3[i],h[i]);
        NTT(dr.first,dr.second,c,-1);
        for(int i=n;i<dr.first;++i) c[i]=0;
    }
    vector<int> vecMul(vector<int> a,vector<int> b,int n)
    {
        pair<int,int> dr=getlen(n);
        for(int i=0;i<a.size();++i) t1[i]=a[i];
        for(int i=0;i<b.size();++i) t2[i]=b[i];
        NTT(dr.first,dr.second,t1,1);
        NTT(dr.first,dr.second,t2,1);
        for(int i=0;i<dr.first;++i) t1[i]=mul(t1[i],t2[i]);
        NTT(dr.first,dr.second,t1,-1);
        vector<int> res;
        for(int i=0;i<=n;++i) res.push_back(t1[i]);
        clear(t1,dr.first);clear(t2,dr.first);
        return res;
    }
    vector<int> DivideMul(int l,int r)
    {
        if(l==r) 
        {
            vector<int> res;
            res.push_back(1),res.push_back(mod-a[l]);
            return res;
        }
        int mid=l+r>>1;
        vector<int> fr=DivideMul(l,mid);
        vector<int> se=DivideMul(mid+1,r);
        return vecMul(fr,se,fr.size()+se.size());
    }
    void deriv(int n,int *a,int *b) {for(int i=1;i<n;++i) b[i-1]=mul(a[i],i);b[n-1]=0;}
    void integ(int n,int *a,int *b) {for(int i=1;i<n;++i) b[i]=mul(a[i-1],qw(i,mod-2));b[0]=0;}
    void getinv(int n,int *a,int *b)
    {
        if(n==1) return b[0]=qw(a[0],mod-2),void();
        getinv(n>>1,a,b);
        int mx=n<<1,tim=round(log(mx)/log(2));
        for(int i=0;i<n;++i) t1[i]=a[i],t2[i]=b[i];
        NTT(mx,tim,t1,1);
        NTT(mx,tim,t2,1);
        for(int i=0;i<mx;++i) t1[i]=mul(t1[i],mul(t2[i],t2[i]));
        NTT(mx,tim,t1,-1);
        for(int i=0;i<n;++i) b[i]=dic(add(b[i],b[i]),t1[i]);
        clear(t1,mx);clear(t2,mx);
    }
    void getln(int n,int *a,int *b)
    {
        int mx,tim;
        for(mx=1,tim=0;mx<=n<<1;mx<<=1,++tim);
        deriv(n,a,f);getinv(n,a,g);
        for(int i=0;i<mx;++i) t1[i]=f[i],t2[i]=g[i];
        NTT(mx,tim,t1,1);
        NTT(mx,tim,t2,1);
        for(int i=0;i<mx;++i) t1[i]=mul(t1[i],t2[i]);
        NTT(mx,tim,t1,-1);
        integ(n,t1,b);
        clear(t1,mx);clear(t2,mx);clear(f,mx);clear(g,mx);
    }
    void getexp(int n,int *a,int *b)
    {
        if(n==1) return b[0]=1,void();
        getexp(n>>1,a,b);
        int mx=n<<1,tim=round(log(mx)/log(2));
        for(int i=0;i<n;++i) A[i]=b[i];
        getln(n,A,B);
        for(int i=0;i<n;++i) B[i]=dic(a[i],B[i]);B[0]=add(B[0],1);
        NTT(mx,tim,A,1);NTT(mx,tim,B,1);
        for(int i=0;i<mx;++i) A[i]=mul(A[i],B[i]);
        NTT(mx,tim,A,-1);
        for(int i=0;i<n;++i) b[i]=A[i];
        clear(A,mx);clear(B,mx);
    }
    void getsqr(int n,int *a,int *b)
    {
        if(n==1) return b[0]=1,void();
        getsqr(n>>1,a,b);
        getinv(n,b,f);
        for(int i=0;i<n;++i) g[i]=a[i];
        pair<int,int>dr=getlen(n);
        NTT(dr.first,dr.second,f,1);NTT(dr.first,dr.second,g,1);
        for(int i=0;i<n<<1;++i) f[i]=mul(f[i],g[i]);
        NTT(dr.first,dr.second,f,-1);
        for(int i=0;i<n;++i) b[i]=mul(add(b[i],f[i]),qw(2,mod-2));
        for(int i=0;i<n<<1;++i) f[i]=g[i]=0;
    }
    void getpow(int n,int k,int *a,int *b)
    {
        getln(n,a,h);
        for(int i=0;i<n;++i) h[i]=mul(h[i],k);
        getexp(n,h,b);
    }
}fntt;