P5110 块速递推
题意
多次询问，求数列
[a_i=egin{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) 和 (1sqrt n、2sqrt n、cdots、sqrt n sqrt n)
显然就可以 (O(1)) 求这个东西了。
但是！询问的数字大小肯定不是在模域范围内的，所以我们需要找循环节。
有一个问题就是，矩阵的循环节并不固定。
但是有一个结论，对角线元素互不相同的下三角矩阵的循环节为 (largemathbf{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
",ans);
}
}
signed main(){
star::work();
return 0;
}