BZOJ4002 [JLOI2015]有意义的字符串 【数学 + 矩乘】

题目链接

BZOJ4002

题解

容易想到\(\frac{b + \sqrt{d}}{2}\)是二次函数\(x^2 - bx + \frac{b^2 - d}{4} = 0\)的其中一根
那么就有
\[x^2 = bx - \frac{b^2 - d}{4}\]
两边乘一个\(x^n\)
\[x^n = bx^{n - 1} - \frac{b^2 - d}{4}x^{n - 2}\]

再观察题目条件，可以发现\(|b^2 - d| < 1\)，所以明显要用到另一个根\(\frac{b - \sqrt{d}}{2}\)
我们设
\[f[i] = (\frac{b + \sqrt{d}}{2})^i + (\frac{b - \sqrt{d}}{2})^i\]
那么就有
\[f[i] = bf[i - 1] - \frac{b^2 - d}{4}f[i - 2]\]
矩乘优化一下就可以算出\(f
\)
\[ans = f
- (\frac{b - \sqrt{d}}{2})^n\]
后面这个玩意是小于\(1\)的，所以我们只需要讨论一下其正负就可以判定出应该向哪边取整了

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define ULL unsigned long long int
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define BUG(s,n) for (int i = 1; i <= (n); i++) cout<<s[i]<<' '; puts("");
using namespace std;
const int maxn = 100005,maxm = 100005,INF = 1000000000;
const ULL P = 7528443412579576937ll;
ULL mul(ULL a,ULL b){
ULL re = 0; b = (b % P + P) % P; a = (a % P + P) % P;
for (; b; b >>= 1,a = (a + a) % P) if (b & 1) re = (re + a) % P;
return re;
}
struct Matrix{
ULL s[2][2]; int n,m;
Matrix(){memset(s,0,sizeof(s)); n = m = 0;}
};
inline Matrix operator *(const Matrix& a,const Matrix& b){
Matrix c;
if (a.m != b.n) return c;
c.n = a.n; c.m = b.m;
for (int i = 0; i < c.n; i++)
for (int j = 0; j < c.m; j++)
for (int k = 0; k < a.m; k++)
c.s[i][j] = ((c.s[i][j] + mul(a.s[i][k],b.s[k][j])) % P + P) % P;
return c;
}
inline Matrix qpow(Matrix a,ULL b){
Matrix c; c.n = c.m = a.n;
for (int i = 0; i < c.n; i++) c.s[i][i] = 1;
for (; b; b >>= 1,a = a * a)
if (b & 1) c = c * a;
return c;
}
int main(){
ULL b,d,n;
cin >> b >> d >> n;
if (n == 0){puts("1"); return 0;}
Matrix A,F,Fn;
A.n = A.m = 2;
A.s[0][0] = (b % P + P) % P; A.s[0][1] = (((d - b * b) / 4 % P) + P) % P;
A.s[1][0] = 1; A.s[1][1] = 0;
F.n = 2; F.m = 1;
F.s[0][0] = (b % P + P) % P; F.s[1][0] = 2;
Fn = qpow(A,n - 1) * F;
if (b * b != d && !(n & 1)) cout << ((Fn.s[0][0] - 1) % P + P) % P << endl;
else cout << Fn.s[0][0] << endl;
return 0;
}