module Main where

import Numeric.AERN.DoubleBasis.RealApprox
import Numeric.AERN.DoubleBasis.MRealApprox

import Control.Monad.ST (runST)

type R = RealApprox

main = do
  putStrLn $ "quadratic 1 1 1 = " ++ show (quadratic 1 1 1)
  putStrLn $ "quadratic 1 2 1 = " ++ show (quadratic 1 2 1)
  putStrLn $ "quadratic 1 5 6 = " ++ show (quadratic 1 5 6)
  putStrLn $ "quadraticInPlace 1 1 1 = " ++ show (quadraticInPlace 1 1 1)
  putStrLn $ "quadraticInPlace 1 2 1 = " ++ show (quadraticInPlace 1 2 1)
  putStrLn $ "quadraticInPlace 1 5 6 = " ++ show (quadraticInPlace 1 5 6)

-- | 
-- Returns list such that the union of its members contains all
-- real roots of each real polynomial u*x^2+v*x+w obtained by
-- choosing coefficients u in a, v in b and w in c. 
quadratic :: R -> R -> R -> [R]
quadratic a b c 
  | certainlyZero discriminant =
    [doubleRoot]
  | certainlyNonnegative discriminant =
    [leftRoot,rightRoot]
  | certainlyNegative discriminant =
    []
  | otherwise =
    [bottom]
  where
  discriminant = b^2-4*a*c
  doubleRoot = -b/(2*a)
  leftRoot = doubleRoot-sqrtDiscriminant/(2*a)  -- (-b-sqrt(b^2-4*a*c))/(2*a)
  rightRoot = doubleRoot+sqrtDiscriminant/(2*a) -- (-b+sqrt(b^2-4*a*c))/(2*a)
  sqrtDiscriminant = sqrt discriminant

-- |
-- In-place quadratic. 
quadraticInPlace :: R -> R -> R -> [R]
quadraticInPlace a b c
  | certainlyZero discriminant =
    [doubleRoot]
  | certainlyNonnegative discriminant =
    [leftRoot,rightRoot]
  | certainlyNegative discriminant =
    []
  | otherwise =
    [bottom]
  where
  discriminant = 
    runST $
      do
      [aM,bM,cM] <- mapM makeMutable [a,b,c]
      bM <^>= 2              -- b^2
      aM <*>|= (4 :: Double) -- 4*a
      aM <*>= cM             -- 4*a*c
      bM <->= aM             -- b^2-4*a*c
      result <- unsafeReadMutable bM
        -- using unsafeReadMutable is safe becasue bM is no longer reachable 
      return result
  doubleRoot =
    runST $
      do
      [aM,bM] <- mapM makeMutable [a,b] 
      bM </>|= (-2 :: Double) -- -b/2 
      bM </>= aM              -- -b/(2*a)
      result <- unsafeReadMutable bM
        -- using unsafeReadMutable is safe becasue bM is no longer reachable 
      return result
  [leftRoot,rightRoot] =
    runST $
      do
      [aM,drM,diM] <- mapM makeMutable [a,doubleRoot,discriminant]
      sqrtOutInPlace diM diM  -- sqrt(b^2-4*a*c)
      diM </>|= (2 :: Double) -- (sqrt(b^2-4*a*c))/2 
      diM </>= aM             -- (sqrt(b^2-4*a*c))/(2*a)
      assignMutable aM drM
      aM <->= diM             -- (-b-sqrt(b^2-4*a*c))/(2*a)
      drM <+>= diM            -- (-b+sqrt(b^2-4*a*c))/(2*a)
      result <- mapM unsafeReadMutable [aM,drM]
        -- using unsafeReadMutable is safe becasue aM and drM  are no longer reachable 
      return result

certainlyZero x =
  case x ==? 0 of
    Just True -> True
    _ -> False 

certainlyNonnegative x =
  case x >=? 0 of
    Just True -> True
    _ -> False 

certainlyNegative x =
  case x <? 0 of
    Just True -> True
    _ -> False