MATLAB Examples

Contents

Fibonacci numbers.

The Fibonacci numbers are an easily accessible area in number theory. As well, they appear in many places in mathematics and in nature , look here too . For those interested in reading more about Fibonacci numbers, you might look at the Fibonacci Quarterly .

I've found many of the Project Euler problems to be interesting and mentally stimulating. Here are a few on Fibonacci numbers: { #2 , #25 , #104 , #137 }.

Less well known than the Fibonacci numbers are the closely related Lucas numbers.

Here is a good reference on the manipulation and derivation of identities relating to the Fibonacci and Lucas numbers

Over the years, I've seen many attempts to generate Fibonacci numbers, and occasionally Lucas numbers. Some of those methods are very good, some very bad, and some are just plain ugly. I think that the Fibonacci sequence is a beautifule one, generated by a simple recursion and accessible to students at every level. As an instructor however, one must remember to teach your students what can be bad about recursion. Why is the factorial function O(N) in its computational complexity when generated recursively, while a recursive method to generate a Fibonacci number has exponential complexity? So when you teach about Fibonacci numbers using recursion, use this as an opportunity to teach them when recursion is a bad thing. Use it as an opportunity to teach techniques like memoization. Next, show what happens when the numbers get large, then teaching about numerical precision issues with floating point numbers. Use Fibnonacci numbers as a jumping off point, then teaching a student a variety of other things that naturally grow from the discussion.

So when it came time for me to write such a function to generate Fibonacci numbers, I decided to think about how best to write one. What must it do? My eventual goals for a Fibonacci tool were:

- It must be able to generate the Fibonacci sequence F(n), starting from n = 1, and extending up to some maximum number n.

- Generate the Lucas numbers for the same values of n.

- Generate Fibonacci and Lucas numbers for a specific value of n, or a disconnected list of values.

- Both sequences are defined for positive, negative, and zero values n. The method employed should be able to handle any integer value of n.

- Return exact integer results for any integer input, recognizing that the resulting number may have many thousands of digits.

- Be as efficient as is reasonably possible.

Along the way, I'll look at how one would NOT want to compute Fibonacci numbers.

First of all, I'll try to make sure that everyone knows what Fibonacci numbers are. Start by defining the first two terms of the Fibonacci sequence as 1. Thus,

F = [1 1];

Next, define all further members of the sequence in terms of the previous elements. This is called a linear recurrence relation . The relation for the Fibonacci sequence is a simple one.

$$F_n = F_{n-1} + F_{n-2}$$

I'll do a few of them for you.

F(3) = F(2) + F(1)
F(4) = F(3) + F(2)
F(5) = F(4) + F(3)
F(6) = F(5) + F(4)
F =
     1     1     2
F =
     1     1     2     3
F =
     1     1     2     3     5
F =
     1     1     2     3     5     8

The Lucas sequence is very similar, using the same fundamental recurrence relation,

$$L_n = L_{n-1} + L_{n-2}$$

merely starting out from a different starting point.

L = [1 3];
L(3) = L(2) + L(1)
L(4) = L(3) + L(2)
L(5) = L(4) + L(3)
L(6) = L(5) + L(4)
L =
     1     3     4
L =
     1     3     4     7
L =
     1     3     4     7    11
L =
     1     3     4     7    11    18

In fact, as you might expect, the two sequences are very closely related.

As an aside, some people might start the Fibonacci sequence from F(0) = 0. Since MATLAB uses an index origin of 1, I'll start at 1 here. In fact, one can define the Fibonacci and Lucas sequences for a negative index too, just run the fundamental recurrences backwards. Useful identities here are:

$$F_{-n} = {-1}^{n+1} F_n$$

$$L_{-n} = {-1}^n L_n$$

Compute the nth Fibonacci number, using a loop

The simplest way to compute the Fibonacci sequence is to use a loop, based on the fundamental recurrences. Suppose we wished to compute F(25)? How might we do so?

First, I'll preallocate a vector for the entire sequence up to 25.

N = 25;
F = zeros(1,N);
F(1:2) = 1;

Now just run the loop. See that the cost, in terms of the number of additions required is only N-2 = 23 additions. On each pass through the loop, I'll need only to do a single add (along with any overhead due to the loop itself.) In fact, this loop is rather efficient, but it does require me to store an array of N elements. Were N to be very large, this might have been a problem. Of course at some point, the size of the number in F(n) might be too large to represent exactly in terms of a floating point integer. Clearly however, if I wish to compute a list of all Fibonacci numbers, this is a very good way to do so.

for n = 3:N
  F(n) = F(n-1) + F(n-2);
end
format long g
F(N)
ans =
       75025
plot(1:N,F,'o')
title('The Fibonacci numbers')
xlabel('N')

How large is F(N)? The largest number that Matlab can represent as an integer exactly when using the double class is 2^53. log2 tells us that the 25th Fibonacci number is only a bit over 2^16.

log2(F(25))
ans =
           16.195083793373

We get the same information from dec2bin. So there is no problem computing F(25). However, had we tried to compute F(100), we might have been unable to compute the exact result in double precision. So Use of large integer utility such as vpi, or that used by the Symbolic toolbox would be essential.

dec2bin(F(N))
numel(dec2bin(F(N)))
ans =
10010010100010001
ans =
    17

Of course, we don't need to save all of those elements if we only wanted to compute the number F(25). So one solution might be to retain only the last few elements of the sequence. This next loop does exactly that. It is slightly more complex, because it must shuffle the results in memory too. This is a frequent tradeoff between many algorithms. A faster algorithm may sometimes be more memory intensive. How much time does the looped algorithm take anyway?

tic
N = 25;
F_nminus2 = 1;
F_nminus1 = 1;
for i = 3:N
  F_n = F_nminus1 + F_nminus2;

  F_nminus2 = F_nminus1;
  F_nminus1 = F_n;
end
toc
F_n
Elapsed time is 0.000082 seconds.
F_n =
       75025

Fibonacci the slow way, using recursion

Next, the Fibonacci sequence is often used as an excuse to teach recursion. In my opinion, while this is an elegant way to write it, recursion can be a computationally horrendous thing to do if done improperly. I'll do it the wrong way first, then compare the time taken and the operation counts.

(By the way, you might notice that I've deliberately left the help on fibrecur and the other functions VERY short. As such, they don't even meet my own standards for help. The reason is because they will be completely listed in this published file, and you should never use them in any serious way. As such, my standard amount of help would be too much clutter. The fibonacci tool that I have included in my vpi toolbox DOES have complete help. It is called fibonacci.m)

dbtype fibrecur
1     function Fn = fibrecur(N)
2     % Computes the Fibonacci number, F(N), using recursion (a poor choice)
3     if N <= 2
4       Fn = 1;
5     else
6       Fn = fibrecur(N-1) + fibrecur(N-2);
7     end
tic,
Fn = fibrecur(25);
toc
Fn
Elapsed time is 1.284648 seconds.
Fn =
       75025

WOW! That was awfully slow. Why is is so terribly slow? Look at what happens. fibrecur(25) calls itself twice, once to evaluate fibrecur(23), and the other time to evaluate fibrecur(24). But then fibrecur(23) and fibrecur(24) each make two calls, each of which will make two more calls. Worse, we end up calling fibrecur(23) twice, one time from fibrecur(25), and the other time from fibrecur(24). You can see how this piles up. In fact, you might guess that the number of recursive calls goes up exponentially with the number being computed. What is the base of the exponential function? Think about it. For example, fibrecur(25) requires calls to both fibrecur(24) and fibrecur(23), each of which require calls to generate the numbers lower than them. So the total number of calls to compute the nth fibonacci number is actually the same as the nth fibonacci number!

We will learn from the Binet formula that the nth fibonacci number grows exponentially, with a base that is essentially

(1+sqrt(5))/2
ans =
          1.61803398874989

To convince ourselves that the time taken is indeed proportional to the n'th Fibonacci number itself, Look at the ratio of times. See that as n grows large, this ratio approaches 1.618...

tic,
Fn = fibrecur(25);
t25 = toc
t25 =
                 1.2507261
tic,
Fn = fibrecur(26);
t26 = toc
t26 =
               2.027842301
t26/t25
ans =
          1.62133204144377

We will show later that a good approximation for the total number of recursive calls for fibrecur(n) is

(1+sqrt(5)/2).^n
----------------
     sqrt(5)

For example, for n = 25, this predicts roughly 75025 recursive calls.

((1+sqrt(5))/2).^25/sqrt(5)
ans =
          75024.9999973343

While 75025 is not yet an immense number, you need not get too large of a number before enponential growth completely overwhelms any computer you will choose to use. For example, with n = 100:

((1+sqrt(5))/2).^100/sqrt(5)
ans =
      3.54224848179263e+20

Recursion is a terrible thing to do to a Fibonacci number. To be honest, I'm not really sure it is a good thing to teach at all, except as an example of what NOT to do. Yes, it is a simple way to teach recursion, for a topic that is easily accessible to beginning students.

Fibonacci numbers, using a faster style of recursion

A nice way to improve the performance of a recursive algorithm is to use Memoization . Why is this code faster?

dbtype fibrecurmemo
1     function [Fn,Fnm1] = fibrecurmemo(N)
2     % Computes the Fibonacci number, F(N), using a memoized recursion
3     if N <= 2
4       Fn = 1;
5       Fnm1 = 1;
6     else
7       [Fnm1,Fnm2] = fibrecurmemo(N-1);
8       Fn = Fnm1 + Fnm2;
9     end

The fibrecurmemo code does not call itself twice for each pass through. Instead, it makes only ONE recursive call, that call returning TWO fibonacci numbers. In effect, we never compute the same lower index Fibonacci number twice. The change from the previous recursive code is a subtle, seemingly minor one, but dramatic in the difference.

Try this version out. While it is recursive in the sense that it calls itself, it is not that slow at all. The recursive function calls involve some system overhead, so it is slower than the looped code. But it is not too bad.

To convince you that this algorithm is basically linear in time, if I should double the size of N, a linear algorithm will take roughly twice as long to execute. Don't try this with the non-memoized version!

tic,
Fn = fibrecurmemo(25);
toc
Fn

tic,
Fn = fibrecurmemo(50)
toc

tic,
Fn = fibrecurmemo(100);
toc

tic,
Fn = fibrecurmemo(200);
toc
Elapsed time is 0.000323 seconds.
Fn =
       75025
Fn =
               12586269025
Elapsed time is 0.000631 seconds.
Elapsed time is 0.000865 seconds.
Elapsed time is 0.001665 seconds.

At some point, even the memoized code will fail, since MATLAB has limits on the recursion depth that it will tolerate. By default, that would have happened at F(502), when the default recursion limit of 500 in MATLAB will be exceeded. Anyway, F(502) is a number with over 100 decimal digits, only a few of which will be correct when done in double precision. If you really want to compute a number that large and do so efficiently, I'll show you how to do it further down. Keep reading.

Another method to compute the Fibonacci numbers

If you recall the plot of the Fibonacci sequence, it was an exponential looking curve.

plot(1:N,F,'o')
title('The Fibonacci numbers')
xlabel('N')

In fact, the Fibonacci sequence is generated by a difference equation , also known as homogeneous recurrence relation. I won't go too deeply into the analysis of a recurrence relation here, but the Fibonacci sequence is generated by the homogeneous relation

F(n) - F(n-1) - F(N-2) = 0

Any such homogeneous linear difference equation has an associated polynomial, the roots of which provide the solution. Here, the associated polynomial is

$$\lambda^2 - \lambda - 1 = 0$$

We can find the roots of the associated polynomial using roots.

roots([1 -1 -1])
ans =
        -0.618033988749895
          1.61803398874989

Symbolically those roots are

$$\lambda_1 = \frac{1-\sqrt{5}}{2}$$

$$\lambda_2 = \frac{1+\sqrt{5}}{2}$$

We should expect the Fibonacci numbers to grow roughly exponentially, as powers of these roots. The actual form, known as Binet's formula , % is given as

$$F(N) = \frac{(\frac{1+\sqrt{5}}{2})^N - (\frac{1-\sqrt{5}}{2})^N}{\sqrt{5}}$$

lambda1 = (1-sqrt(5))/2;
lambda2 = (1+sqrt(5))/2;
binet = @(N) (lambda2.^N - lambda1.^N)/sqrt(5);

See that binet computes the 25th Fibonacci number essentially correctly, but as a floating point number, not as a true integer. As you should expect from an expression that uses real numbers, there is some error down in the least significant bits of the result. As such, this is not an acceptable way to compute the Fibonacci numbers, but it does give us the ability to compute an approximation to their magnitude.

binet(25)
ans =
          75025.0000000001

A Fibonacci determinant?

Not all methods to compute Fibonacci numbers need be efficient. I'll take a quick look at a couple of them before I go on to better ideas.

Mathematical beauty can arise in a variety of ways. Suppose we were to create a tridiagonal matrix, with ones on the main diagonal of the matrix, and i = sqrt(-1) on the first sub and super diagonals? The fibmat function will do exactly that. See that the matrix is indeed nicely tridiagonal.

fibmat = @(N) eye(N) + diag(repmat(sqrt(-1),N-1,1),1) + diag(repmat(sqrt(-1),N-1,1),-1);
spy(fibmat(25))

What does this have to do with the Fibonacci numbers? Look at what happens when I compute the determinants of this sequence of matrices.

for N = 1:10
  det(fibmat(N))
end
ans =
     1
ans =
     2
ans =
     3
ans =
     5
ans =
     8
ans =
    13
ans =
    21
ans =
    34
ans =
    55
ans =
    89

It is initially surprising to many to see the Fibonacci numbers pop out from this expression, det(fibmat(N-1)) yields F(N), the Nth Fibonacci number. Ok, I'll admit that use of a matrix determinant is wildly overkill to compute a Fibonacci number. On the other hand, it is a result that never fails to make me smile. Why should Fibonacci numbers pop out of a determinant? Without deriving it for you, I'll just suggest that an expansion by minors for the determinant (along the first row of the NxN matrix) will actually yield a recurrence relation, relating the desired determinant to the determinants of lower order matrices. In fact, you will end up deriving the Fibonacci recurrence when you do so. Try it. (Please.)

Now, dig a little deeper yet. Recall that the determinant of a matrix is the product of the eigenvalues of that matrix. So next, look at the eigenvalues of this family of tridiagonal matrices. As it turns out, they all have a real part of 1. (At least, approximately so. Remember, we are working in floating point arithmetic in all of this.)

format short g
eig(fibmat(25))
ans =
            1 +     1.9854i
            1 +     1.9419i
            1 +       1.87i
            1 +     1.7709i
            1 +      1.646i
            1 +      1.497i
            1 +     1.3262i
            1 +     1.1361i
            1 +    0.92945i
            1 +    0.70921i
            1 +    0.47863i
            1 +    0.24107i
            1 - 1.5461e-16i
            1 -    0.24107i
            1 -    0.47863i
            1 -     1.9854i
            1 -     1.9419i
            1 -       1.87i
            1 -     1.7709i
            1 -      1.646i
            1 -      1.497i
            1 -     1.3262i
            1 -    0.70921i
            1 -     1.1361i
            1 -    0.92945i

So discard the real parts of the eigenvalues, sort them, and transform them with an arc sine transformation. (I have my devious reasons in all of this.)

plot(asin(sort(imag(eig(fibmat(25))))/2),'o')

Lo and behold, the result is a straight line, so these numbers are equally spaced. Looking at the first differences might help to convince the skeptics.

diff(asin(sort(imag(eig(fibmat(25))))/2))'
ans =
  Columns 1 through 11
      0.12083      0.12083      0.12083      0.12083      0.12083      0.12083      0.12083      0.12083      0.12083      0.12083      0.12083
  Columns 12 through 22
      0.12083      0.12083      0.12083      0.12083      0.12083      0.12083      0.12083      0.12083      0.12083      0.12083      0.12083
  Columns 23 through 24
      0.12083      0.12083

The eigenvalues of an NxN matrix are the roots of an Nth degree characteristic polynomial. A little effort will allow you to transform these roots into the roots of a Chebyshev polynomial of the second kind . When you do the transformation, look at the recurrence relation for these polynomials. It should be no surprise that it turns out to look exactly like that of the Fibonacci numbers. In my opinion, this is yet another of those cases where mathematics shows a true beauty, in that different, seemingly unrelated areas of mathematics touch upon each other. Here we have brought together linear algebra (determinants), special functions (orthogonal polynomials), and number theory (Fibonacci numbers).

N. Cahill, J. D'Errico, J. Spence; "Complex Factorizations of the Fibonacci and Lucas Numbers", Fibonacci Quarterly, V.41, No. 1

Filternacci

As long as I'm looking at alternatives methods for the Fibonacci numbers, try this one liner:

filter([1],[1 -1 -1],[1 1, zeros(1,30)])
ans =
  Columns 1 through 12
           1           2           3           5           8          13          21          34          55          89         144         233
  Columns 13 through 24
         377         610         987        1597        2584        4181        6765       10946       17711       28657       46368       75025
  Columns 25 through 32
      121393      196418      317811      514229      832040     1346269     2178309     3524578

I'll leave it to the student to understand why the filter call works. A closely related idea:

n = 30;
A = eye(n) - diag(ones(n-1,1),-1) - diag(ones(n-2,1),-2);
A(2,1) = 0;
A\[ones(2,1);zeros(n-2,1)]
ans =
           1
           1
           2
           3
           5
           8
          13
          21
          34
          55
          89
         144
         233
         377
         610
         987
        1597
        2584
        4181
        6765
       10946
       17711
       28657
       46368
       75025
      121393
      196418
      317811
      514229
      832040

The flaw with these methods, i.e., Fibonacci determinants, Fibonacci filters and/or linear systems of equations, is all of them are limited in their accuracy. Mathematical and programmatical elegance aside, they fail a main goal of this effort, i.e., to return an exact integer result. They all work in double precision in MATLAB. That will greatly limit us in the size of n if we stick to doubles. In fact, F(78) = 8944394323791464 is the largest Fibonacci number that can be represented exactly in 53 bits of arithmetic precision. The Lucas numbers will exceed that limit above L(77) = 12360848946698171.

Therefore any method used must allow large integers to be represented exactly if the goals is to return exact results for large n. Clearly a loop will be adequate for reasonable values of n, as long as a big integer tool is employed.

tic,
N = 1000;
F_nminus2 = vpi(1);
F_nminus1 = vpi(1);
for i = 3:N
  F_n = F_nminus1 + F_nminus2;

  F_nminus2 = F_nminus1;
  F_nminus1 = F_n;
end
toc
F_n
Elapsed time is 0.333272 seconds.
F_n =
    43466557686937456435688527675040625802564660517371780402481729089536
555417949051890403879840079255169295922593080322634775209689623239873322
471161642996440906533187938298969649928516003704476137795166849228875   

I'd not want to generate F(1000000) that way though.

Large scale Fibonacci numbers

Suppose your goal was to compute a single really large Fibonacci number? Perhaps you wanted to compute F(100000)? Even given the presence of a biginteger toolbox, how might you do this? (I'll use my own vpi tools to do this computation, but the symbolic toolbox will work as well.) First of all, how many digits are there in such a number?

An answer comes from Binet's formula . For any reasonably large n, F(n) will be dominated by one term in the expression.

$$F(N) \approx \frac{(\frac{1+\sqrt{5}}{2})^N}{\sqrt{5}}$$

Therefore we can find the approximate number of decimal digits in the number using a base 10 logarithm.

$$\log_{10}{(F(n))} \approx N \log_{10}{(\frac{1+\sqrt{5}}{2})} - \log_{10}{\sqrt{5}}$$

Thus F(100000) will have just under 21,000 decimal digits.

N = 100000;
N*log10((1+sqrt(5))/2) - log10(sqrt(5))
ans =
        20898

Here are a few identities that will prove useful for really large Fibonacci numbers.

$$F_{2n-1} = {F_{n}}^2 + {F_{n-1}}^2$$

$$F_{2n} = 2F_{n-1}F_n + {F_n}^2$$

For those interested in how one might derive these identities, look at this paper . Start with the F(2*n) identity, then replace the Lucas numbers in that relation.

The trick is to recursively build the nth Fibonacci number. For example, starting with n = 1000000, we first recognize that we con compute it from the second of these identities, IF we but knew the values of F(500000), and F(499999). Next, we can take one more step backwards, since we can use the same identities to find the values of F(500000), and F(499999), given the as yet unknown values of F(250000), and F(249999). Continuously looking downwards, we will eventually end at F(1), a known value.

What happens if n is an odd number? Then n - 1 must be even.

The fibs1 code is fully recursive, but will never generate the same Fibonacci number more than once. As well it requires only O(log2(N)) operations to compute the Nth Fibonacci number. Note that fibs uses the recent release of my vpi toolbox , as found on the MATLAB Central File Exchange .

dbtype fibs1
1     function Fn = fibs(n)
2     % fibs: vpi tool to efficiently compute the n'th and (n-1)'th Fibonacci numbers
3     % usage: Fn = fibs(n)
4     %
5     % A pair of identities are used to cut n by a
6     % factor of 2 on each iteration through. This
7     % makes the call to fibs and O(log2(n)) operation
8     % in time.
9     % 
10    %  F(2*m-1) = F(m)^2 + F(m-1)^2
11    %  F(2*m) = F(m)^2 + 2*F(m-1)*F(m)
12    
13    if (nargin~=1) || (numel(n) ~= 1) || (n <= 0) || (n~=round(n))
14      error('n must be scalar and a positive integer, <= 2^53')
15    elseif isa(n,'vpi')
16      % much faster if n is not a vpi but a double.
17      % also, we don't need to worry about n being
18      % larger than 2^53
19      n = double(n);
20    end
21    
22    % a few special cases to end the recursion,
23    % depending where we came in.
24    if n <= 3
25      if n == 1
26        Fn = vpi([0 1]);
27      elseif n == 2
28        Fn = vpi([1 1]);
29      else
30        % n must be 3
31        Fn = vpi([1 2]);
32      end
33    elseif iseven(n)
34      % N is at least 4, and even
35      
36      % get fibs(n/2)
37      Fnover2 = fibs1(n/2);
38        
39      % combine these results using F(2n), the
40      % 2n Fibonacci relations:
41      % F(2*m-1) = F(m)^2 + F(m-1)^2
42      % F(2*m) = F(m)^2 + 2*F(m-1)*F(m)
43      Fn = Fnover2.*Fnover2;
44      Fn(1) = Fn(1) + Fn(2);
45      Fn(2) = Fn(2) + 2.*Fnover2(1)*Fnover2(2);
46        
47    else
48      % N is at least 5 and odd
49      
50      Fnover2 = fibs1((n-1)/2);
51      
52      % shift up, using the F(2n) relations
53      Fn = Fnover2.*Fnover2;
54      Fn(1) = Fn(1) + Fn(2);
55      Fn(2) = Fn(2) + 2.*Fnover2(1)*Fnover2(2);
56      
57      % and shift up just one step
58      Fn = [Fn(2), sum(Fn)];
59    end
60    
61    function result = iseven(n)
62    % tests if a scalar value is an even integer
63    if isnumeric(n)
64      result = (mod(n,2) == 0);
65    elseif isa(n,'vpi')
66      % must have been a vpi
67      result = (mod(trailingdigit(n,1),2) == 0);
68    else
69      error('n must be either numeric or vpi')
70    end

By way of comparison, recall that the loop above to compute F(1000) took about 1.4 seconds on my machine to do its work.

F_nminus1
F_n
F_nminus1 =
    43466557686937456435688527675040625802564660517371780402481729089536
555417949051890403879840079255169295922593080322634775209689623239873322
471161642996440906533187938298969649928516003704476137795166849228875   
F_n =
    43466557686937456435688527675040625802564660517371780402481729089536
555417949051890403879840079255169295922593080322634775209689623239873322
471161642996440906533187938298969649928516003704476137795166849228875   

Compare the use of fibs1.m to compute F(999) and F(1000),

tic
Fn = fibs1(1000);
toc
Fn
Elapsed time is 0.023941 seconds.
Fn =
 
vpi element: (1,1)
    26863810024485359386146727202142923967616609318986952340123175997617
981700247881689338369654483356564191827856161443356312976673642210350324
634850410377680367334151172899169723197082763985615764450078474174626   
 
vpi element: (1,2)
    43466557686937456435688527675040625802564660517371780402481729089536
555417949051890403879840079255169295922593080322634775209689623239873322
471161642996440906533187938298969649928516003704476137795166849228875   

How about F(100000)? Don't bother using the loop to do so. That would be simply too slow here. The funny thing is, fibs took considerably less time to compute F(100000), a number with 21,000 decimal digits than the poorly written fibrecur took to compute F(25). Even worse, fibrecur was working with the standard double precision arithmetic in MATLAB, so in theory, it had a significant advantage in speed over the vpi based fibs1.m.

tic
Fn = fibs1(100000);
toc
Elapsed time is 0.150936 seconds.

Rather than display these entire numbers, was our estimate using Binet's formula correct?

log10(Fn)
ans =
        20898        20898

Large scale Fibonacci and Lucas numbers

Clearly fibs1 works fairly nicely, to compute isolated Fibonacci numbers, but it does not generate the Lucas sequence. The first 25 numbers in the pair of related sequences are easily generated with a loop.

F = [1;1;zeros(23,1)];
L = [1;3;zeros(23,1)];
for n = 3:25
  F(n) = F(n-1) + F(n-2);
  L(n) = L(n-1) + L(n-2);
end
format long g
disp('           n        F(n)        L(n)')
disp([(1:25)',F,L])
           n        F(n)        L(n)
           1           1           1
           2           1           3
           3           2           4
           4           3           7
           5           5          11
           6           8          18
           7          13          29
           8          21          47
           9          34          76
          10          55         123
          11          89         199
          12         144         322
          13         233         521
          14         377         843
          15         610        1364
          16         987        2207
          17        1597        3571
          18        2584        5778
          19        4181        9349
          20        6765       15127
          21       10946       24476
          22       17711       39603
          23       28657       64079
          24       46368      103682
          25       75025      167761

A fundamental relation exists between the Fibonacci numbers and the Lucas numbers. Much like in trigonometry, where we know the sin and cos functions are related by a fundamental identity, sin(x)^2 + cos(x)^2 == 1, there also exists an identity that links these two integer sequences.

$$L_n^2 - 5F_n^2 = (-1)^n 4$$

Test the identity here. The result should alternate in sign, as +/- 4.

L.^2 - 5*F.^2
ans =
    -4
     4
    -4
     4
    -4
     4
    -4
     4
    -4
     4
    -4
     4
    -4
     4
    -4
     4
    -4
     4
    -4
     4
    -4
     4
    -4
     4
    -4

Again, as with the trig family of functions, where we can derive double angle formulas, there are identities that allow you to compute the value of F(2*n) given only the values of F(n) and L(n).

$$F_{2n} = F_n L_n$$

$$L_{2n} = \frac{5 F_n^2 + L_n^2}{2}$$

A simple verification of these identities seems appropriate.

n = (1:12)';
disp('Compute F(2*n)')
disp('           n         2*n      F(2*n)   F(n)*L(n)')
[n,2*n,F(2*n),F(n).*L(n)]
Compute F(2*n)
           n         2*n      F(2*n)   F(n)*L(n)
ans =
           1           2           1           1
           2           4           3           3
           3           6           8           8
           4           8          21          21
           5          10          55          55
           6          12         144         144
           7          14         377         377
           8          16         987         987
           9          18        2584        2584
          10          20        6765        6765
          11          22       17711       17711
          12          24       46368       46368
disp('Compute L(2*n)')
disp('           n         2*n      L(2*n)  (5*F(n).^2 + L(n).^2)/2')
[n,2*n,L(2*n),(5*F(n).^2 + L(n).^2)/2]
Compute L(2*n)
           n         2*n      L(2*n)  (5*F(n).^2 + L(n).^2)/2
ans =
           1           2           3           3
           2           4           7           7
           3           6          18          18
           4           8          47          47
           5          10         123         123
           6          12         322         322
           7          14         843         843
           8          16        2207        2207
           9          18        5778        5778
          10          20       15127       15127
          11          22       39603       39603
          12          24      103682      103682

The identities above work nicely when n is even. If n is odd, just subtract 1, and it will be even, allowing us to use the 2*n rules we just saw. Afterwards, two more identities allow us to shift back up one step.

$$F_{m+1} = \frac{F_m + L_m}{2}$$

$$L_{m+1} = \frac{5 F_m + L_m}{2}$$

fibs2 uses these identities to efficiently compute the pair F(n) and L(n).

dbtype fibs2
1     function [Fn,Ln] = fibs2(n)
2     % fibs2: vpi tool to efficiently compute the n'th Fibonacci number and the n'th Lucas number
3     % usage: [Fn,Ln] = fibs2(n)
4     %
5     % the "2*n" identities are employed here:
6     %  F(2*m) = F(m)*L(m)
7     %  L(2*m) = (5*F(m)^2 + L(m)^2)/2
8     %
9     % coupled with the addition identities
10    %
11    %  F(m+1) = (F(m) + L(m))/2
12    %  L(m+1) = (5*F(m) + L(m))/2
13    
14    
15    if (nargin~=1) || (numel(n) ~= 1) || (n~=round(n)) || (abs(n) > 2^53)
16      error('n must be scalar and an integer, <= 2^53 in absolute value')
17    end
18    
19    % much faster if n is not a vpi but a double.
20    % also, we don't need to worry about n being
21    % larger than 2^53, as this would have a vast
22    % number of digits.
23    n = double(n);
24    
25    % ensure that n is positive, handle the negative
26    % cases too here.
27    if n < 0
28      n = abs(n);
29      [Fn,Ln] = fibs2(n);
30      if iseven(n)
31        Fn = -Fn;
32      else
33        Ln = -Ln;
34      end
35      return
36    end
37    
38    % a few special cases to end the recursion,
39    % depending where we came in.
40    if n <= 3
41      switch n
42        case 2
43          Fn = vpi(1);
44          Ln = vpi(3);
45        case 3
46          Fn = vpi(2);
47          Ln = vpi(4);
48        case 1
49          Fn = vpi(1);
50          Ln = vpi(1);
51        case 0
52          Fn = vpi(0);
53          Ln = vpi(2);
54      end
55      return
56    end
57    % n is now known to be at least 4.
58    %
59    % Next, we will use several identities relating
60    % F(2*n) and L(2*n) to F(n) and L(n), and between
61    % the Fibonacci and Lucas numbers.
62    
63    % Is n a multiple of 2?
64    if iseven(n)
65      % even, so use 
66      [Fnover2,Lnover2] = fibs2(n./2);
67      
68      % Use the double argument formulas
69      Fn = Fnover2.*Lnover2;
70      Ln = (5 .*Fnover2.*Fnover2 + Lnover2.*Lnover2)./2;
71      
72    else
73      % n must be odd. So if we subtract 1, then
74      % n-1 will be even. I could have added 1
75      % instead, but I want to make n as small
76      % as possible, as quickly as I can.
77      [Fnover2,Lnover2] = fibs2((n-1)./2);
78      
79      % Use the double argument identities
80      Fnminus1 = Fnover2.*Lnover2;
81      Lnminus1 = (5 .*Fnover2.*Fnover2 + Lnover2.*Lnover2)./2;
82      
83      % this has only gotten us to F(n-1) and
84      % L(n-1). use the addition identities to
85      % move up one more step.
86      Fn = (Lnminus1 + Fnminus1)./2;
87      Ln = (5 .*Fnminus1 + Lnminus1)./2;
88    end
89      
90    function result = iseven(n)
91    % tests if a scalar value is an even integer
92    if isnumeric(n)
93      result = (mod(n,2) == 0);
94    elseif isa(n,'vpi')
95      % must have been a vpi
96      result = (mod(trailingdigit(n,1),2) == 0);
97    else
98      error('n must be either numeric or vpi')
99    end
100   
101   

Apply fibs2 to compute F(100000) and L(100000) as a coupled pair.

tic,[Fn,Ln] = fibs2(100000);toc
Elapsed time is 0.599756 seconds.

See that fibs2 was a bit slower than fibs1, but it returns both the last two Fibonacci numbers, as well as the last two Lucas numbers. And, fibs2 is still a recursive function, so there is function call overhead wasted every time fibs2 gets invoked.

fibs2 is not the end to how might optimize the computation of the Fibonacci numbers. There are other steps one can take. In fact, very often careful programming can allow a recursive code to be re-written as an simply looped code. This has been done with fibs3, by working with the binary representation of the original number n. Think about what was done by fibs2. Each pass through, fibs2 divides by 2 if n is even, or subtracts 1 and then divides by 2 if it is odd. Look at the result of such an operation at each step. I'll start with n = 197 as an example:

n = 197
dec2bin(n)
n =
   197
ans =
11000101

n was odd, so subtract 1, then divide by 2.

n = (n - 1)/2,dec2bin(n)
n =
    98
ans =
1100010

n was even, so just divide by 2.

n = n/2,dec2bin(n)
n =
    49
ans =
110001

Repeat each operation as indicated by the next bit in the binary representation.

n = (n - 1)/2,dec2bin(n)
n = n/2,dec2bin(n)
n = n/2,dec2bin(n)
n = n/2,dec2bin(n)
n = (n - 1)/2,dec2bin(n)
n =
    24
ans =
11000
n =
    12
ans =
1100
n =
     6
ans =
110
n =
     3
ans =
11
n =
     1
ans =
1

fibs3 starts with a binary representation of n, then works UPWARDS through this same chain in a simple loop. Also, fibs3 uses a slightly different pair of identities to minimize the number of multiplies and divides at each step.

dbtype fibs3
1     function [Fn,Ln] = fibs3(n)
2     % fibs: vpi tool to efficiently compute the n'th Fibonacci number and the n'th Lucas number
3     % usage: [Fn,Ln] = fibs3(n)
4     %
5     % For efficiency, fibs3 uses a variety of tricks to
6     % maximize speed. While computation of fibonacci numbers
7     % is commonly done recursively, fibs3 does so using a
8     % direct iterative scheme given the binary representation
9     % of n. In addition, several Fibonacci and Lucas number
10    % identities are employed to maximize throughput.
11    %
12    % The methods employed by fibs3 will be O(log2(n)) in time.
13    
14    if (nargin~=1) || (numel(n) ~= 1) || (n~=round(n)) || (abs(n) > 2^53)
15      error('n must be scalar and an integer, <= 2^53 in absolute value')
16    end
17    
18    % much faster if n is not a vpi but a double.
19    % also, we don't need to worry about n being
20    % larger than 2^53, as this would have a vast
21    % number of digits.
22    n = double(n);
23    
24    % ensure that n is positive, handle the negative
25    % cases later.
26    nsign = sign(n);
27    n = abs(n);
28    
29    % get the binary representation of n. Thus
30    % nbin is a character vector, of length
31    % ceil(log2(n)).
32    nbin = dec2bin(n);
33    
34    % get the 4 highest order bits from nbin
35    k = min(numel(nbin),4);
36    
37    % zero is a special case to stop at.
38    if n == 0
39      Fn = vpi(0);
40      Ln = vpi(2);
41      return
42    else
43      % start the sequence from the top
44      % few bits of n.
45      nhigh = bin2dec(nbin(1:k));
46      Fseq = [1 1 2 3  5  8 13 21 34  55  89 144 233 377  610];
47      Lseq = [1 3 4 7 11 18 29 47 76 123 199 322 521 843 1364];
48      
49      Fn = vpi(Fseq(nhigh));
50      Ln = vpi(Lseq(nhigh));
51    end
52    
53    % We need to loop forwards. Essentially, we started
54    % with the highest order bit(s) of the binary representation
55    % for n. Look at each successively lower order bit.
56    % If the next bit is 0, then we are essentially doubling
57    % the index at this step. If the next bit is odd, then
58    % we are moving to 2*n+1.
59    for k = 5:numel(nbin)
60      bit = (nbin(k) == '1');
61      if bit
62        % the next bit was odd. Use
63        % a 2*n+1 rule to step up.
64        F2n = Fn.*Ln;
65        % we want to do this...
66        %   L2n = (5 .*Fn.*Fn + Ln.*Ln)./2;
67        % Instead use the identity that
68        %   5F(n)^2 + L(n)^2 = 2*L(n)^2 + 4*(-1)^(n+1)
69        % to make that expression more efficiently
70        % computed. See that the form used below
71        % has only a single multiplication between a
72        % pair of large integers, whereas the prior
73        % form for L2n would have had several multiples
74        % as well as a divide.
75        L2n = Ln.*Ln + 2*(-1)^(nhigh+1);
76        
77        Fn = (L2n + F2n)./2;
78        Ln = (5 .*F2n + L2n)./2;
79      else
80        % the next bit was even. Use the 2*n
81        % rule to step up.
82        F2n = Fn.*Ln;
83        Ln = Ln.*Ln + 2*(-1)^(nhigh+1);
84        Fn = F2n;
85      end
86      
87      % update the top bits of n
88      nhigh = 2*nhigh + bit;
89    end
90    
91    % if n was negative, then we may need to apply a
92    % sign change to Fn and Ln.
93    if nsign < 0
94      if iseven(n)
95        Fn = (-nsign).*Fn;
96        Ln = nsign.*Ln;
97      else
98        Fn = nsign.*Fn;
99        Ln = (-nsign).*Ln;
100     end
101   end
102   
103   % ==================================
104   % End mainline, begin subfunctions.
105   % ==================================
106   
107   function result = iseven(n)
108   % tests if a scalar value is an even integer
109   if isnumeric(n)
110     result = (mod(n,2) == 0);
111   elseif isa(n,'vpi')
112     % must have been a vpi
113     result = (mod(trailingdigit(n,1),2) == 0);
114   else
115     error('n must be either numeric or vpi')
116   end
117   
118   

fibs3 has become fairly efficient.

tic,
[Fn,Ln] = fibs3(100000);
toc
Elapsed time is 0.093461 seconds.

My final version is embodied in fibonacci.m. This code uses the ideas in fibs3 to compute isolated Fibonacci and Lucas numbers efficiently. In the event that sets of numbers are requested, an addition identity is employed, relating F(m + n) and L(m + n) with lower order numbers.

dbtype fibonacci
1     function [Fn,Ln] = fibonacci(n)
2     % fibonacci: vpi tool to efficiently compute the n'th Fibonacci number and the n'th Lucas number
3     % usage: [Fn,Ln] = fibonacci(n)
4     %
5     % Compute the nth Fibonacci number as well as the nth Lucas
6     % Lucas number. In the event that all members of these
7     % sequences from 1 to n are desired, then a simple,
8     % direct loop would be more efficient, and more direct.
9     %
10    % Both the Fibonacci numbers and the Lucas numbers
11    % are defined by the same basic recursion:
12    %
13    %  F(n) = F(n-1) + F(n-2)
14    %  L(n) = L(n-1) + L(n-2)
15    %
16    % The difference is the starting point. The Fibonacci
17    % numbers start with F(1) = F(2) = 1, whereas the Lucas
18    % sequence starts with L(1) = 1, and L(2) = 3. The first
19    % few members of these sequences are:
20    %
21    %  Fibonacci: [1 1 2 3 5 8 13 21 ... ]
22    %  Lucas:     [1 3 4 7 11 18 29 ... ]
23    %
24    % These sequences are also defined for n = 0 and for
25    % negative values of n.
26    %
27    % For efficiency, fibonacci uses a variety of tricks to
28    % maximize speed. While computation of fibonacci numbers
29    % is commonly done recursively, fibonacci does so using a
30    % direct iterative scheme given the binary representation
31    % of n. In addition, several Fibonacci and Lucas number
32    % identities are employed to maximize throughput.
33    %
34    % The methods employed by fibonacci will be O(log2(n)) in time.
35    %
36    %
37    % Arguments: (input)
38    %  n - any non-negative integer, vpi or numeric.
39    %      n must be less than 2^53 (in theory. Even that
40    %      number would be impossibly huge to compute.
41    %      Don't bother to try it.) In practice, recognize
42    %      that F(1e6) is a number with 208995 digits.
43    %
44    %      When n is a vector or array, fibonacci will
45    %      generate each of the indicated Fibonacci and
46    %      Lucas numbers.
47    %
48    % Arguments: (output)
49    %  Fn, Ln - scalar vpi number, containing the nth
50    %      Fibonacci number and nth Lucas numbers in
51    %      their respective sequences.
52    %
53    % Example:
54    %  [Fn,Ln] = fibonacci(150)
55    %
56    % Fn =
57    %    9969216677189303386214405760200
58    % Ln =
59    %    22291846172619859445381409012498
60    %
61    %  See also: 
62    %
63    %  Author: John D'Errico
64    %  e-mail: woodchips@rochester.rr.com
65    %  Release: 1.0
66    %  Release date: 4/30/09
67    
68    if (nargin ~= 1)
69      error('fibonacci accepts only 1 argument')
70    elseif any(n(:)~=round(n(:))) || any(abs(n) > 2^53)
71      error('n must be an integer, <= 2^53 in absolute value')
72    end
73    
74    % The first 15 Fibonacci and Lucas numbers to
75    % start things off efficiently.
76    Fseq = [0 1 1 2 3  5  8 13 21 34  55  89 144 233 377  610];
77    Lseq = [2 1 3 4 7 11 18 29 47 76 123 199 322 521 843 1364];
78    
79    % intialize Fn and Ln to the proper size,
80    % in case n is a vector or array
81    Fn = repmat(vpi(0),size(n));
82    Ln = Fn;
83    
84    % much faster if n is not a vpi but a double.
85    % also, we don't need to worry about n being
86    % larger than 2^53, as this would have a vast
87    % number of digits.
88    n = double(n(:));
89    
90    % catch any zeros in n first.
91    k = (n(:) == 0);
92    % Fn(k) is already zero
93    if any(k)
94      Ln(k) = vpi(2);
95    end
96    
97    % Negative values for n will be inconvenient
98    % in a loop, so make them all positive. deal
99    % with any signs later.
100   nsign = 2*(n>=0) - 1;
101   if any(n<0)
102     neven = iseven(n);
103     n = abs(n);
104   end
105   
106   if numel(n) > 1
107     % sort them and work upwards
108     [n,ntags] = sort(n);
109     
110     flag = false;
111     for i = 1:numel(n)
112       if n(i) <= 15
113         % small values of n are pre-computed
114         Fn(i) = Fseq(n(i) + 1);
115         Ln(i) = Lseq(n(i) + 1);
116       elseif (i == 1) || ((n(i) - n(i-1)) > 15)
117         % the smallest value of n is > 15,
118         % or the difference from the last value
119         % computed was too large
120         [Fn(i),Ln(i)] = fibonacci(n(i));
121         
122         flag = false;
123       elseif n(i) == n(i-1)
124         % will happen if there were positive and
125         % negative values in the list
126         Fn(i) = Fn(i-1);
127         Ln(i) = Ln(i-1);
128         
129         flag = false;
130       elseif (n(i) == (n(i-1)+1))
131         % we can use an addition formula to
132         % get Fn & Ln. Which one?
133         if flag
134           % the last two values of n were
135           % separated by 1, so just use the
136           % two term recurrence
137           Fn(i) = (Fn(i-1) + Fn(i-2));
138           Ln(i) = (Ln(i-1) + Ln(i-2));
139         else
140           % we must use the addition formula
141           Fn(i) = (Fn(i-1) + Ln(i-1))./2;
142           Ln(i) = (5 .*Fn(i-1) + Ln(i-1))./2;
143           
144           % flag indicates whether the last two
145           % numbers were consecutive
146           flag = true;
147         end
148         
149       else
150         % 1 < n(i) - n(i-1) <= 15
151         m = n(i) - n(i-1);
152         Fm = Fseq(m+1);
153         Lm = Lseq(m+1);
154         
155         % use an addition formula
156         Fn(i) = (Fm.*Ln(i-1) + Lm.*Fn(i-1))./2;
157         Ln(i) = (5 .*Fm.*Fn(i-1) + Lm.*Ln(i-1))./2;
158         
159         flag = false;
160       end
161       
162     end % for i = 1:numel(n)
163     
164     % shuffle Fn and Ln to reflect the sort
165     Fn(ntags) = Fn;
166     Ln(ntags) = Ln;
167     
168   else
169     % For only one value of n, compute it individually.
170     % Uses an efficient iterative (recursive, but not
171     % really so) scheme.
172     
173     % get the binary representation of n. Thus
174     % nbin is a character vector, of length
175     % ceil(log2(n)).
176     nbin = dec2bin(n);
177     
178     % get the 4 highest order bits from nbin
179     k = min(numel(nbin),4);
180     
181     % start the sequence from the top
182     % few bits of n.
183     nhigh = bin2dec(nbin(1:k));
184     
185     Fn = vpi(Fseq(nhigh+1));
186     Ln = vpi(Lseq(nhigh+1));
187     
188     % We need to loop forwards. Essentially, we started
189     % with the highest order bit(s) of the binary representation
190     % for n. Look at each successively lower order bit.
191     % If the next bit is 0, then we are essentially doubling
192     % the index at this step. If the next bit is odd, then
193     % we are moving to 2*n+1.
194     for k = 5:numel(nbin)
195       bit = (nbin(k) == '1');
196       if bit
197         % the next bit was odd. Use
198         % a 2*n+1 rule to step up.
199         F2n = Fn.*Ln;
200         % we want to do this...
201         %   L2n = (5 .*Fn.*Fn + Ln.*Ln)./2;
202         % Instead use the identity that
203         %   5F(n)^2 + L(n)^2 = 2*L(n)^2 + 4*(-1)^(n+1)
204         % to make that expression more efficiently
205         % computed. See that the form used below
206         % has only a single multiplication between a
207         % pair of large integers, whereas the prior
208         % form for L2n would have had several multiples
209         % as well as a divide.
210         L2n = Ln.*Ln + 2*(-1)^(nhigh+1);
211         
212         Fn = (L2n + F2n)./2;
213         Ln = (5 .*F2n + L2n)./2;
214       else
215         % the next bit was even. Use the 2*n
216         % rule to step up.
217         F2n = Fn.*Ln;
218         Ln = Ln.*Ln + 2*(-1)^(nhigh+1);
219         Fn = F2n;
220       end
221       
222       % update the top bits of n
223       nhigh = 2*nhigh + bit;
224     end % for k = 5:numel(nbin)
225     
226   end % if numel(n) > 1
227   
228   % if n was negative, then we may need to apply a
229   % sign change to Fn and Ln.
230   if any(nsign < 0)
231     k = neven & (nsign < 0);
232     Fn(k) = -Fn(k);
233   
234     k = (~neven) & (nsign < 0);
235     Ln(k) = -Ln(k);
236   end
237   
238   % ==================================
239   % End mainline, begin subfunctions.
240   % ==================================
241   
242   function result = iseven(n)
243   % tests if a scalar value is an even integer, works
244   % for either numeric or vpi inputs
245   if isnumeric(n)
246     result = (mod(n,2) == 0);
247   elseif isa(n,'vpi')
248     % must have been a vpi
249     result = (mod(trailingdigit(n,1),2) == 0);
250   else
251     error('n must be either numeric or vpi')
252   end
253   
254   

I'll give up here. Perhaps there will be some ideas gained by some of my readers on how to compute and work with Fibonacci numbers. Even if not, I'll hope that someone has learned some tricks on how to write recursive functions, and how to recognize when recursion may be inappropriate.

tic,
[Fn,Ln] = fibonacci(100000);
toc
Elapsed time is 0.097995 seconds.

F(100000), the 100,000th Fibonacci number

Fn
Fn =
    25974069347221724166155034021275915414880485386517696584724770703952
534543511273686265556772836716744754637587223074432111638399473875091030
965697382188304493052287638531334921353026792789567010512765782716356080
730505322002432331143839865161378272381247774537783372999162146340500546
698603908627509966393664092118901252719601721050603003505868940285581036
751176582513683774386849364134573388343651587754253719124105003321959913
300622043630352137565254218239986908485563740801792517616293917549634585
586163007628199160811098365263529954406942842065710460449038056471363460
330005208522777075544467947237090309790190148604328468198579610159510018
506082649192345873133991501339199323631023018641725364771362664750801339
824312317034314529641817900511879573167668349799016820118499077566864568
450662873924856039140476051995500662888263458771894106803700918793650017
330117100283104739474562560914449328213748555738640805798130282666402703
542944121049199958031318768058991865134251759599115205631553377039969410
355182752749199598022575079020377981030899229849963044962558140455170002
502997643221934621653662108418767454282982613982344783665815880408190033
073829395000821320093747154851310272208173054322648669496309879147143629
255542526240439996153269798768075106468190687921182991679644091782718685
617029181022126792674013626504997849688436809752547001310045741864064482
994858725517447466956518791269169932445648176733222571493149677633458466
238303338202397024368594782876418757885729107101337003000942293335972927
791914092128049015459762627910570552481588840517794181929052167695766087
488155678601288183543542923073978101547857013284386127286201766539534449
930019800629538936985500723286651317181135886613537472684585432548981137
176605194616937916884425342594781263103889520479565943807153019112539648
471126389007133628569101551453423329441284357220996286746119420951661002
309740709965531900508158669911445442647882872642845017253320486483194578
920399848938236367456182203750973485668474338872490493370316338265717607
297788917989136673251906232471180372801739215723908227692280772924566627
505383375006926077210593619421268920302567443565378008318306375933345023
502569729065152853271943677560156660399164048825639676930792905029514886
934137991251748566670747175149389790386533381395346848378086126737554383
821108448976538368483182588363399173104558509056638462025014631311831087
429077292622159430204291594740306101839816855066950261973761508571761199
475875722129872053120607918649803615960923395941041186351688548839119185
179061511562752936158490008721501922265117853150892510275280451512386037
921846921215338292871369243215273327141574788295902601571954853164447945
467502858402360002383447905203451080332820138038807089807348326201227952
633606773669875783326254859449060219173688677862411205621098369850197290
177157801120404586491539351157834995461006366357454485082418882790675313
599505192062229760153765297973085881648731173082370598284894044874039320
535929359764541655607954724778620299692329561389719894679422187273605123
365595211331087787582288795975803204596084790245063851941743126163775104
599211024868794963417068620929088930685252348056925998333775103901013166
178123051145719327066291671254465121517468025481903583516889717075706778
656188008220346836321018130262329960275994035799977740462449521145315883
703579044832931500072461734173558055678321534543411700202585608091662941
986374015145695722728369219632295111877625307534025947814482046574602884
855000628069348113982760168555840795421620575435572915106415375929390228
843561207926437055600623679865443824643739469724719459965557955058380348
255978396827760847315302517889517186307227611036305093600742622617173630
586132915440246954329046162586917746305785076749374879923291817501634840
688134655343709975893536074051729094126976575932951568186247471276364688
365517570183534172746626073065104511957628663499228486787805910851189856
535554349587616640164475880286336297040462890970677362565843002353147494
612339120686321466370878446992104275415694109122465685712047172411333784
898167640969249816334211768571503116710400681753031921154156119580425706
586931272762137106974722260296555246110537155545324997508432752001992143
019105053629960070429632978051030666506387862681576587726837451289768507
963663710593809112254288358391941211547737599813019216509521401333060709
873137329265181692268450634439540567298120315463923249817937804691037934
221694952291007930299492375072993250630509428139027930841344730614116433
556147640931044259184813639305423693789765205264563476483182726333715121
120306292338892864879492097378478618848682608046473195392008403983080088
038690495574197562192939221108257663976813610444900247209483403267967688
376213967440757138872928630798218493143438797780887379588968409461434159
271317578365114578289355818599029235343888888465874521308381377794436361
197628390368945957601203165022798579015453447473527069728514545998614229
027372911314637820455162254475353567736227936485450357102086445412089842
350389087702230398493802147348096874333362254491501174117515707045610508
952740002063804979679604026178186644812485472696308234733772455433905198
413087697812765659167642290229481817630757102557933650081522863836344931
380899717850870708636322058690189383777660630060667577324272729292474212
952650007066467227300099561241914091389846752249557907293984956087504566
942177715511073466304566039441362358884436762152739285970722879373559667
239246138274687032178584599482575147454064364609970593161205968415604732
343966524572316503177928338605903883604176914287327357039868033426046700
717173635730911229813069032861371225979370966057751729645282637574340757
922821807443529086696068540217185978911663338638585897362091142484321786
450394791954242081916260885710691104339948014730131008698488664307212167
624731196181907378207665829682807960794822595490363282665780069948568253
005364366748225346037051345036031521542969439918662368576380623512098844
487411386001711736476321260299614085619255997075668278667787323774194444
622753999092910446977164761511186723272386792081333673061819448493966071
233452718565202536436219641987827529788130600803131418170693144682211892
757849782810943677515407101063505537980038422190455084822393869932969266
592211127426981330623000734656284980936366930494468016285537126334126203
784919194986000972008367278766507868863069334189952257683143908324848863
403189401941610369798438333466086767094316436535384309121578155435128520
777208580989020995864496024794919706872307656871092343807195098248144731
578137800806393584187566550985013218828528401849814076907385073695353777
118803885289353476009303385986916082893354211477229365619072762646037260
272393209911878204070674122722581207667290400719242379303309721323641840
939561029959712917998282900095391473824378027790511120309545825328887211
461701334403859396540478061993332245473178034073409025121302172795957538
631581488103929524754109438805550983826276331276067181261710220113561818
007754002275167341441692164249731756213631285882819780057888324545345815
224349372681334339977105125320814783450671398350383329013139459864818202
723220433419309290119078328965692228783374973543015617228291156273294688
148532819221007523736268276431526857354932230280181014496490090155292486
383388856648930022509743436012008143651536253691994467097111269519667257
800618912154402224875646015546328120919458246535574320476442126507906552
082083379760714651275083204871652715774723258872757611283575921325539344
462894332581050286335836692918285668947362235082502949640657986308096143
416968304675951743553132243626642071976084590242630174733922252912483663
164280065528709750519975049130098594680710136023364401644001791886108532
307649917143720544678235972117604651532001630853363193515896458906817223
728123103202718979179512727996560536940321112428465909945563802154613161
062675216338056643943188812681994940055370686976218552318589211009634410
129335357339184596681975398342846968228894600763520316889220020219313183
697575569620611157743058263055358620156378912460312206729339926173783796
251509999354036487314232088739779689089083699962929953919772177965334212
492919783837514600620549673416628334873410110977705358980664981360113955
715843283087139405825352740560810115039079416880791972129331483030726386
786314110384431282159949368243429981887197686376044963425975242568861886
889789808883158650762626048564650043228968561492550639688114044004295038
942458723822335431010786915173283336047792627277656860761777056168740502
577437499837758301438561354272738385897741335269491654839297215195547935
789238667625027453701046609093824496266269353213037445388924792161611888
897020779104485631995148266308028795495464535838663073444237533197122791
588617072896520901498483054359832007713266534072906620167757064096901837
712013068232453334779666605253254908736019614803782415660712716503835822
572892157082093695109958901328594907243061833257552012080900071750220229
497428018234454137119162984499147222541965946822214682606449618392542496
709031040075814888579716722463228870164384039084638567311643081695373267
903031145836805750211196399056151691547085104597005420985717973180155647
414061723341458471112685479298924430013914682891036791792169786165824890
073220335913767065276765213071439853027609884780562169946596554613791749
856597392273794167264953778019920983554278661791231266993747307777305693
244301668393330115545155426568649374921286870491217542459678311329692484
924667442619990339728256748734602011504422287804661243201830161082321839
086547710423982285313165596856880052265714744288233175394565438819286244
326625033453881995900851052113831244918618026244321955404339857228413412
544094117717221568670862917421240531106205228429861992736294062088347548
536451281232796090972139537753600230767656942082199430346487833485444927
135394502245913343746649377016556057633846970629187257454265058794146301
766397604574743110815567470916527087481252671599137932405273046136939611
698925898083119063225107779285620719994594877006118010022961323045882945
584409524966111583428049086438608807964405577636918577437540258968559272
525145634043852178258905995539546274513854544529167610429692679708935800
562345019185714890304184957674008193599732187119574963570959678251710962
647520688908064076514458931328707674541696071079316927042851680934133110
463535062422098103632167719104207861621842137639381946256972867814136363
896201239769104654189568061973231484142245500716172158513213020306841760
872158927020988791089380810459033972765473264169168454456276007595613671
035845756490944306924525320850030910687831575615198475675691912847846546
925586651115579134612724253360836351313421839051771545112284644551360160
135132289485432715047608393075561009087860966638706122786902748318193316
067014849571630047052622282384062668184487883745481319943803876138301288
598852642019922861882084995886408885213525014576153964826474510259025307
431729568996364996157075518558371659353671254485150893629045677366300355
624573747791009879924991469672240414816012895309440154889426137831400878
043114317418580718261851490511387448313584390672289494082582860216502889
272283874264327861686903819605301558944594518087351972460082215293439808
282541261282571572093509853828007385604729109411840060844852353778335033
068619777245018863640703449733664731006020181287928869918618244184539689
947772594821691371336474704531729798092458443611296189975956962409718455
640205114325895918447249209429303016514887130798021023790655365251547802
980594075294405131458075515377948616358799011581920198088796949671874482
241568364635343261602426329347616344581638901638051238941845239734218414
968892623984896486420934098166814947711551770095626690298501015135375998
012725012419711198715265937474847789354887778151929311714311674447738829
410646150287513277094745047639228748906629898415402593508340351420351361
688192482389980277066669163421334243120545073593886166876911881857761181
357713324839652098820859823912986063868228047543624089565229214108598520
373305446259532613402348646892750605268937551484032985420869912210525970
056285767077025676953009789700464089200098521069802954196998021380532957
981594782899344432454915653278452238405512404452082264354206563133107029
407223715527705042634820739844548895892488613976570791454144276535845729
513297190919476944119109667974742626755909538320391696734942613600322630
774286841050400613510521944137781580950057145268460098103521092490400279
580507364369610212411377397171648695254931148050401265683512688295984139
832226763778045006265072417317573952197968907548251993292596498016270686
656580301788774056151671597319273204793762473755058550528396602945669925
221736008740812120142090710419375985717214313380174251415824918247109050
847159772494170493202541652393232332588515888933370971363108925715314177
619783260337501090262840664158013713593565292780884563059517700814439941
146742918503607488523666547448699280832305168157116029118363741479584921
008605289814695477508123388969431528610212027367470499039304170351713421
269234867005666275062290586369118822289031705103054068820969708755453293
694340639812976964780318254516421783473477164710584232385945801830527562
139101869976043058440686657123468696794560441557421000391797583489799358
827518815246759308789281592434921975453876683056846684207754098217812470
533545231947973989533201759886402810588255576980043971205383124594289573
776960018574973352499650135093689259580218638117259065064368821271568157
510217129007659927503702282839639629159732511734185867210234973177659694
542836255193715560091436803293119628425466284031424443706484323903749064
108113007928489557672434812000903098884572709077508736388732996425550504
738125289759629348228789176199207251383093882882925104168376227582040819
189336036538752841167857037209897188329869219278166296758445801749118091
196630481874341550677908639488314892415043004767045279712834822115222028
370628573142441078237925136450866775666228049772113971406216641163247567
842166129614771090188260946773776864061767214842938949766713801227889413
090265535110961183470125651975408070953840609168639369066737866272094294
342642604029021583173450037274625889926220498771211784055633484924903260
035085690993823927772974984135656148307882623633223683807098223460122742
413790364734517359252157547571609342709351929017239549214264906911152715
233381091240428121028937384881673589539345089306977155229891996989038858
832754090443003219868340034702712200201596993716906503305475770953987485
806700244910455048900617271891680313945280361656339415713346372225504775
474607560550241087643821216888489169403712589019484906853797222445620094
838194915327245022762185891695074057949837598210066044819965193601102615
769471762025717020486849146168940684041408335875621183192108380056321445
620189415059457800253187474719116048406779977654148306221790693308538751
292989830095802775541454350587689849441791365358916200987252220490551835
546037065331831767161107380097866252474886914760776644701471930744763024
116603356717655648744405779905319962716329720091094492492164560306188277
729477507647774464525863289191591074442523200829182095180210837003538813
309832158946086801279542247520719241346483349639150948130975414332442092
999307514810779190023461281223301617994299306188005334145506339321393396
468616164169552202164479954172431711657444713641977332048993650747678441
499295480730258564429423817876415064928783617679786771585107842357026402
133880188756019892340568684232155856285086455252583770106205322242449879
906252634840107743224881725586022333020763999338541520153438477254429178
951306370503204449177977523708719582779767996861136265322911186296311646
851599346606934605575459560631558300336976340002766851512938436388860908
283761411577320035275651587459065670254394379311048385713132944906049265
823631089495350900826731544972263966480886180415739778884728921746189741
897217007700098624496537590127270152276345108749069480122106849520630025
190116559635805524291802055869042596852610474128345184667369385800277002
529653563667216198836724282269339503259303909945831686655422346548570208
755046175205218537215672826799034181355206029998953664701065579005321295
413369244724922124363245230428951884617791223380696742339806948872705875
033892283950951352091231092581590069603951563677360671090505662996035718
764232479207528361608055976977787564767672105212223271848214844466312614
875842260926088757643317310232637688648225946912110323677375581221334705
568059580083101274816739620195835980239674144898672768458698193767837571
679367232130815861910459950589709910646869194634480385741438296295471313
721736698361845581445057486761243224515199433621829161914680260911217930
018647880500613516031443500761892134416024880917410512322903571792054979
279709245024799408426961588184426161637800447594782122408732041244211691
998055726491182436619218357147628914258057718717436880003241130087048193
739622950171430900984769272374988759386399425305953316078916188108635059
824445789427993465149159528848697574880258233535716778648268280511408854
297327881977657369660057277001625924043016886599468629837172705958098087
309018201209310034300587965526947880498092054843054676110346547480672906
743997636125924346377199958438628123919854702024148800768808188480878923
915913694632931132768493297772016466417275872591223547844808134333280500
877588552646861195769621722393086937957571658218524162043419723839899327
348034292623407223381551022091012629492497424232716988420232973032601617
905756731112354658902982983131151236076067739689981538122869996420146098
525797936912460163460887623212862056342159014791886321946596374834825642
916162785329482393132294402310432772887681395502133482663886874532592815
878545038909915619496324788550350902893909737189880039990261320158726786
378730956781096253110080544894188579835659020636806996431650339120299443
277267708693052407184165920700961392864019667257500870122181497331336958
096003697517649513500402859262492033981110149532275336218445007443315624
345324842179861083462613458975912348399707518542232816771872159568272432
459108290198863903697845426225669125427470560975679848571366236790238784
781612014779829390805131502581745237735295101652969345627861222411507835
877553733483727644398380820006672147400344663227769189369676128789834889
420946881023084270364528545049667596973188360444967028531906373969163579
809288657199353977234954867871804164014152814894437850362910715178052858
575839877111454742401564164771941163913549354667555935926088492005463846
854030280809364172505836533680934072253108208447235702268098269514261624
512040407115014487478561999228146645658939384880286438223138498523284523
606670458051136796637510392481633361732745472757756368109773445392758275
605974251607054686896577945305216023159398657809748015154149870977780787
053570580084723768924221897503127585271401731176212798987449584061998439
133656802977212087519349885044997139142851580323248230213406303125860726
245416377652345055220510863182853596585207081733927095664450114040551065
790550374177803933516583609045430477214222818168325396136349825252152322
576909202542164096574526180660517779015929028842405999988827536919575401
169546961522704012808575797661547221929256559639918209488946426575122887
663303021337463674492174493516371047257329808328127264681877593565842183
835947027920136639076897417389622525757826639908097926470114075803678505
993818871845600946958332707751261812820153910417739509182441375619999378
192403624695582359241714787027794484431087519018074141102903707060520851
629757983617542510416422448675773507563380188953792631833898559559565278
572279261555244947393636655339045286562154642883431622829211232904518422
125328881014158840616199391950422300598983499665694635801868167170748188
232158486477343867809115646607551753855522244285240494680336922999893007
839000206901215177406964285739301969105009882785230537976379402579689532
951124361667789105855572133817890899454539479159273749586002682378444868
720372434888346168562900978505324970369333619424398028823643235538082080
038757417109692897254998785662530488670330951505184521269449892515963920
794214526065085160523256148619382824898380008150853515646427617008320964
831179444019717801492133453359033366723767192297220699707660554824522474
169277746375221352017162317221376324456991540223954941582274189305899117
469317737765187358500323180144328839163742437958546956912217740989486115
155640466095650945381155209218637115186845625432750478705300069984231401
801694211091059254935961167194576309623288312712683285017603217716804002
496576741869271132155732700499357099423244163870892424275844076512155726
760379247653418089843126769411103131659514294793776706988812496434219332
874043904855382221608370889075982773901842041381978110258545370885867014
506235785139601099874760525354501004393530620724397099764451467909933814
489946446097809577319536049387349500268605645556932242296918156302939224
876064708734311663842054424896287602136502469918930401125131038350856219
080602708666048735858490017042009239297891939381251167984217881152092591
304355723216356608956035143838839390189531662743556099700156997802892363
62349895374653428746875                                                 

L(100000), the 100,000th Lucas number

Ln
Ln =
    58079784712681363560225068004602387220597828666126993481369551040014
180085634250846413212772503556474496035646834092151621184846502569215643
705785563183985951051271568672996604488909117823563040247200507184170529
281909059775547361728552133368146449785122325937213921078027674153231878
136009746094804504639290393032048794328647653870894017561209528488489920
649217449211806523692022271425674828185637648711204777410882242337120171
417851209726811197241221172404683907723121592431802259702892175490710329
635726903693224180274935452650762143273330747762180869427854869506878118
482886707127382883033106840622621893811836729047677208209383173859033251
736376749304450916136680946052798288241467947657393034883266148541335867
969183039898299755134152500712539394036248506774881062180638682321286095
321769620196006517811296291597400103220436316771196831269599635638657822
346682298074522140153826312847400694082615403941811427422791180940055767
532376630680770041459485325816843995825012690123665997316719367566815658
509949416964831697675245988938147812252657208544162688320158747100887524
548244804565134700402862597366497235484281374100840947550778691568234069
876041972842393561202615919282887522173999067350327241183789577198765682
428379125802686085284760067732899651472601566186816997489553599943772841
519980137033886850593190833871229810535737735309879294772077930750862195
349837037594236521483568603426721416840226455072757754936388527887500596
676439332410021643031670716381114562596118455787907955236282562136604016
442610985074024153099458901832696074462421729255100716028292275989831140
722104259645862184064861607385948928251983305307881264433643931889922362
005016696534582730422002794127275631853713803004364319694520783440528329
430042566874086681037676891483580233282299742113964587863312061011420661
149756796081236943613587130914744388810378742373405947955531436718737838
694401084160294485965960052462548732700221434687091897111612325589571464
620665848090614282492910416238978065993616081371411822921819379800198549
705644869192489850513972400521097797143560669000216782085007923664159276
399273661764092393367196080333758700143834321537822664514611905923118753
067139346549942182589034474440353582430248692502410360520166140513313614
078290643270921680859092690137974590149960174570214151867847426889208909
538371377955035741771724290575481746310985328542443017060506551850784769
613471038744406228920028268753505689803875589199272212896805910818792748
057309562064156556084644643140431379737491693564807884764801151902504208
080165717186380775226092158354587035167677750895524294446593732938442057
480091335194493822484104151392057402835123684669043069780820229494423972
894448201050314337031466591508729775713630407335638311486239364898907727
528154534884453823171394976630172070767742707838208738206021179229938403
154121848196868242396394993806845365063872869276675182125015469328404521
821404730549062815301217319586022796216132654305535029991353855397776837
428623810000723360444415772236634824016460130129066807843728731585148469
048714979470215198675022492898369904199348680598557239227330826215548006
436872401733339890877680041443233302650750346790228311579109072540092084
446029553805086170959661631972151540175703718118276671200287490670759030
891183800393840793292960373185742630136293955718765294694961179761489413
640164036877383885920946783126683814271667627138017976199246168039807836
491900065341316670414652747328764610087476911975831758907828858573257963
998935195102178724032303850644279598011641880163718150376076900610666615
309467654617881092626057920832364327873244387207052148291825322244715580
339548105657544019826504947710364625853350675447699231967320128776602235
800410670190026082820650167315093316724352381119765934044005480676891940
152966438683780241532962615908259078622977449661368993934881835575187578
984719324973630619628676488724345225067950838164641675408949529030883545
514128243690099135196596042174328734542016728461115666188119123137511291
274027240306446184374487609591866382187364705164278956014374908821553382
189485656329094552429994101776987958199714480062035196486347992468649979
763939804301851394645321168329040063146301418275449791446772803785418204
687020541250994379352101316909546445153495694810566261654454758899356532
478901456005196000597757639726639725284178223545411040667685164603938315
557931600912384234283785443695172098026362146409828225114747238012375659
677930840067874085163883998025907224238260818233788848520987368536361040
457333251121328710483693230298027714078353339503220367634623673447521918
563901369304390440842549277395416276396962805384768422559337211539106337
635325035157649688234464291101582545836144918255546020037657350188441439
326334602260505011069699370837935466184934620636366770762689771526226324
820749102619978553382333040460238182411420164651646799940557378857769395
270905709754796559236983238212481172314816291069553131973158399392096049
288962262891138003636062236066555964593064751983276813034836440616657750
803966728887827091431259746986708544639118861856791827975896852484745566
146377979701533769650431772503696515539270389899220303091083736199454540
111568574698769541836459522168286164837427691210909635329492698809409572
317770090561119204948970054492015771438839653800889861993346515700574951
078486297682866761979251420612403994201501942538536370474084183979290862
368936183432827954410938976795968292616881093798343780061438932892593773
193908169762781696864466739035822143056672326754740020958426457601496248
571022717927795218828315492439001598844726531821874172751506847330186041
692605213897937917864555400078507198623245492695372660650957111986844902
199561509650416342188546790750959178584748685728659099965783287710628744
300398755124928314979130660732157817619668188728841492516604448649799319
346762279036565754483987280231082055304109916523296783413615769506011317
970093607085779320859710877836456918401997527549213989714436143822723914
126039486610434635207956525909851863606636869107013044356630126900729601
732218760887891833763536064211067082022341075864142685424495461568177649
347476584804590763393479474737676675779831318978215650904741796026635414
351065159738014786765404978325529031044014478847092788395170278152626631
995220299338265185252782624932936737691934983175257906962500388783756398
623228378448247782575402471635930389353137200793129663129831176653701827
036052080106895468476900553914555874976702231109523682046177002659781714
402093179568077750126233412136413045866072370156881251120575091119960959
404677328783327544970846437607343376010710520344258509183233195638970991
724172518499638154051291941084677984880644195840877115328389648347075606
564859245629609588857789468870910351459536055896455941232961247768825418
041557267747044689970456598266035247967222374671708711981576732613685379
107462846955184489999059877961710295373550550862965671955392684990662072
570485670341844006575494693487497905042959148958358106702402241935063150
155612551495787719474191362983586304957420341439565200336094931504583266
496717052819922472490858335759573191736352915889362803710428745796977380
211427931378737323420803909359299454913050829691822080586600328908771013
986846089904869771378509712695572007842145820382185394125788673960027762
061074649259883176755813661285051551190127338534291350054096844694897001
377509133964235923503276189569144030326086763191410223965309032828186909
424092720284849826618589075702480064746863628015936114607552000736937018
276730243306365421031933439064414870664904383599396286230603906235900920
615295558401350425501303492656169659468555610106552904516113521111592438
799002330074267323320506692592352878043951154021674823035498055446706176
730170980112268699098075668730126949084020506948581580225212592135399375
573060121704149324040654360635977804336561218553584029932051446815203219
782720866863431463905814668042202009618848294879863783152360377492719524
342188106451768615257720864518051170554000130676981933729048026787018408
156966573183483758100567989810000983169145358335434167006350406754397078
537716821812800835289658215394781738671260335046376949633899165359594809
444574552007103202350444959617680199509283924602487571093612132913698343
882054515750629740120304778528467956942823557283329199317438800141827098
036681792696051837838818361203172450037709857281974037186370145125410108
934079272003949679829335153458813376654268565738454676054353343106356887
523268962363772751894246125824247489931114799356629086829236308619896963
076128361661347035852702820602632416012586010970219515560760258746206599
975206148977828912183515172557089131759439404923025559857675834947102279
324210588705903139703458253260091256176800955118365256472874191327922275
617816422443111506938063733344299867068671170982836287955096389068341713
761087898332894037165717743903658129751219415298903479695564391946190430
243037237810901478816494701926223755220937229824188227740076176112014312
969631269345954291292512155397845032794411940711726782203771677681800428
455546673460849328594335348133701337252515411699258454171721765901773493
636804718800172206719473510374127525924801125081641042513709406973347364
982576651711134187701290134161192148849812545279691879111480486268842439
325466181292950966909091487815684058175429165937514603446062073810114631
871723197561569564821808394337216889618776590191960262268033841903834555
193189555274053019342850013231227813783565267155028587938429555570621563
857196608494175776861313014631466634420982255757545583205629576174792283
990176140630987994177128748795283277081747700465522617713336485964944463
439738589064136122713224258306513431700513675825110491470089561676376835
119463845890139584539878717057958367684684991718037208995832218059046987
660755781727764621811445648156680480367361084842311882677522167322334366
779911234541604542786177815758623581409972574032665466404320678891283818
424366739723932744069380952081423357102450469260426621981546456049459764
846732661925808566282504488508881679152424321082904435445838107747577860
375646890069253133130192814603445905299941466262328666257009033930645250
474489404123191887867934436301468137462046160255396306953610491140062788
928617112915287462199954489345265846795319852585638584113478252240395034
451079217867042189392359871437590474827646718914237484409870726508014376
393826091926433537393713900294983409016383517160420012953819291874614098
433953454626451537219109610932173308743640688610140319567052014336473897
558517675030171487653671166667552982339641120486998437326762433281682450
587880194478847160212390644033544601365893500069122990642834389492284386
370272645430195713410881658328586222595084116533305287015833903994910017
081601638428154389812916844662154031041774243315813099814575703497096252
927761751743815466262926010344811253694587189875022740295605414695966950
080662060055851905973302515423728018997208551789524943181565000435400541
575937151079982717421350115717683336450839611001372555247841993017018203
375460525859725291009477911002433554465177369644881375719145166888197082
330530883580382290314090659022141737575623337525351438986211977867918805
937849158248993717322497654401996006425161700786730503156335128770780751
390205460529159304278549644747615724080593818891051602342981665373495450
788555152001497122155460803546173228578041748987471730948032955130602226
149837730784674455021370966260506323837156768155191786701081595015607713
108527989603126897571475747615957836091777929630431700888415807679363556
059413659506723145371990814630915321297038230208320611655019736694499933
870709245019098728487509067257572904351228845163699948019770422700078559
484814130059208729739831024426259294907011087270907950192445823232814867
302690024500848215305897148620004281471016512447883408045690588991762036
922744697679969082916979207359604409419948029431877430159724329967734386
080844470666547925374002949065128189590834446133508826075018175566010534
329388708910440393909926667136817016774583349767321995353863512514694590
059438231868755692792825609445371823248381718907426573027573967762180707
400040568429561980812697977052180575560634380720194540564203126797434257
096967691539240005303791503284349142907183159605667893282058252519380870
437982050970277871912659875799270938834205117769923907049472174943515361
902952032380430735191439405855004192296694013026907033327158186641352035
695867102003849849435946730181866042149600896582026441518186350483956494
866749574992106421461390292803659098702847771659912059862982619146276328
309600283621981650066313448220441677223603954409879575663267081194592177
443908890169107688077004921825907808728519201796494788469184087022802411
449573144099986092226981668295626048483140199026941139851754156926206864
859227210840079642349462369941164157052843353830901778236812642063635361
917304163252139337524156559476802866223382661450439828455085634576192386
653829410251279191235342771044293203594926482982479781791600665838798251
089635124498017143298187390763129298009115507060747781864336470129429921
107173822290212299933785800689644078797833810184290975085049246209671624
934017812289763389007729030435608013166147832925728954889406708182626313
978393392836265904674897593582061827920048232226432704689292941758438861
495608340666579892528889053105132090049865754486578371861368438386889562
817357263149568685941131831142005723481091443875779496487201631426265647
056707415148711722031590615180967611092452665147143837311288575190950870
911906898283178678497132187818150754406646605396432125158045754277665157
726423381453618810233115455288169868427387130373142589985800994506586180
516043608646240355822628142635658814700446440026155436425667685040059034
806883291543706169024453582617237081410356119849012994688241684527778205
930043721903282652008501938277007443477301301365499896583536449918751617
773114475105458180424214628297418907518961504347458744136544445869684810
240743992602882012861626686982677921739234150113577089919595943455644562
277060367576237275370512506596614428640314357513857038510726869428667450
822648645046898995470432818986695617082593950782881990548175825787935868
435826446013079315524032892159825409981122797653525536687164598961784292
412891641776560423581550679541420505143466172513176105522515066403943770
931172238654330365953061224333634366767244664096043795566857040136330946
910524210673247792385915549932003405150158304411738629162956325510543429
128037722896804279651761536823908521802547526248204509946586479912486382
104026919455341573542801018740215292276078693492473394222791203382785220
311854688529843624914337793032581477012263369912587589483470007240336671
062571686601372741729328807053502595796388608502840534500131393578125356
973752888405504752164444837964136792706157540602419170253139695808331428
852803316471545688594852347534106388674010254678097513365296325541394891
670547905134595895244593168897676054717619628445115741176534925984597710
522619428442462910795889936016391105110025177612769906740610097111328867
918095572306323885394674758621199578991908660222726427418879586719674486
013750602333923726571045099507047272446087246700205750087857015926662221
582240732469002475671152431813391916821485285554227210470717734439529145
253954979543299917776671208893204989241281984998680148969359611626255320
565135001656226492316591460096477452801993960973186785893895807034117185
315787481484365158029313381589795077525845632388090706369661744736739296
939686455577101399880163447099706789954013818644438998679618908989040091
050067437539000449222010612973576656672368723628072463323152525638692089
629440821015104941992567553505916467213288338507372603194883595181375412
640612473840933710569491831034457722754177492684757763450863313912811177
059170395958815828267412102729825431309981708118531959367841806906419021
799990746064566988230290674113985053672351504319059464964392637644032713
690071152108203765792275640812030019248014433114748711540552906228060418
114984129272288530971592231506206795548254696311721944069604471858046583
226729285002296896289020240217693783614535339849307676372872211226740490
798913954412204558019666092383204277193383057698696281977218812426999150
431128882561009835078745713487048616117143874512308858951313678540178473
358792272686609586764849857414394943964855334846054949282548120715975467
025336233197081769633759989677893523431002981297102108078632792595443728
306920387267676042454699157464866900325783273795527344686964722354254664
498866398750560343170418543872215959391746238143813683524147766111382778
096590524474807981914314607857362555253812746306047305797572066918061827
167995721275979792872710326917293100816121476459559064595674265120642179
343724625780210014893639448937192058688341704333606205611532660916781648
494504455778612885268415894527358809061224276123247386787247521432123990
701729413936289401899160412954713504909009370646481168368448079340528860
449743905755530050758362195938829169504108490599083435867368132484997708
417317271294005874774613267052138564563147717116538628800175268336586302
917717003537702442837925149808303573742701483968988035105201191264310865
867787221248915547995257084559595180175595276719948780327225563125733090
847763948221660920132167100090213579143072894402293061647132513761013597
504232360945310318566368387926967581691731962421638873726733174490463898
003305712546371870145238828595132196073294863769271340919180466996668725
677944411931796776278531100743353792111418365843092198101392562730178276
860727819138942327922747638289906137716905141015910858529914393931839144
645559129251156928989461017941618779243149740341277122677109869540289710
959871275594883140697060842692462402204621494784231635797451986719481837
052822232196607587137983962156136084590575282163519404201122338249257649
216253924577434707008131768601601431024301332072073759084090924320380827
835605952336307933290611276930308931892699226229612385080521496154691850
312299741059797979456944796376944506738567571200048410926809002423531361
377019967944704418962896172109372221030727191400683232704423369438946975
920562077327528671746585556770068701367390604606850364692827051515448826
020133499018697258139470386968837820833587595204387235814612664173340581
794533545528382490777061666275382424705332849049192858760849016164425106
092945988978040397666515632255986776937776499562564261049791504050847157
011861790066253439542995065638064058383772953374315112952344395146934857
820157646700413469530211294303070678796056293460140058576910428407506227
987459767852623290789817704553129365316941634678399547913743645368318916
481261376601474331558419737496750580395984807814789426303237835781937632
645476104161265510262042770664466394002731078142589963683481209050662186
434915126596953632069492589772947977039248742908572390240397308367971148
948278046937075750271212895911626525641305064579043815205603287295413397
659671220453593036565215736156361743040741498510459069892048183364879220
656893945864497776773527331926807871907491769684294858842807813360485555
101623031540294376279125428187379242128213825589926668352587488239287735
615692612113764291047221185351022390755068497686078024526684690116266840
224575010173423493194290266464410795190969789429286128482052284114501970
350820914076607992443259885075821953664308320748147628303530185047161728
525262223625882373364391967307076097506110796881115838157874036496179431
797254899042919372620019189953164283846253667235383608576835106133813707
671457308859799262199263385776068191759400278416670603957206727889967467
328684291715513607352345472407177842176366554393310690256962298423074590
044140571447948005730532471148216349072097401456895332682070995931553421
736440544770957067256658164538650780364335506798781127759337550940905813
502615211197662886490502696681342107350645560234007104245117342795485590
895490524163454471908919968941239442202760355939084715149597891109017883
520312600017153662094943432384684811614374681376044285178118955091936151
249620184429869700084435298574079086058565368698291770225434941079988835
717155862157376364068529666381511626842107524133965371324878843791020296
347239600263461408548293821660509918289999066265681827296236484843081021
882279608203554758523983034836058667446584335931660888206743750144359854
261841410199753556467357919675941809856290675833410868781938067902189977
540911040685701424377040717482429725926131902120763721100475155672534630
273294581058834261642158781096842237671425827317966991328109342282019147
215255893176136620741277836581862273983222637723020246583040330809856468
000072895117697824013063063874913749127338505842909128238535262612405532
562987692480186169008830464694921241300064653168928299011289107146321577
501911965250511762607684877517000433798401657419013263665448769664623956
114597574734812593620084728934280170836051157192439383443023510264068441
867266186385229308491919684032163205801766426324019511655072244660457801
213172511965312038825663642430040821613523705577387012399443947875161095
578866214408830777341331569057888612466788287625956582244221984888408364
231815786779844768787180610820221612035297648215970593038474024513794778
57726676584685986328127