{-# LANGUAGE ConstraintKinds  #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE PatternSynonyms  #-}
{-# LANGUAGE ViewPatterns     #-}
module Numeric.Basics
( -- * Constants
pattern M_E, pattern M_LOG2E, pattern M_LOG10E, pattern M_LN2, pattern M_LN10
, pattern M_SQRT2, pattern M_SQRT1_2
, pattern M_PI, pattern M_PI_2, pattern M_PI_3, pattern M_PI_4
, pattern M_1_PI, pattern M_2_PI, pattern M_2_SQRTPI
, Epsilon (..), pattern M_EPS
-- * Functions
, RealExtras (..), RealFloatExtras (..)
, negateUnless
-- * Exceptions
, IterativeMethod, TooManyIterations, tooManyIterations
) where

import Data.Int
import Data.Word
import GHC.Exception
import GHC.Stack
import Numeric.PrimBytes

-- | $$e$$
pattern M_E :: (Eq a, Fractional a) => Fractional a => a
pattern $bM_E :: a$mM_E :: forall r a.
(Eq a, Fractional a) =>
a -> (Fractional a => r) -> (Void# -> r) -> r
M_E = 2.718281828459045235360287471352662498

-- | $$\log_2 e$$
pattern M_LOG2E :: (Eq a, Fractional a) => Fractional a => a
pattern $bM_LOG2E :: a$mM_LOG2E :: forall r a.
(Eq a, Fractional a) =>
a -> (Fractional a => r) -> (Void# -> r) -> r
M_LOG2E = 1.442695040888963407359924681001892137

-- | $$\log_{10} e$$
pattern M_LOG10E :: (Eq a, Fractional a) => Fractional a => a
pattern $bM_LOG10E :: a$mM_LOG10E :: forall r a.
(Eq a, Fractional a) =>
a -> (Fractional a => r) -> (Void# -> r) -> r
M_LOG10E = 0.434294481903251827651128918916605082

-- | $$\log_e 2$$
pattern M_LN2 :: (Eq a, Fractional a) => Fractional a => a
pattern $bM_LN2 :: a$mM_LN2 :: forall r a.
(Eq a, Fractional a) =>
a -> (Fractional a => r) -> (Void# -> r) -> r
M_LN2 = 0.693147180559945309417232121458176568

-- | $$\log_e 10$$
pattern M_LN10 :: (Eq a, Fractional a) => Fractional a => a
pattern $bM_LN10 :: a$mM_LN10 :: forall r a.
(Eq a, Fractional a) =>
a -> (Fractional a => r) -> (Void# -> r) -> r
M_LN10 = 2.302585092994045684017991454684364208

-- | $$\sqrt{2}$$
pattern M_SQRT2 :: (Eq a, Fractional a) => Fractional a => a
pattern $bM_SQRT2 :: a$mM_SQRT2 :: forall r a.
(Eq a, Fractional a) =>
a -> (Fractional a => r) -> (Void# -> r) -> r
M_SQRT2 = 1.414213562373095048801688724209698079

-- | $$\frac{1}{\sqrt{2}}$$
pattern M_SQRT1_2 :: (Eq a, Fractional a) => Fractional a => a
pattern $bM_SQRT1_2 :: a$mM_SQRT1_2 :: forall r a.
(Eq a, Fractional a) =>
a -> (Fractional a => r) -> (Void# -> r) -> r
M_SQRT1_2 = 0.707106781186547524400844362104849039

-- | $$\pi$$
pattern M_PI :: (Eq a, Fractional a) => Fractional a => a
pattern $bM_PI :: a$mM_PI :: forall r a.
(Eq a, Fractional a) =>
a -> (Fractional a => r) -> (Void# -> r) -> r
M_PI = 3.1415926535897932384626433832795028841971
-- | $$\frac{\pi}{2}$$
pattern M_PI_2 :: (Eq a, Fractional a) => Fractional a => a
pattern $bM_PI_2 :: a$mM_PI_2 :: forall r a.
(Eq a, Fractional a) =>
a -> (Fractional a => r) -> (Void# -> r) -> r
M_PI_2 = 1.570796326794896619231321691639751442

-- | $$\frac{\pi}{3}$$
pattern M_PI_3 :: (Eq a, Fractional a) => Fractional a => a
pattern $bM_PI_3 :: a$mM_PI_3 :: forall r a.
(Eq a, Fractional a) =>
a -> (Fractional a => r) -> (Void# -> r) -> r
M_PI_3 = 1.047197551196597746154214461093167628

-- | $$\frac{\pi}{4}$$
pattern M_PI_4 :: (Eq a, Fractional a) => Fractional a => a
pattern $bM_PI_4 :: a$mM_PI_4 :: forall r a.
(Eq a, Fractional a) =>
a -> (Fractional a => r) -> (Void# -> r) -> r
M_PI_4 = 0.785398163397448309615660845819875721

-- | $$\frac{1}{\pi}$$
pattern M_1_PI :: (Eq a, Fractional a) => Fractional a => a
pattern $bM_1_PI :: a$mM_1_PI :: forall r a.
(Eq a, Fractional a) =>
a -> (Fractional a => r) -> (Void# -> r) -> r
M_1_PI = 0.318309886183790671537767526745028724

-- | $$\frac{2}{\pi}$$
pattern M_2_PI :: (Eq a, Fractional a) => Fractional a => a
pattern $bM_2_PI :: a$mM_2_PI :: forall r a.
(Eq a, Fractional a) =>
a -> (Fractional a => r) -> (Void# -> r) -> r
M_2_PI = 0.636619772367581343075535053490057448

-- | $$\frac{2}{\sqrt{\pi}}$$
pattern M_2_SQRTPI :: (Eq a, Fractional a) => Fractional a => a
pattern $bM_2_SQRTPI :: a$mM_2_SQRTPI :: forall r a.
(Eq a, Fractional a) =>
a -> (Fractional a => r) -> (Void# -> r) -> r
M_2_SQRTPI = 1.128379167095512573896158903121545172

-- | Define a meaningful precision for complex floating-point operations.
class (Eq a, Floating a) => Epsilon a where
-- | A small positive number that depends on the type of @a@.
--   Compare your values to @epsilon@ to tell if they are near zero.
epsilon :: a

instance Epsilon Double where
epsilon :: Double
epsilon = Double
1e-12

instance Epsilon Float where
epsilon :: Float
epsilon = Float
1e-6

-- | A small positive number that depends on the type of @a@.
pattern M_EPS :: Epsilon a => a
pattern $bM_EPS :: a$mM_EPS :: forall r a. Epsilon a => a -> (Void# -> r) -> (Void# -> r) -> r
M_EPS <- ((epsilon ==) -> True)
where
M_EPS = a
forall a. Epsilon a => a
epsilon

-- | Negate if @False@.
--   This is a useful alternative to signum,
--     when @signum 0 == 0@ causing some troubles.
negateUnless :: Num t => Bool -> t -> t
negateUnless :: Bool -> t -> t
negateUnless Bool
True  = t -> t
forall a. a -> a
id
negateUnless Bool
False = t -> t
forall a. Num a => a -> a
negate
{-# INLINE negateUnless #-}

-- | Extra functions for Real types.
class (Real a, PrimBytes a) => RealExtras a where
-- | @copysign x y@ returns a value with the magnitude of x and the sign of y.
--
--   NB: in future, this function is to reimplemented using primops
--       and should behave the same way as its C analogue.
copysign :: a -> a -> a
copysign a
x a
y
| (a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
0) Bool -> Bool -> Bool
forall a. Eq a => a -> a -> Bool
== (a
y a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
0) = a
x
| Bool
otherwise = a -> a
forall a. Num a => a -> a
negate a
x
{-# INLINE copysign #-}

instance RealExtras Int
instance RealExtras Int8
instance RealExtras Int16
instance RealExtras Int32
instance RealExtras Int64
instance RealExtras Word
instance RealExtras Word8
instance RealExtras Word16
instance RealExtras Word32
instance RealExtras Word64

instance RealExtras Float where
copysign :: Float -> Float -> Float
copysign = Float -> Float -> Float
c'copysignf
{-# INLINE copysign #-}

instance RealExtras Double where
copysign :: Double -> Double -> Double
copysign = Double -> Double -> Double
c'copysignd
{-# INLINE copysign #-}

-- | Extra functions for RealFrac types.
class (Epsilon a, RealExtras a, RealFloat a) => RealFloatExtras a where
-- | $$\sqrt{ x^2 + y^2 }$$.
--
--   NB: in future, this function is to reimplemented using primops
--       and should behave the same way as its C analogue.
hypot :: a -> a -> a
hypot a
x a
y = Int -> a -> a
forall a. RealFloat a => Int -> a -> a
scaleFloat Int
ea (a -> a
forall a. Floating a => a -> a
sqrt (a
ana -> a -> a
forall a. Num a => a -> a -> a
*a
an a -> a -> a
forall a. Num a => a -> a -> a
+ a
bna -> a -> a
forall a. Num a => a -> a -> a
*a
bn))
where
x' :: a
x' = a -> a
forall a. Num a => a -> a
abs a
x
y' :: a
y' = a -> a
forall a. Num a => a -> a
abs a
y
(a
a,a
b) = if a
x' a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
y' then (a
x', a
y') else (a
y', a
x')
(Integer
_, Int
ea) = a -> (Integer, Int)
forall a. RealFloat a => a -> (Integer, Int)
decodeFloat a
a
an :: a
an = Int -> a -> a
forall a. RealFloat a => Int -> a -> a
scaleFloat (Int -> Int
forall a. Num a => a -> a
negate Int
ea) a
a
bn :: a
bn = Int -> a -> a
forall a. RealFloat a => Int -> a -> a
scaleFloat (Int -> Int
forall a. Num a => a -> a
negate Int
ea) a
b
{-# INLINE hypot #-}
-- | Maximum finite number representable by this FP type.
maxFinite :: a

instance RealFloatExtras Float where
hypot :: Float -> Float -> Float
hypot = Float -> Float -> Float
c'hypotf
{-# INLINE hypot #-}
maxFinite :: Float
maxFinite = Float
3.40282347e+38
{-# INLINE maxFinite #-}

instance RealFloatExtras Double where
hypot :: Double -> Double -> Double
hypot = Double -> Double -> Double
c'hypotd
{-# INLINE hypot #-}
maxFinite :: Double
maxFinite = Double
1.7976931348623157e+308
{-# INLINE maxFinite #-}

foreign import ccall unsafe "hypot"  c'hypotd :: Double -> Double -> Double
foreign import ccall unsafe "hypotf" c'hypotf :: Float -> Float -> Float
foreign import ccall unsafe "copysign"  c'copysignd :: Double -> Double -> Double
foreign import ccall unsafe "copysignf" c'copysignf :: Float -> Float -> Float

{- |
This is just an alias for the HasCallStack constraint.
It servers two goals:

1. Document functions that use iterative methods and can fail
(supposedly in extremely rare cases).
2. Provide a call stack in case of failure.

Use tooManyIterations function where necessary if you implement an iterative
algorithm and annotate your implementation function with IterativeMethod constraint.

-}
type IterativeMethod = HasCallStack

-- | Throw a TooManyIterations exception.
tooManyIterations ::
IterativeMethod
=> String  -- ^ Label (e.g. function name)
-> a
tooManyIterations :: String -> a
tooManyIterations String
s
= TooManyIterations -> a
forall a e. Exception e => e -> a
throw TooManyIterations :: String -> CallStack -> TooManyIterations
TooManyIterations
{ tmiUserInfo :: String
tmiUserInfo  = String
s
, tmiCallStack :: CallStack
tmiCallStack = CallStack
HasCallStack => CallStack
callStack
}

{- |
Typically, this exception can occur in an iterative algorithm;
when the number of iterations exceeds a pre-defined limit, but the stop condition
is not fulfilled.

Use the IterativeMethod constraint to annotate functions that throw this exception.
-}
data TooManyIterations
= TooManyIterations
{ TooManyIterations -> String
tmiUserInfo  :: String
-- ^ Short description of the error, e.g. a function name.
, TooManyIterations -> CallStack
tmiCallStack :: CallStack
-- ^ Function call stack.
--   Note, this field is ignored in the Eq and Ord instances.
}

-- | Note, this instance ignores oodCallStack
instance Eq TooManyIterations where
== :: TooManyIterations -> TooManyIterations -> Bool
(==) TooManyIterations
a TooManyIterations
b = TooManyIterations -> String
tmiUserInfo TooManyIterations
a String -> String -> Bool
forall a. Eq a => a -> a -> Bool
== TooManyIterations -> String
tmiUserInfo TooManyIterations
b

-- | Note, this instance ignores oodCallStack
instance Ord TooManyIterations where
compare :: TooManyIterations -> TooManyIterations -> Ordering
compare TooManyIterations
a TooManyIterations
b = String -> String -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (TooManyIterations -> String
tmiUserInfo TooManyIterations
a) (TooManyIterations -> String
tmiUserInfo TooManyIterations
b)

instance Show TooManyIterations where
showsPrec :: Int -> TooManyIterations -> ShowS
showsPrec Int
p TooManyIterations
e = String -> ShowS
errStr
where
s
= let someE :: SomeException
someE = String -> CallStack -> SomeException
errorCallWithCallStackException String
s (TooManyIterations -> CallStack
tmiCallStack TooManyIterations
e)
errc :: ErrorCall
errc :: ErrorCall
errc = case SomeException -> Maybe ErrorCall
forall e. Exception e => SomeException -> Maybe e
fromException SomeException
someE of
Maybe ErrorCall
Nothing -> String -> ErrorCall
ErrorCall String
s
Just ErrorCall
ec -> ErrorCall
ec
in  Int -> ErrorCall -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
p ErrorCall
errc
errStr :: String
errStr = [String] -> String
unlines
[ String
"An iterative method has exceeded the limit for the iteration number."
, String
"User info: " String -> ShowS
forall a. [a] -> [a] -> [a]
++ TooManyIterations -> String
tmiUserInfo TooManyIterations
e
]

instance Exception TooManyIterations