「Luogu P5110」块速递推

发表于 |
字数统计: 867 | 阅读时长 ≈ 4

给定一个数列$a$，满足如下递推式：

$a_n=233a_{n-1}+666a_{n-2}$。

求出这个数列第$n$项$a_n\bmod 10^9+7$的值，一共有$T$组数据，为了减少输出量，请输出这些答案的异或和。

Data Range：$T\leq 5\times 10^7,0\leq n<2^{64}$

11

Solution

首先求$a_n=233a_{n-1}+666a_{n-2}$的通项式。

考虑这个数列的特征方程，为$x^2=233x+666$，于是解出$x_1=\cfrac{233+\sqrt{56953}}{2},x_2=\cfrac{233-\sqrt{56953}}{2}$。

所以可以设这个数列的通项式为$a_n=p(\cfrac{233+\sqrt{56953}}{2})^n+q(\cfrac{233-\sqrt{56953}}{2})^n$。

可别忘了，$a_0=0,a_1=1$，所以解出$p=\cfrac{1}{\sqrt{56953}},q=-\cfrac{1}{\sqrt{56953}}$。

所以$a_n=\cfrac{1}{\sqrt{56953}}[(\cfrac{233+\sqrt{56953}}{2})^n-(\cfrac{233-\sqrt{56953}}{2})^n]$

现在来解决除号以及根号。首先，枚举得$\sqrt{56953}\equiv 188305837\pmod {10^9+7}$。

所以说

所以，最终的通项式就是

。

但是，由于毒瘤出题人的毒瘤给分制度，前面六个点每点$1 pts$，所以用龟速幂只能水$6 pts$。我尝试用$O(1)$的按位龟速乘优化了一下后，还是只有$6 pts$。看来我们只能优化计算上面的两个幂的过程。

很明显可以先优化一下指数，由欧拉定理得：

$a_{n+\varphi(10^9+7)}\equiv a_n\pmod {10^9+7}$

于是可以把指数模$\varphi(10^9+7)=10^9+6$，但是还是只能得$6 pts$。

于是考虑$100 pts$的分块打表大法。由上面的结论可知，可以将指数控制在$10^9+6$以内，所以可以考虑分块打表。

可以知道，$\lceil\sqrt{10^9+6}\rceil=31623$，所以以这个数作为块长分块，将它记为$k$。

首先，预处理出

$94153035^0,94153035^1\cdots 94153035^{k-1}$

和$905847205^0,905847205^1\cdots 905847205^{k-1}$

然后，用龟速幂可算出

$94153035^k\equiv 37348318\pmod {10^9+7}$，$905847205^k\equiv 482464047\pmod {10^9+7}$

最后，预处理出

$94153035^k,94153035^2k\cdots 94153035^{k(k-1)}$

和
$905847205^k,905847205^k\cdots 905847205^{k(k-1)}$

显然，这些步骤的时间复杂度为$O(k)$的。

把询问$m$拆成$ak+b$的形式，其中$0\leq b\leq k-1$。
用数学化的语言来讲，$a=\lfloor \frac{m}{k}\rfloor,b=m\bmod k$，所以可以拆成整块和散块来做到$O(1)$查询啦qwq。

时间复杂度为$O(T+\sqrt{10^9+6})$。

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
#include<bits/stdc++.h>
using namespace std;
typedef unsigned long long int ll;
const ll MAXN=31651,MOD=1e9+7,blockSize=31623;
ll test,num,res,resx,resy;
ll blockx[MAXN],powx[MAXN],blocky[MAXN],powy[MAXN];
inline ll read()
{
    register ll num=0,neg=1;
    register char ch=getchar();
    while(!isdigit(ch)&&ch!='-')
    {
        ch=getchar();
    }
    if(ch=='-')
    {
        neg=-1;
        ch=getchar();
    }
    while(isdigit(ch))
    {
        num=(num<<3)+(num<<1)+(ch-'0');
        ch=getchar();
    }
    return num*neg;
}
inline void setup()
{
    blockx[0]=powx[0]=blocky[0]=powy[0]=1;
    for(register int i=1;i<blockSize;i++)
    {
        blockx[i]=blockx[i-1]*37348318ull%MOD;
        blocky[i]=blocky[i-1]*482464047ull%MOD;
        powx[i]=powx[i-1]*94153035ll%MOD;
        powy[i]=powy[i-1]*905847205ll%MOD;
    }
}
namespace maker{
    ll sa,sb,sc;
    inline void setup()
    {
        sa=read(),sb=read(),sc=read();
    }
    inline ll randInt()
    {
        sa^=sa<<32,sa^=sa>>13,sa^=sa<<1;
        ll t=sa;
        sa=sb,sb=sc,sc^=t^sa;
        return sc;
    }
}
int main()
{
    test=read();
    maker::setup(),setup();
    for(register int i=0;i<test;i++)
    {
        num=maker::randInt()%(MOD-1);
        resx=blockx[num/blockSize]*powx[num%blockSize]%MOD;
        resy=blocky[num/blockSize]*powy[num%blockSize]%MOD;
        res^=(233230706*((resx-resy+MOD)%MOD)%MOD);
    }
    printf("%llu",res);
}