Project Euler 304

Project Euler 304

题目

Primonacci

For any positive integer n the function next-prime(n) returns the smallest prime p such that p>n.

The sequence a(n) is defined by:a(1)=next-prime(1014) and a(n)=next-prime(a(n1)) for n>1.

The fibonacci sequence f(n) is defined by: f(0)=0,f(1)=1 and f(n)=f(n1)+f(n2) for n>1.

The sequence b(n) is defined as f(a(n)).

Find b(n) for 1n100000. Give your answer mod 1234567891011.

解决方案

将斐波那契数的通项公式写成矩阵形式:

[0111][f(i2)f(i1)]=[f(i1)f(i)]

B=[0111],那么通过矩阵快速幂,可以 O(logn) 的时间复杂度内求解出第 n 项的值。

另外,上面的等式可以写成

Bd[f(id1)f(id)]=[f(i1)f(i)]

105 次的询问是独立的,并且之间的间隔 d 也比较小。我们可以基于上式来减少矩阵乘法的运算次数,用上一次的结果继续转移,计算出当前结果。

求一个数的下一个质数,使用 sympynextprime 方法。

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
from sympy import nextprime

M = 10 ** 14
Q = 10 ** 5
mod = 1234567891011


def mul(a: list, b: list):
return [[sum(a[i][k] * b[k][j] for k in range(len(b))) % mod for j in range(len(b[0]))] for i in range(len(a))]


b_pre = [[[0, 1], [1, 1]]]
for i in range(len(bin(M+Q))+2):
b_pre.append(mul(b_pre[-1], b_pre[-1]))


def cal(a, n):
i = 0
while n:
if n & 1:
a = mul(a, b_pre[i])
i += 1
n >>= 1
return a


a = [[0, 1]]
ans = 0
pre = 0
for n in range(Q):
M = nextprime(M)
val = M - pre
a = cal(a, val)
ans = (ans + a[0][0]) % mod
pre = M
print(ans)