我这里和别人不一样的是取模加和乘（我觉得这个东西很好，可能会快一点？
Solution
首先，它让我们算方案数，所以我们设一个数组 (F(x)) 表示完美数恰好为 (x) 的排列数，又由这熟悉的恰好，可以再设一个数组 (G(x)) 表示完美数至少为 (x) 的排列数，那么显然可以得到： (G(m)=sumlimits_{i=m}^ninom imF(i)) ，来一个正常的二项式反演就可以得到 (F(m)=sumlimits_{i=m}^n(-1)^{i-m}inom imG(i)) ，现在考虑如何求出 (G(x)) 。
因为每一位的完美情况只和前后有关系，所以我们可以设 (f_{i,j,0/1,0/1}) 来表示当前位是 (i) ，有 (j) 个完美数， (i) 是/否被选， (i+1) 是/否被选的方案数。
现在考虑状态转移方程：
选择第 (i) 位为完美位，并让 (i-1) 在第 (i) 位，只需要保证 (i-1) 一定没被选。
[f_{i,j,0,0}+=f_{i-1,j-1,0,0},f_{i,j,1,0}+=f_{i-1,j-1,0,1}
]
选择第 (i) 位为完美位，并让 (i+1) 在第 (i) 位，只需要保证 (i+1) 一定没被选。
[f_{i,j,0,1}+=f_{i-1,j-1,0,0}+f_{i-1,j-1,1,0},f_{i,j,1,1}+=f_{i-1,j-1,0,1}+f_{i-1,j-1,1,1}
]
不选择第 (i) 位为完美位
[f_{i,j,0,0}+=f_{i-1,j,0,0}+f_{i-1,j,1,0},f_{i,j,1,0}+=f_{i-1,j,0,1}+f_{i-1,j,1,1}
]
初始状态： (f_{1,0,0,0}=1,f_{1,1,0,1}=1)
因为第 (n) 位只能选择 (n-1) 为完美数，所以需要重新转移或者减去 (n+1) 的部分。我选择的重新转移。
现在就能得到 (G(m)=(f_{n,m,0,0}+f_{n,m,1,0}) imes(n-m)!) ，其中 ((n-m)!) 是剩下 (n-m) 个位置随便排的全排列。
最后反演即可。
代码
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define ll long long
using namespace std;
const int N=2050,mod=1e9+7;
int n,m;
ll f[N][N][2][2],G[N],fac[N],inv[N];
inline int read(){
int x=0,f=1;
char ch=getchar();
while(!isdigit(ch)){if(ch=='-') f=-1;ch=getchar();}
while(isdigit(ch)){x=x*10+(ch^48);ch=getchar();}
return x*f;
}
ll mul(ll a,ll b){return a*b%mod;}
ll add(ll a,ll b){return a+b>=mod?a+b-mod:a+b;}
inline void init(){
fac[0]=inv[0]=inv[1]=1;
for(int i=2;i<=n;i++)
inv[i]=mul(inv[mod%i],mod-mod/i);
for(int i=1;i<=n;i++){
fac[i]=mul(fac[i-1],i);
inv[i]=mul(inv[i],inv[i-1]);
}
}
ll C(int n,int m){
if(n<m) return 0;
return mul(fac[n],mul(inv[n-m],inv[m]));
}
int main(){
n=read();m=read();
init();
f[1][0][0][0]=f[1][1][0][1]=1;
for(int i=2;i<=n;i++){
f[i][0][0][0]=1;
for(int j=1;j<=i;j++){
f[i][j][0][0]=add(f[i-1][j-1][0][0],add(f[i-1][j][0][0],f[i-1][j][1][0]));
f[i][j][0][1]=add(f[i-1][j-1][0][0],f[i-1][j-1][1][0]);
f[i][j][1][0]=add(f[i-1][j-1][0][1],add(f[i-1][j][1][1],f[i-1][j][0][1]));
f[i][j][1][1]=add(f[i-1][j-1][0][1],f[i-1][j-1][1][1]);
}
}
f[n][0][0][0]=1;
for(int i=1;i<=n;i++){
f[n][i][0][0]=add(f[n-1][i-1][0][0],add(f[n-1][i][0][0],f[n-1][i][1][0]));
f[n][i][1][0]=add(f[n-1][i-1][0][1],add(f[n-1][i][0][1],f[n-1][i][1][1]));
}
//上面是重新转移的n
for(int i=0;i<=n;i++)
G[i]=mul(add(f[n][i][0][0],f[n][i][1][0]),fac[n-i]);
ll ans=0;
for(int i=m;i<=n;i++){
ll res=mul(C(i,m),G[i]);
ans=add(ans,(i-m)&1?mod-res:res);
}
printf("%lld
",ans);
return 0;
}