@codeforces - 717A@ Festival Organization

@description@
@solution@
@accepted code@
@details@

@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})\) 的形式。

然后直接把通项公式往里面代：

\[\sum_{i=0}^{n}f_{i}^p = \sum_{i=0}^{n}(a*A^i + b*B^i)^p
\\ = \sum_{i=0}^{n}\sum_{j=0}^{p}C_p^j*(a*A^i)^j*(b*B^j)^{p-j}
\\ = \sum_{j=0}^{p}C_p^j*a^j*b^{p-j}\sum_{i=0}^{n}(A^j*B^{p-j})^i\]

等比数列求和解后面的式子即可。

最后复杂度 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 的勘破神机？

巴特西

@codeforces - 717A@ Festival Organization

@description@

@solution@

@accepted code@

@details@

最新文章

热门文章