又是一道常系数线性递推的题，一下子便能想到矩阵快速幂，但一看数据规模
\[1≤T≤5×10^7 \]我们可以知道，普通的矩阵快速幂肯定是过不了的
显然，本题只能接受\(O(T)\)时间复杂度的方法(LYT:?)
考虑光速幂
对于\(x^y=(x^b)^{\lfloor\frac{y}{b}\rfloor}x^{y\%b}\),,我们预处理出\(x^n(n<b)\),再预处理出\((x^b)^n\)，令\(b=a^k(k\in (0,1))\),则预处理时空复杂度为\(O(b+\frac{y}{b})\),单次查询\(O(1)\)。
更一般的,\(x^y=x^{a_0}(x^b)^{a_1}(x^{b^2})^{a_2}\cdots\Rightarrow y=a_0+a_1b+a_2b^2+a_3b^3\cdots\),就是将\(y\)进行\(b\)进制分解,对每一位进行预处理，令\(b=a^k(k\in (0,1))\) 预处理空间复杂度为\(O(b=a^k)\),查询时间复杂度为\(O(\frac{1}{k})\).
1.矩阵光速幂
但我们可以发现,因为\(a_i< 2^{64}\),所以\(k=0.25\)，对于\(5*10^7\),哪怕处理了光速幂，\(3\)次矩阵乘法(约\(30\)的常数)也是难以通过的
2.特征多项式？
我不会，这里推荐一个大佬的博客：【学习笔记】特征多项式
但\(3\)次多项式乘法(约\(15\)的常数)依然难以通过
3.使用生成函数解出该式通解
先推导一些生成函数对应的函数
\[\frac{1}{1-ax}=g\\ 1=g-axg\\ axg+1=g\\ g_0=1,g_{n+1}=ag_n\\ \] \[\frac{1}{a-x}=g\\ 1=ag-gx\\ g=\frac{1+gx}{a}\\ g_0=\frac{1}{a},g_{n+1}=\frac{g_n}{a} \]进入正题
\[a_n=233a_{n-1}+666a_{n-2},a_0=0,a_1=1\\ 我们设函数g为序列a的生成函数\\ 则有:\\ \begin{aligned} &233xg+666x^2g+x=g\\ \Leftrightarrow &g=\frac{x}{1-233x-666x^2}\\ \Leftrightarrow &g=-x\left(\frac{1}{666x^2+233x-1}\right)\\ \Leftrightarrow &-\frac{g}{x}=\frac{1}{666(x-\frac{-233+\sqrt{56953}}{1332})(x-\frac{-233-\sqrt{56953}}{1332})}\\ \Leftrightarrow &-\frac{666g}{x}=\frac{1}{(x-\frac{-233+\sqrt{56953}}{1332})(x-\frac{-233-\sqrt{56953}}{1332})}\\ &-\frac{666g}{x}=\frac{\frac{1}{x-\frac{-233+\sqrt{56953}}{1332}}-\frac{1}{x-\frac{-233-\sqrt{56953}}{1332}}}{\frac{\sqrt{56953}}{666}}\\ &\frac{g}{x}=\frac{\frac{1}{\frac{-233+\sqrt{56953}}{1332}-x}-\frac{1}{\frac{-233-\sqrt{56953}}{1332}-x}}{\sqrt{56953}}\\ &a_n=\frac{\left(\frac{233+\sqrt{56953}}{2}\right)^n+\left(\frac{233-\sqrt{56953}}{2}\right)^n}{\sqrt{56953}} \end{aligned} \]这里补充一下最后一步的推导
\[令a=\frac{-233+\sqrt{56953}}{1332}\\ \Rightarrow \frac{1}{a}=\frac{233+\sqrt{56953}}{2}\\ \frac{1}{\frac{-233+\sqrt{56953}}{1332}-x}\\ =\sum_{i=0}^{+\infty}\frac{x^i}{a^{i+1}} \]到了这一步，我们就可以通过\(2\)种方式来解决
Ⅰ.扩域
可以看出，上式中出现了无理数，但在模意义下无理数是无意义的，因为模运算是定义在\(L=Z\cap [0,mod=1000000007)\)上的。(如无特殊说明，下文中四则运算皆为模意义下)
我们不妨定义域\(K=(\{a+b\sqrt{56953}|a,b\in L\},+,*,-)\)
证明\(K\)为\(L\)的扩域
- 封闭性
-
存在\(+,*\)
显然
之后就可以像复数运算一样操作了
#include<bits/stdc++.h>
using namespace std;
# define ll long long
# define read read1<ll>()
# define Type template<typename T>
# define fre(k) freopen(k".in","r",stdin);freopen(k".out","w",stdout)
Type T read1(){
T t=0;
char k;
bool vis=0;
do (k=getchar())=='-'&&(vis=1);while('0'>k||k>'9');
while('0'<=k&&k<='9')t=(t<<3)+(t<<1)+(k^'0'),k=getchar();
return vis?-t:t;
}
namespace Mker
{
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;
}
}
# define mod 1000000007
# define M 32767
# define N 32768
# define inv2 500000004ll
struct A{
int x,y;
A& operator *=(const A &b){
int t=y;
y=(1ll*x*b.y+1ll*y*b.x)%mod;
x=(1ll*x*b.x+56953ll*t%mod*b.y)%mod;
return *this;
}
A operator *(const A &b){A t=*this;return t*=b;}
}f[2][N|1],f1[2][N|1],o,p;
int main(){
int T=read,ans=0;
Mker::init();o.x=233*inv2%mod,o.y=inv2;p.x=o.x,p.y=mod-inv2;
f[0][0].x=f1[0][0].x=1;f[1][0].x=f1[1][0].x=1;
for(int i=1;i<=N;++i)f[0][i]=f[0][i-1]*o,f[1][i]=f[1][i-1]*p;
for(int i=1;i<=N;++i)f1[0][i]=f1[0][i-1]*f[0][N],f1[1][i]=f1[1][i-1]*f[1][N];
while(T--){
unsigned ll n=Mker::rand()%(mod-1);
int w=n>>15&M,t=(1ll*f[0][n&M].x*f1[0][w].y+1ll*f[0][n&M].y*f1[0][w].x-1ll*f[1][n&M].x*f1[1][w].y-1ll*f[1][n&M].y*f1[1][w].x)%mod;
if(t<0)t+=mod;ans^=t;
}cout<<ans;
return 0;
}
Ⅱ.二次剩余
我们经过二次剩余可以知道，对于素数\(mod\)，我们可以找到唯一的\(x\)满足\(x^2\equiv 56953 \mod mod\).
经过打表可知，对于本题，\(188305837 \equiv 56953 \mod mod\)
(可能你想知道
Mker::rand()%(mod-1)
是为什么
通过\(扩域\)我们可以知道域\(K\)是满足域\(L\)的运算规则的，自然而然满足费马定理，所以在模意义下\(a_n与a_{n+k(mod-1)}\)相等