P5110 块速递推-光速幂、斐波那契数列通项
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}
\]
然而我用矩阵光速幂水过去了。
之后学了上面的东西之后可能会试着推一下。
代码
#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;
}