Luogu 块速递推 | 光速幂模板

Luogu P5110 块速递推

---

• 题意
  给定一个数列 \(a\) 满足递推式：
  \[a_n=233a_{n-1}+666a_{n-2},a_0=0,a_1=1 \]
  求这个数列第 \(n\) 项 \(mod\) \(10^9+7\) 的值，有 \(T\) 组询问。
  \(n\) 在 unsigned long long 的范围内； \(1 \leq T \leq 5*10^7\)。

---

• 题解
  我们先列出这个递推式的特征方程：
  \[233x+666=x^2 \]
  \[x^2-233x-666=0 \]
  \[x_1=\frac{233+\sqrt{56953}}{2};x_2=\frac{233-\sqrt{56953}}{2} \]
  \[a_n=k_1x_1^n+k_2x_2^n \]
  \[\begin{cases} k_1+k_2=0 \\ k_1x_1+k_2x_2=1 \end{cases}\]
  \[\begin{cases} k_1=\frac{1}{\sqrt{56953}} \\ k_1=\frac{1}{-\sqrt{56953}} \end{cases}\]
  \[\therefore a_n=\frac{1}{\sqrt{56953}}((\frac{233+\sqrt{56953}}{2})^n-(\frac{233-\sqrt{56953}}{2})^n) \]
  \[\because 188305837=\sqrt{56953}(mod　10^9+7) \]
  \[\because a_n=233230706*(94153035^n-905847205^n) \]
  然后我们接下来要做的就是 \(O(1)\) 求这个式子（光速幂）
  • 我们可以设 \(x^n=32769k+r\)
    \[f_1(n)=x^{32768n};f_2(n)=x^n \]
    \[\therefore x^p=f_1(\frac{p}{32768})f_2(p　mod　32768) \]

---

• 接下来就是一点也不激动人心的代码了：

#include<bits/stdc++.h>
/* Forever_chen */
#define RT register
#define int long long
#define P 1000000007
using namespace std;
template<class t> inline t read(t &x){
	char c=getchar();bool f=0;x=0;
	while(!isdigit(c)) f|=c=='-',c=getchar();
	while(isdigit(c))x=(x<<1)+(x<<3)+(c^48),c=getchar();
	if(f)x=-x;return x;
}
template<class t>inline void write(t x){
	if(x<0)putchar('-'),write(-x);
	else{if(x>9)write(x/10);putchar('0'+x%10);}
}
template<class t>inline void writeln(t x){
	write(x);putchar('\n');
	return;
}
template<class t>inline void write_blank(t x){
	write(x);putchar(' ');
	return;
}
namespace Mker{
    unsigned long long SA, SB, SC;
    void init() { 
		scanf("%llu%llu%llu", &SA, &SB, &SC); 
	}
    inline 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;
    }
}
in
int x1,x2,x3,x4;
int f1[65550],f2[65550],f3[65550],f4[65550];
inline int Pow_1(int x) { 
	return 1ll*f3[x>>16]*f1[x&65535]%P; 
}
inline int Pow_2(int x) { 
	return 1ll*f4[x>>16]*f2[x&65535]%P; 
}
int ttt,ans;
signed main(){
	//freopen(".in","r",stdin);
	//freopen(".out","w",stdout);
	f1[0]=1;
	f2[0]=1;
	f3[0]=1;
	f4[0]=1;
	x1=94153035;
	x2=905847205; 
	x3=64353223;
	x4=847809841;
    for(RT int i=1;i<65536;i++){
    	f1[i]=1ll*f1[i-1]*x1%P;
	} 
    for(RT int i=1;i<65536;i++){
    	f2[i]=1ll*f2[i-1]*x2%P;
	} 
    for(RT int i=1;i<65536;i++){
    	f3[i]=1ll*f3[i-1]*x3%P;
	}
    for(RT int i=1;i<36636;i++){
    	f4[i]=1ll*f4[i-1]*x4%P;
	} 
	read(ttt);
	Mker::init();
	while(ttt--){
		unsigned int lop;
		lop=Mker::rand()%(P-1);
		ans^=233230706*(Pow_1(lop)-Pow_2(lop)+P)%P; 
	}
	write(ans);
	//system("pause");
	return 0;
}

注：可能需要卡卡常