import sympy
import math
import time
It is possible to compute the (n+1)-th Fibonacci number $F_{n+1}$ using a matrix
$$ \mathbf Q = \begin{pmatrix} 0 & 1 \\ 1 & 1 \end{pmatrix}, $$
by multplying $\mathbf Q$ with a vector containing the current last 2 terms of the sequence:
$$ \begin{align} \begin{pmatrix} 0 & 1 \\ 1 & 1 \end{pmatrix} \cdot \begin{pmatrix} F_{n-2} \\ F_{n-1} \end{pmatrix} &= \begin{pmatrix} F_{n-1} \\ F_{n-2} + F_{n-1} \end{pmatrix}\\ &= \begin{pmatrix} F_{n-1} \\ F_{n} \end{pmatrix} \end{align} $$
The next element can be computed like this:
$$ \begin{align} \begin{pmatrix} 0 & 1 \\ 1 & 1 \end{pmatrix} \cdot \begin{pmatrix} F_{n-1} \\ F_{n} \end{pmatrix} &= \begin{pmatrix} 0 & 1 \\ 1 & 1 \end{pmatrix}^2 \begin{pmatrix} F_{n-2} \\ F_{n-1} \end{pmatrix} \\ &= \begin{pmatrix} F_{n} \\ F_{n-1} + F_{n} \end{pmatrix} \\ &= \begin{pmatrix} F_{n} \\ F_{n+1} \end{pmatrix} \end{align}. $$
Note that we only needed to square the matrix $\mathbf Q$ and multiply it with $(F_{n-2} \ \ F_{n-1})^\mathsf T$ in order to obtain $F_{n+1}$.
Generally, we find:
$$ \begin{align} \begin{pmatrix} F_{n+k-1} \\ F_{n+k} \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ 1 & 1 \end{pmatrix}^{k+1} \begin{pmatrix} F_{n-2} \\ F_{n-1} \end{pmatrix} , \ k \geq 0. \end{align} $$
If we start with $F_{n-2}=0$ and $F_{n-1}=1$, we have ($n=2$):
$$ \begin{align} \begin{pmatrix} F_{2+k-1} \\ F_{2+k} \end{pmatrix} = \begin{pmatrix} F_{k+1} \\ F_{k+2} \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ 1 & 1 \end{pmatrix}^{k+1} \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \ k \geq 0 \end{align} $$
or equivalently:
$$ \begin{align} \begin{pmatrix} F_{k-1} \\ F_{k} \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ 1 & 1 \end{pmatrix}^{k-1} \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \ k \geq 2 \end{align} $$
If we are only interested in $F_k$, we can simply write:
$$ \begin{align} F_{k} = \begin{pmatrix} 0 & 1 \\ 1 & 1 \end{pmatrix}^{k-1} \Bigg|_{2,2}, \ k \geq 2, \end{align} $$
where the index $2,2$ refers to the element in the bottom-right corner.
Q = sympy.ones(2)
Q[0,0] = 0
Q
For example, we can compute $F_{10}=55$ with:
Q**(10-1)
We find $F_{10}$ in the bottom right corner of the resulting matrix.
We get $F_{100}=354,224,848,179,261,915,075$ with:
Q**(100-1)
Furthermore, we can use a combination of exponential squaring and modular exponentation for fast computation:
\begin{align} Q^2 &\equiv (Q\bmod m) \cdot (Q\bmod m) \pmod{m} \\ Q^4 &\equiv (Q^2\bmod m) \cdot (Q^2\bmod m) \pmod{m} \\ Q^8 &\equiv (Q^4\bmod m) \cdot (Q^4\bmod m) \pmod{m} \\ & \ \ \vdots \end{align}
Additionally, if we have already obtained the power $Q^p$, it is possible to calculate $Q^{p'} = Q^{p+\Delta p'}$:
\begin{align} Q^{p'} &\equiv (Q^p\bmod m) \cdot (Q^{\Delta p'}\bmod m) \pmod{m}, \\ \end{align}
where $(Q^{\Delta p'}\bmod m)$ can be determined based on the exponential squaring rules above.
Now, we have everything to solve our example problem stated earlier.
First we specify our constants:
m = 1110987654321 # modulus
p_min = 10**15 # the first prime should be greater than this value
N = 1000 # We want to compute 1000 Fibonacci numbers F_p for the given primes
q = 10 # We want to use every q-th prime to compute our Fibonacci numbers
1. Collect every q-th prime number larrger than p_min¶
start = time.time()
prime_list = []
p = p_min
for _ in range(N):
p = sympy.nextprime(p, q)
prime_list.append(p)
assert len(prime_list) == N
print("Done computing primes. Time in seconds:", round(time.time() - start, 2))
# This is already the case in this example, but in general the list has to
# be sorted
prime_list = sorted(prime_list)
print("First 10 primes in the list:")
print(prime_list[:10])
Done computing primes. Time in seconds: 4.08 First 10 primes in the list: [1000000000000279, 1000000000000783, 1000000000001209, 1000000000001623, 1000000000002193, 1000000000002467, 1000000000002683, 1000000000003001, 1000000000003309, 1000000000003871]
2. Create a cache for $\mathbf Q$, $\mathbf Q^2$, $\mathbf Q^4$, $\mathbf Q^8$, $\ldots$¶
# 2. Compute powers of the Fibonacci matrix Q^1, Q^2, Q^4, Q^4, etc. and store
# them for later
start = time.time()
i = 1
fibo_power_cache = [] # contains the powers of the fibonacci matrix to 1,2,4,8,16,...
Q = sympy.ones(2)
Q[0,0] = 0
while i < max(prime_list):
fibo_power_cache.append(Q)
Q = (Q @ Q) % m
i = i * 2
print("Done creating cache. Time in seconds:", round(time.time() - start, 2))
Done creating cache. Time in seconds: 0.02
3. Compute the actual Fibonacci Numbers $F_p \bmod m$ for the given primes $p$¶
# 3. Compute the Fibonacci Numbers F_p for the given primes p
# Since we have to compute Q^p, we can use the precomputed powers in the cache.
# E.g., Q^11 = Q^8 @ Q^2 @ Q^1
# and Q^14 = Q^11 @ Q^2 @ Q^1
start = time.time()
last_p = 1
fib_list = []
fib = sympy.eye(2)
for idx, p in enumerate(prime_list):
p_diff = p - last_p
while p_diff > 0:
e = p_diff.bit_length() - 1 # math.floor(math.log2(p_diff))
fib @= fibo_power_cache[e]
fib = fib % m
p_diff -= 2**e
last_p = p
fib_list.append((p, fib[-1,-1]))
print("Done computing Fibonacci numbers. Time in seconds:", round(time.time() - start, 2))
print("Last 10 Fibonacci numbers:")
for i in range(10):
print(f"F_{fib_list[-i-1][0]} mod {m} = {fib_list[-i-1][1]}")
Done computing Fibonacci numbers. Time in seconds: 1.19 Last 10 Fibonacci numbers: F_1000000000349191 mod 1110987654321 = 382892723032 F_1000000000348753 mod 1110987654321 = 1071406309063 F_1000000000348337 mod 1110987654321 = 618378886060 F_1000000000348009 mod 1110987654321 = 71281864435 F_1000000000347763 mod 1110987654321 = 584619443096 F_1000000000347581 mod 1110987654321 = 970898498456 F_1000000000347211 mod 1110987654321 = 35392582421 F_1000000000346819 mod 1110987654321 = 982178757842 F_1000000000346411 mod 1110987654321 = 520718657315 F_1000000000346197 mod 1110987654321 = 894644052794