洛谷P5110 块速递推 题解

题目链接：P5110 块速递推

题意：给定一个数列 $a$ 满足递推式

$$a_0=0,a_1=1 \\a_n = 233a_{n-1}+666a_{n-2}$$

求 $a_n \bmod (10^9+7)$

多组询问

这个题是有个循环节的，正好是 $10^9+6$ ，据出题人说是凑好的

关于循环节怎么求的我也不太清楚，先留个坑，研究好了就补上来

upd.20220721 现在知道怎么求循环节了，详见本文最后

也就是，如果数列 $a_n$ 在模 $M$ 意义下存在循环节 $p$ ，则有 $a_n \equiv a_{n\,\bmod\, p} \bmod M$

本题的解法就是手推通项公式

具体方法如下

对于二阶线性递推数列

$$a_0=A,a_1=B \\a_n = pa_{n-1}+qa_{n-2} ,n\ge 2$$

考虑使用特征方程求解

$a_n=pa_{n-1}+qa_{n-2}$ 的特征方程为

$$x^2=px+q$$

可以求出两个特解（不一定是实数） $x_1,x_2$ ，则

$$a_n=\alpha x_1^{n} + \beta x_2^{n}$$

注意，如果数列从 $a_1$ 开始，这里就是 $a_n=\alpha x_1^{n-1} + \beta x_2^{n-1}$

然后将 $a_0=A,a_1=B$ 代入可得

$$\begin{cases} \alpha + \beta = A \\\alpha x_1+\beta x_2=B \end{cases}$$

解出 $\alpha,\beta$ 即可

本题的通项公式为

$$a_n=\dfrac{1}{\sqrt{56953}}\left(\left(\dfrac{233+\sqrt{56953}}{2}\right)^n-\left(\dfrac{233-\sqrt{56953}}{2}\right)^n\right)$$

关于特征方程法的更多内容，详见 数列通项的求解

注意到这里有个 $\sqrt{56953}$ ，而我们要求它模意义下的值

也就是求出所有的 $x$

$$x^2\equiv 56953 \bmod(10^9+7)$$

考虑二次剩余求解

什么？不会二次剩余？ 那么就打个暴力好了，最多几秒钟就跑出来了

#include <bits/stdc++.h>
using namespace std;
#define int long long
signed main()
{
    for(int i=1; i<=1000000000; i++)
        if(i*i%1000000007==56953)cout << i << endl;
    return 0;
}

然后有两个解 $188305837,811694170$ 取个小点的代入就好了

则有

$$a_n \equiv 233230706 \times(94153035^n-905847205^n)$$

注意到询问有 $10^7$ 个，快速幂过不了

考虑光速幂，$O(\sqrt{n})$ 预处理 $f(x)=x^{65536t},g(x)=x^t$

询问直接查询 $f(x/65536)\times g(n\%65536)$ 即可

这里可以用位运算加速

代码：

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define uint unsigned long long
#define INF 0x3f3f3f3f3f3f3f3f
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 N (int)(7e4+15)
const int p=(int)(1e9+7);
int Q;
uint ans;
uint pw1[2][N],pw2[2][N];
const int a=94153035;
const int b=905847205;
void init()
{
    pw1[0][0]=pw2[0][0]=1;
    for(int i=1; i<1<<16; i++)
        pw1[0][i]=pw1[0][i-1]*a%p;
    pw2[0][1]=pw1[0][(1<<16)-1]*a%p;
    for(int i=2; i<1<<16; i++)
        pw2[0][i]=pw2[0][i-1]*pw2[0][1]%p;

    pw1[1][0]=pw2[1][0]=1;
    for(int i=1; i<1<<16; i++)
        pw1[1][i]=pw1[1][i-1]*b%p;
    pw2[1][1]=pw1[1][(1<<16)-1]*b%p;
    for(int i=2; i<=1<<16; i++)
        pw2[1][i]=pw2[1][i-1]*pw2[1][1]%p;
}
int pow1(int n)
{
    return pw1[0][n&65535]%p*pw2[0][n>>16]%p;
}
int pow2(int n)
{
    return pw1[1][n&65535]%p*pw2[1][n>>16]%p;
}
uint solve(int n)
{
    return 233230706*(pow1(n)-pow2(n)+p)%p;
}
signed main()
{
    init();
    scanf("%lld",&Q);
    Mker::init();
    while(Q--)
        ans^=solve(Mker::rand()%(p-1));
    printf("%llu\n",ans);
    return 0;
}
//233230706×(94153035^n−905847205^n) 

如何求解循环节？

如果数列 $a_n$ 在模 $M$ 意义下存在循环节 $p$ ，则有 $a_n \equiv a_{n\,\bmod\, p} \bmod M$

问了 @zx2017 老师 CCOrz

这个东西其实求起来非常简单，

在已知 $a_n$ 递推式和 $M$ 的情况下

直接暴力去枚举就好了

具体的，这道题

$$a_0=0,a_1=1 \\a_n = 233a_{n-1}+666a_{n-2}$$

注意到 $a_n$ 是从 $a_{n-1}$ 和 $a_{n-2}$ 推来的

那我们枚举 $n$ ，直到出现了这样的情况

$$\begin{cases} a_i\equiv a_2 \\\\a_{i-1}\equiv a_1 \\\\a_{i-2}\equiv a_0 \end{cases} \mod M$$

那么 $i-2$ 就是循环节

是不是很神奇？我太菜了根本没想到 QAQ

当然可能会有更优秀的做法（比起暴力枚举）

不过起码知道怎么简单求就很好了 qwq