**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:

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$.

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[5] = 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[1] = 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.

### Comments

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.

1I took at shot at what I think you're after. If this is inadequate please explain why and I'll try again. – Mr.Wizard – 2013-01-26T02:59:14.970

1

Useful references: http://www.math.toronto.edu/mccann/assignments/199S/FibModK.pdf (very elementary introduction), http://www.fq.math.ca/Scanned/1-4/jarden.pdf (original paper finding the periods modulo powers of $10$), http://en.wikipedia.org/wiki/Pisano_period (the Wikipedia article on these "Pisano periods"), and http://www.math.temple.edu/~renault/fibonacci/fib.html (an elementary but full account of helpful information concerning the periods, based on a PhD thesis: includes a Web applet).

– whuber – 2013-01-27T19:47:59.800