algebraicReal :: Equation -> Flow [Double]
  algebraicReal (Equation lhs rhs) = solve $ simplify $ reduce lhs rhs
    where
      solve :: Expr -> Flow [Double]
      solve (Minus (Constant c) (Variable x)) = return [c]
      solve (Plus (Constant c) (Variable x)) = return [(-c)]
      solve (Plus (Constant c) (Mult (Constant a) (Variable x))) = return [c / a]
      solve (Plus (Constant c)
                  (Plus (Mult (Constant a) (Pow (Variable x1) (Constant 2)))
                        (Mult (Constant b) (Variable x2)))) = quadratic a b c
      solve (Plus (Constant c)
                  (Plus (Pow (Variable x1) (Constant 2))
                        (Mult (Constant b) (Variable x2)))) = quadratic 1 b c
      solve (Plus (Constant c)
                  (Plus (Pow (Variable x1) (Constant 2))
                        (Variable x2))) = quadratic 1 1 c
      solve (Plus (Pow (Variable x1) (Constant 2))
                  (Variable x2)) = quadratic 1 0 1
      solve x = do
        tell ["Don't know how solve " ++ (show x)]
        throwError $ "Fail at " ++ (show x) 

      quadratic a b c =
        let d = b * b - 4 * a * c
        in  if (d >= 0) then 
              return [(-b / 2 * a) + sqrt d, (-b / 2 * a) - sqrt d] 
            else 
              throwError "Complex solution"
algebraicReal :: Equation -> Flow [Double]
algebraicReal (Equation lhs rhs) = solve $ simplify $ reduce lhs rhs
    where
      solve :: Expr -> Flow [Double]
      solve (Minus (Constant c) (Variable x)) = return [c]
      solve (Plus (Constant c) (Variable x)) = return [(-c)]
      solve (Plus (Constant c) (Mult (Constant a) (Variable x))) = return [c / a]
      solve (Plus (Constant c)
                  (Plus (Mult (Constant a) (Pow (Variable x1) (Constant 2)))
                        (Mult (Constant b) (Variable x2)))) = quadratic a b c
      solve (Plus (Constant c)
                  (Plus (Pow (Variable x1) (Constant 2))
                        (Mult (Constant b) (Variable x2)))) = quadratic 1 b c
      solve (Plus (Constant c)
                  (Plus (Pow (Variable x1) (Constant 2))
                        (Variable x2))) = quadratic 1 1 c
      solve (Plus (Pow (Variable x1) (Constant 2))
                  (Variable x2)) = quadratic 1 0 1
      solve x = do
        tell ["Don't know how solve " ++ (show x)]
        throwError $ "Fail at " ++ (show x) 

      quadratic a b c =
        let d = b * b - 4 * a * c
        in  if (d >= 0) then 
              return [(-b / 2 * a) + sqrt d, (-b / 2 * a) - sqrt d] 
            else 
              throwError "Complex solution"
algebraicReal (Equation lhs rhs) = solve $ simplify $ reduce lhs rhs
    where
      solve :: Expr -> Flow [Double]
      solve (Minus (Constant c) (Variable x)) = return [c]
      solve (Plus (Constant c) (Variable x)) = return [(-c)]
      solve (Plus (Constant c) (Mult (Constant a) (Variable x))) = return [c / a]
      solve (Plus (Constant c)
                  (Plus (Mult (Constant a) (Pow (Variable x1) (Constant 2)))
                        (Mult (Constant b) (Variable x2)))) = quadratic a b c
      solve (Plus (Constant c)
                  (Plus (Pow (Variable x1) (Constant 2))
                        (Mult (Constant b) (Variable x2)))) = quadratic 1 b c
      solve (Plus (Constant c)
                  (Plus (Pow (Variable x1) (Constant 2))
                        (Variable x2))) = quadratic 1 1 c
      solve (Plus (Pow (Variable x1) (Constant 2))
                  (Variable x2)) = quadratic 1 0 1
      solve x = do
        tell ["Don't know how solve " ++ (show x)]
        throwError $ "Fail at " ++ (show x) 

      quadratic a b c =
        let d = b * b - 4 * a * c
        in  if (d >= 0) then 
              return [(-b / 2 * a) + sqrt d, (-b / 2 * a) - sqrt d] 
            else 
              throwError "Complex solution"