Description：

输入 $b,d,n$ ，求： $$ \lfloor \left ( \frac{b+\sqrt{d}}{2} \right ) ^n \rfloor \mathrm{mod} \ p $$ $1\leqslant n\leqslant 10^{18},p=7528443412579576937,b\%2=1,d\%4=1,0<b^2\leqslant d<(b+1)^2\leqslant10^{18}$

Solution：

直接拿 $double$ 快速幂精度就挂了，考虑用一些别的方法，观察式子 $\frac{b+\sqrt{d}}{2}$ ，可以发现这个是二次方程 $x^2-bx+\frac{b^2-d}=4$ 的解，移个项就是 $x^2=bx+\frac{d-b^2}4$ ，两边同时乘以 $x^{n-2}$ 就是 $x^n=bx^{n-1}+\frac {d-b^2}4x^{n-2}$ ，如果设 $f[i]=\left(\frac{b+\sqrt{d}}{2}\right)^i+\left(\frac{b-\sqrt{d}}{2}\right)^i$ ，由于 $b\%2=1$ ，所以 $b^2\%4=1$ ，所以 $\frac{d-b^2}4$ 是整数，所以 $f[i]$ 是整数，那就可以把上面那个看成一个递推式然后跑矩阵快速幂，这样就可以得到了 $f[n]$ 的值，因为 $b^2\leqslant d<(b+1)^2$ ，所以 $b\leqslant\sqrt d<b+1$ ，也就是说 $-1<b-\sqrt d\leqslant 0$ ，那么当 $n$ 为奇数且 $b^2 \ne d$ 时 $ans=f[n]$ 否则 $ans=f[n]-1$ ，因为数太大所以要用快速乘，还要用 $unsigned\ long\ long$ 。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

typedef unsigned long long ll;

ll b,d,n;

#define MOD 7528443412579576937

ll mul(ll a,ll b)

{

	ll res = 0;

	while(b > 0)

	{

		if(b & 1)res = (res + a) % MOD;

		a = (a + a) % MOD;

		b = b >> 1;

	}

	return res;

}

struct matrix

{

	ll m[3][3];

	matrix(){memset(m,0,sizeof(m));}

	friend matrix operator * (matrix a,matrix b)

	{

		matrix c;

		for(int i = 1;i <= 2;++i)

			for(int j = 1;j <= 2;++j)

				for(int k = 1;k <= 2;++k)

					c.m[i][j] = (c.m[i][j] + mul(a.m[i][k],b.m[k][j])) % MOD;

		return c;

	}

	friend matrix operator ^ (matrix a,ll b)

	{

		matrix res = a;--b;

		while(b > 0)

		{

			if(b & 1)res = res * a;

			a = a * a;

			b = b >> 1;

		}

		return res;

	}

}f;

int main()

{	

	cin >> b >> d >> n;

	if(n == 0)

	{

		cout << 1 << endl;

		return 0;

	}

	ll b_ = b / 2;

	ll c = ((d / 4 + MOD - mul(b_,b_)) % MOD + MOD - b_) % MOD;

	f.m[1][1] = b;f.m[2][1] = c;f.m[1][2] = 1;f.m[2][2] = 0;

	f = f ^ (n - 1);

	ll ans = (mul(2,f.m[2][1]) + mul(b,f.m[1][1])) % MOD;

	if(n % 2 == 0 && b * b != d)cout << (ans + MOD - 1) % MOD << endl;

	else cout << ans << endl;

	return 0;

}