Easier program for period of Fibonacci sequence modulo p
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.
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[19]
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}
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[2] = 3;
pisanoPeriod[5] = 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] = Algebra`MatrixPowerMod[{{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[1];
PisanoPeriod[RandomInteger[{1, 10^14}]] // AbsoluteTiming
]
{0.00423, 176901917708180}
BlockRandom[
SeedRandom[1];
PisanoPeriod[RandomInteger[{1, 10^30}]] // AbsoluteTiming
]
{0.004367, 789393120074808363816808176}
BlockRandom[
SeedRandom[1];
PisanoPeriod[RandomInteger[{1, 10^50}]] // AbsoluteTiming
]
{0.12467, 161707870628752674660608903834667515301268480}