import Test.QuickCheck
import Erf
import Dif

-- Normal distribution
norm :: Double -> Double
norm x = (2/sqrt pi) * exp(-x*x)

-- Derivative of the erf function
erfDx :: Double -> Double
erfDx = deriv erf

relDiff :: Double -> Double -> Double
relDiff x y = if x+y == 0 then 1e-18 else abs (2*(x-y) / (x+y))

prop_erf_0 :: Bool
prop_erf_0 = erf 0.0 == (0.0::Double)

-- The differential of erf is the normal distribution function.
prop_deriv :: Double -> Bool
prop_deriv x = erfDx x =~= norm x
  -- At large (absolute) values of the argument the precision drops.
  where a =~= b = a == b ||
	          (     if abs x > 25 then True
                   else if abs x > 11 then relDiff a b < 1e-13
                   else if abs x >  3 then relDiff a b < 1e-14
                   else                    relDiff a b < 1e-15)

-- At small values of erf we get cancellation,
-- so ignore those.
prop_erf_erfc :: Double -> Bool
prop_erf_erfc x = erf x =~= (1 - erfc x)
  where a =~= b = abs b < 1e-5 || abs(a-b) < 1e-10

-- At large (absolute) arguments we can't scale the
-- result back with precision, so ignore those.
prop_erfc_erfcx :: Double -> Bool
prop_erfc_erfcx x = (erfcx x / exp(x*x)) =~= erfc x
  where a =~= b = abs x > 6 || abs(a-b) < 1e-14

myTest :: (Testable a) => a -> IO ()
myTest = check config
  where config = defaultConfig {
	    configMaxTest = 10000,
	    configEvery = \ _ _ -> ""
	    }


main :: IO ()
main = do

        putStrLn "Testing erf 0"
	myTest prop_erf_0

        putStrLn "Testing differential"
	myTest prop_deriv

	putStrLn "Testing erf - erfc"
	myTest prop_erf_erfc

	putStrLn "Testing erfc - erfcx"
	myTest prop_erfc_erfcx