{-# 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
addLoc String
errStr
    where
      addLoc :: String -> ShowS
addLoc String
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