考虑DP，我们设g[i]表示第i-m+1天~第i天为第一次出现非法方案的方案数,f[i]为第i天之前合法方案的方案数，h[i]为第i天前非法方案的方案数。
根据定义容易得出：h[i]=h[i-1]*m+g[i]（即第i-1天前的方案数乘以第i天可召唤的种类加上g[i]）
f[i]=mi-h[i]
g[i]=m!*f[i-m]-$\sum _{k=1}^{m-1}k!\ast g\left[ i-k\right]$=m!*mi-m-m!*h[i-m]-$\sum _{k=1}^{m-1}k!\ast g\left[ i-k\right]$
对于g[i]，在m!*f[i-m]个方案中会存在一些不符合g[i]定义的方案，即不是第i-m+1天~第i天才开始出现第一次的非法方案。
如：对于序列.....5 4 1 2 3 4 5......（红色部分为g[i]的范围）
我们发现这个实际上在g[i-2]就出现非法方案了，不是g[i]的非法方案，我们需要减去这部分的方案。
我们发现，在上面这个序列，第i-2-5天到第i-2天是一个排列，第i-5天到第i天也是一个排列。
所以在第i-2-5天到第i-2天的某一个排列里，在第i-2天到第i天必有2!个方案多算的。
所以进一步整理我们可以得出g[i]=m!*f-[i-m]-1!*g[i-1]-2!*g[i-2]-…-(m-1)!*g[i-m+1],也就是上面那个例子。
所以我们得到了递推式：
g[i]=m!*mi-m-m!*h[i-m]-$\sum _limits{k=1}^{m-1}k!\ast g\left[ i-k\right]$
h[i]=h[i-1]*m+g[i]
不过这会超时，但是因为这是递推式，我们可以采用矩阵优化。
这是一个关于i的递推式，我们考虑列出矩阵出来。
| mi-m |
| g[i-m+1] |
| g[i-m+2] |
| …… |
| g[i-1] |
| hi-m |
考虑它要转移到i+1的矩阵
| mi-m+1 |
| g[i-m+2] |
| g[i-m+3] |
| …… |
| g[i] |
| hi-m+1 |
对照以上递推式我们得出转移的系数：
| m | * | mi-m | |||||
| 1 | g[i-m+1] | ||||||
| …… | g[i-m+2] | ||||||
| 1 | …… | ||||||
| m! | -(m-1)! | -(m-2)! | …… | -1! | -m! | g[i-1] | |
| 1 | m | hi-m |
（例：由mi-m到mi-m+1是乘以一个m
因为hi-m+1=hi-m*m+g[i-m+1],所以g[i-m+1]的系数为1，h[i-m]的系数为m，其余同理）
所以用个矩阵快速幂后就可以可以得出答案为mn-hn,即幂运算后的系数矩阵第一列的第一个数减去最后一次数即可。
1 #include <cstring>神奇的代码
2 #include <cstdio>
3 #include <iostream>
4 #include <algorithm>
5 #define mo 1000000007
6 #define N 109
7 using namespace std;
8 long long ans,n,m;
9 struct data{
10 long long ju[N][N];
11 data operator * (const data &a) const{
12 data b;
13 for (long long i=1;i<=m+1;++i)
14 for (long long j=1;j<=m+1;++j){
15 b.ju[i][j]=0;
16 for (long long k=1;k<=m+1;++k)
17 b.ju[i][j]=(b.ju[i][j]+ju[i][k]*a.ju[k][j]%mo+mo)%mo;
18 }
19 return b;
20 }
21 }qwq,tmp;
22 long long jie(long long x){
23 long long s=1;
24 for (long long i=2;i<=x;i++)
25 s=s*i%mo;
26 return s;
27 }
28 void kuai(){
29 long long b=n;
30 for (long long i=1;i<=m+1;++i) tmp.ju[i][i]=1;
31 while (b){
32 if (b&1) tmp=tmp*qwq;
33 qwq=qwq*qwq;
34 b>>=1;
35 }
36 }
37 void init(){
38 scanf("%lld%lld",&n,&m);
39 qwq.ju[1][1]=m;
40 for (long long i=2;i<=m-1;++i)
41 qwq.ju[i][i+1]=1;
42 for (long long i=2;i<=m;++i)
43 qwq.ju[m][i]=-jie(m-i+1);
44 qwq.ju[m][1]=jie(m);
45 qwq.ju[m][m+1]=-qwq.ju[m][1];
46 qwq.ju[m+1][2]=1;
47 qwq.ju[m+1][m+1]=m;
48 }
49 void solve(){
50 ans=tmp.ju[1][1]-tmp.ju[m+1][1];
51 printf("%lld\n",(ans+mo)%mo);
52 }
53 int main(){
54 init();
55 kuai();
56 solve();
57 return 0;
58 }
（人生第一次写矩阵快速幂QAQ）