### Investigating the Fibonacci numbers modulo m

 Leonardo of Pisa

The Fibonacci sequence is, without doubt, one of the most popular sequences in mathematics and in popular culture, named after Italian mathematician Leonardo of Pisa (also known as Fibonacci, Leonardo Bonacci, Leonardo of Pisa, Leonardo Pisano Bigollo, or Leonardo Fibonacci), who first introduced the numbers in Western European with his book Liber Abaci, in 1202.

The sequence is elegantly defined as:

F(0) = 0
F(1) = 1
F(n) = F(n-1) + F(n-2)

where the first 20 terms are:

0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181

In this post, we investigate the Fibonacci numbers modulo a positive integer m. It is well known that for every positive integer m, the modular Fibonacci sequence F(n) mod m is eventually periodic. This period is called the Pisano period.

For example, let's take a look at F(n) mod 4:

0, 1, 1, 2, 3, 1, 0, 1, 1, 2, 3, 1, 0, 1, 1, 2, 3, 1, ...

we observe, as highlighted above in yellow, the Pisano cycle has a length of 6 for m=4, and continues to repeat forever.

### # Efficient computation of Fibonacci numbers

The n-th Fibonacci number can be efficiently computed in O(log_2(n)) steps with the following algorithm:
func fibonacci(n) {

var (f, g) = (0, 1)
var (a, b) = (0, 1)

for k in (0 .. n.ilog2) {
(f, g) = (f*a + g*b, f*b + g*(a+b)) if n.getbit(k)
(a, b) = (a*a + b*b, a*b + b*(a+b))
}

return f
}

This algorithm can be easily derived by noticing that raising the following 2×2 matrix at increasing powers of n, will generate the Fibonacci numbers:

A = [[1, 1], [1,0]]

For example, raising A to the power 12 we have:

A^12 = [[233, 144], [144, 89]]

from which we can extract the 12-th Fibonacci number, which is F(12) = 144.

However, this algorithm is efficient only if we can efficiently raise a matrix to a large power. Fortunately, there is an efficient method for doing this, called exponentiation by squaring, recursively defined as:

A^n = A * (A*A)^((n-1)/2)   if n is odd
A^n = (A*A)^(n/2)        if n is even

with the base-case A^0 = I, where I is the identity matrix.

For better efficiency, we prefer to implement this method iteratively, as:
func matrix_exp(A, n) {
var B = [[1, 0], [0, 1]]

for k in (0 .. n.ilog2) {
B = B.matrix_mul(A) if n.getbit(k)
A = A.matrix_mul(A)
}

return B
}

By slightly modifying the first algorithm, gives us an efficient method for computing the n-th Fibonacci number modulo m in O(log_2(n)) steps, which will be useful later on:
func fibmod(n, m) {

var (f, g) = (0, 1)
var (a, b) = (0, 1)

for k in (0 .. n.ilog2) {
(f, g) = ((f*a + g*b)%m, (f*b + g*(a+b))%m) if n.getbit(k)
(a, b) = ((a*a + b*b)%m, (a*b + b*(a+b))%m)
}

return f
}

### # Fibonacci primality test

For each prime number n > 5, its corresponding congruence is true:

F(n-1) ≡ 0 mod n,   if n ≡ ±1 mod 5
F(n+1) ≡ 0 mod n,   if n ≡ ±2 mod 5

However, there also exist odd composite numbers n that satisfy this congruences, which are called Fibonacci pseudoprimes (A081264), where the first few terms are:

323, 377, 1891, 3827, 4181, 5777, 6601, 6721, 8149, 10877, 11663, 13201, 13981, 15251, ...

Combining this test with a base-2 Fermat primality test, the rate of false-positives is considerably reduced, where the first few pseudoprimes that pass both tests are (A214434):

6601, 13981, 30889, 68101, 219781, 252601, 332949, 399001, 512461, 642001, 721801, ...

Research problem: find a number n ≡ ±2 mod 5 that is simultaneously a base-2 pseudoprime and a Fibonacci pseudoprime. Pomerance, Selfridge, and Wagstaff offer \$620 for the first example or for a proof that no such number exists.

### # Frobenius x^2-x-1 primality test

A stronger primality test that is satisfied by all primes n > 5, uses the following two congruences:

F(n - (n/5)) ≡ 0 mod n
F(n + (n/5)) ≡ 1 mod n

where (n/5) is the Legendre symbol, where the first few pseudoprimes are (A212424):

4181, 5777, 6721, 10877, 13201, 15251, 34561, 51841, 64079, 64681, 67861, 68251, ...

Combining this test with a base-2 Fermat primality test, the smallest pseudoprimes to both tests are:

219781, 252601, 399001, 512461, 722261, 741751, 852841, 1024651, 1193221, ...

### # Divisibility of Fibonacci numbers

Given a positive integer n, find the smallest positive integer k such that F(k) is divisible by n.

In other words, we want to find the least value of k such that F(k) ≡ 0 mod n. (A001177)

Assuming that the prime factorization of n can be easily computed, the following algorithm efficiently finds the least value of k that satisfies this criteria.

Let:

n = \prod_{k=1} p_k^(e_k)

a(p, 0) = 1
a(5, e) = 5*5^(e-1)
a(2, e) = b(2, e-1),     if e > 2
a(p, e) = b(p, e)

where b(p, e) is defined as:

b(p, e) = p^(e-1) *  smallest d|(p-1) such that F(d) ≡ 0 mod p,    if p ≡ ±1 mod 5
b(p, e) = p^(e-1) *  smallest d|(p+1) such that F(d) ≡ 0 mod p,    if p ≡ ±2 mod 5

Then the smallest positive integer k such that F(k) ≡ 0 mod n is given by:

k = lcm(a(p_1, e_1), a(p_2, e_2), ..., a(p_k, e_k))

The relation a(p, e) = p^(e-1)*a(p,1) is called Wall's conjecture and has been verified for primes up to 10^14. Primes for which this relation fails are called Wall-Sun-Sun primes, although no such prime number is currently known.

### # The Pisano period

Given a positive integer n, find the smallest positive integer t such that F(t) ≡ 0 mod n and F(t+1) ≡ 1 mod n.

It is known that a zero can occur in a Pisano cycle either once, twice or four times only.

Using the algorithm described in the previous section, we compute the least value of k such that F(k) ≡ 0 mod n, then we can jump through all the zeros in the Pisano cycle by multiplying k by 1, 2 and 4.

The sequence of numbers n for which the Pisano cycle contains only one zero is given in A053031.

If the cycle contains only one zero, the length of the cycle is k.
If the cycle contains only two zeros, the length of the cycle is 2k.
If the cycle contains four zeros, the length of the cycle is 4k.

In practice, we test the following congruence:

F(v k + 1) ≡ 1 mod n

for each v={1, 2, 4}. The first value v k + 1 that satisfies the congruence, gives us: t = v k, which is the length of the Pisano period modulo n.

Bellow we have a simple implementation of this algorithm, which uses the "fibmod(n, m)" function defined in the first section:
func prime_power_divisor(p, k=1) {
divisors(p - legendre(p, 5)).first_by {|d| fibmod(d, p) == 0 } * p**(k-1)
}

func pisano_period(n) {

return 0 if (n <= 0)
return 1 if (n == 1)

var d = Math.lcm(n.factor_exp.map {|p| prime_power_divisor(p[0], p[1]) }...)

for k in ([1, 2, 4]) {
if (fibmod(k*d + 1, n) == 1) {
return k*d
}
}
}
[Try it Online!]

### # Implementations

The presented algorithms are implemented in the Sidef programming language: