看到这个 \(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;
}