Random Fibonacci Sequence

Author

Eric W. Weisstein
February 9, 2007

This notebook downloaded from http://mathworld.wolfram.com/notebooks/IntegerSequences/RandomFibonacciSequence.nb.

For more information, see Eric's MathWorld entry http://mathworld.wolfram.com/RandomFibonacciSequence.html.

©2007 Wolfram Research, Inc. except for portions noted otherwise

Fibonacci Q-Matrices with Arbitrary Signs

<<Utilities`Typesetting`

With[{max=10^5},
m=#[[1,1]]&/@FoldList[Dot,IdentityMatrix[2],{{0,1},{1,#}}&/@((-1)^Table[Random[Integer],{max}])];
ListPlot[Table[{n,N[Abs[((m[[n]])^(1/n))],10]},{n,1,max,100}],PlotJoined->True,PlotStyle->Red,AxesLabel->TraditionalForm/@{n,Abs[(Subscript[a, n])]^InlineFraction[1,n]},
Epilog->{Blue,Line[{{0,1.13198824},{max,1.13198824}}]}
]
]//Timing

Calculated exactly

With[{n=10^6},
m1=Fold[Dot,IdentityMatrix[2],{{0,1},{1,#}}&/@((-1)^Table[Random[Integer],{n}])];
]//Timing

Calculated numerically

With[{n=10^6},
m2=Fold[Dot,N[IdentityMatrix[2],20],{{0,1},{1,#}}&/@((-1)^Table[Random[Integer],{n}])];
]//Timing

With[{n=10^6},
{
N[((Abs[m1[[1,1]]])^(1/n)),10]//Timing,
N[Exp[Log[Sqrt[Max[Eigenvalues[Transpose[m1].m1]]]]/n],10]//Timing,
Exp[Log[Sqrt[Max[Eigenvalues[N[Transpose[m1].m1,10]]]]]/n]//Timing
}
]

This is hard to compute; you need many many random matrices to do it by brute force.  Its value is known to be 1.13198824…, where the next digit is 7, 8, or 9.

$NumberOfMatrices=10^7;
$BlockSize=10^5;

Generate exact integer matrices

ms={{0,1},{1,#}}&/@((-1)^Table[Random[Integer],{$NumberOfMatrices}]);//Timing

ByteCount[ms]/1024.^2 Megabytes

Split data into chunks of length $BlockSize and Dot each chunk exactly (via Fold; can't do it directly with Dot since it accepts a maximum of 4095 elements).  Easy to parallelize ;)

Length[chunks=Fold[Dot,{{1,0},{0,1}},#]&/@Partition[ms,$BlockSize,$BlockSize,{1,1},{}]]//Timing

Reassemble the chunks and dot them (can do this since matrix multiplication is associative)

prod=Dot@@chunks;//Timing

Compute the constant

v=With[{p=N[prod,10^3]},N[Exp[Log[Sqrt[Max[Eigenvalues[Transpose[p].p]]]]/$NumberOfMatrices],20]]

Compare it with the known value

Abs[1.131988248`9-v]