FFT,NTT目害学笔记
我的YY理解
FFT可以处理形如:
的卷积形式。暴力是\(O(n^2)\)而FFT可以做到\(O(nlogn)\)级别.本质就是多项式相乘
具体用到了复数在二维平面上的表示,形如\(cos\theta+isin\theta\)的形式加减乘,实现点值表示下的\(O(n)\)快速卷积。
以单位根为x代入到多项式中,由于其特殊性质就可以求得系数。但是只能求长度(最大次数)为\(2^n-1\)的数,这也没有关系,只取需要的项即可。
为了计算需要将下标按照奇偶性分到\(1\)~\(n/2\),\(n/2+1\)~\(n\),
这样会发现每个数被分到的位置就是二进制下的数倒转过来表示的数。
例如在\(n=3\)时,\(6(110)\)->\(3(011)\),\(n=4\)时\(5(0101)\)->\(10(1010)\)等。
于是可以预处理出来每个数要分到哪里去,FFT时直接交换即可。
同时我们发现在mod意义下,原根和单位根有着一样的性质,于是NTT就诞生了。
NTT的话就是将FFT中的复数改成mod意义下的数,其他的都一样,不会丢精度,听说速度比FFT快,但是我交到luogu上似乎FFT更快?
ps:NTT中取的膜数必须满足
多项式乘法
是FFT板子题。
放个码:
#include<cstring>
#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>
namespace EMT{
typedef long long ll;typedef double db;//(double)clock() / (double)CLOCKS_PER_SEC;
#define pf printf
#define F(i,a,b) for(register int i=a;i<=b;i++)
#define D(i,a,b) for(register int i=a;i>=b;i--)
inline int read(){int x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}while(ch>='0'&&ch<='9')x=x*10+ch-'0',ch=getchar();return x*f;}
inline void file(){freopen("in.in","r",stdin);freopen("my.out","w",stdout);}
inline int max(int a,int b){return a>b?a:b;}inline int min(int a,int b){return a<b?a:b;}
inline void pi(int x){pf("%d ",x);}inline void pn(){pf("\n");}
const int N=1e7+10;
const db PI=acos(-1.0);
struct comp{
db x,y;
friend comp operator +(comp a,comp b){return (comp){a.x+b.x,a.y+b.y};}
friend comp operator -(comp a,comp b){return (comp){a.x-b.x,a.y-b.y};}
friend comp operator *(comp a,comp b){return (comp){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
}a[N],b[N];
int lim=1,n,m,ans[N],l,r[N];
inline void fft(comp *a,int tp){
F(i,0,lim-1)if(i<r[i])std::swap(a[i],a[r[i]]);
for(int mid=1;mid<lim;mid<<=1){
comp wn=(comp){cos(PI/mid),tp*sin(PI/mid)};
for(int j=0;j<lim;j+=(mid<<1)){
comp w=(comp){1,0};
for(int k=0;k<mid;k++,w=w*wn){
comp x=a[j+k],y=w*a[j+mid+k];
a[j+k]=x+y;
a[j+mid+k]=x-y;
}
}
}
}
inline short main(){
n=read(),m=read();
while(lim<=n+m)lim<<=1,l++;
F(i,0,n)a[i].x=read();
F(i,0,m)b[i].x=read();
F(i,0,lim-1)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
fft(a,1),fft(b,1);
F(i,0,lim)a[i]=a[i]*b[i];
fft(a,-1);
F(i,0,n+m)pi((int)(a[i].x/lim+0.5));
return 0;
}
}
signed main(){return EMT::main();}
快速傅立叶之二
将b数组倒置,发现
因为\(b_{n+1}\)~\(b_{n+k}\)都是0
\(a_{n+1}\)~\(a_{n+k}\)都是0
所以C还可以写成
设
则
其中D满足卷积形式。于是可以用fft求出D下标+n输出即可。
力
\(q_j\)被除掉之后就消掉了。
另外,设\(g_i=\frac{1}{i^2}\),\(qq_i=q_{n-i}\)
所以得到式子:
前一项已经是卷积形式了,考虑后一项:
原式=
设这个为\(r_j\)
再设
所以\(r_j=R_{n-j}\)
于是左右两项都出来了,相减就是答案
ps:最后一定要除以\(lim\)!!!,没发现调了好久QAQ
Normal
无奖竞猜哪个是暴力
暴力都能过,有点迷惑...
这个题就是点分治上利用FFT算出最终树上每个路径长度的个数,
最终答案就是
每个长度上的点被选到的概率是\(\frac{1}{dis+1}\),乘上长度二倍(正反),最后加上单点的,由于单点没正反所以不乘二倍。
ps:清空的时候一定用到多少就清空多少,又差点卡住了。
求和
由于i>j时S(i,j)=0所以可以把j的范围扩大到n
原式就是
由第二类斯特林数性质可得:
移项可得:
当\(j-k=1\)时,\(\frac{\sum\limits_{i=0}^{n}(j-k)^i}{(j-k)!}=n+1\)
当\(j-k=0\)时,\(\frac{\sum\limits_{i=0}^{n}(j-k)^i}{(j-k)!}\)中除了\(i=0\)别的都是零,众所周知,\(0^0=1\)所以原式为1
当\(j-k>1\)时,由等比数列求和公式可得
于是\(\sum\limits_{j=0}^{n}2^j\times (j!)\)用\(O(n)\)可求出,剩下的:
设
所以原式就变成了
剩下的就符合卷积形式,可用NTT求出。于是本题就做完了。
ps:\(0^0=1\)!!!!差点又调了半天。
Triple
这个题用到了生成函数与容斥,
生成函数在本题中的应用是:
设数列\(A=\sum a_ix^i\)
其中x没有实际意义,i就是这个斧头的价值,a_i是出现的方案数,
考虑将A立方,但这样得出来的系数并不是答案,因为有自乘的成分,所以要容斥,
设\(B=\sum a_{i}x^{2i}\),\(C=\sum a_{i}x^{3i}\)则
另外扩展lim时要按照值域扩展,就不是n了。要不是我有数据,差点又栽了。
万径人踪灭
设\(f_i\)表示以i为对称轴有多少个对称点对。字符串中如果\(s_i\)为a则\(a_i=1\)否则\(b_i=1\)
显然
也就是
这样就满足卷积形式了,fft两次求和就能得到f数组,要刨掉连续的对称串,所以可以用马拉车求一下。答案就是:
序列统计
考虑首先考虑在mod意义下的乘法,于是可以想到原根,对于m的原根g来说,\(g^{0,1,2..m-2}%m\)都是不同的值,将问题转化到了原根的几次幂上,再考虑将乘法转化为加法,于是对每个数取膜数意义上的log,就可以满足卷积形式用NTT求解了。
另外由于要算卷积的n次方,可以模拟快速幂的形式将答案求出来。
ps:算出来答案后一定要乘逆元!又迷惑了好久
按位或
依照min-max容斥原理,
于是用FWT求解max即可。
染色
又忘除以lim了!!!淦
考虑容斥
设\(f_i\)表示至少有i种==s的方案数,则:
设\(g_i\)表示恰好有i种==s的方案数,由二项式反演,得:
其中lim是最多有几种能够==s,\(lim=min{m,n/s}\)
将C拆开得:
设
则:
设\({F}'=F_{lim-i}\)
则
设
则满足卷积形式,可以用ntt求解。
最后乘一个\(w_i\)求和即可。
分治FFT
用来处理类似:
的形式,
具体做法是类似CDQ分治,考虑左半边对右半边的贡献,递归求解即可求出\(f\)数组。
多项式乘法逆
问题是给出F,求
中的G,
那么显然有:
设想我们已经得到了H数组,使得:
上下式相减可得:
那么有:
两边平方可得:
两边同时乘上\(F(x)\)并移项可得:
于是递归以此求解H再求得下一个H(G)即可。
显然只有一项时,\(H(x)\)就是\(F(x)\)的逆元,于是边界条件也有了。
【UR #3】链式反应
题目变成人话就是给出一棵满足堆性质的树,对于儿子方向每个结点有两条红色边(中子轰击),c 条黄色边(破坏死光),c 属于集合 A,统计方案数。
设\(f_i\)表示大小为i的子树的方案数。i=1时显然什么都不能做,则\(f_i=1\)
否则:
具体含义就是,我这个i点爆了,剩下i-1个点,从里面选大小是j的子树,再从剩下的选大小是k的子树,并且要符合属于集合A的条件。
由于有可能当前我选的两个子树交换一下就是另外一种方案但这其实不能算,所以要除以2.
暴力转移是\(n^3\)的,考虑优化。
我们把C拆开,设
于是f就能表示成:
移项可得:
观察到g和\(\frac{2f_i}{(i-1)!}\)都满足卷积形式,于是进行分治NTT就能求解了。
由于g是f和f相乘,当l到mid卷mid+1到r时mid+1到r区间的值并没有准确算出来,所以再算mid+1到r时补偿上去,将系数\(\times 2\)即可。(link)
图的价值
考虑枚举所有点的贡献,于是答案就是:
具体意思就是和这个点无关的边随便连,剩下的枚举连几条边,从\(n-1\)条里选,贡献是k次方,每个点互不影响,然后乘n。
前面那两项好说,主要问题就是计算:
首先我们回顾一下第二类斯特林数。\(S_n^m\)表示n个标号不同的小球分到m个相同的盒子里,盒子非空的方案数。
于是有个通项公式:
于是这显然是一个卷积,可以ntt求解。
然后回过头来,如果不要求盒子非空,那么首先显然是\(m^n\),其次还可以枚举几个盒子不空,用\(\sum\limits_{i=0}^{m}P_m^iS_n^i\)表示,其中P是排列数。于是得出一个柿子:
于是原式化为:
移个项转化一下就变成:
显然就可以\(O(k\log n)\)求解了。
拉格朗日插值2
考虑拉格朗日插值的公式则:
然后由于取值连续于是有:
然后我们设:
这样的话有:
于是答案的话可以直接让下标加n把系数乘一下就行了。
残缺的字符串
FFT还能匹配字符串,新技能get啊!
具体的话,由于有通配符,我们把通配符的权值设成0,其他的就是\(c-'a'+1\)
然后A在B中的匹配函数就是:
当\(p(x)=0\)时,说明存在一个以x结尾的和A匹配的子串。
但是这并不能卷,怎么办?老套路,翻转:
这样卷下来,\([m-1,n-1]\)的位置满足\(p(i)=0\)的就是可以匹配的地方。
CF553E Kyoya and Train
由于这是个DAG,我们设\(f_{i,j}\)表示j时间到i点最少需要花费多少到n。
显然n点的出边可以忽略。
我们先跑一遍floyed(话说好像dij就行)求出每个点之间的最短距离,这样:
对于n点\(j<=t\)时\(f_{i,j}=0\)否则\(f_{i,j}=x\)
对于其他点\(j>=t\)时\(f_{i,j}=dis_{i,n}+x\),因为反正也超时了,不如走最短路程。
那么问题就是处理其他点\(j< t\)时的情况。
由于状态转移和时间有关,我们可以从高时间转到低时间,由此可以想到分治FFT,
设\(g_{i,j}\)表示j时间走完第i条边的最小期望时间,那么有:
其中\(w_k\)是边权。那么我们处理好\([mid+1,r]\)之后贡献到\([l,mid]\)里面即可。
任意模数多项式乘法
前置知识很烦。
由于最后乘出来的数一定是\(1e9\times 1e9 \times 1e5<=2e23\)的,所以我们可以自己选择用三个\(1e8/1e9\)级别的质数模数分别做一次NTT,最后用中国剩余定理合并即可。(由于不会中国剩余定理所以我用了EXCRT)
基本的NTT不说了,主要问题是合并,由于值域范围太大了,所以可以选择使用int128来合并,但是实测还是会爆,所以还要用龟速乘来实现合并。
这个不知道有啥好证的。。。EXCRT?
为了避免:
为什么没有证明?读者很疑惑!求博主更更!
于是证一下,刚好之前都没证过。(关于我在FFT专题里证明EXCRT这档事)
对于同余方程:
我们有:
把右边两个等式拿出来并移个项:
我们设\(g\)为\(\gcd(b_1,b_2)\),\(p_1=\frac{b_1}{g},p_2=\frac{b_2}{g}\)。
由裴蜀定理得若\((a_2-a_1)\%g=0\)则一定有解,否则无解。
然后两边同除\(g\)则:
我们珂以用exgcd先求出:
的一组解,则:
然后就有了
两个方程合并成了一个方程:
一直合并下去最后方程即可任取解了。
多项式开根
给定一个n−1次多项式A(x),求一个在\(\bmod\ x^n\)意义下的多项式B(x),使得\(B^2(x) \equiv A(x)\bmod x^n\)。若有多解,请取零次项系数较小的作为答案。
设
第二个式子移项平方得到
然后多项式乘法逆即可。
多项式对数函数
设:
求导可得:
所以多项式求个逆在积分回去就好了~
遗忘的集合
考虑写出每个物品的生成函数,枚举这个物品有几个,也就是\(\sum\limits_{j=0}^{inf}x^{ij}=\frac{1}{1-x^i}\)
然后设每个物品是否存在,即为:\((\frac{1}{1-x^i})^{t_i}\)
则拼起来是x的方案数是这些东西的卷积,即为
这样无法做出,考虑两边同时求ln,则:
然后对右边求导得到:
积分得到:
带回原式得到:
然后多项式求ln之后即可求\(t_i\)。
多项式指数函数
首先有泰勒展开:
\(f^n(0)\)表示对0这个点进行n阶求导
然后对多项式进行牛顿迭代:
对于一个函数\(G(x)\),求满足条件的\(F(x)\)使得\(G(F(x))\equiv 0(\bmod\ x^n)\)
首先当n等于1的时候\(G(F(x))\equiv 0(\bmod x)\)单独求出,
然后假设已经求出了:
然后把\(G(F(x))\)在\(F_0(x)\)处进行泰勒展开:
首先有\(F(x)\)和\(F_0(x)\)至少0到\(\left\lceil\frac{n}{2}\right\rceil\)项相同,那么\((F(x)-F_0(x))^2\)的最低项已经到达了\(2\left\lceil\frac{n}{2}\right\rceil\)次,所以在\(\bmod\ x^n\)意义下没有了。
于是得到:
然后因为\(G(F(x))=0\)所以有:
然后进入正题。
给定\(A(x)\),求\(F(x)\equiv e^{A(x)} (\bmod\ x^n)\)。
首先对两边取ln得\(\ln F(x)-A(x)=0\)
设\(G(F(x))=\ln F(x) -A(x)\)
然后A已知看做常量,求导得到:\(G'(F(x))=\frac{1}{F(x)}\)
带入上面公式得
然后就能进行求解了。
好像最近学了导数之后只做了做板子,还没有练什么题,把各个板子放一下:
MTT
如果不取模由于数太大FTT精度会炸裂,所以可以把数字拆开系数然后再乘。
设要对\(F(x),G(x)\)进行乘法,然后拆出来一个\(\sqrt{n}\)级别的数bas,变成:
然后答案就变成了:
这样的话需要7次DFT即可求解,然后通过一些我不会的黑科技可以变成4次,背板子就完事了
inv,ln,exp,sqrt
namespace poly{
int r[N];
inline int ksm(int a,int b){
int ans=1;
while(b){
if(b&1)ans=1ll*a*ans%mod;
a=1ll*a*a%mod;b>>=1;
}return ans;
}
inline void ntt(int *a,int tp,int lim){
for(int i=0;i<lim;i++)if(i<r[i])std::swap(a[i],a[r[i]]);
for(int mid=1;mid<lim;mid<<=1){
int wn=ksm(tp==1?g:gi,(mod-1)/(mid<<1));
for(int j=0;j<lim;j+=(mid<<1)){
int w=1;
for(int k=0;k<mid;k++,w=1ll*w*wn%mod){
int x=a[j+k],y=1ll*w*a[j+mid+k]%mod;
a[j+k]=(x+y)%mod;
a[j+mid+k]=(x-y+mod)%mod;
}
}
}
if(tp==-1){
int inv=ksm(lim,mod-2);
for(int i=0;i<lim;i++)a[i]=1ll*a[i]*inv%mod;
}
}
inline void getinv(int *a,int *b,int n){
static int A[N],B[N];
b[0]=ksm(a[0],mod-2);int len,lim,l;
for(len=1;len<(n<<1);len<<=1){
lim=1,l=0;
while(lim<=len)lim<<=1,l++;
for(int i=0;i<lim;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
for(int i=0;i<len;i++)A[i]=a[i],B[i]=b[i];
ntt(A,1,lim),ntt(B,1,lim);
for(int i=0;i<lim;i++)b[i]=1ll*(2ll-1ll*A[i]*B[i]%mod+mod)*B[i]%mod;
ntt(b,-1,lim);
for(int i=len;i<lim;i++)b[i]=0;
}
for(int i=0;i<len;i++)A[i]=B[i]=0;
for(int i=n;i<len;i++)b[i]=0;
}
inline void getdao(int *a,int *b,int n){
for(int i=1;i<n;i++)b[i-1]=1ll*i*a[i]%mod;b[n-1]=0;
}
int inv[N];
inline void init(int n){
inv[0]=inv[1]=1;
for(int i=2;i<=n;i++)inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
}
inline void jif(int *a,int *b,int n){
for(int i=1;i<n;i++)b[i]=1ll*inv[i]*a[i-1]%mod;b[0]=0;
}
inline void getln(int *a,int *b,int n){
static int A[N],B[N];
getdao(a,A,n);getinv(a,B,n);
int lim=1,l=0;
while(lim<=(n<<1))lim<<=1,l++;
for(int i=0;i<lim;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
ntt(A,1,lim),ntt(B,1,lim);
for(int i=0;i<lim;i++)A[i]=1ll*A[i]*B[i]%mod;
ntt(A,-1,lim);jif(A,b,n);
}
inline void getsqrt(int *a,int *b,int n){
static int A[N],B[N];
b[0]=1;int len,lim=1,l=0;
for(len=1;len<(n<<1);len<<=1){
lim=1,l=0;
while(lim<=len)lim<<=1,l++;
for(int i=0;i<len;i++)A[i]=a[i];
getinv(b,B,lim>>1);
for(int i=0;i<lim;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
ntt(A,1,lim),ntt(B,1,lim);
for(int i=0;i<lim;i++)A[i]=1ll*A[i]*B[i]%mod;
ntt(A,-1,lim);
for(int i=0;i<len;i++)b[i]=1ll*(b[i]+A[i])*inv2%mod;
for(int i=len;i<lim;i++)b[i]=0;
}for(int i=0;i<len;i++)A[i]=B[i]=0;
for(int i=n;i<len;i++)b[i]=0;
}
inline void getexp(int *a,int *b,int n){
static int A[N];
b[0]=1;int len,lim,l;
for(len=1;len<(n<<1);len<<=1){
lim=1,l=0;getln(b,A,len);
A[0]=(a[0]+1-A[0]+mod)%mod;
for(int i=1;i<len;i++)A[i]=(a[i]-A[i]+mod)%mod;
while(lim<=len)lim<<=1,l++;
for(int i=0;i<lim;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
ntt(A,1,lim),ntt(b,1,lim);
for(int i=0;i<lim;i++)b[i]=1ll*b[i]*A[i]%mod;
ntt(b,-1,lim);
for(int i=len;i<lim;i++)b[i]=A[i]=0;
}
}
}
MTT
namespace poly{
struct comp{
db r,i;
friend comp operator +(comp a,comp b){return {a.r+b.r,a.i+b.i};}
friend comp operator -(comp a,comp b){return {a.r-b.r,a.i-b.i};}
friend comp operator *(comp a,comp b){return {a.r*b.r-a.i*b.i,a.i*b.r+a.r*b.i};}
inline comp conj(){return {r,-i};}
}w[N];
int r[N];
inline void fft(comp *a,int tp,int lim){
for(int i=0;i<lim;i++)w[i]={cos(2.0*i*PI/lim),sin(2.0*i*PI/lim)};
if(tp==-1)for(int i=0;i<lim;i++)w[i]=w[i].conj();
for(int i=0;i<lim;i++)if(i<r[i])std::swap(a[i],a[r[i]]);
for(int mid=1;mid<lim;mid<<=1){
for(int j=0;j<lim;j+=(mid<<1)){
for(int k=0;k<mid;k++){
comp x=a[j+k],y=w[lim/(mid<<1)*k]*a[j+mid+k];
a[j+k]=x+y;
a[j+mid+k]=x-y;
}
}
}
}
inline void mul(int *f,int *g,int *h,int n){
static comp a[N],b[N],c[N],d[N];int lim=n;
for(int i=0;i<n;i++)a[i].r=f[i]>>15,a[i].i=f[i]&32767;
for(int i=0;i<n;i++)c[i].r=g[i]>>15,c[i].i=g[i]&32767;
fft(a,1,lim),fft(c,1,lim);
for(int i=1;i<n;i++)b[i]=a[n-i].conj();b[0]=a[0].conj();
for(int i=1;i<n;i++)d[i]=c[n-i].conj();d[0]=c[0].conj();
for(int i=0;i<n;i++){
comp aa=(a[i]+b[i])*(comp){0.5,0},bb=(a[i]-b[i])*(comp){0,-0.5};
comp cc=(c[i]+d[i])*(comp){0.5,0},dd=(c[i]-d[i])*(comp){0,-0.5};
a[i]=aa*cc+(comp){0,1}*(aa*dd+bb*cc),b[i]=bb*dd;
}fft(a,-1,lim),fft(b,-1,lim);
for(int i=0;i<lim;i++){
int aa=(ll)(a[i].r/n+0.5)%mod;
int bb=(ll)(a[i].i/n+0.5)%mod;
int cc=(ll)(b[i].r/n+0.5)%mod;
h[i]=((1ll*aa*(1<<30)+1ll*bb*(1<<15)+cc)%mod+mod)%mod;
}
}
inline int ksm(int a,int b){
int ans=1;
while(b){
if(b&1)ans=1ll*a*ans%mod;
a=1ll*a*a%mod;b>>=1;
}return ans;
}
int inv[N];
inline void init(int n){
inv[0]=inv[1]=1;
F(i,2,n)inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
int lim=1;while(lim<=n)lim<<=1;
}
inline void getinv(int *a,int *b,int n){
static int fh[N],fhh[N],A[N],B[N];
b[0]=ksm(a[0],mod-2);int len,lim=1,l=0;
for(len=1;len<(n<<1);len<<=1){
lim=1,l=0;
while(lim<=len)lim<<=1,l++;
for(int i=0;i<lim;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
for(int i=0;i<len;i++)A[i]=a[i],B[i]=b[i];
mul(A,B,fh,lim),mul(fh,B,fhh,lim);
for(int i=0;i<len;i++)b[i+len]=0,b[i]=(2ll*b[i]-fhh[i]+mod)%mod;
}
for(int i=0;i<len;i++)A[i]=B[i]=0;
for(int i=n;i<len;i++)b[i]=0;
}
inline void getdao(int *a,int *b,int n){
for(int i=1;i<n;i++)b[i-1]=1ll*i*a[i]%mod;b[n-1]=0;
}
inline void jif(int *a,int *b,int n){
for(int i=1;i<n;i++)b[i]=1ll*inv[i]*a[i-1]%mod;b[0]=0;
}
inline void getln(int *a,int *b,int n){
static int A[N],B[N],C[N];
getdao(a,A,n),getinv(a,B,n);
int lim=1,l=0;while(lim<=n<<1)lim<<=1,l++;
for(int i=0;i<lim;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
mul(A,B,C,lim);
jif(C,b,n);
}
}
生成树计数
转化一下题意就是:
也就是连通块里面的点谁连出去随便选。
然后把问题转化到\(purfer\)序列上,\(c_i\)表示i连通块在序列中出现了几次,有\(d_i=c_i+1\)
则为:
然后把前面的\(\frac{1}{\prod c_i!}\)分配进去有:
然后考虑把前后分别设成生成函数:
然后答案就变成了:
然后对A,B函数分别积分,以B函数为例:
然后把\(j^m\)用第二类斯特林数展开:
发现后面是泰勒展开式:
然后求个导,得到:
然后A函数也类似:
然后设\(\alpha_i(x)=\frac{A_i(x)}{e^{a_ix}},\beta_i(x)=\frac{B_i}{e^{a_ix}}\)
然后把答案转化成:
然后可以把后面的进行分治,暴力通分即可。通分后发现分母正好和前面的消掉了,所以只用分子然后乘上系数即可。
青蕈(\(\grave{xun}\))领主
有几点结论:
- 两个连续区间不可能有交,只能包含与被包含
证明就是由于连续区间值域连续所以如果有交就可以合并成一个大区间,所以不合法。 - 如果把一段连续区间看成一个点,包含它的最短区间看成父亲,那么这些区间可以看做一棵树
然后先考虑不合法情况,如果\(L_1\not=1|L_n\not=n\)或者区间有交就无解。
然后考虑求出答案。
设\(f_i\)表示长度为i+1的排列中\(\forall i\in[1,n]\ L_i=1\)的排列数量,
那么设\(c_i\)表示i有几个直接儿子,由于区间互不相交,显然答案就是\(\prod f(c_i)\)
\(c_i\)可以很容易地通过单调栈求出,然后考虑如何求f。
考虑从长度为\(n\)的排列转移到\(n+1\)的排列,也就是原排列整体加一,把\(1\)插到排列\([2,n+1]\)中。
- 原排列本来就合法
那么可以插到除了2左右的部分,即\(f(n-1)\times (n-1)\)
合法性?考虑如果插完之后使得\([l,r]\)成为\([1,x]\)并且\(x\)不是最大值,那么原序列\([l,r-1]\)一定是\([1,x-1]\),与原排列合法相矛盾。 - 原排列不合法,并且原排列刚好有一个不经过最大值的连续区间(因为多于一个就不可能让它变得合法了)。
设这个原排列中非法区间值域是[x,x+l-1],其中\(x,l>=2,x+l-1<n\to 2\leq x\leq n-l,2\leq l \leq n-2\),那么x的取值有\(n-l-1\)种,再考虑这样的非法区间的方案数,考虑插完1之后这个区间就合法了,那么有\(f(l)\)的方案数,然后其它区间和它并起来同样合法,有\(f(n-l)\)的方案数,那么这个贡献就是\(\sum\limits_{j=2}^{n-2}(j-1)f(j)f(n-j)\)
那么有式子:
然后可以进行分治求解。
但是这里有“自己卷自己的形式”,本来\([l,mid]\)卷\([mid+1,r]\)无法在递归左区间时产生贡献,只能在递归到右半边时再计算贡献,然后贡献是\((j-1)f_jf_{i-j}+(i-j-1)f_{i-j}f_j=(i-2)f_{i-j}f_{j}\),那么只需要卷一次乘上\(i-2\)就好了。
排列计数
考虑设\(f(k)\)为至少k个升高的方案数,它的生成函数是:\(\sum\limits_{i=1}\frac{x^i}{i!}=e^x-1\)然后拼成m=n-k段,就是\((e^x-1)^m\)
然后根据二项式定理得到:
然后可以进行一次卷积出f。
然后二项式反演再卷积一次得到g即可。
多项式快速幂
求\(G(x)\equiv F^k(x)\pmod{x^n}\)
两边同时求ln可得\(\ln G(x)\equiv k\ln F(x)\pmod{x^n}\)
那么只需要对F求ln之后所有系数和k相乘对mod取模即可。
为什么可以直接对mod取模?
我不知道因为k到外面之后乘以一个多项式的表现就是对每个系数乘k,由于系数对mod取模所以可以这样做。
如果第一项是1可以直接这样做,
否则可以找到第一个不为0的项之后把所有系数都除以它最后再乘回去即可,注意这样也需要平移,把握好边界。