BZOJ-4002. [JLOI2015]有意义的字符串

给定正整数 $b, d, n$，计算

\[ \left \lfloor \left (\frac{b + \sqrt d}{2}\right )^n \right \rfloor \bmod 7528443412579576937 \]

其中 $0 < b^2 \leq d < (b + 1)^2 \leq 10^{18}, n \leq 10^{18}$，并且 $b \bmod 2 = 1, d \bmod 4 = 1$

题目连接：BZOJ-4002（这里的题面不完全，我上面打出来的是完整的）

这题的样例中 $b = 1, d = 5, n = 9$，带入之后会发现要求的是 $ \left (\frac{1 + \sqrt 5}{2}\right )^9$

如果对 Fibonacci 数列比较熟悉的话就会发现它正是其通项公式中的一部分

而且，根据特征根法，二阶常系数线性递推数列的通项基本都是长成下面这样子

\[c_1 (A + \sqrt B)^n + c_2 (A - \sqrt B)^n\]

然后再看题目要求的式子，可以猜想它是某个这种线性递推数列通项的一部分，如果真的是这样，那么另外一部分就是 $\left(\frac{b - \sqrt d}{2}\right)^n$， 根据题目中给出的条件可以知道 $\frac{b - \sqrt d}{2} \in (-0.5, 0]$，也就是如果我们找到一个数列 $a_n = \left (\frac{b + \sqrt d}{2}\right )^n + \left (\frac{b - \sqrt d}{2}\right )^n $ 的递推式，就可以利用矩阵快速幂来求出第 $n$ 项，如果 $a_n$ 每一项都是整数，那么答案与 $a_n$ 的误差不会超过 $1$，现在来求这个数列的递推式

设 ${a_n}$ 的特征方程是 $Ax^2 + Bx + C = 0$

根据其通项可以知道

\[\begin{eqnarray*} \left\{ \begin{aligned} B^2-4AC&=d\\ -B&=b \\ 2A&=2 \end{aligned} \right. \end{eqnarray*}\]

最后可以得到 $A=1, B = -b, C = \frac{b^2 - d}{4}$，因此这个数列是

\[a_n = ba_{n-1}-\frac{b^2-d}{4}a_{n-2}, a_0 = 2, a_1 = b\]

再根据题目的条件可以知道 $\frac{b^2-d}{4} \in \mathbb Z$，这样 $a_n \in \mathbb Z$，所以刚刚所有的要求都满足了

现在来考虑什么时候会产生误差，因为

\[\left (\frac{b + \sqrt d}{2}\right )^n = a_n - \left (\frac{b - \sqrt d}{2}\right )^n\]

并且 $b - \sqrt d < 0$，所以当 $n$ 是偶数的时候，答案是 $a_n - 1$，否则答案是 $a_n$。

/* BZOJ-4002: [JLOI2015]有意义的字符串
 *   线性递推+矩阵快速幂 */
#include <cstdio>
#include <algorithm>

typedef long long int64;
typedef unsigned long long uint64;
const int64 mod_v = 7528443412579576937ll;

int64 mod_add(int64 a, int64 b)
{
	uint64 z = uint64(a) + uint64(b);
	if(z >= uint64(mod_v)) return z - mod_v;
	return z;
}

int64 mod_mul(int64 a, int64 b)
{
	if(a < b) std::swap(a, b);
	int64 x = 0;
	while(b)
	{
		if(b & 1) x = mod_add(x, a);
		a = mod_add(a, a);
		b >>= 1;
	}
	return x;
}

struct matrix_t
{
	int64 m[2][2];
	matrix_t operator * (const matrix_t& x) const
	{
		matrix_t z;
		for(int i = 0; i != 2; ++i)
		{
			for(int j = 0; j != 2; ++j)
			{
				int64 v = 0;
				for(int k = 0; k != 2; ++k)
					v = mod_add(v, mod_mul(m[i][k], x.m[k][j]));
				z.m[i][j] = v;
			}
		}

		return z;
	}
};

matrix_t matrix_pow(matrix_t a, int64 p)
{
	matrix_t v;
	v.m[0][0] = v.m[1][1] = 1;
	v.m[1][0] = v.m[0][1] = 0;
	while(p)
	{
		if(p & 1) v = a * v;
		a = a * a;
		p >>= 1;
	}
	return v;
}

int main()
{
	int64 b, d, n;
	std::scanf("%lld %lld %lld", &b, &d, &n);
	if(n == 0) 
	{
		std::puts("1");
		return 0;
	}

	matrix_t m;
	m.m[0][0] = b;
	m.m[0][1] = mod_v - ((b * b - d) >> 2);
	m.m[1][0] = 1;
	m.m[1][1] = 0;
	m = matrix_pow(m, n - 1);
	int64 ans = mod_add(mod_mul(m.m[0][0], b), mod_mul(m.m[0][1], 2));
	if(n % 2 == 0) ans -= 1;
	std::printf("%lld\n", ans);
	return 0;
}

（你们可以试试看把模数 $7528443412579576937$ 直接解释成字符串输出看看会发生什么）

#include <cstdio>

char str[10];
int main()
{
	*reinterpret_cast<long long*>(str) = 7528443412579576937ll;
	std::puts(str);
	return 0;
}