令$b_n=6a_n+5$，则有$b_1=11, b_{n+1}=b_n^2-2$。
若将$b_1$写作$k+1/k$的形式，则有$k^2+1/k^2=(k+1/k)^2-2$，从而
$$b_n=k^{2^{n-1}}+\dfrac{1}{k^{2^{n-1}}}.$$
$k+1/k=11$的解是$(11\pm3\sqrt{13})/2$，因此若13是p的二次剩余，则可以直接在有限域$\mathbb Z_p^*$中解出$k$，否则可在域$\mathbb Z_p(\sqrt{13})$中解出$k$，这是因为$\mathbb Z_p^*(\sqrt{13})$中的数都具有$a+b\sqrt{13}$的形式，结果中$k^m+k^{-m}$保证了$\sqrt{13}$前的系数必定为$0$。
代入计算即可得最终结果$242586962923928$。
注：因使用sympy作$\mathbb Z_p$中的相关计算，以下代码为Python 3。
import sympy
start=10**9
interval=10**7
n=10**15
f=lambda ax,ay,ux,uy:((ax*ux+ay*uy*13)%p,(ax*uy+ay*ux)%p)
res=0
for p in sympy.primerange(start,start+interval):
inv2=(p+1)//2
if pow(13,(p-1)//2,p)==1:
x=pow((3*sympy.sqrt_mod(13,p)+11)*inv2%p,pow(2,n-1,p-1),p)
s=x+sympy.mod_inverse(x,p)
else:
ux,uy=11*inv2%p,3*inv2%p
m=pow(2,n-1,p*p-1)
ax,ay=1,0
while m:
if m%2:
ax,ay=f(ax,ay,ux,uy)
ux,uy=f(ux,uy,ux,uy)
m//=2
s=2*ax
res+=(s-5)*sympy.mod_inverse(6,p)%p
print(res)
|
