The Fibonacci sequence modulo p 2 -- An investigation by computer for p< 10 14
aa r X i v : . [ m a t h . N T ] J un THE FIBONACCI SEQUENCE MODULO p –AN INVESTIGATION BY COMPUTER FOR p < ANDREAS-STEPHAN ELSENHANS AND J ¨ORG JAHNEL
Abstract.
We show that for primes p < the period length κ ( p ) of theFibonacci sequence modulo p is never equal to its period length modulo p .The investigation involves an extensive search by computer. As an application,we establish the general formula κ ( p n ) = κ ( p ) · p n − for all primes less than 10 . Introduction
The Fibonacci sequence { F k } k ≥ is defined recursively by F = 0, F = 1,and F k = F k − + F k − for k ≥
2. Modulo some integer l ≥
2, it must ultimatelybecome periodic as there are only l different pairs of residues modulo l . There is nopre-period since the recursion may be reversed to F k − = F k − F k − . The minimalperiod κ ( l ) of the Fibonacci sequence modulo l is often called the Wall number asits main properties were discovered by D. D. Wall [Wa].Wall’s results may be summarized by the theorem below. It shows, in particular,that κ ( l ) is in general a lot smaller than l . In fact, one always has κ ( l ) ≤ l whereas equality holds if and only if l = 2 · n for some n ≥ Theorem 1.2 (Wall) . a) If gcd( l , l ) = 1 then κ ( l l ) = lcm( κ ( l ) , κ ( l )) .In particular, if l = Q Ni =1 p n i i where the p i are pairwise different prime numbersthen κ ( l ) = lcm( κ ( p n ) , . . . , κ ( p n N N )) .It is therefore sufficient to understand κ on prime powers. b) κ (2) = 3 and κ (5) = 20 . Otherwise, • if p is a prime such that p ≡ ± then κ ( p ) | ( p − . • If p is a prime such that p ≡ ± then κ ( p ) | (2 p + 2) but κ ( p ) ∤ ( p + 1) . Mathematics Subject Classification.
Primary 11-04, 11B39. Secondary 11Y55, 11A41.
Key words and phrases.
Fibonacci sequence. Wall number. Period length. Prime power. Mont-gomery representation. Wieferich problem.While this work was done the first author was supported in part by a Doctoral Fellowship ofthe Deutsche Forschungsgemeinschaft (DFG).The computer part of this work was executed on the Linux PCs of the Gauß Laboratory forScientific Computing at the G¨ottingen Mathematical Institute. Both authors are grateful toProf. Y. Tschinkel for the permission to use these machines as well as to the system adminis-trators for their support.
ANDREAS-STEPHAN ELSENHANS AND J ¨ORG JAHNEL c) If l ≥ then κ ( l ) is even. d) If p is prime, e ≥ , and p e | F κ ( p ) but p e +1 ∤ F κ ( p ) then (1) κ ( p n ) = (cid:26) κ ( p ) for n ≤ e,κ ( p ) · p n − e for n > e. The Open Problems
The Period Length Modulo a Prime.2.1.1.
It is quite surprising that the Fibonacci sequence still keeps secrets. But thereare at least two of them.
Problem 2.1.2.
The first open problem is “What is the exact value of κ ( p )?”.Equivalently, one should understand precisely the behaviour of the quotient Q given by Q ( p ) := p − κ ( p ) for p ≡ ± Q ( p ) := p +1) κ ( p ) for p ≡ ± Q ( p ) in terms of p but, may be, that istoo optimistic. It is known that Q is unbounded. This is an elementary result due to D. Jar-den [Ja, Theorem 3].On the other hand, Q does not at all tend to infinity. If fact, in his unpublishedPh.D. thesis [G¨o], G. G¨ottsch computes a certain average value of Q . To be moreprecise, under the assumption of the Generalized Riemann Hypothesis, he proves X p ≡± p ≤ x, p prime Q ( p ) = C x log x + O (cid:16) x log log x log x (cid:17) where C = Q p prime (1 − pp − ) ≈ .
331 055 98. The proof shows as well that thedensity of { p prime | Q ( p ) = 1 , p ≡ ± } within the set of all primes isequal to C = Q p prime (1 − p ( p − ) ≈ .
265 705 45.Not assuming any hypothesis, it is still possible to verify that the right handside constitutes an upper bound. For that, the error term needs to be weakenedto O ( x log log log x log x log log x ).For the case p ≡ ± X p ≡± p ≡ p ≤ x, p prime Q ( p ) ≤ C x log x + O (cid:16) x log log log x log x · log log x (cid:17) where C = Q p prime ,p =2 , (1 − pp − ) ≈ .
210 055 99. The density of the set { p prime | Q ( p ) = 1 , p ≡ ± , p ≡ } within the set of all primesis at most C = Q p prime ,p =2 , (1 − p ( p − ) ≈ .
196 818 85.
HE FIBONACCI SEQUENCE MODULO p It seems, however, that the inequalities could well be equalities. In addition,the restriction to primes satisfying p ≡ p < · by computer. Up tothat bound, there are 317 687 prime numbers such that p ≡ ± p ≡ Q ( p ) = 1 exactly 250 246 times which is a relativefrequency of 0 .
787 712 434 . . . = 4 · .
196 928 108 . . . .On the other hand, there are 317 747 primes p satisfying p ≡ ± p ≡ Q ( p ) = 1 occurs 250 353 times which is basically thesame frequency as in the case p ≡ The Period Length Modulo a Prime Power.Problem 2.2.1.
There is another open problem. In fact, there is one questionwhich was left open in the formulation of Theorem 1.2: What is the exact valueof e in dependence of p ? Experiments for small p show that e = 1. Is this alwaysthe case? In other words, does one always have(2) κ ( p n ) = κ ( p ) · p n − similarly to the famous formula for Euler’s ϕ function?This is the most perplexing point in D. D. Wall’s whole study of the Fibonaccisequence modulo m . For p < , it was investigated by help of an electroniccomputer by Wall in 1960, already. We continued Wall’s investigation concerning Problem 2.2.1 on today’s ma-chines. Our main result is Theorem 4.4 below. The purpose of the present article isto give a description of our approach, particularly of the various algorithms devel-oped and optimizations used.
Definition 2.2.3.
We call a prime number p exceptional if equation (2) is wrongfor some n ≥ Fundamental Lemma 2.2.4.
Assume p = 2 and l to be a multiple of κ ( p ) .Then, p e | F l is sufficient for l being a period of { F k mod p e } k ≥ , i.e. for κ ( p e ) | l . Proof.
The claim is that, in our situation, F l +1 ≡ p e ) is automatic.For that, we note that there is the standard formula F l +1 F l − − F l = ( − l = 1which we explain in (5) below. Here, by assumption, F l ≡ p e ) and, by virtueof the recursion, F l − ≡ F l +1 (mod p e ). Therefore, F l +1 ≡ p e ).On the other hand, the condition κ ( p ) | l implies that F l +1 ≡ p ).As p = 2, Hensel’s lemma says that the lift is unique. This shows F l +1 ≡ p e )which is our claim. (cid:3) Proposition 2.2.5.
Let p be a prime number. Then, the following assertionsare equivalent. ANDREAS-STEPHAN ELSENHANS AND J ¨ORG JAHNEL i) p is exceptional, ii) F κ ( p ) is divisible by p . Proof. “ i) = ⇒ ii)” Assume, to the contrary, that p ∤ F κ ( p ) . By definition of κ ( p ),we know for sure that nevertheless p | F κ ( p ) . Together, these statements mean, The-orem 1.2.d) may be applied for e = 1 showing κ ( p n ) = κ ( p ) · p n − for every n ∈ N .This contradicts i).“ ii) = ⇒ i)” We choose the maximal e ∈ N such that p e | F κ ( p ) . By assumption, e ≥ κ ( p ) = κ ( p ) which shows equation (2) to be wrongfor n = 2. p is exceptional. (cid:3) Proposition 2.2.6.
Let p = 2 , be a prime number. I. If p ≡ ± then the following assertions are equivalent. i) p is exceptional, ii) F p − is divisible by p , iii) For r = √ ∈ Z /p Z one has r p − = 1 . II. If p ≡ ± then the following assertions are equivalent. i) p is exceptional, ii) F p +2 is divisible by p , iii) F p +1 is divisible by p . iv) In R p := Z /p Z [ r ] / ( r − r − one has r p +1 = − . Proof.
I. “ i) = ⇒ iii)” We put s = −√ and use formula (3) below. By Proposi-tion 2.2.5, F κ ( p ) is divisible by p . Therefore, r κ ( p ) = s κ ( p ) = ( r κ ( p ) ) − ∈ ( Z /p Z ) ∗ ,i.e. ( r κ ( p ) ) = 1. Since κ ( p ) | ( p − r p − ) = 1 from this.On the other hand, we know r p − ≡ p ) by Fermat’s Theorem. Unique-ness of Hensel’s lift implies r p − = 1.“ iii) = ⇒ ii)” We have r p − s p − = ( − p − = 1. Thus, r p − = 1 implies s p − = 1.Consequently, ( F p − mod p ) = r p − − s p − √ = 0 and F p − is divisible by p .“ ii) = ⇒ i)” As ( p −
1) is a multiple of κ ( p ), Lemma 2.2.4 may be applied. It shows κ ( p ) | ( p − n = 2. p is exceptional.II. “ i) = ⇒ ii)” By Proposition 2.2.5, F κ ( p ) is divisible by p . In that situation,Lemma 2.2.4 implies that κ ( p ) is actually a period of { F k mod p } k ≥ . By conse-quence, (2 p + 2) is a period of { F k mod p } k ≥ , too. This shows p | F p +2 .“ ii) = ⇒ iii)” Since F p +2 = F p +1 V p +1 , all we need is p ∤ V p +1 . This, however, is clearas V p +1 = r p +1 + r p +1 ≡ − p ).“ iii) = ⇒ iv)” The assumption implies r p +1 = r p +1 , i.e. r p +1 ∈ Z /p Z . As rr = − r p +1 ) = r p +1 r p +1 = ( − p +1 = 1 from this. Hensel’s lemmaimplies r p +1 = − r p +1 ≡ − p ) is known. HE FIBONACCI SEQUENCE MODULO p “ iv) = ⇒ i)” r p +1 = − F p +2 mod p ) = ( r p +1 ) − ( r p +1 ) √ = 0.As (2 p +2) is a multiple of κ ( p ), Lemma 2.2.4 may be applied. It shows κ ( p ) | (2 p +2).This contradicts equation (2) for n = 2. p is exceptional. (cid:3) Remark 2.2.7.
By Proposition 2.2.6, the problem of finding exceptional primes isin perfect analogy to the problem of finding
Wieferich primes .In the Wieferich case, one knows 2 p − ≡ p ) and would like to understandthe set of all primes for which even 2 p − ≡ p ) is valid. Here, we know F κ ( p ) ≡ p ) and look for the primes which fulfill F κ ( p ) ≡ p ).At least in the case p ≡ ± r . Remark 2.2.8.
One might want to put the concept of an exceptional prime intothe wider context of algebraic number theory. We work in the number field Q (cid:0) √ (cid:1) in which r = √ is a fundamental unit.By analogy, we could say that an odd prime number p is exceptional for the real qua-dratic number field K = Q (cid:0) √ d (cid:1) if, for ε a fundamental unit in K , ε p − ≡ p )when (cid:0) dp (cid:1) = 1, ε p +2 ≡ p ) when (cid:0) dp (cid:1) = −
1, or ε p ( p − ≡ p ) in theramified case d | p . A congruence modulo p is, of course, supposed to mean equalityin O K / ( p ) which, as p = 2, is isomorphic to Z [ √ d ] / ( p ) = Z /p Z [ X ] / ( X − d ).Note that there is no ambiguity coming from the choice of ε since all exponentsare even.For many real quadratic number fields it does not require sophisticated program-ming to find a few exceptional primes. Below, we give the complete list of all ex-ceptional primes p < for the fields Q (cid:0) √ d (cid:1) where d is square-free and up to 101.Thereby, primes put in parentheses are those such that Q (cid:0) √ d (cid:1) is ramified at p . d exceptional primes p d exceptional primes p d exceptional primes p2 13, 31, 1 546 463 17 34 37, 547, 4 7333 103 19 79, 1 271 731, 13 599 893, 31 352 389 35 23, 577, 1 325 6635 21 46 179 311 37 7, 89, 257, 6316 (3), 7, 523 22 43, 73, 409, 28 477 38 57 23 7, 733 39 5, 7, 37, 163 409, 795 490 66710 191, 643, 134 339, 25 233 137 26 2 683, 3 967, 18 587 41 29, 53, 7 21111 29 3, 11 42 (3), 5, 43, 7113 241 30 43 3, 47914 6 707 879, 93 140 353 31 157, 261 687 119 46 (23)15 (3), 181, 1 039, 2 917, 2 401 457 33 (3), 29, 37, 6 713 797 47 5 762 437d exceptional primes p d exceptional primes p d exceptional primes p51 (3), 5, 37, 4 831 67 3, 11, 953, 57 301 83 3, 19 699, 2 417 37753 5 69 (3), 5, 17, 52 469 057 85 3, 204 520 55955 571 70 (5), 59, 20 411 86 1 231, 5 77957 59, 28 927, 1 726 079, 7 480 159 71 67, 2 953, 8 863, 522 647 821 87 (3), 17, 757, 1 12358 3, 23, 4 639, 172 721, 16 557 419 73 5, 7, 41, 3 947, 6 079 89 5, 7, 13, 5959 1 559, 17 385 737 74 3, 7, 1 171 91 (13), 1 218 69161 77 3, 418 270 987 93 (3), 1362 3, 5, 263, 388 897 78 (3), 19, 62 591 94 7365 1 327, 8 831, 569 831 79 3, 113, 4 049, 6 199 95 6 257, 10 93766 21 023, 106 107 779 82 3, 5, 11, 769, 3 256 531, 624 451 181 97 17, 3 331101 7, 19 301 ANDREAS-STEPHAN ELSENHANS AND J ¨ORG JAHNEL
Among these 158 exceptional primes, there are exactly nine for which even thestronger congruence modulo p is true. These are p = 3 for d = 29, 42, 67, and 74, p = 5 for d = 62, 73, and 89, p = 17 for d = 69, and p = 29 for d = 41. We do notobserve a congruence of the type above modulo p .3. Background
Part a) of Wall’s theorem is trivial.For the proof of b), Binet’s formula(3) F k = r k − s k √ , where r = √ and s = −√ , is of fundamental importance. It is easily establishedby induction. If p ≡ ± p and,therefore, ±√ ∈ F p . Fermat states their order is a divisor of p − ±√ ∈ F p are elements of norm ( − N : F ∗ p → F ∗ p is surjective, its kernel is a group of order p − p − = p +1 and N − ( { , − } ) = 2 p +2.As F ∗ p is cyclic, we see that N − ( { , − } ) is even a cyclic group of order 2 p + 2. N ( r ) = N ( s ) = − r and s are not contained in its subgroup ofindex two. Therefore,(4) r p +1 ≡ s p +1 ≡ − p ) . From this, we find F p +2 ≡ r p +2 − s p +2 √ ≡ − r + s √ ≡ − F ≡ − p ) which shows p + 1is not a period of { F k } k ≥ modulo p .c) In the case p ≡ ± k ∈ N , one has F k +1 F k − − F k = r k + s k − r k +1 s k − − r k − s k +1 − r k + s k − r k s k − ( − k − ( r + s ) + 2( − k − k as rs = − r + s = 3. On the other hand, F κ ( l )+1 F κ ( l ) − − F κ ( l ) ≡ · − ≡ l ) . As l ≥ κ ( l ) is even.For d), it is best to establish the following p -uplication formula first. HE FIBONACCI SEQUENCE MODULO p Lemma 3.2 (Wall) . One has (6) F pk = 12 p − p X j =1 j odd (cid:18) pj (cid:19) j − F jk V p − jk . Here, { V k } k ≥ is the Lucas sequence given by V = 2 , V = 1 , and V k = V k − + V k − for k ≥ . Proof.
Induction shows V k = r k + s k . Having that in mind, it is easy to calculateas follows. F pk = ( r k ) p − ( s k ) p √ V k + √ F k ) p − ( V k −√ F k ) p √ . The assertion follows from the Binomial Theorem. (cid:3)
The fundamental Lemma 2.2.4 allows us to prove d) for p = 2 in a somewhatsimpler manner than D. D. Wall did it in [Wa].First, we note that for n ≤ e , 2.2.4 implies κ ( p n ) | κ ( p ). However, divisibility theother way round is obvious.For n ≥ e , by Lemma 2.2.4, it is sufficient to prove ν p ( F κ ( p ) · p n − e ) = n , i.e. that p n | F κ ( p ) · p n − e but p n +1 ∤ F κ ( p ) · p n − e . Indeed, the first divisibility implies κ ( p n ) | κ ( p ) · p n − e while the second, applied for n − n , yields κ ( p n ) ∤ κ ( p ) · p n − e − . The resultfollows as κ ( p ) | κ ( p n ).For ν p ( F κ ( p ) · p n − e ) = n , we proceed by induction, the case n = e being known byassumption. One has F κ ( p ) · p n − e +1 = p − pF κ ( p ) · p n − e V p − κ ( p ) · p n − e + p − p P j =3 j odd (cid:0) pj (cid:1) j − F jκ ( p ) · p n − e V p − jκ ( p ) · p n − e . In the second term, every summand is divisible by F κ ( p ) · p n − e , i.e. by p n . The claimwould follow if we knew p ∤ V κ ( p ) · p n − e . This, however, is easy as there is the formula(7) V l = F l − + F l +1 which implies V l ≡ p ) for l any multiple of κ ( p ). For p = 2, as always, things are a bit more complicated. We still have κ (2 n ) = 3 · n − . However, for n ≥
2, one has 2 n +1 | F · n − for which there is noanalogue in the p = 2 case. On the other hand, ν ( F · n − +1 −
1) = n which issufficient for our assertion.The duplication formula provided by Lemma 3.2 is(8) F k = F k V k = F k ( F k − + F k +1 ) = F k + 2 F k F k − . As F = 8, a repeated application of this formula shows 2 n +1 | F · n − for every n ≥ ANDREAS-STEPHAN ELSENHANS AND J ¨ORG JAHNEL
We further claim F k +1 = F k + F k +1 . Indeed, this is true for k = 0 as 1 = 0 + 1 and we proceed by induction as follows: F k +3 = F k +1 + F k +2 = F k + F k +1 + F k +1 + 2 F k +1 F k == F k +1 + ( F k + F k +1 ) = F k +1 + F k +2 . (9)The assertion ν ( F · n − +1 −
1) = n is now easily established by induction. We notethat F = 13 ≡ F · n +1 = F · n − + F · n − +1 where the first summand is even divisible by 2 n +2 .The second one is congruent to 1 modulo 2 n +1 , but not modulo 2 n +2 , by consequenceof the induction hypothesis.4. A heuristic argument
We expect that there are infinitely many exceptional primes for Q (cid:0) √ (cid:1) .Our reasoning for this is as follows. p | F κ ( p ) is known by definition of κ ( p ). Thus,for any individual prime p , ( F κ ( p ) mod p ) is one residue out of p possibilities. If wewere allowed to assume equidistribution then we could conclude that p | F κ ( p ) shouldoccur with a “probability” of p . Further, by [RS, Theorem 5],log log N + A −
12 log N ≤ X p prime p ≤ N p ≤ log log N + A + 12 log N , at least for N ≥ A ∈ R is Mertens’ constant which is given by A = γ + X p prime (cid:20) p + log (cid:16) − p (cid:17)(cid:21) = 0 .
261 497 212 847 642 783 755 . . . whereas γ denotes the Euler-Mascheroni constant.This means that one should expect around log log N + A exceptional primes lessthan N . On the other hand, p | F κ ( p ) should occur only a few times or even not at all.Indeed, if we assume equidistribution again, then for any individual prime p , p | F κ ( p ) should happen with a “probability” of p . However, ∞ X p =2 p prime p = 0 .
452 247 420 041 065 498 506 . . . . is a convergent series.
Remark 4.3.
It is, may be, of interest that, for any exponent n ≥
2, one has theequality P p prime 1 p n = P ∞ k =1 µ ( k ) k log ζ ( nk ) where the right hand converges a lot fasterand may be used for evaluation. This equation results from the Moebius inversionformula and Euler’s formula log ζ ( nk ) = P p prime − log(1 − p nk ) = P ∞ j =1 1 j P p prime 1 p jnk . HE FIBONACCI SEQUENCE MODULO p We carried out an extensive search for exceptional primes but, unfortunately,we had no success and our result is negative.
Theorem.
There are no exceptional primes p < .Down the earth, this means that one has κ ( p n ) = κ ( p ) · p n − for every n ∈ N andall primes p < . Algorithms
We worked with two principally different types of algorithms. First, in the p ≡ ± r p − mod p ). A second and morecomplete approach is to compute ( F p − mod p ) in the p ≡ ± F p +2 mod p ) or ( F p +1 mod p ) in the case p ≡ ± Remark 5.0.2.
In the case p ≡ ± p = 2, exceptionality is equiva-lent to r p +2 ≡ p ). Unfortunately, an approach based on that observationturns out to be impractical as it involves the calculation of a modular power in R p = Z /p Z (cid:2) √ (cid:3) = Z (cid:2) √ (cid:3) / ( p ) in a situation where √ Z /p Z . In comparisonwith Z /p Z , multiplication in R p is a lot slower, at least in our (naive) implemen-tations. This puts a modular powering operation in R p out of competition with adirect approach to compute F p +2 (or F p +1 ) modulo p .5.1. Algorithms based on the computation of √ .5.1.1. If p ≡ ± r p − mod p ). Thealgorithm should consist of four steps.i) Compute the square root of 5 in Z /p Z .ii) Take the Hensel’s lift of this root to Z /p Z .iii) Calculate the golden ratio r := √ ∈ Z /p Z .iv) Use a modular powering operation to find ( r p − mod p ).We call algorithms which follow this strategy algorithms powering the golden ratio. Here, the final steps iii) and iv) are not critical at all. For iii), it is obvious thatthis is a simple calculation while for iv), carefully optimized modular poweringoperations are available. Further, ii) can be effectively done as r ≡ p )implies w := r − r − p · ( r mod p ) · p is a square root of 5 modulo p . Thus, themost expensive operation should be a run of Euclid’s extended algorithm in orderto find ( r mod p ).In fact, there is a way to avoid even this. We first calculate ∈ F p . This iseasier than an arbitrary division in residues modulo p . We may put := p +15 if p ≡ := p +15 if p ≡ − v of ( r mod p ) can be computed as v = r · . We get away with one integer divisionand one multiplication. Thus, the most interesting point is i), the computation of √ ∈ F p . In gen-eral, there is a beautiful algorithm to find square roots modulo a prime number dueto Shanks [Co, Algorithm 1.5.1]. We implemented this algorithm but let it finallyrun only in the p ≡ p p ≡ w := (5 p +14 mod p ) to find a square rootof 5 by one modular powering operation.If p ≡ w := (5 p +38 mod p )as long as 5 p − ≡ p ) and(11) w := (10 · p − mod p )if 5 p − ≡ − p ). Note that 5 is a quadratic residue modulo p . Hence, wealways have 5 p − ≡ ± p ).For sure, (5 p − mod p ) can be computed using a modular powering operation.In fact, we implemented an algorithm doing that and let it run through the in-terval [10 , · ].However, (5 p − mod p ) is nothing but a quartic residue symbol. For that reason,there is an actually faster algorithm which we obtained by an approach using thelaw of biquadratic reciprocity. Proposition 5.1.3.
Let p be a prime number such that p ≡ and p ≡ ± and let p = a + b be its (essentially unique) decomposition intoa sum of two squares. a) Then, a and b may be normalized such that a ≡ and b is even. b) Assume a and b are normalized as described in a). Then, there are only thefollowing eight possibilities. i) a ≡ , , , or
19 (mod 20) and b ≡
10 (mod 20) .In this case, p − ≡ p ) , i.e. is a quartic residue modulo p . ii) a ≡
15 (mod 20) and b ≡ , , , or
18 (mod 20) .Here, p − ≡ − p ) , i.e. is a quadratic but not a quartic residue modulo p . Proof. a) As p is odd, among the integers a and b there must be an even and anodd one. We choose b to be even and force a ≡ a by ( − a ),if necessary.b) We first observe that a ≡ b ≡ b ≡ a and b must be divisible by 5. Indeed,otherwise we had a , b ≡ ± a + b ≡ ± a and b cannot be both divisible by 5. HE FIBONACCI SEQUENCE MODULO p If a is divisible by 5 then a ≡ a ≡
15 (mod 20). b ≡ b not divisible by 5 yield the four possibilities stated. On the other hand, if b isdivisible by 5 then b ≡ b ≡
10 (mod 20). a ≡ a not divisible by 5 show there are precisely the four possibilities listed.For the remaining assertions, we first note that (5 p − mod p ) tests whether x ≡ p ) has a solution x ∈ Z , i.e. whether 5 is a quartic residue modulo p .By [IR, Lemma 9.10.1], we know(5 p − mod p ) = χ a + bi (5)where χ denotes the quartic residue symbol. The law of biquadratic reciprocity[IR, Theorem 9.2] asserts χ a + bi (5) = χ ( a + bi ) . For that, we note explicitly that a + bi ≡ i (mod 4), 5 ≡ N (5) − = 6 is even. Let us now compute χ ( a + bi ): χ ( a + bi ) = χ − i ( a + bi ) · χ − − i ( a + bi )= χ − − i ( a − bi ) · χ − − i ( a + bi )= (cid:16) a + b (cid:17) · χ − − i ( a − bi ) · χ − − i ( a + bi )= (cid:16) a + b (cid:17) · χ − − i ( p ) . Here, the first equation is the definition of the quartic residue symbol for compositeelements while the second is [IR, Proposition 9.8.3.c)].For the third equation, we observe that χ − − i ( a − bi ) is either ± ± i . By sim-ply omitting the complex conjugation, we would make a sign error if and onlyif χ − − i ( a − bi ) = ± i . By [IR, Lemma 9.10.1], this means exactly that a − bi defines, under the identification 2 i = −
1, not even a quadratic residue modulo 5.Therefore, the correction factor is ( a + b ). The final equation follows from [IR, Propo-sition 9.8.3.b)].We note that, by virtue of [IR, Lemma 9.10.1], χ − − i ( p ) tests whether p is a quarticresidue modulo 5 or not. As p is for sure a quadratic residue, we may write χ − − i ( p ) = (cid:26) p ≡ , − p ≡ − χ − − i ( p ) = ( p mod 5).The eight possibilities could now be inspected one after the other. A more conceptualargument works as follows. In case i), we have (cid:16) a + b (cid:17) = (cid:16) a (cid:17) = ( a mod 5) = ( a + b mod 5) = ( p mod 5) . Therefore, (5 p − mod p ) = 1. On the other hand, in case ii), (cid:16) a + b (cid:17) = (cid:16) b (cid:17) = (cid:16) b (cid:17) = ( − b mod 5) = ( − a − b mod 5) == − ( p mod 5) . Hence, (5 p − mod p ) = − (cid:3) Although we are not going to make use of it, let us state the complementaryresult for p ≡ Proposition.
Let p be a prime such that p ≡ and p ≡ ± andlet p = a + b be its (essentially unique) decomposition into a sum of two squares. a) Then, a and b may be normalized such that a ≡ and b is even. b) Assume a and b are normalized as described in a). Then, there are only thefollowing eight possibilities. i) a ≡ , , , or
17 (mod 20) and b ≡ .In this case, p − ≡ p ) , i.e. is a quartic residue modulo p . ii) a ≡ and b ≡ , , , or
16 (mod 20) .Here, p − ≡ − p ) , i.e. is a quadratic but not a quartic residue modulo p . Proof. a) As p is odd, among the integers a and b there must be an even and anodd one. We choose b to be even and force a ≡ a by ( − a ),if necessary.b) We first observe that a ≡ b ≡ | b . Then, werealize that one of the two numbers a and b must be divisible by 5. Indeed, otherwisewe had a , b ≡ ± a + b ≡ ± a and b cannot be both divisible by 5.If a is divisible by 5 then a ≡ a ≡ | b and b notdivisible by 5 yield the four possibilities stated. On the other hand, if b is divisibleby 5 then 4 | b implies b ≡ a ≡ a not divisible by 5 showthere are precisely the four possibilities listed.The proof of the remaining assertions works exactly in the same way as the proofof Proposition 5.1.3 above. We note explicitly that a + bi ≡ (cid:3) As the transformation a
7→ − a does not affect any of the three statementsbelow, we may formulate the following theorem. Actually, this is the result we needfor the application. HE FIBONACCI SEQUENCE MODULO p Theorem.
Let p be a prime number such that p ≡ and p ≡ ± and let p = a + b be its decomposition into a sum of two squares. We normalize a and b such that a is odd and b is even. Then, the following three statementsare equivalent. i) 5 is a quartic residue modulo p . ii) b is divisible by . iii) a is not divisible by . Remark 5.1.6.
We note that the restrictions on p exclude only trivial cases.If p
6≡ ± p . If p ≡ Algorithm 5.1.7.
The square sum sieve algorithm for prime numbers p such that p ≡ ,
29 (mod 40) runs as follows.We investigate a rectangle [ N , N ] × [ M , M ] of numbers. We will go through therectangle row-by-row in the same way as the electron beam goes through a screen.a) We add 0, 1, 2, or 3 to M to make sure M ≡ b go from M to M in steps of length four.b) For a fixed b we sieve the odd numbers in the interval [ N , N ].Except for the odd case that l | a, b which we decided to ignore as the density of thesepairs is not too high, l | a + b implies that ( −
1) is a quadratic residue modulo l ,i.e. we need to sieve only by the primes l ≡ l which is below a certain limit we cross out all those a such that a ≡ ± v l b (mod l ). Here, v l is a square root of ( −
1) modulo l , i.e. v l ≡ − l ).For practical application, this requires that the square roots of ( −
1) modulo therelevant primes have to be pre-computed and stored in an array once and for all.c) For the remaining pairs ( a, b ), we compute p = a + b and do steps i) through iv)from 5.1.1. In step i), if b is divisible by 5 then we use formula (10) to compute thesquare root of 5 modulo p . Otherwise, we use formula (11). In practice, we ran the square sum sieve algorithm on the rectangles[0 , × [1 580 000 , , × [0 , p ∈ [5 · , . · ] such that p ≡ ,
29 (mod 40) plusseveral others.In fact, on the second rectangle we ran a modified version, the inverted square sumsieve, where the two outer loops are reversed. That means, we let a go through theodd numbers in [ N , N ] in the very outer loop. This has some advantage in speedas longer intervals are sieved at once. In other words, we go through the rectanglecolumn-by-column.We implemented the square sum sieve algorithms in C using the mpz functions ofGNU’s GMP package for arithmetic on long integers. On a single 1211 MHz Athlon processor, the computations for the first rectangle took around 22 days of CPU time.The computations for the smaller second rectangle were finished after nine days. For primes p such that p ≡ p ≡ ± w := (5 p +14 mod p ) for the square root of 5 makes things a lot easier. Instead ofthe square sum sieve we implemented the sieve of Eratosthenes. Caused by thelimitations of main memory in today’s PCs, we could actually sieve intervals of onlyabout 250 000 000 numbers at once. For each such interval the remainders of itsstarting point have to be computed (painfully) by explicit divisions. Algorithm 5.1.10.
More precisely, the algorithm powering the golden ratio forprimes p ≡ ,
19 (mod 20) runs as follows.We investigate an interval [ N , N ]. We assume that N − N is divisible by 5 · and that N is divisible by 20.a) We let an integer variable i count from 0 to N − N · − i we work on the interval I = [ N + 5 · · i, N + 5 · · ( i + 1)].For each prime l which is below a certain limit, we compute ( N + 5 · · i mod l ).Then, we cross out all p ∈ I , p ≡
11 (or 19) mod 20 which are divisible by l .c) For the remaining p ∈ I , p ≡
11 (or 19) mod 20 we do steps i) through iv)from 5.1.1. In step i), we use the formula w := (5 p +14 mod p ) to compute the squareroot of 5 modulo p . In practice, we ran this algorithm in order to test all prime numbers p ∈ [10 , · ] such that p ≡
11 (mod 20) or p ≡
19 (mod 20). It was im-plemented in C using the mpz functions of the GMP package.Later, when testing primes above 10 , we used the low level mpn functions forlong natural numbers. In particular, we implemented a modular powering functionwhich is hand-tailored for numbers of the considered size. It uses the left-rightbase 2 powering algorithm [Co, Algorithm 1.2.3] and the sliding window improve-ment from mpz powm.Having done all these optimizations, work on the test interval [4 · , · +5 · ]of 250 000 000 numbers p such that p ≡
11 (mod 20), among them 19 955 067 primes,lasted 7:50 Minutes CPU time on a 1211 MHz Athlon processor.
Similarly, for prime numbers p satisfying the simultaneous congruences p ≡ p ≡ ± p . HE FIBONACCI SEQUENCE MODULO p Algorithm 5.1.13.
More precisely, the algorithm powering the golden ratio forprimes p ≡ , N , N ]. We assume that N − N is divisible by 10 and that N is divisible by 40.a) We let an integer variable i count from 0 to N − N − i we work on the interval I = [ N + 10 · i, N + 10 · ( i + 1)]. For eachprime l which is below a certain limit, we compute (( N + 10 · i ) mod l ). Then,we cross out all p ∈ I , p ≡ l .c) For the remaining p ∈ I , p ≡ p . We ran this algorithm on the interval [10 , · ]. After all optimiza-tions, the test interval [4 · , · + 10 ] of 250 000 000 numbers p such that p ≡ p ≡
11 (mod 20)or p ≡
19 (mod 20). The difference comes entirely from the more complicatedprocedure to compute √ ∈ F p . Remark 5.1.15.
At a certain moment, such a running time was no longer found rea-sonable. A direct computation of the Fibonacci numbers could be done as well. Afterseveral optimizations of the code of the direct methods, it turned out that only the3 mod 4 case could still compete with them. We discuss the direct methods in thesubsection below.5.2.
Algorithms for a direct computation of Fibonacci numbers.Algorithm 5.2.1.
A nice algorithm for the fast computation of a Fibonacci numberis presented in O. Forster’s book [Fo]. It is based on the formulae F k − = F k + F k − ,F k = F k + 2 F k F k − . (12)and works in the spirit of the left-right binary powering algorithm using bits.Our adaption uses modular operations modulo p instead of integer operations.An implementation in O. Forster’s Pascal-style multi precision interpreter languageARIBAS looks like this. (*------------------------------------------------------------------*)(*** Schnelle Berechnung der Fibonacci-Zahlen mittels der Formeln** fib(2*k-1) = fib(k)**2 + fib(k-1)**2** fib(2*k) = fib(k)**2 + 2*fib(k)*fib(k-1)**** Dabei werden alle Berechnungen mod m durchgef¨uhrt*)function fib(k,m : integer): integer;var b, x, y, xx, temp: integer;beginif k <= 1 then return k end;x := 1; y := 0;for b := bit_length(k)-2 to 0 by -1 doxx := x*x mod m;x := (xx + 2*x*y) mod m;y := (xx + y*y) mod m;if bit_test(k,b) thentemp := x;x := (x + y) mod m;y := temp;end;end;return x;end.(** ein systematischer Versuch**)function test() : integervar p,r,r1 : integer;ptest : boolean;beginfor p := 90000000001 to 95000000001 by 2 doif (p mod 10000) = 1 thenwriteln("getestete Zahl: ", p);end;ptest := rab_primetest(p);if (ptest = true) thenif ((p mod 5 = 2) or (p mod 5 = 3)) thenr := fib(2*p+2,p*p);elser := fib(p-1,p*p);end;if (r <= 30000000000000000) thenr1 := r div p;writeln(p," ist eine interessante Primzahl. Quotient ", r1);end;end;end;return(0);end. A call to fib(k,m) computes ( F k mod m ). test is the main function. test() executes an outer loop which contains a Rabin-Miller composedness test. For apseudo prime p , it uses the function fib to compute ( F p − mod p ) or ( F p +2 mod p ).As these are divisible by p we output the quotient instead. Note that in order to limitthe output size we actually write an output only when the quotient is rather small. ARIBAS is fast enough to ensure that this algorithm could be run from p = 7up to 10 . We worked on ten PCs in parallel for five days. That was our first biggercomputing project concerning this problem. It showed that no exceptional primes p < do exist, thereby a establishing a lightweight version of Theorem 4.4. HE FIBONACCI SEQUENCE MODULO p The running time made it clear that we had approached to the limits of aninterpreter language. For a systematic test of larger prime numbers, the algorithmwas ported to C. For the arithmetic on long integers we used the mpz functionsof GMP. After only one further optimization, the integration of a version of thesieve of Eratosthenes, the interval [10 , ] could be attacked. A test intervalof 250 000 000 numbers was dealt with on a 1211 MHz Athlon processor in around40 Minutes CPU time. Again, we did parallel computing on ten PCs. The searchthrough [10 , ] was finished in less than five days. For the interval [10 , ], the methods which compute √ ∈ F p and squarethe golden ratio were introduced as they were faster than our implementation ofO. Forster’s algorithm at that time. For this reason, only the case p ≡ ± Optimizations
Sieving.6.1.1.
Near 10 , one of about 32 numbers is prime. We work in a fixed primeresidue class modulo 10, 20, or 40 but still, only one of about 13 numbers is prime.We feel that the computations of ( F p ± mod p ) should take the main part of therunning time of our programs. Our goal is, therefore, to rapidly exclude (most of)the non-primes from the list and then to spend most of the time on the remain-ing numbers.There are various methods to generate the list of all primes within an interval. Un-fortunately, this section of our code is not as harmless as one could hope for. In fact,for an individual number p , one might have the idea to decide whether it is prob-ably prime by computing ( F p ± mod p ). That is the Fibonacci composedness test.It would, unfortunately, not reduce our computational load a lot as it is almostas complex as the main computation. This clearly indicates the problem that thestandard “pseudo primality tests” which are designed to test individual numbersare not well suited for our purposes. In this subsection, we will explain what we didinstead in order to speed up this part of the program. Our first programs in ARIBAS in fact used the internal primality test tocheck each number in the interval individually. At the ARIBAS level, this is optimalbecause it involves only one instruction for the interpreter.When we migrated our programs to C, using the GMP library, we first tried the same.We used the function mpz probab prime with one repetition for every number tobe tested. It turned out that this program was enormously inefficient. It tookabout 50 per cent of the running time for primality testing and 50 per cent for thecomputation of Fibonacci numbers. However, it could easily be tuned by a naiveimplementation of the sieve of Eratosthenes in intervals of length 1 000 000.
We first combined sieving by small primes and the mpz probab prime functionbecause sieving by huge primes is slow. This made sure that the computa-tion of Fibonacci numbers took the major part of the running time. However,mpz probab prime is not at all intended to be combined with a sieve. In fact, itchecks divisibility by small primes once more. Thus, an optimization of the code forthe Fibonacci numbers reversed the relation again. It became necessary to carry outa further optimization of the generation of the list of primes. We decided to aban-don all pseudo primality tests. Further, we enlarged the length of the array of up to250 000 000 numbers to minimize the number of initializations.In principle, the sieve works as follows. Recall that we used different algorithmsfor the computation of the Fibonacci numbers, depending on the residue class of p modulo 10, 20, or 40. This leads to a sieve in which the number S ( i ) := starting point + residue + modulus · i is represented by array position i . Since all our moduli are divisible by 2 and 5 wedo no longer sieve by these two numbers.Such a sieve is still easy to use. Given a prime p = 2 ,
5, one has to compute thearray index i of the first number which is divisible by p . Then, one can cross outthe numbers at the indices i , i + p, i + 2 p, . . . until the end of the sieve is reached. An array of the size above fitsinto the memory of today’s PCs but it does not fit into the cache. Thus, the speed-limiting part is the transfer between CPU and memory. Sieving by big primes is likea random access to single bytes. The memory manager has to transfer one block tothe cache memory, change one byte, and then transfer the whole block back to thememory. This is the limiting bottleneck.To avoid this problem as far as possible, we built a two stage sieve.
In the first stage, we sieve by the first 25 000, the “small”, primes. For that, we dividethe sieve further into segments of length 30 000. These two constants were found tobe optimal in practical tests. They are heavily machine dependent.The first stage is now easily explained. In a first step, we sieve the first segment by allsmall primes. Then, we sieve the second segment by all small primes. We continuein that way until the end of the sieve is reached.In the second stage, we work with all relevant “big” primes on the complete sieve,as usual.The result of this strategy is a sieve whose segments fit into the machine’s cache.Thus, the speed of the first sieve stage is the speed of the cache, not the speed ofthe memory. The speed of the second stage is limited by the initialization.On our machines the two stage sieve is twice as fast as the ordinary sieve.
HE FIBONACCI SEQUENCE MODULO p The choice of the prime limit for sieving is a point of interest, too. As wesearch for one very particular example, it would do no harm if, from to time, wetest a composite number p for p | F p ± . When the computer would tell us p divides F p ± which, in fact, it never did then it would be easy to do a reliable primality test.As long as we sieve by small primes, it is clear that lots of numbers will be crossedout in a short time and this will reduce the running time as it reduces the numberof times the actual computation of ( F p ± mod p ) is called. Afterwards, when wesieve by larger primes, the situation is no longer that clear. We will often crossout a number repeatedly which was crossed out already before. This means, it canhappen that further sieving costs actually more time than it saves.Our tests show nevertheless that it is best to sieve almost till to the square rootof the numbers to be tested. We introduced an automatic choice of the variableprime limit as √ p log √ p which means we sieve by the first [ √ p log √ p ] primes. Here, p meansthe first prime of the interval we want to go through. Another optimization was done by looking at the prime three. Every thirdnumber is crossed out when sieving by this prime and, when sieving by a big-ger prime, every third step hits a number which is divisible by three and alreadycrossed out.Thus, we can work more efficiently as follows. Let p be a prime bigger than three andcoprime to the modulus. We compute i , the first index of a number divisible by p .Then, we calculate the remainder of the corresponding number modulo three. If itis zero then we skip i and continue with i := i + p . Now, i corresponds to the firstnumber in the sieve which is divisible by p but not by three. Thus, we must cross out i , i + p, i + 3 p, i + 4 p, i + 6 p, . . . or i , i + 2 p, i + 3 p, i + 5 p, i + 6 p, . . . depending on whether i + 2 p corresponds to a number which is divisible by threeor not.6.2. The Montgomery Representation.6.2.1.
The algorithms for the computation of Fibonacci numbers modulo m ex-plained so far spend the lion’s share of their running time on the divisions by m whichoccur as the final steps of modular operations such as x := (xx + 2*x*y) mod m .Unfortunately, on today’s PC processors, divisions are by far slower than multipli-cations or even additions.An ingenious method to avoid most of the divisions in a modular powering operationis due to P. L. Montgomery [Mo]. We use an adaption of Montgomery’s method toO. Forster’s algorithm which works as follows.Let R be the smallest positive integer which does not fit into one machine word.That will normally be a power of two. On our machines, R = 2 . Recall that all op-erations on unsigned integers in C are automatically modular operations modulo R . We choose some exponent n such that the modulus m = p fulfills m ≤ R n . In oursituation p < , therefore m = p < < , such that n = 3 will be sufficient.Instead of the variables x, y, . . . ∈ Z /m Z , we work with their Montgomery rep-resentations x M , y M , . . . ∈ Z . These numbers are not entirely unique but boundto be integers from the interval [0 , R n ) fulfilling x M ≡ R n x (mod m ). This meansthat modular divisions still have to be done in some initialization step, one for eachvariable that is initially there, but these turn out to be the only divisions we aregoing to execute!A modular operation, for example x := (( x + 2 xy ) mod m ), is translated into its Montgomery counterpart.
In the example this is x M := (cid:16) x M + 2 x M y M R n mod m (cid:17) . We see here that x M , y M < R n implies x M + 2 x M y M < · R n . An inspection ofO. Forster’s algorithm shows that we always have to compute ( AR n mod m ) forsome A < · R n = R n . A (cid:16) AR n mod m (cid:17) is Montgomery’s REDC function. It occurs everywhere in the algorithm wherenormally a reduction modulo m , i.e. A ( A mod m ), would be done.This looks as if we had not won anything. But, in fact, we won a lot as for computerhardware it is much easier to compute ( AR n mod m ), which is a “reduction frombelow”, than ( A mod m ) which is a “reduction from above” and usually involvestrial divisions.Indeed, A fits into 2 n machine words. It has 2 n so-called limbs . The rightmost,i.e. the least significant, n of those have to be transformed into zero by adding somesuitable multiple of m . Then, these n limbs may simply be omitted.Which multiple of m is the suitable one that erases the rightmost limb A of A ?Well, q · m for q := ( − A · m − mod R ) will do. This operation is in fact anordinary multiplication of unsigned integers in C as ( − A ) on unsigned integersmeans ( R − A ) and multiplication is automatically modulo R . We add q · m to A andremove the last limb. This procedure of transforming the rightmost machine wordof A into zero and removing it has to be repeated n times.Still, m needs to be inverted modulo R = 2 . The naive approach for this would beto use Euclid’s extended algorithm which, unfortunately, involves quite a numberof divisions. At least, we observe that it is necessary to do this only once, not n times although there are n iterations. However, for the purpose of inverting anodd number modulo 2 , there exists a very elegant and highly efficient C macroin GMP, named modlimb invert. It uses a table of the modular inverses of all oddintegers modulo 2 and then executes two Hensel’s lifts in a row. Note that, if HE FIBONACCI SEQUENCE MODULO p i · n ≡ N ) then (2 i − i · n ) · n ≡ N ). We observe that, in thisparticular case, we need no division for the Hensel’s lift.What is the size of the representative of ( AR n mod m ) found? We have A < R n .We add to that less than R n m and divide by R n . Thus, the representative is less than R n + R n mR n = R n m. We want REDC( A ) < R n , the same inequality we have for all variables in Mont-gomery representation. To reach that, we may now simply subtract m in the casewe found an outcome ≥ R n . (This is the point where we use m ≤ R n .)Our version of REDC looks as follows. In order to optimize for speed, we designedit as a C macro, not as a function. It is typically invoked as
REDC (m, REDC BREITE, invm, ?); , with various vari-ables in the place of the ?, after invm is set by modlimb invert (invm, m[0]); and invm = -invm; . Up to now, we always had REDC BREITE = 3.At the very end of our algorithm we find ( F k ) M , the desired Fibonacci number in itsMontgomery representation. To convert back, we just need one more call to REDC.Indeed,( F k mod m ) = (cid:18) F k R n R n mod m (cid:19) = (cid:18) ( F k ) M R n mod m (cid:19) = REDC(( F k ) M ) . Further, ( F k ) M < R n impliesREDC(( F k ) M ) < R n + R n · mR n = 15 + m, i.e. REDC(( F k ) M ) ≤ m .We note explicitly that there is quite a dangerous trap at this point. The residue 0,the one we are in fact looking for, will not be reported as 0 but as m . We workaround this by outputting residues of small absolute value. If ( r mod m ) is foundand r is not below a certain output limit then m − r is computed and comparedwith that limit. Remark 6.2.2.
The integration of the Montgomery representation into our algo-rithm allowed us to avoid practically all the divisions. This caused a stunningreduction of the running time to about one third of its original value.6.3.
Other Optimizations.6.3.1.
We introduced several other optimizations. One, which is worth a mention,is the integration of a pre-computation for the first seven binary digits of p . Note, ifwe let p go linearly through a large interval then its first seven digits will change veryslowly. This means, as a study of our algorithm for the computation of ( F p mod p )shows, that the same first seven steps will be done again and again. We avoid thisand do these steps once, as a pre-computation. As 10 consists of 47 binary digitsthis saves about 14 per cent of the running time.Of course, p is not a constant for the outer loop of our program and its first sevenbinary digits are only almost constant. One needs to watch out for the momentwhen the seventh digit of p changes. Another improvement by a few per cent was obtained through the switchto a different algorithm for the computation of the Fibonacci numbers. Our hand-tailored approach computes the k -th Fibonacci number F k simultaneously with the k -th Lucas number V k . It is based on the formulae F k = F k V k ,V k = V k + 2( − k +1 ,F k +1 = F k V k + V k − k +1 ,V k +1 = F k +1 + 2 F k V k . (13)This is faster than the algorithm explained above as it involves only one multipli-cation and one squaring operation instead of one multiplication and two squaringoperations. It seems here that the number of multiplications and the number ofsquaring operations determine the running time. Multiplications by two are notcounted as multiplications as they are simple bit shifts. Bit shifts and additions area lot faster than multiplications while a squaring operation costs about two thirdsof what a multiplication costs.From that point of view there should exist an even better algorithm. One can makeuse of the formulae F k +1 = 4 F k − F k − + 2( − k ,F k − = F k + F k − ,F k = F k +1 − F k − (14)which we found in the GMP source code. If we meet a bit which is set then wecontinue with F k +1 and F k . Otherwise, with F k and F k − . HE FIBONACCI SEQUENCE MODULO p Here, there are only two squaring operations involved and no multiplications, at all.This should be very hard to beat. Our tests, however, unearthed that the programmade from (14) ran approximately ten per cent slower than the program madefrom (13). For that reason, we worked finally with (13). Nevertheless, we expectthat for larger numbers p , in a situation where additions and bit shifts contributeeven less proportion to the running time, an algorithm using (14) should actuallyrun faster. It is possible that this is the case from the moment on that p > doesno longer fit into three limbs but occupies four. Some other optimizations are of a more practical nature. For example, in-stead of GMP’s mpz functions we used the low level mpn functions for long naturalnumbers. Further, we employed some internal GMP low level functions althoughthis is not recommended by the GMP documentation.The point is that the size of the numbers appearing in our calculations is a-prioriknown to us and basically always the same. When, for example, we multiply twonumbers, then it does not make sense always to check whether the base case multi-plication, the Karatsuba scheme, or the FFT algorithm will be fastest. In our case,mpn mul basecase is always the fastest of the three, therefore we call it directly.6.4.
The Performance Finally Achieved.6.4.1.
As a consequence of all the optimizations described, the CPU time it took ourprogram to test the interval [4 · , · + 2 . · ] of 250 000 000 numbers p suchthat p ≡ In a project of somewhat largerscale, we ran the optimized algorithm on all primes p in the interval [10 , ] suchthat p ≡ ± √ ∈ F p are no longer faster, we ran it, too, on all prime numbers p ∈ [4 · , ]such that p ≡ ± p ∈ [1 . · , · ] such that p ≡ p ≡ ± , ]. To dothis took us around 820 days of CPU time. The computational work was done inparallel on up to 14 PCs from July till October 2004. Output Data
Neither our earlier computations for p < nor themore recent ones for the interval [10 , ] detected any exceptional primes. Aswe covered the intervals systematically and tested each individual prime, this estab-lishes the fact that for all prime numbers p < one has p ∤ F κ ( p ) . There are noexceptional primes below that limit. Theorem 4.4 is verified.
We do never find ( F p ± mod p ) = 0. Does thatmean, we have found some evidence that our assumption, the residues ( F p ± mod p )should be equidistributed in { , p, p, . . . , p − p } , is wrong? Actually, it does not.Besides the fact that the value zero does not occur, all other reasonable statisticalquantities seem to be well within the expected range.Indeed, a typical piece of our output data looks as follows. Durchsuche Fenster mit Nummer 34304.Beginne sieben.Restklassen berechnet.Beginne sieben mit kleinen Primzahlen.Sieben mit kleinen Primzahlen fertig.Fertig mit sieben.Initialisierex mit 110560307156090817237632754212345,y mit 247220362414275519277277821571239und vorz mit 1.10786 Quotient 1912354 mit p := 85760594147971.10787 Quotient 1072750 mit p := 85760627258851.10788 Quotient -1617348 mit p := 85760847493241.10789 Quotient -3142103 mit p := 85761104075891.Initialisierex mit 178890334785183168257455287891792,y mit 400010949097364802732720796316482und vorz mit -1.10790 Quotient -9341211 mit p := 85761921174961.Fertig mit Fenster mit Nummer 34304.Durchsuche Fenster mit Nummer 34305.Beginne sieben.Restklassen berechnet.Beginne sieben mit kleinen Primzahlen.Sieben mit kleinen Primzahlen fertig.Fertig mit sieben.Initialisierex mit 178890334785183168257455287891792,y mit 400010949097364802732720796316482und vorz mit -1.10791 Quotient 3971074 mit p := 85763512710481.10792 Quotient 2441663 mit p := 85764391244491.Fertig mit Fenster mit Nummer 34305.
To make the output easier to understand we do not print ( F p ± mod p ) which isautomatically divisible by p but R ( p ) := ( F p ± mod p ) /p ∈ Z /p Z . Such a quotientmay be as large as p . We output only those which fall into ( − , ) which is verysmall in comparison to p .The data above were generated by a process which had started at 8 · and workedon the primes p ≡ . · , it found 10 792 primes p such that R ( p ) = ( F p − mod p ) /p ∈ ( − , ). HE FIBONACCI SEQUENCE MODULO p On the other hand, assuming equidistribution we would have predicted to find sucha particularly small quotient for around . · X p =8 · p ≡ p prime · − p ≈ (2 · − · ϕ (10) · (log(log(8 . · )) − log(log(8 · )))= 2 · − · (log(log(8 . · )) − log(log(8 · ))) ≈
10 856 . R (82 789 107 950 701) = −
42. We find 1 074 quotients of absolute value lessthan 1 000 000, 98 quotients of absolute value less than 100 000, and 10 of absolutevalue less than 10 000. These are, besides the one above, R (80 114 543 961 461) = − ,R (80 607 583 847 341) = − ,R (80 870 523 194 401) = − ,R (81 232 564 906 631) = 3579 ,R (81 916 669 933 751) = − ,R (83 575 544 636 251) = − ,R (84 688 857 018 011) = − ,R (84 771 692 838 421) = 2281 ,R (85 325 902 236 661) = − . There have been 5 235 positive and 5 557 negative quotients detected.
Remarks 7.3. a) We note explicitly that this is not at all a constructed example.One may basically consider every interval which is not too small and will observethe same phenomena.b) Being very sceptical one might raise the objection that the computations donein our program do not really prove that the 10 792 numbers p which appear in thedata are indeed prime.It is, however, very unlikely that one of them is composite as they all passed two tests.First, they passed the sieve which in this case makes sure they have no prime divi-sor ≤ p | F p − .It is easy to check primality for all of them by a separate program. A more spectacular interval is [0 , ]. One mayexpect a lot more small quotients as all small prime numbers are taken into consid-eration.Here, we may do some statistical analysis on the small positive values of the quo-tient R ′ which is given by R ′ ( p ) := ( F p − mod p ) /p for p ≡ ± R ′ ( p ) := ( F p +2 mod p ) /p for p ≡ ± Except for 0 and 9, all one-digit numbers do appear.Further, the counts are again well within the expected range. For example, con-sider one-digit numbers. R (3) and R (7) are automatically one-digit. Therefore, theexpected count is2 + X p =10 p prime p ≈ · (log(log 10 ) − log(log 10)) ≈ .
849 066which is surprisingly close the 30 one-digit quotients which were actually found.We note that already for two-digit quotients, it is no longer true that they ap-pear only within the subinterval [0 , ]. In fact, there are twelve prime num-bers p ∈ [10 , ] such that R ′ ( p ) < HE FIBONACCI SEQUENCE MODULO p Once again, we may compare this to the expected count which is here X p =10 p prime p ≈ · (log(log 10 ) − log(log 10 )) ≈ .
701 137 73 . References [Co] Cohen, H.: A course in computational algebraic number theory,
Springer,
Graduate TextsMath. 138, Berlin 1993[Fo] Forster, O.: Algorithmische Zahlentheorie (Algorithmic number theory),
Vieweg , Braun-schweig 1996[G¨o] G¨ottsch, G.: ¨Uber die mittlere Periodenl¨ange der Fibonacci-Folgen modulo p (On the averageperiod length of the Fibonacci sequences modulo p ), Dissertation,
Fakult¨at f¨ur Math. undNat.-Wiss., Hannover 1982[IR] Ireland, K., Rosen, M.: A classical introduction to modern number theory, Second edition,
Springer,
Graduate Texts Math. 84, New York 1990[Ja] Jarden, D.: Two theorems on Fibonacci’s sequence,
Amer. Math. Monthly
Math. Comp.