记$p(n)=n^4+4n^3+2n^2+5n$，首先注意到
$$p(-1)+p(1)=2b+2,\quad p(-2)+p(2)=8b+32,$$
而$\gcd(2b+2,8b+32)=\gcd(2b+2,24)$，因此m必定是24的约数（进一步计算其它$p(-n)+p(n)$的值并不能得到更强的结论）。
经过简单的分析可以得到
\begin{align*}
2 | m &\Leftrightarrow a+b+c=1 \pmod 2 \\
3 | m &\Leftrightarrow b=2\pmod 3, \text{且 } a+c=0 \pmod 3 \\
4 | m &\Leftrightarrow a=c=0\pmod 2, \text{且 } a+b+c=3 \pmod 4 \\
8 | m &\Leftrightarrow a=c=2\pmod 4, \text{且 } a+b+c=7 \pmod 8
\end{align*}
因此可以首先按a,b,c模24的不同值预处理找出对应的$M(a,b,c)$。
接下来由斐波那契数列的定义可得其模任意整数所得的序列必定都有周期性，再进一步分块求和。
最后，$10^9=2^9\times 5^9$，可分别计算$\sum_{k=2}^{1234567890123}S(F_k)$模$2^9$和模$5^9$的结果，再用中国剩余定理加以组合，最终结果是$356019862$。
注：中国剩余定理的计算中需求乘法逆，此处用sympy实现，因此代码为Python 3。此外，$5^9$为千万级，因此运行需要若干分钟，但未打印进度。
import sympy
target=1234567890123
block=24
z=[[0]*block,[0]*block,[0]*block,[0]*block]
for a in range(1,block+1):
for b in range(1,block+1):
for c in range(1,block+1):
d=1
if b%3==2 and (a+c)%3==0:
d*=3
if a%4==2 and c%4==2 and (a+b+c)%8==7:
d*=8
elif a%2==0 and c%2==0 and (a+b+c)%4==3:
d*=4
elif (a+b+c)%2==1:
d*=2
for j in range(block):
z[0][j]+=d
z[1][j]+=d*((j>=a)+(j>=b)+(j>=c))
z[2][j]+=d*((j>=max(a,b))+(j>=max(a,c))+(j>=max(b,c)))
if j>=a and j>=b and j>=c:
z[3][j]+=d
x,y=1 << 9,5**9
res=[]
for p in (x,y):
m=p*block
q=target//m
r=target%m
a,b=0,1
s=t=0
for i in range(m):
u=(b//block)%p
v=b%block
k=(z[0][v]*pow(u,3,p)+z[1][v]*pow(u,2,p)+z[2][v]*u+z[3][v])%p
s=(s+k)%p
if i < r:
t=(t+k)%p
a,b=b,(a+b)%(m)
print((q*s+t-z[3][1])%p)
res.append((q*s+t-z[3][1])%p)
print((sympy.mod_inverse(y%x,x)*res[0]*y+sympy.mod_inverse(x%y,y)*res[1]*x)%(10**9))
|
