Error finding the number of positive integers, x, such that x(t-x)>d in haskell

Question

With the input I gave to numWays in main I am getting an incorrect answer, but I can't see any fault in my algorithm,

isInteger :: RealFrac a => a -> Bool
isInteger x = x == fromIntegral (round x)

big :: Float -> Int
big x | x <= 0 = error "big x"
    | otherwise = if isInteger x then round (x+1) else ceiling x

small :: Float -> Int
small x | x <= 0 = error "small x"
    | otherwise = if isInteger x then round (x-1) else floor x

numWays :: Float -> Float -> Int
numWays t d = small ((t + sqrt (t^2 - 4*d))/2) - big ((t- sqrt (t^2-4*d))/2) + 1

main :: IO ()
main = do
    print $ numWays 58819676.0 434104122191218.0

the output from this is:

ghci> :r  
[1 of 2] Compiling Main             ( d6.hs, interpreted ) [Source file changed]
Ok, one module loaded.
ghci> main
41513106

The correct answer is 41513103.

I've tried asking chatGPT and a friend. ChatGPT says this is because of rounding errors, but I don't know how to get around that.

This is day 6, part 2 of advent of code. I have ran my code on the example and I get the right answer (7503)

Answer 1

Caveat lector: I have not read the description of Advent of Code problem 6. There may be a cleaner, more direct solution. I am basing this answer only off of the solution posted in this StackOverflow question, showing how to get an exact calculation without modifying the core algorithmic idea.

One way to avoid rounding is to never round -- or at least, to carefully control exactly where rounding happens. I observe that we sometimes need to know whether a number has a fractional part, but never need any of the details of the exact value of that fractional part. Therefore I propose the following fixed-point type:

data Fraction = Fraction
    { wholePart :: Integer
    , fractionalPart :: Bool
    } deriving (Eq, Ord, Read, Show)

A quick contrast between floating-point types like Float and Double: for small numbers, the floating-point types keep track of much more information about the fractional part, but for large numbers, floating-point types give up on tracking the fractional part at all and in fact may not even store the correct whole part. So, if we want to carefully control our rounding to make sure we know the correct value of both parts at all times, we will need to carefully avoid Double and Float.

The trickiest bit of avoiding them is finding an analog of sqrt. There are various techniques for implementing it, such as Newton's method, that I'll let you read up on and try yourself if you're interested. But for the purposes of this answer, I'll simply use the integer-roots package, which exports an integerSquareRoot function. We'll write a small adapter that uses this to inhabit our Fraction type:

import Math.NumberTheory.Roots

isqrt :: Integer -> Fraction
isqrt n = Fraction
    { wholePart = w
    , fractionalPart = w*w<n
    } where
    w = integerSquareRoot n

We'll also need some basic arithmetic: addition, subtraction, and division by integers. I'll put a . on the side that takes a Fraction.

(+.) :: Integer -> Fraction -> Fraction
n +. f = f { wholePart = wholePart f + n }

(-.) :: Integer -> Fraction -> Fraction
n -. f = f
    { wholePart = n - wholePart f - if fractionalPart f then 1 else 0
    }

(./) :: Fraction -> Integer -> Fraction
f ./ n = Fraction
    { wholePart = wholePart f `div` n
    , fractionalPart = fractionalPart f || (wholePart f`mod`n > 0)
    }

Finally, we can implement your big and small. I find your implementations a bit odd; are you really sure you got them right? See the clear asymmetry in the analogous Fraction versions to see what I mean:

big :: Fraction -> Integer
big = succ . wholePart

small :: Fraction -> Integer
small f = wholePart f - if fractionalPart f then 0 else 1

In any case, with this machinery in place, we can now transliterate your numWays:

numWays :: Integer -> Integer -> Integer
numWays t d = small ((t +. isqrt (t^2 - 4*d)) ./ 2) - big ((t -. isqrt (t^2 - 4*d)) ./ 2) + 1

Because there are no floating point types at any point, we can now handle even larger inputs than the Double-y version proposed in the other answer. For example, smashing the keyboard for a while to get a big number, it's not hard to cause rounding to happen once again even with the larger precision of Double:

>               235956230947^2 - 4*13918835730678500119202 :: Integer
40001
> numWays       235956230947       13918835730678500119202
200
> numWaysDouble 235956230947       13918835730678500119202
2896

(Note that 40001 = 200^2+1.)