CF392C & Yet Another Number Sequence 题解

纪念一下完全凭借自己想出来的的第一道黑题。

思路

首先看到数据范围 $n\le 10^{17}$ 就知道本题是矩阵快速幂或者推式子。

首先发现对于不同的 $i,j$ $f_i\times i^k,f_j\times j^k$ 非常难合并，除非用出斐波那契通项公式，所以放弃这条道路，走矩阵快速幂道路。

考虑一个最简单的初始矩阵：

$$\left[\begin{array}{cc}{f}_i\times i^k& f_{i+1}\times (i+1{)}^k\end{array}\right]$$

我们要转移到的目标矩阵：

$$\left[\begin{array}{cc}{f}_{i+1}\times {(i+1)}^k& f_{i+2}\times (i+2{)}^k\end{array}\right]$$

我们发现这是不可行的。因为 $(i+2{)}^k$ 不是 $(i+1{)}^k$ 且增广矩阵中的量必须为常数的缘故，转移到第一个数十分简单，但是第二个数就很困难了。

于是我们尝试着倒退，将 $f_{i+2}\times (i+2{)}^k$ 拆成与 $f_i\times i^k{f}_{i+1}\times (i+1{)}^k$ 有关的数，并尝试将其加入初始矩阵与增广矩阵中。

$$f_{i+2}\times (i+2{)}^k$$

$$=f_i\times (i+2{)}^k+f_{i+1}\times (i+2{)}^k$$

由二项式定理：

$$(i+2{)}^k=((i+1)+1{)}^k=\sum _{j=0}^k\left(\genfrac{}{}{0ex}{}{k}{j}\right)(i+1{)}^k$$

$$(i+2{)}^k=\sum _{j=0}^k(\genfrac{}{}{0ex}{}{k}{j})i^j\times 2^{j-k}$$

将上式子带回原式：

$$=f_i\times \sum _{j=0}^k\left(\genfrac{}{}{0ex}{}{k}{j}\right)i^j\times 2^{j-k}+f_{i+1}\times \sum _{j=0}^k\left(\genfrac{}{}{0ex}{}{k}{j}\right)(i+1{)}^k$$

将外面的乘进去：

$$=\sum _{j=0}^k{f}_i\times i^j\left(\genfrac{}{}{0ex}{}{k}{j}\right)2^{j-k}+\sum _{j=0}^k{f}_{i+1}\times (i+1{)}^j(\genfrac{}{}{0ex}{}{k}{j})$$

你会惊奇的发现：加和中的 $\left(\genfrac{}{}{0ex}{}{k}{j}\right)$ 以及 $\left(\genfrac{}{}{0ex}{}{k}{j}\right)2^{j-k}$ 均为可以提前预处理的常数。结合 $k$ 又非常小，可以将 $p\in [0,k] $ $f_i\times i^p$ 全部预处理出来，塞进初始矩阵里就可以了。

矩阵大小为 $40\times 40$，一次矩阵乘法复杂度为 $40^3$，所以总时间为 $6400\times {\log }_2(10^{17})<57\ast 6400$ 可过

还有，题目中要求求的是和，你还得在初始矩阵中额外添加一个和的数据，

代码

#include <bits/stdc++.h>
using namespace std;

#define MAXN 40
#define int long long 

const int mod = 1e9 + 7;

struct matrix {
	int n, m;
	int a[(MAXN<<1)+5][(MAXN<<1)+5];
	matrix() { memset(a, 0, sizeof a); }
} ans, p;

void debug(matrix x) {
	cout << x.n << ' ' << x.m << ":\n";
	for (int i = 1; i <= x.n; ++i) {
		for (int j = 1; j <= x.m; ++j) cout << x.a[i][j] <<' ';
		cout << endl;
	}
}

matrix operator * (matrix x, matrix y) {
	matrix z;
	z.n = x.n, z.m = x.m;
	for (int i = 1; i <= x.n; ++i) 
		for (int j = 1; j <= y.m; ++j)
			for (int k = 1; k <= x.m; ++k)
				(z.a[i][j] += x.a[i][k] * y.a[k][j] % mod) %= mod;
	return z;
}

matrix mpow(matrix x, int y) {
	matrix ans;
	ans.n = x.n, ans.m = x.m;
	for (int i = 1; i <= min(ans.n, ans.m); ++i) ans.a[i][i] = 1; 
	while (y) {
		if (y & 1) ans = ans * x;
		x = x * x;
		y >>= 1;
	}
	return ans;
}

int mpow(int x, int y) {
	if (y < 0) return 0;
	int s = 1;
	while (y) {
		if (y & 1) s = s * x % mod;
		x = x * x % mod;
		y >>= 1;
	}
	return s;
}

int C(int a, int b) {
	if (a < b) return 0;
	if (b > a / 2) b = a - b;
	int s = 1;
	for (int i = a - b + 1; i <= a; ++i) s = s * i % mod;
	for (int i = 1; i <= b; ++i) s = s * mpow(i, mod - 2) % mod;
	return s;
}

int n, k;

signed main() {
	cin >> n >> k;
	if (n == 1) {
		puts("1");
		return 0;
	} 
	ans.n = 1, ans.m = 2 * k + 3;
	for (int i = 0; i <= k; ++i) 
		ans.a[1][i + 1] = 1, ans.a[1][i + k + 2] = 2 * mpow(2, i) % mod;
	ans.a[1][2 * k + 3] = (1 + mpow(2, k + 1)) % mod;
	p.n = p.m = 2 * k + 3;
	for (int j = 1; j <= 2 * k + 3; ++j) {
		if (j <= k + 1) 
			p.a[j + k + 1][j] = 1;
		else if (j <= 2 * k + 2) {
		 	// F[i] * i ^ p
		 	// F[i + 1] * (i + 1) ^ p
		 	// ↓
			// F[i + 1] * (i + 2) ^ p + F[i] * (i + 2) ^ p 
			for (int i = k + 2; i <= 2 * k + 2; ++i) // F[i + 1]
				p.a[i][j] = C(j - k - 2, i - k - 2);
			for (int i = 1; i <= k + 1; ++i) 
				p.a[i][j] = mpow(2, j - k - 1 - i) * C(j - k - 2, i - 1) % mod;
		}
		else {
			for (int i = k + 2; i <= 2 * k + 2; ++i) // F[i + 1]
				p.a[i][j] = C(k, i - k - 2);
			for (int i = 1; i <= k + 1; ++i) 
				p.a[i][j] = mpow(2, k + 1 - i) * C(k, i - 1) % mod;
			p.a[2 * k + 3][j] = 1;
		}
	}
//	debug(ans), debug(p);
	ans = ans * mpow(p, n - 2);
//	debug(ans);
	cout << ans.a[1][2 * k + 3] << endl;
	return 0;
}