@codeforces - 717A@ Festival Organization
@description@
一个长度为 n 的 01 序列是好的，当且仅当该序列任意两个 0 不相邻。
求从 [l, r] 中选出 k 个长度相等的 01 序列的方案数。
1 ≤ k ≤ 200, 1 ≤ l ≤ r ≤ 10^18。
@solution@
定义 f[i] 表示合法 01 串数量，对最后一个是 1 还是 0 进行讨论得到 f[i] = f[i-1] + f[i-2]。
递推式有斐波那契数列形式，但注意 f[1] = 2 所以不完全是斐波那契，不过可以通过移位变成斐波那契。
以下都把它当作斐波那契数列来解决。
则问题相当于问 \(\sum_{i=l}^{r}C_{f_i}^k = \sum_{i=0}^{r}C_{f_i}^k - \sum_{i=0}^{l-1}C_{f_i}^k\)。
可以通过第一类斯特林数化成 \(\sum_{i=0}^{n}C_{f_i}^k = \frac{1}{k!}\sum_{i=0}^{n}(\sum_{j=0}^{k}S_{k}^{j}*f_{i}^j*(-1)^{k-j})\)。
问题可以等价于求 \(\sum_{i=0}^{n}f_{i}^p\)。当 p = 1, 2 的时候都可以搜到相应的式子。是否 p 增大也有相应的式子呢？
除了把递推式写成矩阵以外，斐波那契还有一个常见解决方法：通项公式。众所周知，斐波那契的通项公式可以写作 \(f[i] = a*A^i + b*B^i\)。
注意通项公式中含有 \(\sqrt{5}\)，而该模数下 5 没有二次剩余。所以需要类似复数，将数存成 \((a + b*\sqrt{5})\) 的形式。
然后直接把通项公式往里面代：
等比数列求和解后面的式子即可。
最后复杂度 O(k^2logA)（还有等比数列快速幂的复杂度）。
@accepted code@
#include <cstdio>
#include <algorithm>
using namespace std;
typedef long long ll;
const int MAXK = 200;
const int MOD = int(1E9) + 7;
struct mint{
int x, y;
mint(int _x=0, int _y=0) : x(_x), y(_y) {}
friend mint operator + (mint a, mint b) {
int x = (a.x + b.x >= MOD ? a.x + b.x - MOD : a.x + b.x);
int y = (a.y + b.y >= MOD ? a.y + b.y - MOD : a.y + b.y);
return mint(x, y);
}
friend mint operator - (mint a, mint b) {
int x = (a.x - b.x < 0 ? a.x - b.x + MOD : a.x - b.x);
int y = (a.y - b.y < 0 ? a.y - b.y + MOD : a.y - b.y);
return mint(x, y);
}
friend mint operator * (mint a, mint b) {
int x = (1LL*a.x*b.x + 5LL*a.y*b.y) % MOD;
int y = (1LL*a.x*b.y + 1LL*a.y*b.x) % MOD;
return mint(x, y);
}
friend mint mpow(mint a, ll k) {
mint r = 1;
while( k ) {
if( k & 1 ) r = r * a;
a = a * a;
k >>= 1;
}
return r;
}
friend mint operator / (mint a, mint b) {
mint p = a * mint(b.x, (MOD - b.y)%MOD), q = mint((1LL*b.x*b.x%MOD + MOD - 5LL*b.y*b.y%MOD)%MOD, 0);
return p * mpow(q, MOD-2);
}
}s[MAXK + 5][MAXK + 5], c[MAXK + 5][MAXK + 5];
const mint A1 = 1/mint(0,1);
const mint A2 = (1 + mint(0,1))/2;
const mint B1 = 1/mint(0,MOD-1);
const mint B2 = (1 - mint(0,1))/2;
void init() {
for(int i=0;i<=MAXK;i++) {
c[i][0] = 1;
for(int j=1;j<=i;j++)
c[i][j] = c[i-1][j] + c[i-1][j-1];
}
for(int i=0;i<=MAXK;i++)
for(int j=0;j<=i;j++)
c[i][j] = c[i][j] * mpow(A1, j) * mpow(B1, i-j);
s[0][0] = 1;
for(int i=1;i<=MAXK;i++)
for(int j=1;j<=i;j++)
s[i][j] = s[i-1][j-1] + (i-1)*s[i-1][j];
for(int i=0;i<=MAXK;i++)
for(int j=0;j<=MAXK;j++)
if( (i + j) & 1 ) s[i][j] = 0 - s[i][j];
}
int k;
mint get(mint a, ll m) {
if( a.x == 1 && a.y == 0 ) return m%MOD + 1;
else return (mpow(a, m + 1) - 1) / (a - 1);
}
mint solve(ll m) {
mint ans = 0;
for(int j=0;j<=k;j++) {
mint del = 0;
for(int p=0;p<=j;p++) {
mint tmp = mpow(A2, p) * mpow(B2, j-p);
del = del + c[j][p] * get(tmp, m);
}
ans = ans + s[k][j] * del;
}
return ans;
}
int main() {
ll l, r; init();
scanf("%d%lld%lld", &k, &l, &r), l += 2, r += 2;
mint p = 1; for(int i=1;i<=k;i++) p = p*i;
printf("%d\n", ((solve(r) - solve(l-1))/p).x);
}
@details@
注意等比数列特判公比为 1 ！！！（老是记不住）
貌似是 BJOI2019 的勘破神机？
