Project Euler Problem 66 Solution

Question

Consider quadratic Diophantine equations of the form:

x2Dy2=1\displaystyle \begin{aligned} x^2-Dy^2=1 \end{aligned}

For example, when D=13D=13, the minimal solution in xx is 6492131802=16492 - 131802 = 1.

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

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

322×22=1223×12=1925×42=1526×22=1827×32=1\displaystyle \begin{aligned} 3^2 - 2 \times 2^2 &= 1 \\ 2^2 - 3 \times 1^2 &= 1 \\ 9^2 - 5 \times 4^2 &= 1 \\ 5^2 - 6 \times 2^2 &= 1 \\ 8^2 - 7 \times 3^2 &= 1 \\ \end{aligned}

Hence, by considering minimal solutions in xx for D7D \leq 7, the largest xx is obtained when D=5D=5.

Find the value of D1000D \leq 1000 in minimal solutions of xx for which the largest value of xx is obtained.

Haskell

import Data.Function (on)
import Data.List (maximumBy)

isInteger :: Double -> Bool
isInteger f = f - (fromIntegral (floor f)) < 0.0000001

isSquare :: Integer -> Bool
isSquare n = isInteger (sqrt $ fromIntegral n)

rationalize :: [Integer] -> (Integer, Integer)
rationalize = foldr (\x (n, d) -> (x*n + d, n)) (1, 0)

convergents :: Integer -> [Integer]
convergents s | isSquare s = []
              | isSquare (s-1) = (a 0) : repeat (fromIntegral $ 2 * (floor $ sqrt $ fromIntegral (s-1)))
              | otherwise = map a [0..]
              where m = (map m' [0..] !!)
                    m' 0 = 0
                    m' n = (d (n-1))*(a (n-1)) - (m (n-1))
                    d = (map d' [0..] !!)
                    d' 0 = 1
                    d' n = (s - (m n)^2) `quot` (d (n-1))
                    a = (map a' [0..] !!)
                    a' 0 = floor $ sqrt $ fromIntegral s
                    a' n = floor $ (fromIntegral ((a 0) + (m n))) / (fromIntegral (d n))

solve :: Integer -> Integer
solve d = head $ [x | n <- [1..], let (x, y) = rationalize $ take n $ convergents d, x^2 - d*y^2 == 1]

main :: IO ()
main = print $ maximumBy (compare `on` solve) $ [1..1000]
$ ghc -O2 -o diophantine diophantine.hs
$ time ./diophantine
real   0m0.312s
user   0m0.308s
sys    0m0.000s

Questions? Comments? Send me an email: z@chdenton.com