# Easier program for period of Fibonacci sequence modulo p

For a little project I need to calculate the period of a Fibonacci sequence modulo p, for which p is a prime number. For example, the Fibonacci sequence modulo 19 would be: $$0, 1, 1, 2, 3, 5, 8, 13, 2, 15, 17, 13, 11, 5, 16, 2, 18, 1, 0, 1, 1, 2...$$ As you can see, the sequence repeats after 18 numbers.

The formula for a Fibonacci sequence is given by: $f(n)=f(n-1)+f(n-2)$. In my attempt at making a program, I came up with the following:

newrow[a_, b_, c_] := Module[{f, g},
f[n_] := f[n] = f[n - 1] + f[n - 2];
f = a;
f = b;
g[n_] := Mod[f[n], Prime[c]];
L = Table[g[n], {n, 0, Prime[c]^2} ];
number1 = Flatten[Position[L, g]];
number2 = Flatten[Position[L, g]];
number3 = Select[number1, (MemberQ[number2, # + 1] ) &];
period = number3[] - 1;
Return[
period
]
]


However, if I want to calculate it for quite some prime numbers at once like this:

Table[newrow[0, 1, i], {i, 1, 50}]


the computer can't handle it. Might there be an easier algorithm for this? I also tried this, but this seems to be even harder for the computer:

anderrijtje[a_, b_, c_] := Module[{f, g},
f[n_] := f[n] = f[n - 1] + f[n - 2];
f = a;
f = b;
g[n_] := Mod[f[n], c];
Return[
For[
k = 1, Not[{g, g} == {g[k], g[k + 1]}], k++, {}
];
k
]
]


Thanks in advance for any help.

Update: now with whuber's suggested improvement.

fibMod[0, b_, m_] := i MultiplicativeOrder[b, m]
fibMod[a_, b_, m_] := fibMod[i++; b, Mod[a + b, m], m]
fibMod[m_] := Block[{i = 1}, fibMod[1, 1, m]]


Test:

fibMod

18


For really long periods you'll need to increase $IterationLimit: $IterationLimit = 1*^6;

fibMod[Prime@15000]

163840


Here is a compiled alternative for greater speed, though if taken to an extreme the terms could exceed machine size integers.

cf =
Compile[{{m, _Integer}},
Module[{a, b, i},
For[a = b = i = 1,
a != 0,
i++,
{a, b} = {b, Mod[a + b, m]}
];
i * MultiplicativeOrder[b, m]
]
];

cf[Prime@145000] // Timing

{0.281, 3885496}

• +1. You can speed up some cases by changing the first line to fibMod[0, n_, p_] := i MultiplicativeOrder[n, p]. This algorithm stops the first time a $0$ reappears. Suppose it's followed by $k$. This means the remainder of the sequence begins $0, k, \ldots$--but that will be exactly $k$ times the sequence beginning $0, 1, \ldots$, because the recursion is linear. Thus we find the sequence breaks down into $q$ length-$i$ subsequences $0, 1, \ldots, 0, k, \ldots, 0, k^2, \ldots, 0, 1$ where $k^q=1$: that's precisely what MultiplicativeOrder computes. The speedup is by a factor of $q$. Jan 26 '13 at 21:48
• @whuber Great observation! I'll amend my answer. Jan 27 '13 at 3:22
• @Mr.Wizard Thanks so much, this is wonderful. I do not fully understand the code you wrote, but it works perfect. I wonder how I can change the starting positions of my sequence, like in my attempt. For example, I start with $2, 2, 4, 6...$. Jan 27 '13 at 15:24
• Jeroen, my answer begins with some analytical remarks that partially address that generalization. For instance, the sequence $2,2,4,6,\ldots$ will be term for term twice the sequence $1,1,2,3,\ldots$, which in turn is just the usual Fibonacci sequence $0,1,1,2,3,\ldots$ shifted one place to the left. As long as you are not working modulo $2$, then, the two sequences will have exactly the same periods. Jan 27 '13 at 17:57

We can substantially speed up the calculation for large primes by making some elementary observations about the Fibonacci series. There are two motivating ideas behind them:

1. Almost all the properties of the Fibonacci series rely only on the field axioms, so that the theory of Fibonacci series (and linear recurrences in general) holds practically without change in (almost) any field, not just the rational or real numbers. This include the fields of integers modulo primes $p$.

2. Solutions of linear recurrences look like geometric series. (The simplest such, of order $1$, is determined by a starting value $a$ and a common ratio $r$, with the recursive definition $f_0 = a$ and $f_{n+1} = r f_n$. A closed formula is $f_n = a r^n$.) Thus we can expect to obtain closed-form formulas even when working modulo $p$ for primes $p$. These can be exploited to limit the search for the point where the series starts repeating.

### Analysis

We work first with coefficients in an arbitrary commutative ring $R$, and later specialize to the finite fields $\mathbb{F}_p$ (which are the rings of integers modulo $p$).

First, because the Fibonacci series $(F_n)$ is generated by a reversible recurrence--that is, we can recover $F_{n-1} = F_{n+1} - F_n$ from its succeeding terms--it is completely determined by any ordered pair $(F_n, F_{n+1})$ occurring within it, whence such ordered pairs all determine the same period. For instance, modulo $7$ the period for the series beginning $\color{red}{5, 1,} 6, 0, 6, 6, 5, 4, \ldots$ must be the same as the period for the Fibonacci series $0, 1, 1, 2, 3, \color{red}{5, 1,} \ldots$ because $5,1$ occurs in both.

Second, because the recurrence is linear, the period determined by an ordered pair $(f, g)$ is the same as that determined by $(\lambda f, \lambda g)$ for any invertible $\lambda$. For instance, modulo $3$ the series $0, 2, 2, 1, \ldots$ is, term-by-term, twice the series $0, 1, 1, 2, \ldots$ because its starting values $(0,2)$ are twice that of the Fibonacci starting values $(0,1)$.

Third, let $\varphi_+$ and $\varphi_{-}$ be the two roots of the polynomial $x^2-x-1$ in $R$ or some extension thereof (in case the roots are "imaginary"). (Notice this asserts two things about the roots: they sum to $1$ and their product is $-1$.) Then, exactly as always, we can obtain a closed form formula for the $n$th Fibonacci number in $R$,

$$F_n = \frac{\varphi_+^n - \varphi_-^n}{\varphi_+ - \varphi_-},$$

assuming that the denominator $\varphi_+ - \varphi_-$ is invertible.

For example, we can see this inductively by verifying that (a) when $n=0$ and $n=1$ are plugged into the formula we recover $0$ and $1$, respectively; and (b) assuming the formula works for $n$, we can show it must work for $n+1$ because

\eqalign{ F_{n+1} = \varphi_+^{n+1} - \varphi_-^{n+1} &= \varphi_+^2(\varphi_+^{n-1}) - \varphi_-^2(\varphi_-^{n-1}) \\ &= (\varphi_+ + 1)(\varphi_+^{n-1}) - (\varphi_{-}+1)(\varphi_-^{n-1}) \\ &= (\varphi_+^{n} - \varphi_-^{n}) + (\varphi_+^{n-1} - \varphi_-^{n-1}) \\ &= F_n + F_{n-1}. }

The crux was the expansion (on the second line) of $\varphi^2 = \varphi + 1$ for each root $\varphi$ of $x^2 - x - 1$.

From now on assume $R$ is actually a field: that is, every nonzero element has a reciprocal. This guarantees that $\varphi_+ - \varphi_-$ is invertible unless both roots are the same. (It is easy to show that the roots are the same only when $p=5$, for only then does the derivative $2x-1$ divide $x^2-x-1$.)

Fourth, this key polynomial equation $x^2 - x - 1= 0$ can be rewritten (by completing the square and multiplying by $4$) as

$$(2x-1)^2 = 5.$$

(The multiplication by $4$ is not invertible in fields of characteristic $2$, so for now on let's assume $p \ne 2$.) Whence

• Both roots lie in $R$ when $5$ is a square (and could be suggestively written $(\pm\sqrt{5}+1)/2$).

• Otherwise, neither root lies in $R$ and instead we need to work within the extension of $R$ generated by these roots.

Quadratic Reciprocity informs us that $5$ is a square modulo $p \ne 5$ if and only if $p$ is a square modulo $5$. The squares (of primes) mod $5$ can be exhaustively enumerated: they are $1$ and $4=-1$.

For example, because $7 = 2$ mod $5$ is not a square, when we work modulo $7$ the roots will lie in the extension field $\mathbb{F}_{7^2}$: they are "imaginary." Because $11 = 1$ mod $5$ is a square, the roots lie in $\mathbb{F}_{11}$. Indeed, $x^2 - x - 1 = (x-4)(x-8)$ mod $11$, whence we can take $\varphi_+=8$ and $\varphi_{-}=4$.

### Application

We are going to write somewhat efficient code, order[p], to find the period (or "order") of the Fibonacci sequence modulo $p$.

We previously identified two special primes--$2$ and $5$, for which the periods are $3$ and $20$, respectively--and otherwise there are two general cases to consider, depending on whether $p$ is in $\{1,-1\}$ or $\{2,-2\}$ modulo $5$. It turns out that the algorithms will work fine for $2$, so we only have to code one special case:

order = 20


The case where both roots exist

The trick is to spot when the sequence is repeating. To save effort, rather than looking for the next occurrence of $0,1$ in the sequence, we only have to find the next $0$. Specifically, suppose $F_n=0$ but all of $F_2, F_3,\ldots, F_{n-1}$ are nonzero. Then the effect of going $n$ steps from the initial pair $(F_0, F_1) = (0,1)$ is to produce $(F_n, F_{n+1}) = (0, F_{n+1})$. This multiplies the initial pair by $F_{n+1}$. By virtue of the linearity, after the next $n$ steps we will find that $F_{2n+1} = F_{n+1}^2$ and in general $F_{kn+1} = F_{n+1}^k$. Thus, the series will next exhibit the pair $\ldots, 0, 1,\ldots$ at the first value of $k$, $k=1, 2, \ldots$, for which $F_{n+1}^k = 1$. This value of $k$ is computed by Multiplicative Order in Mathematica. The order of the series is then $n k$.

This reduces the problem to finding $n$; that is, we need to find the smallest positive $n$ for which

$$0 = F_n = \frac{\varphi_+^n - \varphi_-^n}{\varphi_+ - \varphi_-}.$$

By clearing the denominator this is equivalent to the equation $\varphi_+^n - \varphi_-^n=0$. Dividing by $\varphi_-^n$ and remembering the product of the roots is $-1$, we obtain

$$(-\varphi_+^2)^n = 1.$$

Once again, the solution $n$ is found using MultiplicativeOrder. Here is the full code for this case (using variables s and t for $\varphi_+$ and $\varphi_-$):

order[p_?PrimeQ] /; Abs[Mod[p, 5, -2]] == 1 := Module[{x, s, t, a, n},
{s, t} = x /. Solve[x^2 - x - 1 == 0, x,  Modulus -> p];
n = MultiplicativeOrder[-s^2, p] ;
{a} = x /. Solve[PowerMod[s, n + 1, p] - PowerMod[t, n + 1, p] == (s - t) x, x, Modulus -> p];
n MultiplicativeOrder[a, p]
]


The case where the roots are "imaginary"

The extension field $\mathbb{F}_{p^2}$ consists of all formal linear combinations $i + j\varphi$ where $i$ and $j$ are in $\mathbb{F}_p$ and $\varphi$ can be taken to be either of $\varphi_+$ or $\varphi_-$. To multiply two such elements, use the defining property $\varphi^2 = \varphi + 1$ to compute

$$(a + b\varphi)(c + d\varphi) = a c + (a d + b c)\varphi + b d \varphi^2 = a c + (a d + b c)\varphi + b d (\varphi +1)\\ = (a c + b d) + (a d + b c + b d)\varphi.$$

Watch what happens when we take powers of $\varphi$ itself. I will show the case $p=3$. Starting with the first power we obtain the sequence

$$\varphi, 1+\varphi, 1 + 2\varphi, 2, 2\varphi, 2 + 2\varphi, 2 + \varphi, 1, \varphi, \ldots.$$

It's clear--and easily proven--that the $n$th term is precisely $F_n + F_{n+1}\varphi$. The order of the series therefore is the smallest value of $n$ for which $\varphi^n=1$.

Unfortunately, Mathematica does not appear to implement an efficient method to compute this value. (It is the analog of MultiplicativeOrder for $\mathbb{F}_{p^2}$. Documentation for the Finite Fields package suggests this can be obtained with FieldInd, but it appears that the values have to be precalculated for every element of the field, which is horribly slow.) We will have to write our own version. At this juncture it helps to use a deep (but elementary) theorem from finite field theory: the multiplicative group of any finite field is cyclic. Because this group includes all but $0$, its order in $\mathbb{F}_{p^2}$ is $p^2-1$. Therefore, $n$ must be a divisor of $p^2-1$. By searching only among those divisors we obtain an efficient algorithm (at least for large $p$ and any $p$ where $p^2-1$ does not have a lot of divisors). Here is the code for this case. The principal challenge was to exploit the computation of powers of small divisors to expedite their computation for larger divisors; this is done by means of a memoized recursion which relies on Divisors returning a sorted list of values.

order[p_?PrimeQ] := Module[{power, divisors,  gf},
Times[gf[a_, b_], gf[c_, d_]] ^:= gf @@ Mod[{a c + b d, b c + a d + b d}, p];
power = gf[0, 1];
divisors = Divisors[p^2 - 1];
power[n_] := power[n] = With[{m = Max[Select[divisors, # < n &]]}, power[m] power[n - m]];
For[i = 2, i <= Length[divisors], i++,
If[power[divisors[[i]]] == gf[1, 0], Return[divisors[[i]]]]];
]


### Results

As a check of the analysis, I compared the values of order for the first $500$ primes to the values of Mr.Wizard's fibMod, which conducts a brute-force search (and therefore appears to be quite trustworthy, as well as remarkably speedy in its own right). No differences were found. For the record, here is the version of fibMod I am using for comparison (it has not been compiled or otherwise optimized):

fibMod[0, n_, p_] := i
fibMod[a_, b_, m_] := fibMod[i++; b, Mod[a + b, m], m]
fibMod[m_] := Block[{i = 1}, fibMod[1, 1, m]]


But was this mathematical work worth while? Well, the demands on iteration and RAM appear less--but that wasn't a terrible restriction in the first place. Let's consider timing. First the direct search:

Timing[Table[fibMod[p], {p, Prime[Range[501, 1000]]}];]


{7.909, Null}

Now our refined implementation:

Timing[Table[order[p], {p, Prime[Range[501, 1000]]}];]


{2.464, Null}

OK, that's a little faster. But these are smallish primes: let's compare the timings for larger primes.

fibMod[Prime[10^6]] // AbsoluteTiming


{71.3233253, 30971728}

order[Prime[10^6]] // AbsoluteTiming


{0.1404003, 30971728}

We have accomplished something: order can compute values that are inaccessible to fibMod. Here, for instance, are ten calculations for ten consecutive seven-digit primes:

AbsoluteTiming[order[#]] & /@ Prime[Range[10^7, 10^7 + 9]]  // TableForm


The output times range from $0$ to $1.83$ seconds each, totaling less than $2.5$ seconds.

The speed of order varies tremendously depending on whether $5$ is a square modulo $p$ and, if not, on the number of divisors of $p^2-1$. Sometimes it can be astonishingly fast while for primes of similar magnitude it may take a long time.

My implementation of the version of MultiplicativeOrder for $\mathbb{F}_{p^2}$ is inefficient, because it requires factoring $p^2-1$ and a rather unintelligent search. Of course we can obtain that factorization from those of $p-1$ and $p+1$ and perhaps the search could be made smarter. There may be more efficient algorithms known by specialists (especially cryptographers).

It is straightforward to generalize this solution to a general second-order linear recurrence of the form $F_{n+1} = \alpha F_n + \beta F_{n-1}$ for integers $\alpha$ and $\beta$. Applying the Chinese Remainder Theorem shows that the order modulo a product of distinct primes is the least common multiple of the orders modulo the primes separately. However, extending the analysis to prime powers (and thence to any modulus at all) requires some additional ideas.

• Oh my... bookmarked for a later read! Jan 27 '13 at 18:09
• Yes, I am sorry about the length @Mr.Wizard. I did not think it would be fair just to post the code, though, and hope that people would believe an unsupported claim that it works! Jan 27 '13 at 18:12
• Please don't be sorry; I love these "dissertations" that a few of you can write. I just need more time than I have now to read it. :-) Jan 27 '13 at 18:17
• Fun fact: among the first $10^5$ primes, $p=1158551$ has the shortest order ($170$) relative to its maximum possible value ($p-1$) among primes for which $5$ is a square (enabling very fast computation of the order): that's about half of all such primes. The calculation took $23$ seconds on four cores. Jan 27 '13 at 19:00
• Very nice. Info/word ratio is high, so length is fine. The thought/word demand on the reader is high, too -- very delightful solution. Thanks. Jan 29 '13 at 13:32

We can get a fast solution assuming no Wall-Sun-Sun primes exist.

Following the outlines here and here:

PisanoPeriod[n_Integer?Positive] := LCM @@ pisanoPeriod @@@ FactorInteger[n]

pisanoPeriod[p_, k_] := pisanoPeriod[p]p^(k-1)

pisanoPeriod = 3;
pisanoPeriod = 20;

pisanoPeriod[p_] /; Abs[Mod[p, 10, -1]] == 1 :=
MultiplicativeOrder[PowerMod[2, -1, p](1 + PowerMod[5, 1/2, p]), p]

(* brute force on a small number of guesses *)
pisanoPeriod[p_] :=
Catch @ Block[{i = 1, divs = Prepend[4 Divisors[(p+1)/2], 1], fibPower},

fibPower[n_] := fibPower[n] = AlgebraMatrixPowerMod[{{1, 1}, {1, 0}}, n, p];

Fold[
With[{m = Mod[#1.fibPower[#2], p]},
i++;
If[m === {{1, 0}, {0, 1}},
Throw[divs[[i]]]
];
m
]&,
{{1, 1}, {1, 0}},
Differences[divs]
];

]


### Tests from other posts:

PisanoPeriod[Prime[10^6]] // AbsoluteTiming

{0.00099, 30971728}

AbsoluteTiming[PisanoPeriod[#]] & /@ Prime[Range[10^7, 10^7 + 9]] // TableForm

0.001150    358849348
0.000079    179424690
0.000153    358849396
0.000050    89712359
0.000044    179424730
0.000636    358849488
0.000048    179424778
0.000151    358849576
0.000251    358849588
0.000455    358849596


### Some large numbers:

BlockRandom[
SeedRandom;
PisanoPeriod[RandomInteger[{1, 10^14}]] // AbsoluteTiming
]

{0.00423, 176901917708180}

BlockRandom[
SeedRandom;
PisanoPeriod[RandomInteger[{1, 10^30}]] // AbsoluteTiming
]

{0.004367, 789393120074808363816808176}

BlockRandom[
SeedRandom;
PisanoPeriod[RandomInteger[{1, 10^50}]] // AbsoluteTiming
]

{0.12467, 161707870628752674660608903834667515301268480}
`