题面

一个由自然数组成的数列按下式定义：

对于 \(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#includetypedef long long valueType;typedef std::vector ValueVector;valueType MOD_;valueType const &MOD = MOD_;class Matrix {public:    typedef long long valueType;    typedef valueType &reference;    typedef size_t sizeType;    typedef std::vector Row;    typedef std::vector 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;}