Question
Consider quadratic Diophantine equations of the form:
\begin{aligned} x^2-Dy^2=1 \end{aligned}
For example, when D=13, the minimal solution in x is 6492 - 131802 = 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:
\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 x for D \leq 7, the largest x is obtained when D=5.
Find the value of D \leq 1000 in minimal solutions of x for which the largest value of x 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.201s
user 0m0.201s
sys 0m0.000s