The Fibonacci Sequence and
Binomial Coefficients mod n

April 25, 2004
NKS 2004
Boston

Eric Rowland
Rutgers University
erowland@math.rutgers.edu

Introduction

Sequences of natural numbers play a large role in mathematics.  We can ask 
the simple question:  What do certain sequences look like when we reduce each term a[n] modulo n?

When a[n] is a polynomial, the sequence {Mod[a[n], n]} is quite simple—it is constant for n greater than the constant term:

Table[Mod[(n^2) + 3  n + 7, n], {n, 20}]

Indeed, polynomials are quite "simple" generators of sequences of natural 
numbers.  Let's look at something more complicated:  an exponential function.

ListPlot[Table[Mod[(5^n), n], {n, 1000}], PlotJoined -> True, 

        PlotRange -> All];

Here is the factorial function.

Table[Mod[(n - 1)!, n], {n, 100}]

ListPlot[Table[Mod[(n - 1)!, n], {n, 150}], 

        PlotJoined -> True, PlotRange -> All];

There seems to be some regularity here.  In fact, Wilson's theorem tells 
us Mod[(n - 1)!, 

                        n] == -1 if and only if n is prime.

In general, which sequences are "simple" when we reduce them termwise modulo n?

By "simple" we could mean one of several things.


 •  Can we find patterns and make them explicit by  exhibiting an alternative way to compute Mod[a[n], n] that makes use of these patterns?

•  Can we compute  Mod[a[n], n] without using the Mod function?

•  Is there a computationally reduced way  of computing Mod[a[n], n] (with or without the use of Mod)?

Our focus will be on the first question.  We will see fairly good evidence 
that the binomial coefficients are simple by this criterion (although the 
approach requires incidental use of Mod).  The Fibonacci numbers likewise appear to be simple in this way, but we 
will only examine this case briefly before considering how the problem can be 
reduced to that for binomial coefficients.

Fibonacci Numbers modulo n

The Fibonacci sequence Fibonacci[n] is defined by the initial conditions Fibonacci[1] = Fibonacci[2] = 1 and the recurrence relation  Fibonacci[n] = 

        Fibonacci[n - 1] + Fibonacci[n - 2] for n ≥ 3.

Array[Fibonacci, 20]

Table[Mod[Fibonacci[n], n], {n, 50}]

ListPlot[Table[Mod[Fibonacci[n], n], {n, 300}], 

        PlotJoined -> True, PlotRange -> All];

For which n do we have Mod[Fibonacci[n], n] == 0?  
(In other words, when does n divide Fibonacci[n]?)

Select[Range[1000], Mod[Fibonacci[#], #] == 0 &]

ListPlot[Select[Range[10000], Mod[Fibonacci[#], #] == 0 &], 

        PlotJoined -> True];

In NKS (p. 891), Stephen Wolfram conjectures that if Mod[Fibonacci[n], n] == 0, then n is of the form (5^k) or 12  r, where r is composite or one of {1, 2, 3, 5}.  Experimentally this appears to hold for large n, but only parts of the conjecture have been proven.

A theorem of Fulton and Morris (Acta Arithmetica v.16, 1969) implies that Mod[Fibonacci[24 (5^k)], 24 (5^k)] == 0 for nonnegative integers k.

Table[Mod[Fibonacci[24 (5^k)], 24 (5^k)], {k, 0, 9}]

There are a handful of other results like this.  I proved the converse of 
the conjecture for n = (5^k):  Mod[Fibonacci[(5^k)], (5^k)] == 0 for all integers k ≥ 1.  However, the converse for the case n = 12  r does not hold; there are (many) examples of composite r such that Mod[Fibonacci[12  r], 12  r] ≠ 0, and no description of these is evident from the data.

({#, Mod[Fibonacci[12  #], 12  #]} &) /@ 

      Select[Range[6, 50], ! PrimeQ[#] &]

ListPlot[({#, Mod[Fibonacci[12  #], 12  #]} &) /@ 

          Select[Range[6, 300], ! PrimeQ[#] &], 

        PlotJoined -> True];

For which n do we have Mod[Fibonacci[n], n] == 1?

Select[Range[500], Mod[Fibonacci[#], #] == 1 &]

ListPlot[Select[Range[10000], Mod[Fibonacci[#], #] == 1 &], 

        PlotJoined -> True];

Although it occurs more frequently, this case has been investigated much 
less in the literature.  In the data, one notices the following pattern for 
prime p ≠ 5:

Mod[Fibonacci[p], p] == p - 1 if p ≡ 2, 3 mod 5, and
Mod[Fibonacci[p], p] == 1 if p ≡ 1, 4 mod 5.

(That is, Fibonacci[p] ≡ 

        JacobiSymbol[p, 5] mod p.)  When we search for composite solutions to Mod[Fibonacci[n], n] == 1, we find that the vast majority are even.  The odd solutions form an 
interesting sequence.

Odd solutions to Mod[Fibonacci[n], n] == 1:

Select[Range[100000], 

      OddQ[#] && ! PrimeQ[#] && 

          Mod[Fibonacci[#], #] == 1 &]

ListPlot[{323, 2737, 4181, 6479, 6721, 7743, 11663, 13201, 15251, 

          18407, 19043, 23407, 27071, 34561, 34943, 35207, 39203, 44099, 

          47519, 51841, 51983, 53663, 54839, 64079, 64681, 65471, 67861, 

          68251, 72831, 78089, 79547, 82983, 86063, 90061, 94667, 96049, 

          97921}, PlotJoined -> True];

This sequence shares many terms with two other, related sequences:  
composite n such that Mod[Fibonacci[n - 1], n] == 0 (i.e., n divides Fibonacci[n - 1]), and composite n such that Mod[Fibonacci[n + 1], n] == 0 (i.e., n divides Fibonacci[n + 1]).  In fact, most of the odd composite n satisfying Mod[Fibonacci[n], n] == 1 belong to exactly one of these two sequences.  The exceptions are 7743, 
27071, 54839, 72831, 217257, 388367, 417601, 575599, 670879, 691447, 701569, 
809999, 850541, 881011, ...

This sequence appears to grow roughly linearly but doesn't seem to be given by a simple, closed-form expression.

ListPlot[{7743, 27071, 54839, 72831, 217257, 388367, 417601, 575599, 

          670879, 691447, 701569, 809999, 850541, 881011}, 

        PlotJoined -> True, PlotRange -> All];

I have not been able to find anything more about this sequence directly, but 
I did find relationships to Lucas numbers quite similar to those related to 
Fibonacci numbers.  Lucas numbers satisfy the Fibonacci recurrence relation 
but have different initial conditions:

Lucas[1] = 1;
Lucas[2] = 3;
Lucas[n_] := Lucas[n - 1] + Lucas[n - 2]

Array[Lucas, 20]

Conjecture. 

Let n ∈ ℕ be composite.  Then n is odd and n divides Fibonacci[n] - 1 if and only if exactly one of the following holds:

n divides Fibonacci[n - 1] and n divides Lucas[
            n] - 1, 
n divides Fibonacci[
            n + 1] and n divides Lucas[n] + 1,  or
n ∈ {7743, 27071, 54839, 72831, 217257, 
            388367, 417601, 575599, 670879, 691447, 701569, 809999, 850541, 
            881011,  ... }.

It is easy to prove that if n divides both Fibonacci[n - 1] and Lucas[n] - 1, or that if n divides both Fibonacci[n + 1] and Lucas[n] + 1, then n divides Fibonacci[n] - 1.  However, I have been unable to make progress in the other direction 
because of the difficulty in understanding the sequence 7743, 27071, 54839, 72831,  ...

A proof of this conjecture would tell us which odd n satisfy Mod[Fibonacci[n], n] == 1, but the even case still remains unsolved.

Another Approach:  Binomial Coefficients modulo n

Rather than searching for a formula for Mod[Fibonacci[n], n] directly, we might try to compute it in terms of something else that we 
know.  The Fibonacci numbers satisfy the identity \(TraditionalForm\`Fibonacci[

          n] == \((\((1 + \@5)\)\^n - \((1 - \@5)\)\^n)\)/\((\(2\^n\) 
\@5)\)\).  When we expand an expression of the form (((x + y))^n), we obtain the binomial coefficients Binomial[n, 

          m] = n!/(m! (n - m)!), which comprise Pascal's triangle:

With[{m = 12}, 

      DisplayForm[

        GridBox[Table[

            If[j ≥ m - i && j ≤ m + i && 

                Mod[i + j - m, 2] == 0, 

              Binomial[i, (1/2) (i + j - m)], ""], {i, 0, m}, {j, 0, 

              2 m}], ColumnsEqual -> True, ColumnSpacings -> 0, 

          RowSpacings -> 1]]]

By expanding (((1 + (Sqrt[5])))^n) and (((1 - (Sqrt[5])))^n) in \(TraditionalForm\`Fibonacci[

          n] == \((\((1 + \@5)\)\^n - \((1 - \@5)\)\^n)\)/\((\(2\^n\) 
\@5)\)\), we obtain the identity Fibonacci[n] == (2^(1 - n)) 

          Sum[Binomial[n, m] (5^((m - 1)/2)), {m, 1, n, 2}].  Thus, if we can find the patterns of Mod[Binomial[n, m], n], then we have (at least indirectly) found a way of computing Mod[Fibonacci[n], n].  Unfortunately, only partial results are known.

For prime p, we have the "freshman's dream" identity:  (((x + y))^p) ≡ (x^p) + (y^p) mod p.  This is because Binomial[p, m]  is divisible by p unless m = 0 or m = p.

For composite n, we do not find such simplicity.

ListDensityPlot[

        Reverse[Table[Mod[Binomial[n, m], n], {n, 1, 35}, {m, 0, 35}]], 

        Mesh -> False, Frame -> False];

NiceRasterGraphics[state_, colors_ : 2, size_ : 1] := 

      With[{dim = Reverse[Dimensions[state]]}, 

        Graphics[Raster[Reverse[1 - (state/(colors - 1))]], 

          AspectRatio -> Automatic, 

          PlotRange -> {{0, 

                dim[[1]]}, {0, 

                dim[[2]]}}, 

          ImageSize -> size dim + 1]]

data = 

        Table[Mod[Binomial[n, i + Floor[(n/2)]], n], {n, 

            (2^10)}, {i, -(2^9), (2^9)}];

This is Mod[Binomial[n, m], n].

Show[NiceRasterGraphics[data]];

The image of Mod[Binomial[n, i], n] is a somewhat like a superposition of the images of Mod[Binomial[n, i], m], in which the mth line of Pascal's triangle mod m is roughly included.  (That is, for example, every second line of Mod[Binomial[n, i], 2] (shown below) appears.)  This isn't exactly what happens, however.

Show[

        NiceRasterGraphics[

          Table[Mod[Binomial[n, i + Floor[(n/2)]], 2], {n, 0, 

              (2^8) - 1}, {i, -(2^7), (2^7)}]]];

Show[

        NiceRasterGraphics[

          Table[Mod[Binomial[n, i + Floor[(n/2)]], 3], {n, 0, 

              (2^8) - 1}, {i, -(2^7), (2^7)}]]];

Show[

        NiceRasterGraphics[

          Table[Mod[Binomial[n, i + Floor[(n/2)]], 4], {n, 0, 

              (2^8) - 1}, {i, -(2^7), (2^7)}]]];

The coefficients Mod[Binomial[n, m], n] appear to depend on those of Mod[Binomial[d, m], d] for factors d of n, but in highly nontrivial ways.  Even the case of m = 2  p for p prime seems somewhat complicated:

Show[

        NiceRasterGraphics[Table[

                Mod[Binomial[#, i + #/2], #], {i, -Prime[100], 

                  Prime[100], 2}] & /@ (2  Prime[Range[100]]), 2, 

          1]];

I have experimentally determined a "formula" for Mod[Binomial[n, m], n] for m = 2  p, but this "formula" is not obviously expressible in the notation of 
traditional mathematics; in fact it is a program.  To demonstrate how this 
program works, we will first consider the case p = 3 and then move to the general expression.

Analyzing Mod[Binomial[n, 2  p], n]

p = 3:

Table[Mod[Binomial[n, 6], n], {n, 100}]

Analyzing Mod[Binomial[n, 2  p], n]

We can better see what is going on by partitioning this sequence into strings of 72 terms and transposing.

Transpose[

      Partition[Table[Mod[Binomial[n, 6], n], {n, 300}], 72]]

Clearly there are patterns here.  In fact, each of the nonzero 
subsequences here is linear, and we have the following values of Mod[Binomial[n, 6], n]:

({{n/2, if n ≡ 0 mod 8 and n ≢ 0, 

            6 mod 9}, {n/3, if n ≡ 6 mod 8 and n ≢ 0, 

            6 mod 9}, {n/3, if n ≡ 0 mod 9 and n ≢ 0, 

            6 mod 8}, {2  n/3, if n ≡ 6 mod 9 and n ≢ 

              0, 6 mod 8}, {5  n/6, if n ≡ 0 mod 9 and n ≡ 0, 

            6 mod 8}, {n/6


            0, if n ≡ 6 mod 9 and n ≡ 0, 

            6 mod 8


            otherwise.}})

This can be simplified to the following expression, which appears to hold 
for all n.

(1/6) n Mod[

        If[Mod[n, 8] == 0 || Mod[n, 8] == 6, 3, 0] + 

          If[Mod[n, 9] == 0, 2, 0] + If[Mod[n, 9] == 6, 4, 0], 

        6]

Analyzing Mod[Binomial[n, 2  p], n]

By examining many more cases in precisely this way, I arrived at a general 
expression for Mod[Binomial[n, 2  p], n]:

ModBin2p[p_, n_] := 

      Module[{j = Floor[Log[2, p]] + 2}, 

        (n/(2  p)) Mod[(∑_(k = 2)^p) If[

                  Mod[n - k p, (p^2)] == 0, (k - 1) (p + 1), 

                  0] + If[MemberQ[Residues[p], Mod[n - 1, (2^j)] + 1], p, 0], 

            2 p]]

({#, 

            Table[ModBin2p[#, i], {i, 100}] == 

              Table[Mod[Binomial[i, 2 #], i], {i, 100}]} &) /@ 

      Prime[Range[45]]

The above function requires one additional component, however.  The 
following, rather intriguing function finds the "exceptional" residues 
modulo 2  p, i.e. those numbers for which Mod[Binomial[n, 2  p], n] == n/2.  These cases appear to be distinct from the others, for they have no 
natural place in the sum in ModBin2p.

Residues[p_] := Residues[p] = 

        Module[{j = Floor[Log[2, p]] + 2}, 

          Prepend[Select[Range[(2^j)], 

              Mod[Binomial[#, 2  p], #] == #/2 &], 2  p]]

Analyzing Mod[Binomial[n, 2  p], n]

TableForm[{#, Residues[#]/2} & /@ Prime[Range[2, 11]], 

      TableDepth -> 2]

If we can understand Residues, then we will have succeeded in capturing the patterns of Mod[Binomial[n, 2  p], n].

MaxArithmeticSequence gives the length of the longest arithmetic subsequence of numbers at the 
beginning of a list.

MaxArithmeticSequence[list_] := 

      Module[{diflist = Rest[list - RotateRight[list]]}, 

        If[Length[Split[diflist]] > 

            1, Flatten[

              Position[

                diflist, Split[diflist][[2, 

                  1]], 2, 

                1]][[1]], 

          Length[list]]]

TableForm[{#, Residues[#]/2, 

            MaxArithmeticSequence[Residues[#]/2]} & /@ Prime[Range[2, 11]], 

      TableDepth -> 2]

Analyzing Mod[Binomial[n, 2  p], n]

We can write these sequences of exceptional residues in two-dimensional 
arrays in which each row consists of a sequence of consecutive numbers.  The 
differences between consecutive rows then completely determine the 
sequence.

3	({{3, 4}})

5	({{5, 6, 7, 8}})

7	({{7, 8}})

11	({{11, 12}, {15, 16}})	difference  =  ({{4}})

13	({{13, 14, 15, 16}})

17	({{17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32}})

19	({{19, 20}, {23, 24}, {27, 28}, {31, 32}})	differences  =  ({{4}, {4}, {4}})

Analyzing Mod[Binomial[n, 2  p], n]

Dimension finds the dimensions of the two-dimensional array representations of a 
sequence of residues, and Spacing determines the differences between successive rows.

Dimension[list_] := 
      Module[{dif = 
            If[Length[list] > 2, MaxArithmeticSequence[list], Length[list]]}, 
        Partition[list, dif]]

Spacing[list_] := 
      Module[{dif = MaxArithmeticSequence[list]}, 
        If[dif < Length[list], 
          Table[list[[i]] - 
              list[[i - dif]], {i, 
              dif + 1, Length[list], dif}], {}]]

Analyzing Mod[Binomial[n, 2  p], n]

\(TableForm[\(({#[[1]]\/2, 

              2\^Floor[Log[2, #[[1]]]], 

              If[MaxArithmeticSequence[#] < Length[#], 

                ToString[Length[Spacing[#\/2]] + 1] <> " × " <> 

                  ToString[

                    Length[#]\/\(Length[Spacing[#\/2]] + 1\)] <> " = " <> 

                  ToString[Length[#]], "1 × " <> 

                  ToString[Length[#]] <> " = " <> ToString[Length[#]]], 

              Spacing[#\/2]} &)\) /@ \(Residues /@ Prime[Range[2, 50]]\), 

      TableDepth -> 2]\)

In general, how do we find the spacing sequence?

In general, how do we find the spacing sequence?

It turns out that the sequence of row differences for a prime p can be computed in the following way:  Set a = (2^(Floor[Log[2, p]] + 1)) - p + 1 and b = Inv[a - (2^IntegerExponent)[a, 2]], where Inv is defined below.

Inv[n_] := 

      With[{max = (If[

                      Length[#] > 

                        0, #[[1]], 

                      Length[IntegerDigits[n, 2]] + 1] &)[

                Flatten[Position[IntegerDigits[n, 2], 0, 1, 1]]] - 3}, 

        n - (∑_(i = 0)^max) (2^(Floor[Log[2, n]] - i))]

Inv /@ Range[100]

Let c = 3 (2^(Floor[Log[2, b/3]] + 1)) - b and d = (c/(2^IntegerExponent))[c, 2].  If d == 1, the spacing sequence is (roughly) composed of the number (2^IntegerExponent)[b, 2].  If d > 1, we look at the reversed base-2 representation of Floor[d/2] - (2^(Floor[Log[2, d]] - 1)) and treat it as a string of 0's and 1's.  This string acts as input to 
the function Modify, which alters a sequence in one of two ways:  "0" means "double the 
value of last element of the list", while "1" means "append twice the 
value of the last element of the list".  Modify is recursively executed (via Fold) on the sequence {4, 8} and the string of instructions to generate the numbers that appear in the 
spacing sequence.
The length of the spacing sequence in either case is given 
by \(TraditionalForm\`2\^\(2 + Floor[Log[2, NonTwoPart[a - TwoPart[a]]]] - 
Length[seq]\) - 1\), where seq is the unrepeated, fundamental sequence.

In general, how do we find the spacing sequence?

Here is the explicit code:

TwoPart[n_] := (2^IntegerExponent)[n, 2]
NonTwoPart[n_] := (n/TwoPart)[n]
Modify[list_, instruction_] := 

      If[instruction == 0, 

        Append[Most[list], 

          2  list[[-1]]], 

        Append[list, {2 

               list[[-1, 1]], 

            list[[-1, 

              2]]}]]
SpacingSequenceComponents[d_] := 

      Fold[Modify, {{4, 1}, {8, 2}}, 

        Reverse[IntegerDigits[Floor[(d/2)] - (2^(Floor[Log[2, d]] - 1)), 2, 

            Floor[Log[2, d]] - 1]]]
ComposeSequence[seq_, len_] := 

      Total[(Table[

                If[Mod[i, #[[2]]] 
== 0, #[[1]], 0], {i, len}] & /@ 

            seq)]
SpacingSequence[twop_] := 

      Module[{p = (twop/2), a, b, c, d, pow, seq}, 

        a = (2^(Floor[Log[2, p]] + 1)) - p + 1; b = Inv[a - TwoPart[a]]; 

        c = 3 (2^(Floor[Log[2, (b/3)]] + 1)) - b; d = NonTwoPart[c]; 

        pow = Log[2, ((b + c)/3)]; 

        If[d > 1, 

          seq = ({(1/4) TwoPart[

                        c] #[[1]], #
[[2]]} &) /@ 

              SpacingSequenceComponents[d]; 

          seq = Append[

              seq, (((2^(pow - Floor[Log[2, d]])) 
{seq[[-1, 1]], (1/2) seq
[[-1, 2]]})/TwoPart)[c]], 

          seq = If[

              TwoPart[b] == b, {{TwoPart[b], 1}}, {{TwoPart[b], 1}, 

                (2^pow) {1, (1/TwoPart)[b]}}]]; 

        ComposeSequence[seq, 

          (2^(2 + Floor[Log[2, NonTwoPart[a - TwoPart[a]]]] - Length[seq])) 
- 1]]
ConstructResidues[p_] := ConstructResidues[p] = 

        Module[{w = TwoPart[(2^(Floor[Log[2, p]] + 1)) - p + 1], 

            seq = SpacingSequence[2 p]}, 

          Flatten[2 Table[

                p + i + Total[Take[seq, j]], {j, 0, Length[seq]}, {i, 0, 

                  w - 1}]]]

ConstructResidues agrees with Residues but is completely constructive (i.e., there are no "search and select" procedures involved).

ConstructResidues /@ Prime[Range[2, 100]] == 

      Residues /@ Prime[Range[2, 100]]

Needless to say, this "formula" for Mod[Binomial[n, 2  p], n] is not of the kinds that have appeared in traditional mathematics.  Yet  it succeeds in explaining the patterns that are clear from the data.