Google Code Jam 2008: Round 1A Question 3

At Google Code Jam 2008 round 1A, there is problem:

Calculate last three digits before the decimal point for the number (3+sqrt(5))^n

n can be big number up to 1000000.
For example: if n = 2 then (3+sqrt(5))^2 = 27.4164079... answer is 027.
For n = 3: (3+sqrt(5))^3 = 3935.73982... answer is 935.

One of the solution is to create matrix M 2x2 : [[0, 1], [-4, 6]] than calculate matrix P = M^n, Where calculation preformed by modulo 1000. and the result is (6*P[0,0] + 28*P[0,1] - 1) mod 1000.

Who can explain me this solution?

I'll present a method to solve this problem without even understanding the solution.

Assuming that you are familiar with the fibonacci numbers:

ghci> let fib = 0 : 1 : zipWith (+) fib (tail fib)
ghci> take 16 fib
[0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610]

And are also familiar with its closed form expression:

ghci> let calcFib i = round (((1 + sqrt 5) / 2) ^ i / sqrt 5)
ghci> map calcFib [0..15]
[0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610]

And you notice the similarity of ((1 + sqrt 5) / 2)ⁿ and (3 + sqrt 5)ⁿ.

From here one can guess that there is probably a series similar to fibonacci to calculate this.

But what series? So you calculate the first few items:

ghci> let calcThing i = floor ((3 + sqrt 5) ^ i)
ghci> map calcThing [0..5]
[1,5,27,143,751,3935]

Guessing that the formula is of the form:

thing_n = a*thing_n-1 + b*thing_n-2

We have:

27 = a*5 + b*1

143 = a*27 + b*5

We solve the linear equations set and get:

thing_n = 4*thing_n-1 + 7*thing_n-2 (a = 4, b = 7)

We check:

ghci> let thing = 1 : 5 : zipWith (+) (map (* 4) (tail thing)) (map (* 7) thing)
ghci> take 10 thing
[1,5,27,143,761,4045,21507,114343,607921,3232085]
ghci> map calcThing [0..9]
[1,5,27,143,751,3935,20607,107903,564991,2958335]

Then we find out that sadly this does not compute our function. But then we get cheered by the fact that it gets the right-most digit right. Not understanding why, but encouraged by this fact, we try to something similar. To find the parameters for a modified formula:

thing_n = a*thing_n-1 + b*thing_n-2 + c

We then arrive at:

thing_n = 6*thing_n-1 - 4*thing_n-2 + 1

We check it:

ghci> let thing =
        1 : 5 : map (+1) (zipWith (+)
          (map (*6) (tail thing))
          (map (* negate 4) thing))
ghci> take 16 thing == map calcThing [0..15]
True

Just to give an answer to a very old question:

Thanks to yairchu i've got the idea to reread the prove of Binet's formula on the wikipedia page. It's there not really that clear, but we can work with it.

We see on the wikipedia page there is a closed form with 'computation by rounding': F_n = ⌊φ/√5⌋ⁿ. If we could replace the φ/√5 with 3 + √5 (call the latter x). We could compute the floor of xⁿ fairly easily, especially mod 1000, by finding the nth term in our freshly constructed sequence (this is the analogon of F (later we will call this analogon U)).

What sequence are we looking for? Well, we'll try following the prove for the Binet's formula. We need a quadratic equation with x as a root. Let's say x² = 6 x-4 this one has roots x and y := 3 - √5. The handy part is now:

Define U_n (for every a and b) such:

U_n = a xⁿ + b yⁿ

by definition of x and y you can see that

U_n = 6 U_n-1 - 4 U_n-2

Now we can choose a and b freely. We need U_n to be integers so I propose choosing a=b=1. Now is U₀ = 2, U₁ = 6, U₂ = 28...

We still need to get our 'computation by rounding'. You can see that yⁿ < 1 for every n (because y ≅ 0.76 < 1) so U_n = xⁿ + yⁿ = ⌈xⁿ⌉.

If we can compute U_n we can find ⌊xⁿ⌋, just subtract 1.

We could compute U_n by it's recursive formula but that would require O(n) computation time. We can do better!

For computing such a recursive formula we can use matrices:

⌈ 0 1⌉ ⌈ U(n-1) ⌉     ⌈ U(n) ⌉
⌊-4 6⌋ ⌊  U(n)  ⌋  =  ⌊U(n+1)⌋

Call this matrix M. Now does M*(U(1), U(2)) compute (U(2), U(3)). Now we can compute P = M^n-1 (notice that I use one less than n, you can see that this is right if you test the small cases: n=0, n=1, n=2) P*(6,28) gives us now the nth and (n+1)th term of our sequence so:

(P*(6,28))₀ - 1 = ⌊xⁿ⌋

Now we can take everything mod 1000 (this is simplifying the calculations (a lot)) and we get the desired result in computation time O(log(n)) (or even better with the computational wonders of powers of matrices (over a cyclic finite field)). This explains the very weird looking solution, I guess.

I don't know how to explain that, but the auther of the problem have compose this analysis.