題面

一個由自然數組成的數列按下式定義：

對于 \(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;}