Project Euler — Problem 66 Solution


Con­sid­er qua­drat­ic Dio­phan­tine equa­tions of the form:

x2 – Dy2 = 1

For exam­ple, when D=13, the min­i­mal solu­tion in x is 6492 – 13x1802 = 1.

It can be assumed that there are no solu­tions in pos­i­tive inte­gers when D is square.

By find­ing min­i­mal solu­tions in x for D = {2, 3, 5, 6, 7}, we obtain the fol­low­ing:

32 – 2x22 = 1

22 – 3x12 = 1

92 – 5x42 = 1

52 – 6x22 = 1

82 – 7x32 = 1

Hence, by con­sid­er­ing min­i­mal solu­tions in x for D <= 7, the largest x is obtained when D=5.

Find the val­ue of D <= 1000 in min­i­mal solu­tions of x for which the largest val­ue of x is obtained.


// define function to get the continued fractions, a0,a1,a2,
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
    |> (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)
    |> (fun (pn, qn, pn', qn') -> (pn, qn))
    |> Seq.append [(p0, q0)]

let answer =
    |> 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)
        |> fst
        |> Seq.head)

After spend­ing much time read­ing the infor­ma­tion on Pell’s equa­tion on math­world and wikipedia, it turns out that the best way to solve the equa­tion for a giv­en D is to work out the con­tin­ued frac­tion con­ver­gents to the square root of D.

For exam­ple, for D = 7, the equa­tion becomes:


The fun­da­men­tal solu­tion (the first pair of x, y which sat­is­fies the equa­tion) can be obtained by going through the sequence of con­ver­gents for the square root of 7:

image image
image image
image image
image image
image image

So how do we work­out the val­ues of the con­ver­gents? The con­tin­ued frac­tion for a num­ber is usu­al­ly writ­ten as


The first ten val­ues gen­er­at­ed by the con­tin­ued frac­tion 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 val­ues are called con­tin­ued frac­tion con­ver­gents, get­ting clos­er to the true val­ue of square root 2 with each step..

So, to find the con­ver­gents to the square root of D, we first need to find the val­ues of a:





Once we have the sequence of val­ues for a, we can then find out the con­ver­gent val­ues:



The rest is sim­ple, for each val­ue of D between 1 and 1000, go through the con­ver­gents to find the fun­da­men­tal solu­tion and return the max val­ue of x in all the fun­da­men­tal solu­tions.