題意：有一張 $(L+1)\times n$ 個點(diǎn)的有向圖，每個結(jié)點(diǎn)有二元組 $(x,y)(0\le x\le L,1\le y\le n)$ 表示。對于所有 $(u_1,v_1),(u_2,v_2)$，若 $u_1<u_2$,則前者向后者連 $w_{v_1,v_2}$ 條邊。對所有 $0\le t<k$，你需要求出從 $(0,x)$ 走模 $k$ 余 $t$ 步到達(dá)任意一個第二維為 $y$ 的點(diǎn)的方案數(shù)。答案對 $p$ 取模。

$L\le 10^8,n\le 3,k\le 2^{16},10^8<p<2^{30},p\in prime,k\mid p?1$

神奇的題

顯然是單位根反演，所以先想辦法搞出走到位置 $i$ 的生成函數(shù)。

顯然是個矩陣的形式。看成只有 $n$ 個點(diǎn)在上面跑，鄰接矩陣為 $M$

那么走一步的生成函數(shù)為

$$s(x)=\sum _{i=1}^{\infty }M{x}^i$$

答案的生成函數(shù)為

$$f(x)=\sum _{i=0}^{\infty }\left(\sum _{j=0}^L[x^j]s^i(x)\right)x^i$$

看上去很不可做，但實(shí)際上中間那坨看起來形式很簡單。冷靜分析，實(shí)際上 $M$ 的次數(shù)就是 $i$，前面的系數(shù)就是選 $i$ 個正整數(shù)和不超過 $L$,等價于選 $i+1$ 個正整數(shù)不超過 $L+1$，即

$$\begin{gathered} f(x)=\sum _{i=0}^L\left(\genfrac{}{}{0ex}{}{L}{i}\right)M^i{x}^i \\ =(Mx+1{)}^L \end{gathered}$$

然后就可以單位根反演了

$$an{s}_t=\frac{1}{k}\sum _{i=0}^{k?1}{\omega }_k^{?it}(M{\omega }_k^i+1{)}^L$$

我們發(fā)現(xiàn)我們只求這個東西的 $x$ 行 $y$ 列，所以可以設(shè) $a_i={\left[(M{\omega }_k^i+1{)}^L\right]}_{x,y}$

然后

$$an{s}_t=\frac{1}{k}\sum _{i=0}^{k?1}{\omega }_k^{?it}{a}_i$$

這是個 CZT 的形式，拆開就可以了

$$an{s}_t=\frac{1}{k}\sum _{i=0}^{k?1}{\omega }_k^{\left(\genfrac{}{}{0ex}{}{i}{2}\right)+\left(\genfrac{}{}{0ex}{}{t}{2}\right)?\left(\genfrac{}{}{0ex}{}{i+t}{2}\right)}{a}_i$$

$$an{s}_t=\frac{{\omega }_k^{\left(\genfrac{}{}{0ex}{}{t}{2}\right)}}{k}\sum _{i=0}^{k?1}{\omega }_k^{\left(\genfrac{}{}{0ex}{}{i}{2}\right)}{a}_i{\omega }_k^{?\left(\genfrac{}{}{0ex}{}{i+t}{2}\right)}$$

直接做減法卷積就可以了。

求個原根出來求單位元，然后 MTT 即可。

因?yàn)橐獙懙臇|西太多了，很多地方都寫的很浪，僅供參考。

#include <iostream> #include <cstdio> #include <cstring> #include <cctype> #include <cmath> #include <algorithm> #define MAXN (1<<18)+5 #define double long double using namespace std; typedef long long ll; inline int read() {int ans=0;char c=getchar();while (!isdigit(c)) c=getchar();while (isdigit(c)) ans=(ans<<3)+(ans<<1)+(c^48),c=getchar();return ans; } int n,k,L,x,y,M,P,g; inline int qpow(int a,int p) {int ans=1;while (p){if (p&1) ans=(ll)ans*a%P;a=(ll)a*a%P,p>>=1;}return ans; } const double Pi=acos(-1.0); struct complex{double x,y;inline complex(const double& x=0,const double& y=0):x(x),y(y){}}; inline complex operator +(const complex& a,const complex& b){return complex(a.x+b.x,a.y+b.y);} inline complex operator -(const complex& a,const complex& b){return complex(a.x-b.x,a.y-b.y);} inline complex operator *(const complex& a,const complex& b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);} inline complex adj(const complex& a){return complex(a.x,-a.y);} int l,lim,r[MAXN]; inline void init(){lim=1<<l;for (int i=0;i<lim;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));} complex rt[2][24]; void fft(complex* a,int type) {for (int i=0;i<lim;i++) if (i<r[i]) swap(a[i],a[r[i]]);for (int L=0;L<l;L++){int mid=1<<L,len=mid<<1;complex Wn=rt[type][L+1];for (int s=0;s<lim;s+=len){complex w(1,0);for (int k=0;k<mid;k++,w=w*Wn){complex x=a[s+k],y=w*a[s+mid+k];a[s+k]=x+y,a[s+mid+k]=x-y;}}}if (type) for (int i=0;i<lim;i++) a[i].x/=lim,a[i].y/=lim; } complex A[MAXN],B[MAXN],C[MAXN],D[MAXN]; void fmul(int* F,int* G,int n,int m) {M=sqrt(P);for (l=0;(1<<l)<n+m;l++);init();for (int i=0;i<n;i++) A[i]=complex(F[i]/M,F[i]%M);for (int i=0;i<m;i++) C[i]=complex(G[i]/M,G[i]%M);fft(A,0),fft(C,0);B[0]=adj(A[0]),D[0]=adj(C[0]);for (int i=1;i<lim;i++) B[i]=adj(A[lim-i]),D[i]=adj(C[lim-i]);for (int i=0;i<lim;i++){complex a=A[i],b=B[i],c=C[i],d=D[i];A[i]=(a+b)*complex(0.5,0);B[i]=(a-b)*complex(0,-0.5);C[i]=(c+d)*complex(0.5,0);D[i]=(c-d)*complex(0,-0.5);}for (int i=0;i<lim;i++){complex a=A[i],b=B[i],c=C[i],d=D[i];A[i]=a*c+a*d*complex(0,1);B[i]=b*c+b*d*complex(0,1);}fft(A,1),fft(B,1);for (int i=0;i<n+m-1;i++){ll a=A[i].x+0.5,b=A[i].y+B[i].x+0.5,c=B[i].y+0.5;a%=P,b%=P,c%=P;F[i]=(a*M%P*M%P+b*M%P+c)%P;} } struct mat {int e[3][3];inline mat(const int& v=0){memset(e,0,sizeof(e));for (int i=0;i<3;i++) e[i][i]=v; } inline int* operator [](const int& i){return e[i];}inline const int* operator [](const int& i)const{return e[i];} }W; inline mat operator +(const mat& a,const mat& b) {mat c;for (int i=0;i<3;i++)for (int j=0;j<3;j++)c[i][j]=(a[i][j]+b[i][j])%P;return c; } inline mat operator *(const mat& a,const mat& b) {mat c;for (int i=0;i<3;i++)for (int k=0;k<3;k++)for (int j=0;j<3;j++)c[i][j]=(c[i][j]+(ll)a[i][k]*b[k][j])%P;return c; } inline mat qpow(mat a,int p) {mat ans(1);while (p){if (p&1) ans=ans*a;a=a*a,p>>=1;}return ans; } inline bool check(int g) {int x=P-1;for (int i=2;i*i<=x;i++)if (x%i==0){if (qpow(g,(P-1)/i)==1) return false;while (x%i==0) x/=i;}if (x>1&&qpow(g,(P-1)/x)==1) return false; return true; } inline int sum(int x){return (ll)x*(x-1)/2%k;} int F[MAXN],G[MAXN]; int main() {for (int i=0;i<20;i++){double a=2*Pi/(1<<i);rt[1][i]=adj(rt[0][i]=complex(cos(a),sin(a)));}n=read(),k=read(),L=read(),x=read()-1,y=read()-1,P=read();while (!check(++g));for (int i=0;i<n;i++) for (int j=0;j<n;j++) W[i][j]=read();int Wn=qpow(g,(P-1)/k);for (int i=0,w=1;i<k;i++,w=(ll)w*Wn%P) F[i]=(ll)qpow(W*mat(w)+mat(1),L)[x][y]*qpow(Wn,sum(i))%P;for (int i=0;i<2*k;i++)G[i]=qpow(Wn,k-sum(i)); // for (int t=0;t<k;t++) // { // int ans=0; // for (int i=0;i<k;i++) ans=(ans+(ll)F[i]*G[i+t])%P; // printf("%d ",ans); // } // puts("");reverse(F,F+k); fmul(F,G,k,2*k); // for (int i=k-1;i<2*k-1;i++) printf("%d ",F[i]); // puts("");int t=qpow(k,P-2);for (int i=k-1;i<2*k-1;i++) printf("%lld\n",(ll)F[i]*qpow(Wn,sum(i-k+1))%P*t%P);return 0; }

總結(jié)

以上是生活随笔為你收集整理的【HNOI2019】白兔之舞【组合数学】【矩阵快速幂】【单位根反演】【Chirp Z-Transform】【原根】【MTT】的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

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