Investigating the Fibonacci numbers modulo m

By

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.

Leonardo of Pisa

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 factorization method

Surprisingly, the Fibonacci numbers can also be used as an integer factorization method, similar in flavor to Pollard's p-1 and Williams' p+1 methods.

The method uses the smallest divisor `d` of `p - (5/p)` for which `F(d) ≡ 0 mod p`, where `(5/p)` is the Legendre symbol.

By selecting a small bound `B`, we compute:

`k = \prod_{p<=B} p^floor(log(B)/log(p))`

which is equivalent with:

`k = lcm(1..B)`

If `p - (5/p)` is B-smooth, then `k` is most likely a multiple of `d`, which gives us the prime factor `p` of `n` as:

`p = gcd(F(k) mod n, n)`

For example, let:

`n = 257221 * 470783`

It can be verified that `(5/470783) = -1` and `470783+1` is 700-smooth, which allows us to factorize `n` as:

`k = lcm(1..700)`
`F(k) ≡ 101758333101 mod n`

Then, by computing the greatest common divisor, we find a non-trivial factor of `n`:

`gcd(n, F(k) mod n) = 470783`

# 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 - (5/n)) ≡ 0 mod n`
`F(n + (5/n)) ≡ 1 mod n`

where `(5/n)` 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, ...`

# Constructing Fibonacci pseudoprimes

Given a list of prime numbers, we can generate a Fibonacci pseudoprime if two or more primes `p` in this list share a divisor `d` of `p - (5/p)`, for which `F(d) ≡ 0 mod p`, where `(5/n)` is the Legendre symbol.

For example, let:

`p = 1410364649`
`q = 8462187889`

The Legendre symbol for both primes is `1`, and we have:

`p-1 = 2^3 * 7^2 * 11 * 327079`
`q-1 = 2^4 * 3 * 7^2 * 11 * 327079`

Because they share many prime factors, there is a good chance that a divisor `d` of `p-1` also divides `q-1` and:

`F(d) ≡ 0 mod p`
`F(d) ≡ 0 mod q`

If that is the case, then this gives us a new Fibonacci pseudoprime: `n = p * q`, which satisfies:

`F(n - (5/n)) ≡ 0 mod n`

In our example, one such divisor is `4579106 = 2 * 7 * 327079`:

`F(4579106) ≡ 0 mod 1410364649`
`F(4579106) ≡ 0 mod 8462187889`

[Try it Online!]

This method works for any number of primes and can be used in constructing arbitrary large Fibonacci pseudoprimes.

# Twin primes and Fibonacci pseudoprimes

Let `p` and `p+2` be prime numbers, with the condition that:

`p ≡ 2 mod 5`

Then `n = p * (p+2)` is a Fibonacci pseudoprime, for which the Legendre symbol `(5/n) = -1`, meaning that `F(n + 1)` is divisible by `n`.

Also, since `p` has the form `p = 5x + 2`, for some odd integer `x`, then `(5x + 2) * (5x + 4) = (5x + 3)^2-1`, from which results that `p * (p+2) + 1` is a square.

This observation gives us an easy way to generate arbitrary large Fibonacci pseudoprimes (which are also Lucas pseudoprimes) that are congruent to `3 mod 5`, as illustrated bellow in an example:

`p = 15868511662415791787`
`q = p + 2`
`n = p q`

`F(n + 1) ≡ 0 mod n`

[Try it Online!]

More generally, we can generate random values of `k` and testing if `(5*(2k+1) + 3)^2 - 1` is a Fibonacci pseudoprime, which is very likely to happen quite often, even for large `k`.

# 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` for which `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.

Additionally, it is well-known that `gcd(F(n), F(k*n)) = F(n)` for any `k >= 1`.

From this, we can ask the question: what is the smallest integer `r > 1` such that `gcd(F(n), F(r)) != 1`?

If `n` is a prime number greater than `2`, the smallest integer `r > 1` that satisfies the above inequality is `r = n`.

Furthermore, for an odd composite number `n>=3`, the smallest number `r > 1` is the smallest prime factor of `n`.

More generally, if `d` divides `n`, then `F(d)` divides `F(n)`, which is useful in factoring Fibonacci numbers at composite indices.

# 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:
See also:

Comments

Post a Comment