题目链接
题解
容易想到(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[n])
[ans = f[n] - (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;
}