根据OEIS A028452中的公式有
\begin{align*}
a(n)&=679a(n-1)-76177a(n-2)+3519127a(n-3)-85911555a(n-4)+1235863045a(n-5)-11123194131a(n-6) \\
& +65256474997a(n-7)-257866595482a(n-8)+705239311926a(n-9)-1363115167354a(n-10)+1888426032982a(n-11) \\
& -1888426032982a(n-12)+1363115167354a(n-13)-705239311926a(n-14)+257866595482a(n-15)-65256474997a(n-16) \\
& +11123194131a(n-17)-1235863045a(n-18)+85911555a(n-19)-3519127a(n-20)+76177a(n-21)-679a(n-22)+a(n-23)
\end{align*}
其中$a(n)=f(2n)$。将这一线性递推式写成矩阵形式,再用矩阵快速幂(例如第258题的方法)即可得结果$96972774$。
n=(10**10000)>>1
m=100000007
a=[1, 229, 117805, 64647289, 69563725, 96149360, 40041351, 13625499, 49444743, 7047816, 63444149, 71774943, 18904145, 72062, 58262981, 68125412, 64015488, 12253197, 9101639, 63312159, 47925805, 75412110, 81862504]
p=[1, 99999328, 76177, 96480880, 85911555, 64137046, 23193354, 43529574, 66577436, 60737445, 15071937, 74099213, 25900794, 84928070, 39262562, 33422571, 56470433, 76806653, 35862961, 14088452, 3519127, 99923830, 679]
def cal(x,y):
z=[0]*46
for i in range(23):
for j in range(23):
z[i+j]=(z[i+j]+x[i]*y[j])%m
for i in range(45,22,-1):
for j in range(23):
z[i+j-23]=(z[i+j-23]+z[i]*p[j])%m
return z
v=[0]*23
b=[0]*23
v[0]=b[1]=1
while n:
if (n & 1)==1:
v=cal(v,b)
b=cal(b,b)
n>>=1
print(sum(v[i]*a[i]%m for i in range(23))%m)
|