Problem 66 -- Haskell, Continued Fractions, and Precision

In theory, problem 66 shouldn't have been very hard for me. I am currently taking a class on Continued Fractions, and Problem 66 involves finding solutions to a Pell's Equation, which involves computing convergents of Continuted Fractions. And indeed, implementing the algorithm for computing continued fraction expansions was not all that hard for me, as I had done it in Python as a way to simplify my work for a problem set just a couple days earlier.

However, for this problem, I found that in order to apply the algorithm accurately, more than the standard amounts of Double and Int precision were necessary. So, I had to use a FixedPoint package to give me access to 128-bit Integers and 256-bit fixed fractions. The use of these custom data types increased the runtime to about 2 minutes on my machine, but it's still fairly reasonable.

And now that I have "freed up Python", a solution I wrote a while ago will be able to cause the first forward momentum on problems in a long time.

import Data.List
import Data.FixedPoint
import Debug.Trace
-- Replace Double with Rational or with FixedPoint
cfr :: FixedPoint256256 -> (Int -> Int128)

cfr = repindex . cfr_helper

repindex :: ([Int128], Int) -> Int -> Int128

repindex (coeffs, rep) index
    | index < (length coeffs) = coeffs !! index
    | otherwise = coeffs !! (rep + ((index - (rep)) `mod` ((length coeffs) - rep)))

cfr_helper :: FixedPoint256256 -> ([Int128], Int)

cfr_helper alpha = let a0 = floor alpha in
                     let alpha1 = alpha - (fromIntegral a0) in
                       cfr_rec alpha1 [a0, floor (1.0/alpha1)] [alpha1]

cfr_rec :: FixedPoint256256 -> [Int128] -> [FixedPoint256256] -> ([Int128], Int)

gauss_map :: FixedPoint256256 -> FixedPoint256256

gauss_map alpha = (1.0 / alpha) - (fromIntegral (floor(1.0 / alpha)))

fuzzyIndex :: FixedPoint256256 -> [FixedPoint256256] -> Maybe Int

fuzzyIndex _ [] = Nothing
fuzzyIndex x (y:ys)
    | abs(x - y) < 1e-5 = Just 0
    | otherwise = case (fuzzyIndex x ys) of
        Just i  -> Just (i + 1)
        Nothing -> Nothing

cfr_rec alpha0 as alphas = let alpha1 = gauss_map(alpha0) in
    case fuzzyIndex alpha1 alphas of
        Just i  -> (as,i+1)
        Nothing -> cfr_rec alpha1 (as ++ [floor(1.0 / alpha1)]) (alphas ++ [alpha1])

min_soln :: Int -> (Int128,Int)
min_soln d = let coeffs = ((cfr . sqrt' . fromIntegral) d) in
        min_soln_rec d (coeffs) (coeffs 0) 1 1 0 1

min_soln_rec :: Int -> (Int -> Int128) -> Int128 -> Int128 -> Int128 -> Int128 -> Int -> (Int128,Int)

min_soln_rec d coeffs p1 q1 p0 q0 n = let an = coeffs n in
    let (pn,qn) = (p1 * an + p0, q1 * an + q0) in
        if (pn^2 - (fromIntegral d)*qn^2) == 1 then (pn,d)
        else min_soln_rec d coeffs pn qn p1 q1 (n+1)

maxSoln :: (Int128,Int) -> (Int128,Int) -> (Int128,Int)
maxSoln (x1,d1) (x2,d2)
    | x1 > x2 = (x1,d1)
    | otherwise = (x2,d2)

ans :: Int -> Int

isNotSquare :: Int -> Bool

isNotSquare x = let xx = fromIntegral x in (floor (sqrt xx))^2 /= (floor xx)

ans b = snd . foldr maxSoln (0,0) $ map min_soln $ filter isNotSquare [1..b]

main = putStrLn $ show $ ans 1000