[洛谷P5110] 块速递推

前言

已经吐了，感觉在用计算机做数学题

题目

洛谷

讲解

由于推导过程太小了，所以我适当调整了一下字体(于是就变得这么丑了QAQ)

part1 求通项公式

待定系数法求通项公式！(这是唯一一个好懂一点的方法了)

我们的条件有：

\(\Large{a_n=233a_{n-1}+666a_{n-2}\pmod {10^9+7},a_0=0,a_1=1}\)

设\(\Large{a_{n+2}+pa_{n+1}=q(a_{n+1}+pa_n)}\)

显然，我们可以得到：

\[\Large{\begin{cases} q-p=233\\pq=666\end{cases}} \]

\[\Large{p=q-233,q(q-233)=666} \]

即 \(\Large{q^2-233q-666=0}\)

解得 \(\Large{q=\dfrac{233±\sqrt{56953}}{2}}\)

通过循环暴力查找，我们发现 \(\Large{188305837^2}≡56953\pmod {10^9+7}\)

然后我们就可以用 \(\Large{188305837}\) 替换 \(\Large{\sqrt{56953}}\) 了

我们令 \(\Large{q_1=\dfrac{233+\sqrt{56953}}{2}}=94153035\)

则 \(\Large{p_1=q_1-233}=94152802\)

令 \(\Large{b_n=a_{n+1}+p_1a_n}\)，显然 \(\Large{b_{n+1}=q_1b_n}\)，\(\Large{b}\) 序列为等比数列

其中 \(\Large{b_1=233+p_1=q_1}\)，所以 \(\Large{b_n=q^n}\)

\(\Large{a_{n+1}+p_1a_n=q_1^n ① }\)

我们令 \(\Large{q_2=\dfrac{233-\sqrt{56953}}{2}}=905847205\)

则 \(\Large{p_2=q_2-233=905846972}\) ，同理得到：

\(\Large{a_{n+1}+p_2a_n=q_2^n ② }\)

\(\Large{①-②}\) 得

\[\Large{\begin{aligned} a_n&=\dfrac{q_1^n-q_2^n}{p_1-p_2}\\ &=\dfrac{94153035^n-905847205^n}{188305837}\\ &=233230706(94153035^n-905847205^n)\\ \end{aligned}}\]

part2 优化

通过 通项公式+费马小定理/打表 可知循环节为 \(\Large{10^9+6}\)

对于每个给出的 \(\Large{n}\) ，我们可以对其取模，使其在 \(\Large{[0,10^9+6)}\)

part3 光速幂

我们发现如果是直接通过上述的通项公式计算答案，复杂度仍然带 \(\log_2\)，所以我们要使用一些奇技淫巧

显然 \(\Large{a^n=(a^{{n/b}})^b*a^{n\%b}},n\in[0,10^9+6)\)

定义函数 \(\Large{f(x)=(a^{x/b})^b=(a^b)^{x/b},g(x)=a^{x\%b}}\) (这里的\为整除)

于是分块预处理，令 \(\Large{b=\sqrt{10^9+6}}\approx31623\) 即可 (用2的幂不香吗？)

\(\Large{a^n=f(n)g(n)}\)

时间 \(\Large{O(1)}\)

代码

一份没有特别去卡常的代码

const int MAXN = 31623 + 5;
const int MOD = 1e9 + 7;
int T,n,ans;
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;
	}
}
using namespace Mker;

int qpow(int x,int y)
{
	int ret = 1;
	while(y){if(y & 1) ret = 1ll * ret * x % MOD;x = 1ll * x * x % MOD;y >>= 1;}
	return ret;
}

int b = 31623;
int f[2][MAXN],g[2][MAXN],num[2] = {94153035,905847205};
void pre()
{
	f[0][0] = g[0][0] = f[1][0] = g[1][0] = 1;
	for(int i = 0;i <= 1;++ i)
		f[i][1] = qpow(num[i],b),g[i][1] = num[i];
	for(int i = 2;i < b;++ i)
	{
		for(int j = 0;j <= 1;++ j)
		{
			f[j][i] = 1ll * f[j][i-1] * f[j][1] % MOD;
			g[j][i] = 1ll * g[j][i-1] * g[j][1] % MOD;
		}
	}
}
int update(int x)
{
	if(x < 0) x += MOD;
	return x;
}

int main()
{
//	freopen(".in","r",stdin);
//	freopen(".out","w",stdout);
	T = Read();
	init();
	pre();
	while(T--)
	{
		n = Rand() % (MOD-1);
		if(!n) continue;
		int A = n/b,B = n%b;
		ans ^= update((1ll * f[0][A] * g[0][B] - 1ll * f[1][A] * g[1][B]) % MOD * 233230706 % MOD);
	}
	Put(ans,'\n');
	return 0;
}