Google Code JamThe "Numbers" problem in Round 1A is considerably more devious than either Minimum Scalar Product, or Milkshakes problems. The problem is stated simply but the solution requires some innovative thinking to get around a precision problem. The fact that the exponent n can be up 2 x 109 and √5 is irrational means that sufficient precision will be lost resulting in incorrect answers, or a stack overflow.

How do we know √5 is irrational? It can easily be proved using Proof by Contradiction.... Assume √5 is rational. In this case, by definition, it can be expressed as a quotient of 2 integers:





 

From number theory, we know that the fundamental theorem of arithmetic says that every positive composite integer has a unique prime factorization.  This theorem also implies that you can't have both an even and odd number of prime factors. It is for this very reason that the theorem is also referred to as the Unique Factorization Theorem.

If a and b are in lowest terms (as supposed), their squares would each have an even number of prime factors. 5b² has one more prime factor than b², meaning it would have an odd number of prime factors. This contradiction forces the supposition wrong, so √5 cannot be rational. It is, therefore, irrational.

The irrationality of √5 is the heart of the precision problem. If it were an integer value then we could easily apply a modulus 1000 to interim calculations to truncate unnecessary precision and avoid overflow. So one strategy which helps remove the irrational component is the use of conjugates. Consider the following:





From (3) we note that  ω = r - ω'.  Now the problem becomes:



From this we can draw several observations. The first is that r is an integer, which is to be expected - that's why we introduced the conjugate in the first place. To see this you only need to examine the binomal expansion of the terms and you'll note that the even exponents of both √5 and -√5 produce integers, and the odd exponents of √5 and -√5 produce irrational components that cancel each other out since there is both a positive and negative value for each term.





The second observation is that (3 - √5) = (3 - 2.2306) < 1  which implies that ω' < 1. Since this term is less than 1 (and gets smaller as n gets larger) we can effectively ignore it since we are only interested in the first three digits in front of the decimal place. We have now reduced our problem to an integer calculation that won't run into any precision issues:



What is now needed is a way to calculate r:



We know the result is an integer but the irrational √5 component is still in the mix. The clue to solving this is given above in the binomial expansion of r. Looking at that it should be clear that the expansion of the first term in r, will produce a result in the following form:



Here the subscript n is used to denote the resulting co-efficients of the rational and irrational components when the original term, (3 + √5), is raised to the power n. When n equals 0 or 1 the values of the coefficients, x and y, are easy to determine. For n=0, the result must equal 1 so clearly x0=1 and y0=0. Since we know x0, y0 it might be worthwhile deriving a recursive solution for higher values of n. To do this we can multiply the nth power by the original term and examine the results:






Equation (5) tells us that the co-efficients for the (n+1)-th power are a linear combination of the co-efficients for the n-th power. Since we know x0 and y0 we can solve for any value of n by using recursion. That is, repeatedly substituting values into (5) to get the desired xn and yn values. Here is the general equation in matrix form:




The astute reader will notice that we've only dealt with the first term in equation (4) above. However, the second term in (4) merely involves swapping the sign of the yn component and in fact when the 2 terms are added together the yn coefficients cancel out the irrational component leaving only the rational component (as shown earlier). Consequently,



So an elegant way to solve do this matrix multiplication problem quickly is to use fast exponentiation which requires far fewer (log n) operations than simple multiplication repeated n times. The difference between these two approaches is considerable when when n is very large.

Finally, here's some C# code to do this:

using System;
using System.IO;

namespace GoogleCodeJam.Numbers
{
    class Program
    {
        static void Main()
        {
            string basePath = @"";
            StreamReader infile = new System.IO.StreamReader(basePath + "large.txt");
            StreamWriter outfile = new System.IO.StreamWriter(basePath + "output.txt");

            int testCases = Int32.Parse(infile.ReadLine());
            for (int caseNo = 1; caseNo <= testCases; caseNo++) // note this is 1-based
            {
                // read in the input data
                uint n = UInt32.Parse(infile.ReadLine());

                int[,] coeffs = new int[2, 2];
                coeffs[0, 0] = 3;
                coeffs[0, 1] = 5;
                coeffs[1, 0] = 1;
                coeffs[1, 1] = 3;

                int[,] result = MatrixPow(coeffs, (int)n); ;

                double answer = (2 * result[0,0] + 999) % 1000;

                // write out the results
                outfile.WriteLine("Case #{0}: {1}", caseNo, answer.ToString("000"));
            }

            infile.Close();
            outfile.Flush();
            outfile.Close();
        }

        // sqare matrix exponentiation (using fast exponentiation)
        static int[,] MatrixPow(int[,] a, int n)
        {
            if (n == 1)
                return a;

            if (n % 2 == 0)
                return MatrixPow(MatrixMultiply(2, a, a), n / 2);
            else
                return MatrixMultiply(2, MatrixPow(MatrixMultiply(2, a, a), n / 2), a);
        }

        // square matrix multiplication
        static int[,] MatrixMultiply(int size, int[,] m1, int[,] m2)
        {
            int[,] result = new int[size, size];

            for (int i = 0; i < size; i++)
            {
                for (int j = 0; j < size; j++)
                {
                    result[i, j] = 0;
                    for (int k = 0; k < size; k++)
                    {
                        checked
                        {
                            result[i, j] += (m1[i, k] * m2[k, j]) % 1000;
                        }
                    }
                }
            }
            return result;
        }
    }
}
Bookmark and Share