洛谷P5110 块速递推题解

q779

已于 2022-06-17 10:34:31 修改

阅读量106

点赞数

分类专栏： OI 文章标签：算法数学

于 2022-05-13 10:27:33 首次发布

本文链接：https://blog.csdn.net/qq_50332374/article/details/124746821

版权

OI 专栏收录该内容

240 篇文章 4 订阅

订阅专栏

洛谷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$ ，据出题人说是凑好的

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

也就是，如果数列 $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)