P5110 块速递推
题意
多次询问，求数列 \[
a_i=\begin{cases}233a_{i-1}+666a_{i-2} & i>1\\
0 & i=0\\
1 & i=1\\
\end{cases}
\] 的第 \(n\) 项在 \(\mod 1e9+7\) 意义下的值的异或和。
思路
首先这个数列是一个广义斐波那契数列。对于广义斐波那契数列，我们一般是用矩阵快速幂求的。
但是，这个题的询问次数是 \(5e7\)
。
所以我们就必须用 \(O(1)\)
的方法处理询问。于是，一个自诩光速幂的东西登场了。
实际上，光速幂就是在 \(\sqrt n\)
的时间复杂度内预处理，然后 \(O(1)\)
查询。具体来讲，我们可以预处理出转移矩阵的 \(1、2、\cdots、\sqrt n\) 和 \(1\sqrt n、2\sqrt n、\cdots、\sqrt n \sqrt
n\)
显然就可以 \(O(1)\)
求这个东西了。
但是！询问的数字大小肯定不是在模域范围内的，所以我们需要找循环节。
有一个问题就是，矩阵的循环节并不固定。
但是有一个结论，对角线元素互不相同的下三角矩阵的循环节为
\(\large\mathbf{\varphi_{mod}}\)
。但是笔者并不会证。
所以这题的正解并不是矩阵光速幂QAQ
我们可以用生成函数或者特征方程或者待定系数法来推出通项公式。具体推导过程与斐波那契数列的推导类似，然后用二次剩余将在根号下的项化成模域下的数，然后我们就得出了数列的通项公式：
\[
a_n=233230706(94153035^n−905847205^n)\pmod{10^9}
\] 然而我用矩阵光速幂水过去了。
之后学了上面的东西之后可能会试着推一下。
代码
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<iostream>
#include<cstdio>
#include<cstring>
#include<cctype>
#include<algorithm>
#define int unsigned
using namespace std;
inline int read(){
int w=0,x=0;char c=getchar();
while(!isdigit(c))w|=c=='-',c=getchar();
while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
return w?-x:x;
}
namespace star
{
const int mod=1e9+7,ring=1e9+6,siz=31623;
struct mat{
int a[2][2];
mat(){memset(a,0,sizeof a);}
inline void set(){a[0][0]=a[1][1]=1;}
inline int* operator [] (const int x){return a[x];}
inline const int* operator [] (const int x) const {return a[x];}
inline mat operator * (const mat &b)const{
mat ans;
for(int i=0;i<2;i++)
for(int j=0;j<2;j++)
for(int k=0;k<2;k++)
(ans[i][j]+=1ll*a[i][k]*b[k][j]%mod)>=mod&&(ans[i][j]-=mod);
return ans;
}
}now,pow[siz+1],Pow[siz+1];
unsigned long long SA,SB,SC;
void init(){scanf("%llu%llu%llu",&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;
}
inline void work(){
now[0][1]=0,now[0][0]=1,pow[1][0][0]=233,pow[1][1][0]=666,pow[1][0][1]=1;
pow[0].set();
Pow[0].set();
for(int i=2;i<=siz;i++)
pow[i]=pow[i-1]*pow[1];
Pow[1]=pow[siz];
for(int i=2;i<=siz;i++)
Pow[i]=Pow[i-1]*Pow[1];
int T=read();
init();
unsigned ans=0;
while(T--){
int zp=rand()%ring;
int x=zp/siz,y=zp%siz;
int res;
ans^=(res=(1ll*Pow[x][0][0]*pow[y][0][1]%mod+1ll*Pow[x][0][1]*pow[y][1][1]%mod))>=mod?res-=mod:res;
}
printf("%u\n",ans);
}
}
signed main(){
star::work();
return 0;
}
|
