[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)modp$。

$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$ 的矩阵，可以列出初始矩阵：

$$\left[\begin{array}{ccccc}{a}_k& a_{k-1}& \cdots & a_2& a_1\end{array}\right]$$

其下一个转移则为：

$$\left[\begin{array}{ccccc}{a}_{k+1}& a_k& \cdots & a_3& a_2\end{array}\right]$$

根据递推式可以列出转移矩阵：

$$\left[\begin{array}{ccccc}{c}_1& 1& 0& \cdots & 0\\ c_2& 0& 1& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & 0\\ c_n& 0& 0& \cdots & 1\end{array}\right]$$

这样我们就可以在 ${\displaystyle \mathcal{O}(K^{2}logN)}$ 的时间里递推出 $a_n$ 的值。但是我们回顾题意要求求出的值：

$$G(n,m)=\left(\sum _{i=m}^n{a}_i\right)modp$$

我们可以设 ${\displaystyle Sum(n)=\sum _{i=1}^n{a}_i}$ ，可以发现：

$$G(n,m)=Sum(n)-Sum(m-1)$$

所以我们可以在矩阵加速的时候一起处理出来 $Sum(n)$，令我们的初始矩阵扩充为：

$$\left[\begin{array}{cccccc}{a}_k& a_{k-1}& \cdots & a_2& a_1& Sum(k-1)\end{array}\right]$$

其下一个转移则为：

$$\left[\begin{array}{cccccc}{a}_{k+1}& a_k& \cdots & a_3& a_2& Sum(k)\end{array}\right]$$

考虑到 $Sum(n)=Sum(n-1)+a_n$，可以得到扩充后的转移矩阵为：

$$\left[\begin{array}{cccccc}{c}_1& 1& 0& \cdots & 0& 1\\ c_2& 0& 1& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ c_{n-1}& 0& 0& \cdots & 0& 0\\ c_n& 0& 0& \cdots & 1& 1\end{array}\right]$$

这样我们就可以在 ${\displaystyle \mathcal{O}(K^{2}logN)}$ 的时间里解决这道题。

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;
}