題目大意

　　對于一個長度為 \(n\) 的排列 \(p\)，我們稱一個區(qū)間 \([l,r]\) 是連續(xù)的當(dāng)且僅當(dāng) \((\max_{l\leq i\leq r}a_i)-(\min_{l\leq i\leq r}a_i)=r-l\)。

　　對于兩個排列 \(p_1,p_2\)，我們稱這兩個排列是等價的，當(dāng)且僅當(dāng)他們的長度相同且連續(xù)區(qū)間的集合相同。

　　給你 \(n\)，對于 \(1\leq i\leq n\)，所有求 \(i\) 階排列形成的等價類個數(shù)對 \(p\) 取模的值。

　　\(n\leq 100000,p=k\times2^{18}+1,k\in \mathbb {N},p\) 是質(zhì)數(shù)。

題解

　　對于一個區(qū)間 \([l,r]\)，如果這個區(qū)間的所有子區(qū)間都是連續(xù)的，就稱這個區(qū)間為強(qiáng)連續(xù)區(qū)間。

　　每次我們找一個最大的強(qiáng)連續(xù)區(qū)間，把這個區(qū)間內(nèi)所有點縮到一起。

　　如果找不到，就找一個最小的連續(xù)區(qū)間，把這個區(qū)間內(nèi)所有點縮到一起。

　　這個過程就形成了一個樹結(jié)構(gòu)。

• 總共有 \(n\) 個葉子。
• 對于一個代表強(qiáng)連續(xù)區(qū)間的點，這個點的兒子個數(shù) \(\geq 2\)。
• 對于一個代表連續(xù)區(qū)間的點，這個點的兒子個數(shù) \(\geq 4\)。（可以發(fā)現(xiàn)，長度為 \(\leq 3\) 的序列是一定有強(qiáng)連續(xù)區(qū)間的。）

　　每一個排列 \(p\) 對應(yīng)了一棵樹。

　　對于一棵樹，可以構(gòu)造出一個排列 \(p\)。

　　那么就只用對這棵樹計數(shù)了。

　　顯然可以 \(O(n^2)\) DP。

　　記 \(f_i\) 為長度為 \(i\) 的排列的方案數(shù)，\(F(x)=\sum_{i\geq 0}f_ix^i\)，那么就有：
\[ F(x)=\frac{(F(x)^2+F(x)^4)}{1-F(x)}+x \]
　　用牛頓迭代法解這個方程即可。
\[ (F(x)^2+F(x)^4)\frac{1}{1-F(x)}-F(x)+x=0\\ G(y)=(y^2+y^4)\frac{1}{1-y}-y+x\\ G'(y)=(2y+4y^3)\frac{1}{1-y}+(y^2+y^4)(-\frac{1}{(1-y)^2})(-1)-1\\ G(F(x))\equiv 0\pmod {x^n}\\ G(F(x))\equiv G(F_0(x))+G'(F_0(x))(F(x)-F_0(x)) \pmod {x^n}\\ F(x)\equiv F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))}\pmod {x^n}\\ \]
　　時間復(fù)雜度：\(O(n\log n)\)

代碼

#include<cstdio> #include<cstring> #include<algorithm> #include<cstdlib> #include<ctime> #include<functional> #include<cmath> #include<vector> #include<assert.h> //using namespace std; using std::min; using std::max; using std::swap; using std::sort; using std::reverse; using std::random_shuffle; using std::lower_bound; using std::upper_bound; using std::unique; using std::vector; typedef long long ll; typedef unsigned long long ull; typedef double db; typedef std::pair<int,int> pii; typedef std::pair<ll,ll> pll; void open(const char *s){ #ifndef ONLINE_JUDGEchar str[100];sprintf(str,"%s.in",s);freopen(str,"r",stdin);sprintf(str,"%s.out",s);freopen(str,"w",stdout); #endif } void open2(const char *s){ #ifdef DEBUGchar str[100];sprintf(str,"%s.in",s);freopen(str,"r",stdin);sprintf(str,"%s.out",s);freopen(str,"w",stdout); #endif } int rd(){int s=0,c,b=0;while(((c=getchar())<'0'||c>'9')&&c!='-');if(c=='-'){c=getchar();b=1;}do{s=s*10+c-'0';}while((c=getchar())>='0'&&c<='9');return b?-s:s;} void put(int x){if(!x){putchar('0');return;}static int c[20];int t=0;while(x){c[++t]=x%10;x/=10;}while(t)putchar(c[t--]+'0');} int upmin(int &a,int b){if(b<a){a=b;return 1;}return 0;} int upmax(int &a,int b){if(b>a){a=b;return 1;}return 0;} const int N=270000; ll p,g; ll fp(ll a,ll b) {ll s=1;for(;b;b>>=1,a=a*a%p)if(b&1)s=s*a%p;return s; } namespace ntt {const int W=262144;int rev[N];int *w[20];void init(){ll s=fp(g,(p-1)/W);w[18]=new int[1<<17];w[18][0]=1;for(int i=1;i<W/2;i++)w[18][i]=w[18][i-1]*s%p;for(int i=17;i>=1;i--){w[i]=new int[1<<(i-1)];for(int j=0;j<1<<(i-1);j++)w[i][j]=w[i+1][j<<1];}}void ntt(ll *a,int n,int t){for(int i=1;i<n;i++){rev[i]=(rev[i>>1]>>1)|(i&1?n>>1:0);if(rev[i]>i)swap(a[i],a[rev[i]]);}for(int i=2,l=1;i<=n;i<<=1,l++)for(int j=0;j<n;j+=i)for(int k=0;k<i/2;k++){ll u=a[j+k];ll v=a[j+k+i/2]*w[l][k];a[j+k]=(u+v)%p;a[j+k+i/2]=(u-v)%p;}if(t==-1){reverse(a+1,a+n);ll inv=fp(n,p-2);for(int i=0;i<n;i++)a[i]=a[i]*inv%p;}}void add(ll *a,ll *b,ll *c,int n,int m,int l){static ll a1[N];int k=max(n,m);for(int i=0;i<=k;i++)a1[i]=0;for(int i=0;i<=n;i++)a1[i]=(a1[i]+a[i])%p;for(int i=0;i<=m;i++)a1[i]=(a1[i]+b[i])%p;for(int i=0;i<=l;i++)c[i]=a1[i];}void mul(ll *a,ll *b,ll *c,int n,int m,int l){static ll a1[N],a2[N];int k=1;while(k<=n+m)k<<=1;for(int i=0;i<k;i++)a1[i]=a2[i]=0;for(int i=0;i<=n;i++)a1[i]=a[i];for(int i=0;i<=m;i++)a2[i]=b[i];ntt(a1,k,1);ntt(a2,k,1);for(int i=0;i<k;i++)a1[i]=a1[i]*a2[i]%p;ntt(a1,k,-1);for(int i=0;i<=l;i++)c[i]=(a1[i]+p)%p;}void getinv(ll *a,ll *b,int n){if(n==1){b[0]=fp(a[0],p-2);return;}getinv(a,b,n>>1);static ll a1[N],a2[N];for(int i=0;i<n<<1;i++)a1[i]=a2[i]=0;for(int i=0;i<n;i++)a1[i]=a[i];for(int i=0;i<n>>1;i++)a2[i]=b[i];ntt(a1,n<<1,1);ntt(a2,n<<1,1);for(int i=0;i<n<<1;i++)a1[i]=a2[i]*(2-a1[i]*a2[i]%p)%p;ntt(a1,n<<1,-1);for(int i=0;i<n;i++)b[i]=(a1[i]+p)%p;} } int e[50],t; int check(int x) {for(int i=1;i<=t;i++)if(fp(x,(p-1)/e[i])==1)return 0;return 1; } void getg() {int _=p-1;for(int i=2;i*i<=_;i++)if(_%i==0){e[++t]=i;while(_%i==0)_/=i;}if(_!=1)e[++t]=_;for(int i=1;;i++)if(check(i)){g=i;return;} } void solve(ll *a,int n) {if(n==1)return;solve(a,n>>1);static ll f1[N],f2[N],f3[N],f4[N],f5[N],a1[N],a2[N],a3[N],a4[N];for(int i=0;i<n;i++)a1[i]=(p-a[i])%p;a1[0]++;for(int i=0;i<n;i++)f1[i]=a[i];ntt::ntt(f1,n<<1,1);for(int i=0;i<n<<1;i++)f2[i]=f1[i]*f1[i]%p;ntt::ntt(f2,n<<1,-1);for(int i=0;i<n;i++)f4[i]=f2[i];ntt::ntt(f4,n<<1,1);for(int i=0;i<n<<1;i++)f3[i]=f4[i]*f1[i]%p;ntt::ntt(f3,n<<1,-1);for(int i=0;i<n<<1;i++)f4[i]=f4[i]*f4[i]%p;ntt::ntt(f4,n<<1,-1);for(int i=0;i<n;i++)f5[i]=f4[i];ntt::ntt(f5,n<<1,1);for(int i=0;i<n<<1;i++)f5[i]=f5[i]*f1[i]%p;ntt::ntt(f5,n<<1,-1);for(int i=0;i<n;i++)f1[i]=a[i];for(int i=0;i<n;i++)a1[i]=(-f1[i]+3*f2[i]-2*f3[i]+f4[i]-f5[i])%p;for(int i=1;i<n;i++)a1[i]=(a1[i]-2*f1[i-1]+f2[i-1])%p;a1[1]++;for(int i=0;i<n;i++)a2[i]=(4*f1[i]-2*f2[i]+4*f3[i]-3*f4[i])%p;a2[0]--;ntt::getinv(a2,a3,n);ntt::mul(a1,a3,a4,n-1,n-1,n-1);for(int i=0;i<n;i++)a[i]=(a[i]-a4[i]+p)%p; } int n; ll s[N]; int main() {open("b");scanf("%d%lld",&n,&p);getg();ntt::init();int k=1;while(k<=n)k<<=1;solve(s,k);for(int i=1;i<=n;i++){s[i]=(s[i]+p)%p;printf("%lld\n",s[i]);}return 0; }

轉(zhuǎn)載于:https://www.cnblogs.com/ywwyww/p/10193748.html

總結(jié)

以上是生活随笔為你收集整理的【XSY3344】连续段 DP 牛顿迭代 NTT的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。