Search code examples
mathhaskellpattern-matchingcomplex-numbers

Haskell floating point error


So I have finished creating my own complex number data type in haskell.

I've also, thanks to another question on here, got a function that will solve a quadratic equation.

The only problem now is that the code generates a parsing error in hugs, when trying to solve a quadratic with complex roots.

i.e. In hugs...

Main> solve (Q 1 2 1)
(-1.0,-1.0)

Main> solve (Q 1 2 0)
(0.0,-2.0)

Main> solve (Q 1 2 2)
(
Program error: pattern match failure: v1618_v1655 (C -1.#IND -1.#IND)

It looks to my like its a problem after the square-root has been applied, but I'm really not sure. Any help trying to pick up what is going wrong or any indications as to what this error means would be brilliant.

Thanks,

Thomas

The Code:

-- A complex number z = (re +im.i) is represented as a pair of Floats

data Complex = C {
re :: Float,
im :: Float
} deriving Eq

-- Display complex numbers in the normal way

instance Show Complex where
    show (C r i)
        | i == 0            = show r
        | r == 0            = show i++"i"
        | r < 0 && i < 0    = show r ++ " - "++ show (C 0 (i*(-1)))
        | r < 0 && i > 0    = show r ++ " + "++ show (C 0 i)
        | r > 0 && i < 0    = show r ++ " - "++ show (C 0 (i*(-1)))
        | r > 0 && i > 0    = show r ++ " + "++ show (C 0 i)


-- Define algebraic operations on complex numbers
instance Num Complex where
    fromInteger n       = C (fromInteger n) 0 -- tech reasons
    (C a b) + (C x y)   = C (a+x) (b+y)
    (C a b) * (C x y)   = C (a*x - b*y) (b*x + b*y)
    negate (C a b)      = C (-a) (-b)

instance Fractional Complex where
    fromRational r      = C (fromRational r) 0 -- tech reasons
    recip (C a b)       = C (a/((a^2)+(b^2))) (b/((a^2)+(b^2)))


root :: Complex -> Complex
root (C x y)
    | y == 0 && x == 0  = C 0 0
    | y == 0 && x > 0   = C (sqrt ( ( x + sqrt ( (x^2) + 0 ) ) / 2 ) )  0
    | otherwise         = C (sqrt ( ( x + sqrt ( (x^2) + (y^2) ) ) / 2 ) ) ((y/(2*(sqrt ( ( x + sqrt ( (x^2) + (y^2) ) ) / 2 ) ) ) ) )


-- quadratic polynomial : a.x^2 + b.x + c
data Quad = Q {
    aCoeff, bCoeff, cCoeff :: Complex
} deriving Eq


instance Show Quad where
    show (Q a b c) = show a ++ "x^2 + " ++ show b ++ "x + " ++ show c

solve :: Quad -> (Complex, Complex)
solve (Q a b c) = ( sol (+), sol (-) )
  where sol op = (op (negate b) $ root $ b*b - 4*a*c) / (2 * a)

Solution

  • Your numbers seem denormalized in your error :

    (C -1.#IND -1.#IND)
    

    In this case, you can't assume that any comparison on float are valid anymore. This is in the definition of floating point numbers. Then your definition of show

    show (C r i)
        | i == 0                        = show r
        | r == 0                        = show i++"i"
        | r < 0 && i < 0        = show r ++ " - "++ show (C 0 (i*(-1)))
        | r < 0 && i > 0        = show r ++ " + "++ show (C 0 i)
        | r > 0 && i < 0        = show r ++ " - "++ show (C 0 (i*(-1)))
        | r > 0 && i > 0        = show r ++ " + "++ show (C 0 i)
    

    leave opportunity for a pattern failure, because of denormalized numbers. You can add the following condition

        | otherwise = show r ++ "i" ++ show i"
    

    Now for the why is it like that, when you evaluate

    b * b - 4 * a * c
    

    with Q 1 2 2, you obtain -4, and then in root, you fall in your last case, and in the second equation :

                  y
    -----------------------------
                ________________
               /        _______
              /        / 2    2
             /   x + \/ x  + y
     2 * \  /   ----------------
          \/           2
    

    -4 + sqrt( (-4) ^2) == 0, from there, you're doomed, division by 0, followed by a "NaN" (not a number), screwing everything else