題面欺詐系列...

因為一個點最多只能連到前k個點，所以只有當前的連續k個點的連通情況是對接下來的求解有用的

那麽就可以計算k個點的所有連通情況，dfs以下發現k=5的時候有52種。

我們把它們用類似於並查集的方式表達（比如12132代表點1和點3連通，2和5連通，3自己），然後再壓縮一下。

但要註意的是，12132和23213這兩種實際對應的是一種連通情況，我們只要把它都化成字典序最小的那種就可以了

然後考慮增加一個新點以後狀態的轉移，可以枚舉這個點與前k個點（始狀態S）的連邊情況，其中有一些是不合法的：

　　1.連到了兩個本來就連通的點上（導致成環）

　　2.在1號點不與其他點連通的情況下，沒有連到1號點（導致不連通）

然後再根據連邊情況得到終狀態E，將trans[S][E]++。最後trans[i][j]表示的就是加入一個點後由i狀態到j狀態的方案數

那麽就可以得到遞推式$f[i][j]=\sum^{總狀態數}_{k=1}{f[i-1][k]*trans[k][j]}$，其中f[i][j]表示以i為結尾的k個點狀態是j的方案數

那麽答案就是f[N][1]（假設1是都連通的狀態）

然後初值就應該是f[K][i]=i狀態的方案數，其中i狀態的方案數為它其中每個聯通塊方案數的乘積。

那麽聯通塊的方案數怎麽算呢？其實題面已經說了...n個點的方案數就是$n^{n-2)}$

然後就可以愉快地dp啦...

誒？$N<=10^{15}$？

大概能卡過去吧（大霧）

容易發現遞推的形式其實和矩陣乘法是相同的，把f[i]和trans看作矩陣，就是$f[N]=f[K]*trans^{N-K}$

然後就可以倍增做矩陣快速冪了。方法和整數的快速冪是一樣的。

復雜度：

　　K=5的狀態數接近50，$log(10^{15})$接近50。

　　所以基本上是$50^3*50=6250000$的。

代碼寫的很蛋疼...

  1 #include<cstdio>
  2 #include<cstring>
  3 #include<algorithm>
  4 #include<map>
  5 #include<vector>
  6
 
 #define LL long long 
  7 #define inf 0x3f3f3f3f
  8 using namespace std;
  9 const LL maxn=1e15+5,maxs=55,p65=50000,mod=65521;
 10 
 11 LL rd(){
 12     LL x=0;char c=getchar();int neg=1;
 13     while(c<‘0‘||c>‘9‘){if(c==‘-‘) neg=-1;c=getchar();}
 14     while(c>=‘0‘&&c<=‘9‘) x=x*10+c-‘0‘,c=getchar();
 15     return x*neg;
 16 }
 17 
 18 LL N;int K,sct,num[6]={0,1,1,3,16,125};
 19 LL f[maxs];
 20 LL trans[2][maxs][maxs],out[2][maxs][maxs];
 21 struct ST{
 22     int s[6];
 23     ST(int x=0){memset(s,x,sizeof(s));}
 24 }sta[maxs];
 25 struct Mat{
 26     LL m[maxs][maxs];
 27     Mat(){memset(m,0,sizeof(m));}
 28 };
 29 int mp[p65];
 30 
 31 void print(ST s){
 32     for(int i=1;i<=K;i++) printf("%d",s.s[i]);
 33 }
 34 int STtoInt(ST s){
 35     int x=0;for(int i=1;i<=K;i++) x=x*6+s.s[i];return x;
 36 }
 37 
 38 ST toLeast(ST s){
 39     ST flag=ST(-1),re=ST();int n=0;
 40     for(int i=1;i<=K;i++){
 41         if(flag.s[s.s[i]]==-1) flag.s[s.s[i]]=n++;
 42         re.s[i]=flag.s[s.s[i]];
 43     }return re;
 44 }
 45 
 46 void dfs1(int x,ST s,int m){
 47     if(x>=K+1){
 48         sta[++sct]=s;mp[STtoInt(s)]=sct;return;
 49     }
 50     for(int i=0;i<=m+1;i++){
 51         s.s[x]=i;
 52         dfs1(x+1,s,max(m,i));
 53     }
 54 }
 55 
 56 void dfs2(int x,ST s,ST from){
 57     if(x>=K+1){
 58         bool flag[6],bb=s.s[1];ST to=ST();int mi=5;
 59         memset(flag,0,sizeof(flag));
 60         //print(from);printf("\t");print(s);printf("\n");
 61         for(int i=1;i<=K;i++){
 62             if(!s.s[i]) continue;
 63             if(flag[from.s[i]]) return;
 64             flag[from.s[i]]=1;mi=min(mi,from.s[i]);
 65         }
 66         for(int i=2;i<=K;i++){
 67             to.s[i-1]=flag[from.s[i]]?mi:from.s[i];
 68             if(from.s[1]==to.s[i-1]) bb=1;
 69         }to.s[K]=mi;
 70         if(!bb) return;
 71         to=toLeast(to);
 72         //print(from);printf("\t");print(s);printf("\t");print(to);printf("\n");
 73         trans[0][mp[STtoInt(from)]][mp[STtoInt(to)]]++;
 74         
 75     }else{
 76         dfs2(x+1,s,from);s.s[x]=1;dfs2(x+1,s,from);
 77     }
 78 }
 79 
 80 void sets(){
 81     for(int i=1;i<=sct;i++){
 82         dfs2(1,ST(),sta[i]);
 83         //for(int j=1;j<=sct;j++) printf("%d ",trans[0][i][j]);printf("\n");
 84         ST cnt=ST();
 85         for(int j=1;j<=K;j++) cnt.s[sta[i].s[j]]++;
 86         f[i]=1;for(int j=0;j<K;j++) if(num[cnt.s[j]])f[i]*=num[cnt.s[j]];
 87     }
 88 }
 89 
 90 void solve(){
 91     LL p=N-K;bool b=1,c=1;int i,j,k;
 92     for(i=1;i<=sct;i++) out[0][i][i]=1;
 93     while(p){
 94         if(p&1){
 95             memset(out[c],0,sizeof(out[c]));
 96             for(i=1;i<=sct;i++){
 97                 for(j=1;j<=sct;j++){
 98                     for(k=1;k<=sct;k++){
 99                         out[c][i][j]=(out[c][i][j]+out[c^1][i][k]*trans[b^1][k][j])%mod;
100                     }
101                 }
102             }c^=1;
103         }
104         memset(trans[b],0,sizeof(trans[b]));
105         for(i=1;i<=sct;i++){
106             for(j=1;j<=sct;j++){
107                 for(k=1;k<=sct;k++){
108                     trans[b][i][j]=(trans[b][i][j]+trans[b^1][i][k]*trans[b^1][k][j])%mod;
109                 }
110             }
111         }p>>=1;b^=1;
112     }LL ans=0;
113     for(i=1;i<=sct;i++){
114         ans=(ans+f[i]*out[c^1][i][1])%mod;
115     }printf("%d\n",ans);
116 }
117 
118 int main(){
119     int i,j,k;
120     K=rd(),N=rd();
121     dfs1(2,ST(),0);sets();
122     solve();
123     return 0;
124 }

bzoj1494 生成樹計數 (dp+矩陣快速冪)