[SDOI2008] 遞迴數列

題面

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

對於 \(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\) 的矩陣，可以列出初始矩陣：

其下一個轉移則為：

根據遞推式可以列出轉移矩陣：

這樣我們就可以在 \(\displaystyle \mathcal{O}(K^2logN)\) 的時間裡遞推出 \(a_n\) 的值。但是我們回顧題意要求求出的值：

我們可以設 \(\displaystyle Sum(n) = \sum_{i = 1}^{n}a_i\) ，可以發現：

所以我們可以在矩陣加速的時候一起處理出來 \(Sum(n)\)，令我們的初始矩陣擴充為：

其下一個轉移則為：

考慮到 \(Sum(n) = Sum(n - 1) + a_n\)，可以得到擴充後的轉移矩陣為：

這樣我們就可以在 \(\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;
}