Project Euler – Problem 66 Solution

Yan Cui

I help clients go faster for less using serverless technologies.

This article is brought to you by

Is your CI build step taking too long? Try Depot for free today and experience up to 40x faster build speed!

Unlock faster CI for FREE

Problem

Consider quadratic Diophantine equations of the form:

x2 – Dy2 = 1

For example, when D=13, the minimal solution in x is 6492 – 13×1802 = 1.

It can be assumed that there are no solutions in positive integers when D is square.

By finding minimal solutions in x for D = {2, 3, 5, 6, 7}, we obtain the following:

32 – 2×22 = 1

22 – 3×12 = 1

92 – 5×42 = 1

52 – 6×22 = 1

82 – 7×32 = 1

Hence, by considering minimal solutions in x for D <= 7, the largest x is obtained when D=5.

Find the value of D <= 1000 in minimal solutions of x for which the largest value of x is obtained.

Solution

// define function to get the continued fractions, a0,a1,a2,...an
let continuedFraction D =
    // define functions for working out the nth term of P, Q and a sequence
    let getPn an' qn' pn' = an' * qn' - pn'
    let getQn D pn qn' = (D - pown pn 2) / qn'
    let getAn a0 pn qn = (a0 + pn) / qn

    // work out the initial terms
    let a0 = bigint(sqrt(double(D)))
    let p1, q1 = a0, D-a0*a0
    let a1 = getAn a0 p1 q1
    let initial = (p1, q1, a1)

    Seq.unfold (fun (pn', qn', an') ->
        let pn = getPn an' qn' pn'
        let qn = getQn D pn qn'
        let an = getAn a0 pn qn
        Some((pn', qn', an'), (pn, qn, an))) initial
    |> Seq.map (fun (pn, qn, an) -> an)
    |> Seq.append [a0]

// define function to get the continued fraction convergents
// e.g. for D = 7: 2/1, 3/1, 5/2, 8/3, ...
let continuedFractionConvergents D =
    let getN an n' n'' = an * n' + n''

    // work out the initial terms
    let fractions = continuedFraction D
    let a0 = Seq.head fractions
    let p0, p1 = a0, a0 * (Seq.nth 1 fractions) + 1I
    let q0, q1 = 1I, Seq.nth 1 fractions
    let initial = (p1, q1, p0, q0)

    Seq.scan (fun (pn', qn', pn'', qn'') an ->
        let pn = getN an pn' pn''
        let qn = getN an qn' qn''
        (pn, qn, pn', qn')) initial (fractions |> Seq.skip 2)
    |> Seq.map (fun (pn, qn, pn', qn') -> (pn, qn))
    |> Seq.append [(p0, q0)]

let answer =
    [1I..1000I]
    |> List.filter (fun d -> sqrt(double(d)) % 1.0 <> 0.0)
    |> List.maxBy (fun d ->
        continuedFractionConvergents d
        |> Seq.filter (fun (x, y) -> x*x - d*y*y = 1I)
        |> Seq.map fst
        |> Seq.head)

After spending much time reading the information on Pell’s equation on mathworld and wikipedia, it turns out that the best way to solve the equation for a given D is to work out the continued fraction convergents to the square root of D.

For example, for D = 7, the equation becomes:

image

The fundamental solution (the first pair of x, y which satisfies the equation) can be obtained by going through the sequence of convergents for the square root of 7:

image image
image image
image image
image image
image image

So how do we workout the values of the convergents? The continued fraction for a number is usually written as

image

The first ten values generated by the continued fraction for the square root of 2 are:

1, 3/2, 7/5, 17/12, 41/29, 99/70, 239/169, 577/408, 1393/985, 3363/2378, …

These values are called continued fraction convergents, getting closer to the true value of square root 2 with each step..

So, to find the convergents to the square root of D, we first need to find the values of a:

clip_image002[8]

clip_image004

clip_image002[10]

clip_image002[6]

Once we have the sequence of values for a, we can then find out the convergent values:

clip_image002[12]

clip_image004[6]

The rest is simple, for each value of D between 1 and 1000, go through the convergents to find the fundamental solution and return the max value of x in all the fundamental solutions.

Whenever you’re ready, here are 4 ways I can help you:

  1. Production-Ready Serverless: Join 20+ AWS Heroes & Community Builders and 1000+ other students in levelling up your serverless game. This is your one-stop shop for quickly levelling up your serverless skills.
  2. Do you want to know how to test serverless architectures with a fast dev & test loop? Check out my latest course, Testing Serverless Architectures and learn the smart way to test serverless.
  3. I help clients launch product ideas, improve their development processes and upskill their teams. If you’d like to work together, then let’s get in touch.
  4. Join my community on Discord, ask questions, and join the discussion on all things AWS and Serverless.

1 thought on “Project Euler – Problem 66 Solution”

  1. Pingback: Can you help me solve this math questions? | Round Cutting Board

Leave a Comment

Your email address will not be published. Required fields are marked *