题解 P5110 【块速递推】

· · 题解

\Large\texttt{My Blog}

Description

题目链接：Luogu 5110

给定一个数列 a 满足递推式

a_n=233\times a_{n-1}+666\times a_{n-2},a_{0}=0,a_{1}=1

求这个数列第 n 项 a_n\bmod (10^9+7) 的值，一共有 T 组询问。为了减少你的输出量，你只需要输出所有询问答案的异或和。

数据范围：0\le n<2^{64}，1\le T\le 5\times 10^7

Solution

先考虑朴素的做法，我们可以直接用矩阵乘法和矩阵快速幂在 O(2^3T\log n) 的时间内求解。矩乘的式子如下：

\begin{aligned}\left[\begin{matrix}a_n\\a_{n-1}\end{matrix}\right]&=\left[\begin{matrix}233\times a_{n-1}+666\times a_{n-2}\\1\times a_{n-1}+0\times a_{n-2}\end{matrix}\right]\\&=\left[\begin{matrix}233&666\\1&0\end{matrix}\right]\times\left[\begin{matrix}a_{n-1}\\a_{n-2}\end{matrix}\right]\\&=\left[\begin{matrix}233&666\\1&0\end{matrix}\right]^{n-1}\times\left[\begin{matrix}a_1\\a_0\end{matrix}\right]\end{aligned}

接下来我们考虑如何优化。由于 T 的范围很大，我们需要保证每次询问的复杂度为 O(1)。首先有一个非常显然的性质 a^b\equiv a^{b\bmod \varphi(p)}\pmod p，因此我们把式子转化为：

\begin{aligned}\left[\begin{matrix}a_n\\a_{n-1}\end{matrix}\right]&=\left[\begin{matrix}233 & 666 \\1 & 0\end{matrix}\right]^{(n-1)\bmod \varphi(10^9+7)}\times\left[\begin{matrix}a_1\\a_0\end{matrix}\right] \\&=\left[\begin{matrix}233 & 666 \\1 & 0\end{matrix}\right]^{(n-1)\bmod (10^9+6)}\times\left[\begin{matrix}a_1\\a_0\end{matrix}\right]\end{aligned}

虽然有这样一个转化，但是我的复杂度瓶颈却在矩阵快速幂的 O(\log n) 上，考虑如何优化掉这个 \log：由于询问很多，我们需要在常数时间内回答，考虑预处理矩阵的幂。

我们令矩阵 A=\left[\begin{matrix}233 & 666\\ 1 & 0\end{matrix}\right]，我们预处理出两部分矩阵的幂。

首先我们令 k=\left\lceil\sqrt{10^9+6}\right\rceil（10^9+6 也就是 n 的最大值）。此处注意一定要上取整，不然预处理的幂可能不充分。

  1. 预处理 A^1,A^2,A^3,\cdots,A^k，将矩阵 A^i 记为 X(i)。
  2. 预处理 A^k,A^{2k},A^{2k},\cdots,A^{k^2}，将矩阵 A^{ik} 记为 Y(i)。

很显然以上两部分的复杂度均为 O(k)。

我们把询问 m（此处的 m 为 (n-1)\bmod \varphi(10^9+7)）拆成 m=ak+b，其中 0\le b<k。转化为数学语言为 a=\left\lfloor\frac{m}{k}\right\rfloor,b=m\bmod k，于是 A^m 转化为 (A^k)^a\times A^b，根据我们上面预处理的部分，我们可以发现，(A^k)^a 和 A^b 都是可以直接查表得到的。所以 A^m=(A^k)^a\times A^b=X(b)\times Y(a)。

这样我们的询问复杂度就降低为 O(1)（其实应该是矩阵乘法复杂度 O(2^3)），加上常数优化即可通过本题！

吐槽：此题卡常！出题人毒瘤！

时间复杂度：O(T+\sqrt{10^9+6})（其中 \sqrt{10^9+6} 约为 32000 左右）

Code

#include <cstdio>
#include <cstring>

const int N=32000,mod=1e9+7;
struct Matrix {
    int mat[2][2];
    Matrix operator * (const Matrix &b) const {
        Matrix ans;
        ans.mat[0][0]=(1LL*mat[0][0]*b.mat[0][0]+1LL*mat[0][1]*b.mat[1][0])%mod;
        ans.mat[0][1]=(1LL*mat[0][0]*b.mat[0][1]+1LL*mat[0][1]*b.mat[1][1])%mod;
        ans.mat[1][0]=(1LL*mat[1][0]*b.mat[0][0]+1LL*mat[1][1]*b.mat[1][0])%mod;
        ans.mat[1][1]=(1LL*mat[1][0]*b.mat[0][1]+1LL*mat[1][1]*b.mat[1][1])%mod;
        return ans;
    }
} p1[N+5],p2[N+5];
unsigned long long SA,SB,SC;

unsigned long long rand() {
    SA^=SA<<32,SA^=SA>>13,SA^=SA<<1;
    unsigned long long t=SA;
    SA=SB,SB=SC,SC^=t^SA;
    return SC;
}
int main() {
    Matrix s;
    s.mat[0][0]=233,s.mat[0][1]=666,s.mat[1][0]=1,s.mat[1][1]=0;
    p1[0].mat[0][0]=p1[0].mat[1][1]=1;
    for(int i=1;i<=N;++i) p1[i]=p1[i-1]*s;
    p2[0].mat[0][0]=p2[0].mat[1][1]=1,s=p1[N];
    for(int i=1;i<=N;++i) p2[i]=p2[i-1]*s;
    int T;
    scanf("%d%llu%llu%llu",&T,&SA,&SB,&SC);
    int ans=0;
    for(int i=1;i<=T;++i) {
        int n=(rand()-1)%(mod-1);
        ans^=(p2[n/N]*p1[n%N]).mat[0][0];
    }
    printf("%d\n",ans);
    return 0;
}