「总结」多项式生成函数相关(1)
实在是太毒瘤了。
大纲。
多项式生成函数相关
默认前置:微积分,各种数和各种反演,FFT,NTT,各种卷积,基本和式变换。
主要内容:
泰勒展开,级数求和,牛顿迭代,主定理。 //例题:在美妙的数学王国中畅游,礼物
多项式全家桶:乘法,求逆,求导,积分,分治,ln,exp,fwt,MTT。 //城市规划,图的价值,染色,遗失的答案,按位或,随机游走。
生成函数:普通型生成函数,指数型生成函数计数原理。 //猎人杀,遗忘的集合,生成树计数
例题。
一、泰勒展开和级数求和
1.泰勒展开
即对于任何函数\(f(x)\),如果在\(x_0\)处\(n\)阶可导,那么满足如下公式:
这里当\(x_0\)为0的时候被称作麦克劳林公式。
先推导麦克劳林公式
即:
这里只证明多项式函数的正确性(其实是因为任意函数太难证了吧)。
设多项式函数:
那么:
而我们的\(x_0=0\)
也就是说每阶导数都只有常数项。
即:
所以
那如果\(x_0!=0\)呢?
这时候就是真正的泰勒展开了。
我们照旧设多项式函数:
将之更换为关于\(w=x-x_0\)的多项式。
用麦克劳林公式展开:
也就是说:$$b_i=\frac{g^{(i)}}{i!}$$
那么:
所以:
这样的话就是说:
证毕。
2.一些级数求和的公式
其实个人认为泰勒展开是级数求和的逆运算。
二、牛顿迭代
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\)处泰勒展开。
得到:
我们只保留其线性部分。
个人猜测牛顿这样做的原因是因为这样迭代一次之后\(x_1\)仅有一个解,如果保留二次项部分,虽然会让迭代的次数减少,不过此时有两个解,剩余部分的迭代增加量变成了指数级别的。
那么也就是说我们需要求:
此时\(x_1=x_0-\frac{f(x_0)}{f'(x_0)}\)
不断迭代下去就可以不断逼近解\(x_z\)
考虑为什么。
由于二阶导符号不变,那么也就是说函数的凹凸性不变,这样的话一阶导的就是单调的了。
其实就是函数切线的斜率单调,我们求出的\(x_1\)就是切线的解,这样就可以不断的逼近\(x_z\)。
2.多项式牛顿迭代
本来这个东西使用来求连续函数的。
现在我们把它搞到离散数学里面来,利用相似的思想。
求解如下方程。
按照套路来。
我们假设已经求出了\(G'(x)\)
满足:
对函数\(F(G(x))\)在\(G'(x)\)处泰勒展开。
可以发现:
所以:
那么:
这样也就是说,不准确的部分只有一阶导部分。
也就是:
每次迭代精度翻倍了。
可以递归求解。
复杂度是:\(T(n)=nlogn+T(\frac{n}{2})=nlogn\)。
可能要疑问这个是干什么用的。
下面的全家桶就知道了。
讲个题:https://www.cnblogs.com/Lrefrain/p/11921361.html.
三、主定理
主定理用来证明一些分治递归算法的复杂度。
即。
对于如下的复杂度:
设\(f(n)=n^d\)
有:
偷一张算导的图:
是递归示意图了。
那么可以发现总的复杂度可以这样来表示:
讨论一下。
1.\(a<b^d\)后面那个东西是小数,所以说是收敛的,是一个常数,所以\(T(n)=n^d\)
2.\(a=b^d\)后面全是1,所以\(T(n)=n^dlog_bn\)
3.\(a>b^d\)后面东西,对于最后一项来说,前面的全是常数,所以舍去前面的项。
推一下:
四、多项式全家桶
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)\)是满足:
假设我们已经求出了一个\(B'(x)\)使得:
那么:
因为他俩前\(\frac{n}{2}\)项一样。
那么:
也就是:
两边同乘\(A(x)\)
这样我们就得到了一个递归算法,使得每次回溯精度翻倍,而在最底层只有常数项,直接费马小定理即可。
注意这里我们求出的逆可能不仅有\(n\)项,但是我们只保留\(n\)项,这样虽然不精确,不过更高的项是没有用的。
4.\(ln\)
对于一个多项式\(A(x)\)(必须保证常数项为1),我们要求的\(B(x)\)是满足:
两侧求导。
这样可通过求逆和求导得到\(B'(x)\)
再对\(B'(x)\)不定积分即可得到\(B(x)\)
同样省略大于\(n\)的高次项。
5.开方\(sqr\)
必须保证\(A(0)=1\)
首先,
对于任何一个多项式运算我们均可以设计一种多项式函数,使得运算结果\(B(x)\)为方程:
的解。
开方就是这样的。
我们要设计的函数就是:
那么我们要求的就是:
这样延续套路,假设我们求出了:
我们对\(F(B(x))\)在\(F(B_0(x))\)处泰勒展开。
由于\(B_0(x)\)和\(B(x)\)的前\(x^{\frac{n}{2}}\)项是相同的,那么包括2次及以上的项均为0。
那么只剩下常数项和一阶导了。
消序展开得到:
这种方法可以求解任何多项式函数。
也被称作多项式牛顿迭代。
在开方中就是:
这样的到了递归算法,回溯一次精度翻倍,递归到最后一层常数项为1。
6.\(exp\)
对于一个函数\(A(x)\)(常数项为0)。
求:
构造$$F(B(x))=ln(B(x))-A(x)$$
那么要求:
利用牛顿迭代公式:
同样是递归算法,递归一次精度翻倍,递归到最后一层常数项为1。
7.快速幂
不是直接倍增的。
要求:
这样我们求个\(ln\),然后给他乘上\(k\),然后再\(exp\)回去即可。
8.分治\(FFT\)
对于这样一个东西:
求全部的\(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*g)\)其实就被分开了:
这样我们对于$$A1,A2,B1,B2$$分别做\(4\)次\(FFT\),对于上面三项分别求出,再做\(3\)次\(FFT\)。
这样每一项的大小都不会超过\(10^{13}\),再分别乘上相对应的\(M\)即可。
虽然乘上\(M\)还是爆了,不过这个时候是可以直接取模的。
另外:
对于如上的全部递归算法(求逆,exp,开方)复杂度均为\(O(nlogn)\)
证明:
我们设复杂度为:\(T(n)\)
那么:
根据主定理得到:
给个板子(乘法,求导,积分,求逆,\(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;