Project Euler problem 66 and continued fractions

June 19th, 2009 at 2:49 pm

Problem 66 is one of those problems that make Project Euler lots of fun. It doesn’t have a brute-force solution, and to solve it one actually has to implement a non-trivial mathematical algorithm and get exposed to several interesting techniques.

I will not post the solution or the full code for the problem here, just a couple of hints.

After a very short bout of Googling, you’ll discover that the Diophantine equation:

$$x^2-Dy^2=1$$

Is quite famous and is called Pell’s equation. From here, further web searches and Wikipedia-reading will bring you to at least two methods for finding the fundamental solution, which is the pair of x and y with minimal x solving it.

One of the methods involves computing the continued-fraction representation of the square root of D. This page is a must read on this topic, and will help you with other Euler problems as well.

I want to post here a code snippet that implements the continued-fraction computation described in that link. Its steps follow the Algebraic algoritm given there:

def CF_of_sqrt(n):
    """ Compute the continued fraction representation of the
        square root of N.

        The first element in the returned array is the whole
        part of the fraction. The others are the denominators
        up to (and not including) the point where it starts
        repeating.

        Uses the algorithm explained here:
        http://www.mcs.surrey.ac.uk/Personal/R.Knott/Fibonacci/cfINTRO.html
        In the section named: "Methods of finding continued
        fractions for square roots"
    """
    if is_square(n):
        return [int(math.sqrt(n))]

    ans = []

    step1_num = 0
    step1_denom = 1

    while True:
        nextn = int((math.floor(math.sqrt(n)) + step1_num) / step1_denom)
        ans.append(int(nextn))

        step2_num = step1_denom
        step2_denom = step1_num - step1_denom * nextn

        step3_denom = (n - step2_denom ** 2) / step2_num
        step3_num = -step2_denom

        if step3_denom == 1:
            ans.append(ans[0] * 2)
            break

        step1_num, step1_denom = step3_num, step3_denom

    return ans

As I said, this still isn’t enough to solve the problem, but with this code in hand, the solution isn’t too far. Read some more about Pell’s equation and you’ll discover how to use this code to reach a solution.

It took my program ~30 milliseconds to find an answer to the problem, by the way. It took less than a second to solve a 10-times larger problem (for D <= 10000), so I believe it to be a pretty good implementation.

Related posts:

  1. Project Euler problem 69
  2. Project Euler problem 12
  3. Project Euler problem 83 – how creeps can help
  4. Project Euler problem 26
  5. the rule of 72

4 Responses to “Project Euler problem 66 and continued fractions”

  1. ShayNo Gravatar Says:

    Hey Eli,

    My name is Shay Levy. I got hooked to the Euler project after being introduced to it a few weeks ago. In my initial research regarding problem 66 I found your website and was really impressed. It seems like I will spend some time in your website for other purposes also.
    I do have a question regarding problem 66. I implemented this problem using java. I wrote the method that implements Pell’s equation and it seems to work fine only when the solution for x is relatively small. For D=778 i got x that equals to 2.178548422E9 and this was the largest x that I found but the answer D=778 is wrong. I read in some forum that the x solution for D=991 was 379516400906811930638014896080.
    All my variables are of type double, Do you think this is what might cause the problem?

    Hope I was clear.

    Thanks,
    Shay

  2. elibenNo Gravatar Says:

    @Shay: yes, this can cause the problem. The precision of floats and doubles is limited. When numbers become large, the exponent has to be increased and thus precision is lost. Try to re-implement it with some big-integer library and compare the numbers you get. You can also find out where the answers of doubles start deviating from the integer (hint: this can discover the length of the mantissa in the floating point representation of double on your Java installation).

  3. John BullNo Gravatar Says:

    What language is that, or is it pseudo-code?

  4. elibenNo Gravatar Says:

    @John: it’s Python. Yes, it indeed looks like pseudocode sometimes :-)

Leave a Reply

To post code with preserved formatting, enclose it in `backticks` (even multiple lines)