题面
一个由自然数组成的数列按下式定义：
对于 \(i \le k\)：\(a_{i}= b_{i}\)。
对于 \(i > k\)：\(\displaystyle a_{i}= \sum_{j=1}^{k}{c_{j} \times a_{i-j}}\)。
其中 \(b_{1\dots k}\) 和 \(c_{1\dots k}\) 是给定的自然数。
写一个程序，给定自然数 \(m \le n\)，计算 \(\left( \sum_{i=m}^{n}{a_{i}} \right) \bmod p\)。
\(1 \le k \le 15\)，\(1 \le m \le n \le 10^{18}\)，\(0 \le b_{i},c_{i} \le 10^{9}\)，\(p \le 10^{8}\)。
题解
因为 \(k\) 很小，\(n, m\) 很大，不难想到矩阵加速递推。
根据 \(\displaystyle a_{i}= \sum_{j=1}^{k}{c_{j} \times a_{i-j}}\)，所以我们的矩阵应该至少是一个 \(1 \times k\) 的矩阵，可以列出初始矩阵：
\[\begin{bmatrix} a_k & a_{k - 1} & \cdots & a_2 & a_1 \end{bmatrix}\]其下一个转移则为：
\[\begin{bmatrix} a_{k + 1} & a_{k} & \cdots & a_3 & a_2 \end{bmatrix}\]根据递推式可以列出转移矩阵：
\[\begin{bmatrix} c_1 & 1 & 0 & \cdots & 0\\ c_2 & 0 & 1 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & 0\\ c_n & 0 & 0 & \cdots & 1 \end{bmatrix}\]这样我们就可以在 \(\displaystyle \mathcal{O}(K^2logN)\) 的时间里递推出 \(a_n\) 的值。但是我们回顾题意要求求出的值：
\[G(n, m) = \left( \sum_{i=m}^{n}{a_{i}} \right) \bmod p \]我们可以设 \(\displaystyle Sum(n) = \sum_{i = 1}^{n}a_i\) ，可以发现：
\[G(n, m) = Sum(n) - Sum(m - 1) \]所以我们可以在矩阵加速的时候一起处理出来 \(Sum(n)\)，令我们的初始矩阵扩充为：
\[\begin{bmatrix} a_k & a_{k - 1} & \cdots & a_2 & a_1 & Sum(k - 1) \end{bmatrix}\]其下一个转移则为：
\[\begin{bmatrix} a_{k + 1} & a_{k} & \cdots & a_3 & a_2 & Sum(k) \end{bmatrix}\]考虑到 \(Sum(n) = Sum(n - 1) + a_n\)，可以得到扩充后的转移矩阵为：
\[\begin{bmatrix} c_1 & 1 & 0 & \cdots & 0 & 1\\ c_2 & 0 & 1 & \cdots & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ c_{n - 1} & 0 & 0 & \cdots & 0 & 0\\ c_n & 0 & 0 & \cdots & 1 & 1 \end{bmatrix}\]这样我们就可以在 \(\displaystyle \mathcal{O}(K^2logN)\) 的时间里解决这道题。
Code
//Luogu - P2461
#include<bits/stdc++.h>
typedef long long valueType;
typedef std::vector<valueType> ValueVector;
valueType MOD_;
valueType const &MOD = MOD_;
class Matrix {
public:
typedef long long valueType;
typedef valueType &reference;
typedef size_t sizeType;
typedef std::vector<valueType> Row;
typedef std::vector<Row> Container;
valueType MOD = ::MOD;
enum TYPE : int {
EMPTY = 0, UNIT = 1
};
protected:
sizeType _row_, _column_;
Container data;
public:
Matrix(sizeType row, sizeType column) : _row_(row), _column_(column), data(_row_) {
for (auto &iter: data)
iter.resize(column, 0);
};
sizeType row() const {
return _row_;
}
sizeType column() const {
return _column_;
}
void set(TYPE type) {
for (auto &iter: data) {
std::fill(iter.begin(), iter.end(), 0);
}
if (type == EMPTY)
return;
if (type == UNIT)
for (sizeType i = 0, end = std::min(_row_, _column_); i < end; ++i)
data[i][i] = 1;
}
reference operator()(sizeType i, sizeType j) {
if (i > this->_row_ || j > this->_column_)
throw std::out_of_range("Too Large.");
if (i == 0 || j == 0)
throw std::out_of_range("0 index access.");
return std::ref(data[i - 1][j - 1]);
}
Matrix operator+(const Matrix &T) const {
if (this->_row_ != T._row_ || this->_column_ != T._column_)
throw std::range_error("Illegal operation.");
Matrix result(this->_row_, this->_column_);
for (sizeType i = 0; i < this->_row_; ++i)
for (sizeType j = 0; j < this->_column_; ++j)
result.data[i][j] = (this->data[i][j] + T.data[i][j]) % MOD;
return result;
}
Matrix operator*(const Matrix &T) const {
if (this->_column_ != T._row_)
throw std::range_error("Illegal operation.");
Matrix result(this->_row_, T._column_);
for (sizeType i = 0; i < this->_row_; ++i) {
for (sizeType k = 0; k < this->_column_; ++k) {
valueType r = this->data[i][k];
for (sizeType j = 0; j < T._column_; ++j)
result.data[i][j] = (result.data[i][j] + T.data[k][j] * r) % MOD;
}
}
return result;
}
Matrix operator^(valueType x) const {
if (x < 1)
throw std::range_error("Illegal operation.");
Matrix result(this->_row_, this->_column_);
Matrix base = *this;
result.set(UNIT);
while (x) {
if (x & 1) result = result * base;
base = base * base;
x = x >> 1;
}
return result;
}
friend std::ostream &operator<<(std::ostream &os, const Matrix &T) {
for (sizeType i = 0; i < T._row_; ++i)
for (sizeType j = 0; j < T._column_; ++j)
os << T.data[i][j] << " \n"[j == T._column_ - 1];
return os;
}
friend std::istream &operator>>(std::istream &os, Matrix &T) {
for (sizeType i = 0; i < T._row_; ++i)
for (sizeType j = 0; j < T._column_; ++j)
os >> T.data[i][j];
return os;
}
};
int main() {
valueType K, M, N;
std::cin >> K;
ValueVector B(K + 30, 0), C(K + 30, 0);
for(int i = 1; i <= K; ++i)
std::cin >> B[i];
for(int i = 1; i <= K; ++i)
std::cin >> C[i];
std::cin >> M >> N >> MOD_;
for(int i = 1; i <= K; ++i) {
B[i] %= MOD;
C[i] %= MOD;
}
Matrix ans(1, K + 1), base(K + 1, K + 1);
ans.set(Matrix::EMPTY);
base.set(Matrix::EMPTY);
for(int i = 1; i <= K; ++i)
base(i, 1) = C[i];
for(int i = 2; i <= K; ++i)
base(i - 1, i) = 1;
base(1, K + 1) = base(K + 1, K + 1) = 1;
for(int i = 1; i <= K; ++i)
ans(1, K + 1 - i) = B[i];
ans(1, K + 1) = std::accumulate(B.begin() + 1, B.begin() + K, 0) % MOD;
valueType resultN = 0, resultM = 0;
++N;
if(N > K) {
Matrix MatrixN = ans * (base ^ (N - K));
resultN = MatrixN(1, K + 1);
} else {
resultN = std::accumulate(B.begin() + 1, B.begin() + N, 0);
}
if(M > K) {
Matrix MatrixM = ans * (base ^ (M - K));
resultM = MatrixM(1, K + 1);
} else {
resultM = std::accumulate(B.begin() + 1, B.begin() + M, 0);
}
valueType result = resultN - resultM;
result = (result % MOD + MOD) % MOD;
std::cout << result << std::flush;
return 0;
}