题解 P3873【[TJOI2010]天气预报】

已知 \(a_{1\cdots n}\) 和 \(w_{1\cdots n}\)，定义 \(w_i=\sum\limits_{j=1}^na_jw_{i-j}\bmod 4147\ (i>n)\)，求 \(w_m\)。

看着感觉很像广义斐波那契那种东西，考虑对递推式进行转化，构造矩阵出来用快速幂加速。

发现求出 \(w_i\) 需要 \(w_{i-1\cdots i-n}\)，这些东西肯定要放到行向量中，那么行向量就是：

\[W_{i-1} = \begin{pmatrix} w_{i-1} & w_{i-2} & \cdots & w_{i-n} \end{pmatrix} \]

那么我们需要构造 \(A\) 使得 \(W_{i-1}A=W_i\)，对着上面递推式和矩阵乘法的运算规则，可以写出下式：

\[\begin{pmatrix} w_{i-1} & w_{i-2} & \cdots & w_{i-n} \end{pmatrix} \begin{pmatrix} a_1 & 1 & 0 & \cdots & 0 & 0 \\ a_2 & 0 & 1 & \cdots & 0 & 0 \\ a_3 & 0 & 0 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ a_{n-1} & 0 & 0 & \cdots & 0 & 1\\ a_n & 0 & 0 & \cdots & 0 & 0 \end{pmatrix} = \begin{pmatrix} w_i & w_{i-1} & \cdots w_{i-n+1} \end{pmatrix} \]

也就是 \(W_{i-1}A=W_i\)。我们已知 \(W_n\)，所以 \(W_m=W_nA^{m-n}\)，套矩阵快速幂即可。

// Problem: P3873 [TJOI2010]天气预报
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P3873
// Memory Limit: 125 MB
// Time Limit: 1500 ms
// 
// Powered by CP Editor (https://cpeditor.org)

//By: OIer rui_er
#include <bits/stdc++.h>
#define rep(x,y,z) for(int x=(y);x<=(z);x++)
#define per(x,y,z) for(int x=(y);x>=(z);x--)
#define debug printf("Running %s on line %d...\n",__FUNCTION__,__LINE__)
#define fileIO(s) do{freopen(s".in","r",stdin);freopen(s".out","w",stdout);}while(false)
using namespace std;
typedef long long ll;
const int N = 105, mod = 4147;

int n, m;
template<typename T> void chkmin(T& x, T y) {if(x > y) x = y;}
template<typename T> void chkmax(T& x, T y) {if(x < y) x = y;}
struct Matrix {
	int n, m, a[N][N];
	Matrix(int x=0, int y=0) : n(x), m(y) {
		memset(a, 0, sizeof(a));
	}
	void e() {
		assert(n == m);
		rep(i, 1, n) a[i][i] = 1;
	}
	friend Matrix operator * (const Matrix& a, const Matrix& b) {
		assert(a.m == b.n);
		Matrix c(a.n, b.m);
		rep(i, 1, a.n) {
			rep(j, 1, b.m) {
				rep(k, 1, a.m) {
					c.a[i][j] += a.a[i][k] * b.a[k][j] % mod;
					c.a[i][j] %= mod;
				}
			}
		}
		return c;
	}
	friend Matrix operator ^ (Matrix a, int k) {
		assert(a.n == a.m);
		Matrix ans(a.n, a.m);
		ans.e();
		for(;k;k>>=1,a=a*a) if(k & 1) ans = ans * a;
		return ans;
	}
};

int main() {
	scanf("%d%d", &n, &m);
	Matrix w(1, n), a(n, n);
	rep(i, 1, n) scanf("%d", &w.a[1][i]);
	rep(i, 1, n) scanf("%d", &a.a[i][1]);
	rep(i, 2, n) a.a[i-1][i] = 1;
	a = a ^ (m - n);
	w = w * a;
	printf("%d\n", w.a[1][1]);
	return 0;
}