[JLOI2015] 有意义的字符串 题解

看到这个 \(7\times 10^{18}\) 的模数已经可以摆烂了。

不是，你告诉我这东西跟矩阵快速幂优化 DP 有关系？？

观察到这个题显然不能硬做，因为你显然不能直接算小数部分，而且还有个取模很难受。

所以我们希望把一切的计算都基于整数。

这个时候我们就要思考，有什么东西可以把根号转化为整数。

直接平方再开方显然是不行的，因为你没办法开方取模。

于是，很神奇的东西出现了，你发现这个东西很像求根公式。

即方程 \(x^2-b\times x+\frac{b^2-d}{4}=0\) 的两个根分别为 \(\frac{b+\sqrt d}{2}\), \(\frac{b-\sqrt d}{2}\)。

这个东西显然可以拿来降幂，即 \(x^n=b\times x^{n-1}+\frac{d-b^2}{4}\times x^{n-2}\)。

所以你就可以直接通过矩阵快速幂求得 \(x^n=a\times x^0+b\times x^1\) 的 \(a,b\)。但是有个问题是 \(x^1\) 并不是整数，而前面还有个取了模的系数 \(b\) 完全没办法向下取整。

发现一个有趣的事实是，对于这个方程的两个根，他们的和为 \(b\)。

故我们考虑对于 \(f_n=(\frac{b+\sqrt d}{2})^n+(\frac{b-\sqrt d}{2})^n\) 进行降幂，即 \(f_i=b\times f_{i-1}+\frac{d-b^2}{4}\times f_{i-2}\)，然后有 \(f_0=2,f_1=b\)。可以轻松求得 \(f_n\)。

另外一个有趣的事实是这个奇怪的数据范围 \(b^2\le d \le (b+1)^2\)。可以发现，当 \(n\) 为偶数，有 \(0\le (\frac{b-\sqrt d}{2})^n \le 1\)，否则为奇数时，有 \(-1\le (\frac{b-\sqrt d}{2})^n \le 0\)。用人话说就是当 \(n\) 为偶数并且 \(b^2\not=d\) 时，答案就是 \(f_n-1\)，否则就是 \(f_n\)。

#include <bits/stdc++.h>
using namespace std;
#define ll long long
const ll mod = 7528443412579576937ll;
ll add(ll x, ll y){return (1ull * x + 1ull * y) % mod;}
ll mul(ll x, ll y)
{
	ll sum = 0;
	while(y)
	{
		if(y & 1) sum = add(sum, x);
		y >>= 1, x = add(x, x);
	}
	return sum;
}
struct matrix
{
	vector<vector<ll> > a;
	matrix operator *(const matrix &c)const
	{
		matrix now;
		int n = a.size(), m = c.a.size(), l = c.a[0].size();
		now.a.resize(n);
		for (int i = 0; i < n; ++i) now.a[i].resize(l); 
		for (int i = 0; i < n; ++i)
		{
			for (int j = 0; j < l; ++j)
			{
				for (int k = 0; k < m; ++k)
				now.a[i][j] = add(now.a[i][j], mul(a[i][k], c.a[k][j]));
			}
		}
		return now;
	}
	matrix(){}
	matrix(int len1, int len2, int op = 1)
	{
		a.resize(len1);
		for (int i = 0; i < len1; ++i) a[i].resize(len2);
		if(op) for (int i = 0; i < len1; ++i) a[i][i] = 1;
	}
};
ll b, d, n;
matrix ksm(matrix x, ll y)
{
	matrix sum(2, 2);
	while(y)
	{
		if(y & 1) sum = sum * x;
		x = x * x, y >>= 1;
	}
	return sum;
}
int main()
{
	cin >> b >> d >> n;
	if(!n)
	{
		puts("1");
		return 0; 
	}
	matrix now(2, 2), a(1, 2);
	now.a[0][0] = b, now.a[1][0] = (d - b * b) / 4, now.a[0][1] = 1, now.a[1][1] = 0;
	a.a[0][0] = b, a.a[0][1] = 2;
	now = ksm(now, n - 1);
	a = a * now;
	if(b * b != d && n % 2 == 0) a.a[0][0]--;
	cout << a.a[0][0] << endl;
}