【BZOJ】4002: [JLOI2015]有意义的字符串

题意

求$\left \lfloor \left( \frac{b+\sqrt{d}}{2} \right)^n \right \rfloor \pmod {7528443412579576937} \(，\)\left( 0 \le n \le 10^{18}, 0 < b^2 \le d < (b+1)^2 \le 10^{18}, b \mbox{ mod } 2 = 1, d \mbox{ mod } 4=1 \right) $

分析

发现这个并不好算，而如果是\(\left( \frac{b-\sqrt{d}}{2} \right)^n\)那么就好算了。于是又想到数列的特征方程得到的解\(a_n = c_1 x_1^n + c_2 x_2^n\)，于是我们搞搞。直接将\(c_1 = c_2 = 1\)，则变成\(a_n = x_1^n + x_2^n\)，而我们知道\(x_1、x_2\)是特征方程的两个解，和上面那个形式极为相似，于是我们继续假设。即\(x_1 = \left( \frac{b+\sqrt{d}}{2} \right), x_2 = \left( \frac{b-\sqrt{d}}{2} \right)\)。则\(a_{n+2} = pa_{n+1} + qa_{n}\)中，\(p = x_1 + x_2, q = - x_1 x_2\)，因此得到\(a_{n+2} = ba_{n+1} - \frac{b^2-d}{4}a_{n}\)

题解

根据上面这个递推式，我们容易算出其中两项，容易得到\(a_1 = b, a_2 = \frac{b^2+d}{2}\)。而发现通项求出来的是整数，因此我们用矩阵乘法求出\(a_n\)即可。最后再根据条件特判一下\(\left( \frac{b-\sqrt{d}}{2} \right)^n\)即可，即\(ans = a_n - [b^2 \neq d \land n是偶数]\)
注意n=0要特判...

#include <bits/stdc++.h>
using namespace std;
typedef unsigned long long ll;
typedef ll mtx[2][2];
const ll mo=7528443412579576937ull, Lim=1e9;
inline void CK(ll &c) {
	if(c>=mo)
		c-=mo;
}
inline ll mul(ll a, ll b) {
	if(a<=Lim && b<=Lim) {
		return a*b;
	}
	if(a<b) {
		swap(a, b);
	}
	ll c=0;
	for(; b; b>>=1, CK(a<<=1)) {
		if(b&1) {
			CK(c+=a);
		}
	}
	return c;
}
void mul(mtx a, mtx b, mtx c, int la, int lb, int lc) {
	static mtx t;
	memset(t, 0, sizeof t);
	for(int i=0; i<la; ++i) {
		for(int j=0; j<lc; ++j) {
			for(int k=0; k<lb; ++k) {
				CK(t[i][j]+=mul(a[i][k], b[k][j]));
			}
		}
	}
	memcpy(c, t, sizeof t);
}
ll b, d, n;
bool spj(ll n) {
	if(n==1) {
		printf("%lld\n", (ll)((((double)b+sqrt(d))/2.0)));
	}
	else if(n==2) {
		printf("%lld\n", (b*b+d)/2);
	}
	return n<=2;
}
mtx a, c;
int main() {
	scanf("%lld%lld%lld", &b, &d, &n);
	if(spj(n)) {
		return 0;
	}
	ll t1=b, t2=(d-b*b)/4;
	CK(t1),	CK(t2);
	a[0][0]=t1, a[0][1]=1;
	a[1][0]=t2, a[1][1]=0;
	c[0][0]=c[1][1]=1;
	for(ll tt=n-2; tt; tt>>=1, mul(a, a, a, 2, 2, 2)) {
		if(tt&1) {
			mul(c, a, c, 2, 2, 2);
		}
	}
	ll a2=(b*b+d)/2, a1=b;
	CK(a1), CK(a2);
	ll ans;
	CK(ans=mul(a2, c[0][0])+mul(a1, c[1][0]));
	if(b*b!=d && (n&1)==0) {
		if(ans==0) {
			ans=mo-1;
		}
		else {
			ans--;
		}
	}
	printf("%llu\n", ans);
	return 0;
}