Google Code 2008 Cork: Round 1A Question 3

Question

In the Google Code Jam 2008 round 1A, there is a problem :

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

n may be a large number up to 1,000,000.
For example: if n = 2, then (3 + sqrt (5)) ^ 2 = 27 .4164079 ... the answer is 027.
For n = 3: (3 + sqrt (5)) ^ 3 = 3 935 .73982 ... the answer is 935.

One solution is to create a matrix M 2x2: [[0, 1], [-4, 6]] than to calculate the matrix P = M ^ n, where the calculation was performed using module 1000. and the result (6*P[0,0] + 28*P[0,1] - 1) mod 1000.

Who can explain this decision to me?

Answer 1

I will introduce a method to solve this problem without even understanding the solution.

Assuming you are familiar with 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 also familiar with the expression of a closed form:

 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 will notice the similarities between ((1 + sqrt 5) / 2) ⁿ and (3 + sqrt 5) ⁿ .

From here one can guess that there probably exists a series similar to the fibonacci for calculating this.

But what series? So, you are calculating the first few elements:

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

Assuming that the formula is:

item _n = a * item _n-1 + b * item _n-2

We have:

27 = a * 5 + b * 1

143 = a * 27 + b * 5

We solve a set of linear equations and get:

item _n = 4 * item _n-1 + 7 * item _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 that, unfortunately, this does not compute our function. But then we are pleased with the fact that he has the right to the correct number. Not understanding why, but encouraged by this fact, we are trying something like that. To find parameters for a modified formula:

thing _n = a * thing _n-1 + b * thing _n-2 + c

Then we come to:

item _n = 6 * item _n-1 - 4 * item _n-2 + 1

We check:

 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 

Answer 2

I do not know how to explain this, but the author of the problem compiled this analysis .

Answer 3

Just to answer a very old question:

Thanks to yairchu , I got the idea to re-read the confirmation of Formula Binet on the wikipedia page. This is not so clear, but we can work with him.

We see a closed form on the wikipedia page with "rounding off": F _n = ⌊φ / √5⌋ ⁿ . If we could replace φ / √5 with 3 + √5 (let us call the last x). We could easily calculate the floor x ⁿ especially mod 1000 by finding the nth term in our freshly prepared sequence (this is an analogue of F (we will call this analogue U later).

What sequence are we looking for? Well, we try to follow the proof of the Binet formula. We need a quadratic equation with x as the root. Let say that x ² = 6 x-4, this one has roots x and y: = 3 - √5. Now the convenient part:

Define U _n (for any a and b) such:

U _n = ax ⁿ + by ⁿ

by the 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 suggest choosing a = b = 1. Now U ₀ = 2, U ₁ = 6, U ₂ = 28 ...

We still need to get our "rounding calculations". You can see that y ⁿ 1 for every n (because y ≅ 0.76 <1), therefore U _n = x ⁿ + y ⁿ = ⌈x ⁿ ⌉.

If we can calculate U _n , we can find ⌊x ⁿ ⌋, just subtract 1.

We could calculate U _{n by a} recursive formula, but this would take O (n) to compute. We can do better!

To calculate 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 M * (U (1), U (2)) compute (U (2), U (3)). Now we can calculate P = M ^n-1 (note that I use one less than n, you can see that this is correct if you check small cases: n = 0, n = 1, n = 2) P * (6 , 28) gives us the nth and (n + 1) -th term of our sequence as follows:

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

Now we can take all mod 1000 (this simplifies the calculations (a lot)), and we get the desired result while calculating O (log (n)) (or even better with the computational wonders of matrix degrees (over a cyclic finite field)). This explains a very strange promising solution, I think.