-- SPDX-License-Identifier: MPL-2.0

{- | Scalar utilities. Generic Num arithmetic follows the chosen numeric type.
Checked real operations require finite inputs and report domain failures and
finite-result overflow explicitly. Exact combinatorics use 'Integer'.
expm1/log1p use convergent series near zero to avoid cancellation; their target
is a few ulps for ordinary IEEE binary32/64 inputs, not correct rounding.
-}
module HMath.Scalar (
    add,
    subtractScalar,
    multiply,
    divide,
    quotientRemainder,
    euclideanDivision,
    absolute,
    sign,
    minimumScalar,
    maximumScalar,
    clampChecked,
    lerpChecked,
    inverseLerp,
    remap,
    step,
    smoothstep,
    smootherstep,
    power,
    integerPower,
    modularPower,
    sqrtChecked,
    cubeRoot,
    nthRoot,
    reciprocalSqrt,
    expChecked,
    logChecked,
    expm1,
    log1p,
    hypot,
    floorChecked,
    ceilChecked,
    roundChecked,
    truncateChecked,
    bankersRound,
    roundDigits,
    integerFraction,
    positiveModulo,
    absoluteError,
    relativeError,
    factorialChecked,
    doubleFactorial,
    risingFactorial,
    fallingFactorial,
    fibonacciChecked,
    lucas,
    permutations,
    combinations,
    multinomial,
    catalan,
    bernoulli,
    harmonic,
    cumulativeSum,
    cumulativeProduct,
    prefixSum,
    fusedMultiplyAdd,
    saturatingAdd,
    saturatingSubtract,
    saturatingMultiply,
    checkedAdd,
    checkedSubtract,
    checkedMultiply,
    checkedFloatMultiply,
) where

import qualified Data.List as L
import Data.Ratio ((%))
import HMath.Error (MathError (..))

add, subtractScalar, multiply :: (Num a) => a -> a -> a

{- | Add two values using their Num instance. No overflow check is added; bounded integral
types retain their usual arithmetic behavior. Example: add 2 3 == 5.
-}
add :: forall a. Num a => a -> a -> a
add = a -> a -> a
forall a. Num a => a -> a -> a
(+)

{- | Subtract the second value from the first using Num arithmetic. Example: subtractScalar
5 2 == 3; use checkedSubtract to reject bounded integer overflow.
-}
subtractScalar :: forall a. Num a => a -> a -> a
subtractScalar = (-)

{- | Multiply two values using Num arithmetic. Example: multiply 3 4 == 12; use
checkedMultiply or checkedFloatMultiply for explicit overflow checks.
-}
multiply :: forall a. Num a => a -> a -> a
multiply = a -> a -> a
forall a. Num a => a -> a -> a
(*)

finite :: (RealFloat a) => a -> Bool
finite :: forall a. RealFloat a => a -> Bool
finite a
x = Bool -> Bool
not (a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
x Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
x)

result :: (RealFloat a) => a -> Either MathError a
result :: forall a. RealFloat a => a -> Either MathError a
result a
x
    | a -> Bool
forall a. RealFloat a => a -> Bool
finite a
x = a -> Either MathError a
forall a b. b -> Either a b
Right a
x
    | Bool
otherwise = MathError -> Either MathError a
forall a b. a -> Either a b
Left MathError
NonFiniteResult

inputs :: (RealFloat a) => [a] -> Either MathError ()
inputs :: forall a. RealFloat a => [a] -> Either MathError ()
inputs [a]
xs
    | (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all a -> Bool
forall a. RealFloat a => a -> Bool
finite [a]
xs = () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    | Bool
otherwise = MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"non-finite scalar input")

{- | Divide the first finite real value by the second. Return Left DomainError for
non-finite input or a zero divisor, and Left NonFiniteResult for overflow. Example:
divide 6 2 == Right 3.
-}
divide :: (RealFloat a) => a -> a -> Either MathError a
divide :: forall a. RealFloat a => a -> a -> Either MathError a
divide a
a a
b = do
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
a, a
b]
    if a
b a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 then MathError -> Either MathError a
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"division by zero") else a -> Either MathError a
forall a. RealFloat a => a -> Either MathError a
result (a
a a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
b)

quotientRemainder, euclideanDivision :: Integer -> Integer -> Either MathError (Integer, Integer)

{- | Return integer quotient and remainder using quotRem (quotient truncates toward zero). A
zero divisor returns Left DomainError. Example: quotientRemainder (-7) 3 == Right
(-2,-1).
-}
quotientRemainder :: Integer -> Integer -> Either MathError (Integer, Integer)
quotientRemainder Integer
_ Integer
0 = MathError -> Either MathError (Integer, Integer)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"division by zero")
quotientRemainder Integer
a Integer
b = (Integer, Integer) -> Either MathError (Integer, Integer)
forall a b. b -> Either a b
Right (Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
quotRem Integer
a Integer
b)

-- | Euclidean remainder satisfies 0 <= r < abs divisor, for either divisor sign.
euclideanDivision :: Integer -> Integer -> Either MathError (Integer, Integer)
euclideanDivision Integer
_ Integer
0 = MathError -> Either MathError (Integer, Integer)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"division by zero")
euclideanDivision Integer
a Integer
b = let r :: Integer
r = Integer
a Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer -> Integer
forall a. Num a => a -> a
abs Integer
b in (Integer, Integer) -> Either MathError (Integer, Integer)
forall a b. b -> Either a b
Right ((Integer
a Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
r) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`quot` Integer
b, Integer
r)

absolute, sign :: (Num a) => a -> a

{- | Return abs using the supplied Num instance, including that instance's overflow
behavior. Example: absolute (-3) == 3.
-}
absolute :: forall a. Num a => a -> a
absolute = a -> a
forall a. Num a => a -> a
abs

-- | Return signum using the supplied Num instance. Example: sign (-3) == -1; sign 0 == 0.
sign :: forall a. Num a => a -> a
sign = a -> a
forall a. Num a => a -> a
signum

minimumScalar, maximumScalar :: (Ord a) => a -> a -> a

{- | Return the smaller of two values under their Ord instance. Example: minimumScalar 2 5
== 2.
-}
minimumScalar :: forall a. Ord a => a -> a -> a
minimumScalar = a -> a -> a
forall a. Ord a => a -> a -> a
min

{- | Return the larger of two values under their Ord instance. Example: maximumScalar 2 5 ==
5.
-}
maximumScalar :: forall a. Ord a => a -> a -> a
maximumScalar = a -> a -> a
forall a. Ord a => a -> a -> a
max

{- | Given lower bound, upper bound, then value, clamp to the closed interval. Non-finite
inputs or reversed bounds return Left DomainError. Example: clampChecked 0 1 2 == Right
1.
-}
clampChecked :: (RealFloat a) => a -> a -> a -> Either MathError a
clampChecked :: forall a. RealFloat a => a -> a -> a -> Either MathError a
clampChecked a
lo a
hi a
x = do
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
lo, a
hi, a
x]
    if a
lo a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
hi then MathError -> Either MathError a
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"reversed clamp bounds") else a -> Either MathError a
forall a b. b -> Either a b
Right (a -> a -> a
forall a. Ord a => a -> a -> a
max a
lo (a -> a -> a
forall a. Ord a => a -> a -> a
min a
hi a
x))

{- | Uses exact rational interpolation followed by one rounding. Extrapolation
is allowed; finite overflow is reported instead of silently producing infinity.
-}
lerpChecked :: (RealFloat a) => a -> a -> a -> Either MathError a
lerpChecked :: forall a. RealFloat a => a -> a -> a -> Either MathError a
lerpChecked a
a a
b a
t = do
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
a, a
b, a
t]
    a -> Either MathError a
forall a. RealFloat a => a -> Either MathError a
result (Rational -> a
forall a. Fractional a => Rational -> a
fromRational (a -> Rational
forall a. Real a => a -> Rational
toRational a
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ a -> Rational
forall a. Real a => a -> Rational
toRational a
t Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* (a -> Rational
forall a. Real a => a -> Rational
toRational a
b Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- a -> Rational
forall a. Real a => a -> Rational
toRational a
a)))

{- | Given endpoints a,b then value x, return (x-a)/(b-a) using rational intermediates and
one conversion. Equal endpoints/non-finite inputs fail; extrapolation is allowed.
Example: inverseLerp 10 20 15 == Right 0.5.
-}
inverseLerp :: (RealFloat a) => a -> a -> a -> Either MathError a
inverseLerp :: forall a. RealFloat a => a -> a -> a -> Either MathError a
inverseLerp a
a a
b a
x = do
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
a, a
b, a
x]
    if a
a a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
b
        then MathError -> Either MathError a
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"equal interpolation endpoints")
        else a -> Either MathError a
forall a. RealFloat a => a -> Either MathError a
result (Rational -> a
forall a. Fractional a => Rational -> a
fromRational ((a -> Rational
forall a. Real a => a -> Rational
toRational a
x Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- a -> Rational
forall a. Real a => a -> Rational
toRational a
a) Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ (a -> Rational
forall a. Real a => a -> Rational
toRational a
b Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- a -> Rational
forall a. Real a => a -> Rational
toRational a
a)))

{- | Given source endpoints a,b, destination endpoints c,d, then x, map x from the source
interval to the destination interval. Extrapolation is allowed. Fails on equal source
endpoints, non-finite input/result. Example: remap 0 10 0 100 5 == Right 50.
-}
remap :: (RealFloat a) => a -> a -> a -> a -> a -> Either MathError a
remap :: forall a.
RealFloat a =>
a -> a -> a -> a -> a -> Either MathError a
remap a
a a
b a
c a
d a
x = a -> a -> a -> Either MathError a
forall a. RealFloat a => a -> a -> a -> Either MathError a
inverseLerp a
a a
b a
x Either MathError a
-> (a -> Either MathError a) -> Either MathError a
forall a b.
Either MathError a
-> (a -> Either MathError b) -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= a -> a -> a -> Either MathError a
forall a. RealFloat a => a -> a -> a -> Either MathError a
lerpChecked a
c a
d

{- | Given edge then x, return 0 below the edge and 1 at/above it. Non-finite inputs return
Left DomainError. Example: step 2 2 == Right 1.
-}
step :: (RealFloat a) => a -> a -> Either MathError a
step :: forall a. RealFloat a => a -> a -> Either MathError a
step a
edge a
x = [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
edge, a
x] Either MathError () -> Either MathError a -> Either MathError a
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> a -> Either MathError a
forall a b. b -> Either a b
Right (if a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
edge then a
0 else a
1)

smoothstep, smootherstep :: (RealFloat a) => a -> a -> a -> Either MathError a

{- | Given strictly increasing edges then x, clamp t=(x-a)/(b-a) to [0,1] and return
t*t*(3-2*t). Invalid edges/non-finite values fail. Example: smoothstep 0 1 0.5 == Right
0.5.
-}
smoothstep :: forall a. RealFloat a => a -> a -> a -> Either MathError a
smoothstep a
a a
b a
x = do
    t <- a -> a -> a -> Either MathError a
forall a. RealFloat a => a -> a -> a -> Either MathError a
boundedParameter a
a a
b a
x
    result (t * t * (3 - 2 * t))

{- | Given increasing edges then x, apply the quintic 6t^5-15t^4+10t^3 to the clamped
interpolation parameter. Invalid edges/non-finite values fail. Example: smootherstep 0
1 0.5 == Right 0.5.
-}
smootherstep :: forall a. RealFloat a => a -> a -> a -> Either MathError a
smootherstep a
a a
b a
x = do
    t <- a -> a -> a -> Either MathError a
forall a. RealFloat a => a -> a -> a -> Either MathError a
boundedParameter a
a a
b a
x
    result (t * t * t * (t * (t * 6 - 15) + 10))

boundedParameter :: (RealFloat a) => a -> a -> a -> Either MathError a
boundedParameter :: forall a. RealFloat a => a -> a -> a -> Either MathError a
boundedParameter a
a a
b a
x
    | a
a a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
b = MathError -> Either MathError a
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"smoothstep requires increasing edges")
    | Bool
otherwise = a -> a -> a -> Either MathError a
forall a. RealFloat a => a -> a -> a -> Either MathError a
inverseLerp a
a a
b a
x Either MathError a
-> (a -> Either MathError a) -> Either MathError a
forall a b.
Either MathError a
-> (a -> Either MathError b) -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= a -> a -> a -> Either MathError a
forall a. RealFloat a => a -> a -> a -> Either MathError a
clampChecked a
0 a
1

{- | Raise a finite real base to a finite real exponent using (**). Non-real/domain results
and overflow return Left NonFiniteResult; non-finite input fails. Example: power 2 3 ==
Right 8.
-}
power :: (RealFloat a) => a -> a -> Either MathError a
power :: forall a. RealFloat a => a -> a -> Either MathError a
power a
a a
b = [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
a, a
b] Either MathError () -> Either MathError a -> Either MathError a
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> a -> Either MathError a
forall a. RealFloat a => a -> Either MathError a
result (a
a a -> a -> a
forall a. Floating a => a -> a -> a
** a
b)

{- | Raise a finite real base to an Integer exponent, allowing negative powers except for a
zero base. Invalid zero reciprocal/non-finite input fails; overflow is NonFiniteResult.
Example: integerPower 2 (-2) == Right 0.25.
-}
integerPower :: (RealFloat a) => a -> Integer -> Either MathError a
integerPower :: forall a. RealFloat a => a -> Integer -> Either MathError a
integerPower a
a Integer
n = do
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
a]
    if a
a a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 Bool -> Bool -> Bool
&& Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 then MathError -> Either MathError a
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"negative power of zero") else a -> Either MathError a
forall a. RealFloat a => a -> Either MathError a
result (a
a a -> Integer -> a
forall a b. (Fractional a, Integral b) => a -> b -> a
^^ Integer
n)

{- | Given base, nonnegative exponent, positive modulus, compute modular exponentiation by
squaring. Invalid exponent/modulus returns Left DomainError. Example: modularPower 2 10
1000 == Right 24.
-}
modularPower :: Integer -> Integer -> Integer -> Either MathError Integer
modularPower :: Integer -> Integer -> Integer -> Either MathError Integer
modularPower Integer
a Integer
n Integer
m
    | Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 Bool -> Bool -> Bool
|| Integer
m Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
0 =
        MathError -> Either MathError Integer
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"modular power needs nonnegative exponent and positive modulus")
    | Bool
otherwise = Integer -> Either MathError Integer
forall a b. b -> Either a b
Right (Integer -> Integer -> Integer -> Integer
forall {a}. Integral a => Integer -> a -> Integer -> Integer
go (Integer
a Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
m) Integer
n (Integer
1 Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
m))
  where
    go :: Integer -> a -> Integer -> Integer
go Integer
_ a
0 Integer
acc = Integer
acc
    go Integer
b a
k Integer
acc = Integer -> a -> Integer -> Integer
go ((Integer
b Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
b) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
m) (a
k a -> a -> a
forall a. Integral a => a -> a -> a
`quot` a
2) (if a -> Bool
forall a. Integral a => a -> Bool
odd a
k then (Integer
acc Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
b) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
m else Integer
acc)

sqrtChecked, cubeRoot, reciprocalSqrt :: (RealFloat a) => a -> Either MathError a

{- | Return the nonnegative square root of a finite nonnegative input; reject
negatives/non-finite values. Example: sqrtChecked 9 == Right 3.
-}
sqrtChecked :: forall a. RealFloat a => a -> Either MathError a
sqrtChecked a
x = [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
x] Either MathError () -> Either MathError a -> Either MathError a
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> if a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 then MathError -> Either MathError a
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"negative square root") else a -> Either MathError a
forall a. RealFloat a => a -> Either MathError a
result (a -> a
forall a. Floating a => a -> a
sqrt a
x)

{- | Return the signed real cube root of a finite input, including negative values, through
nthRoot 3. Example: cubeRoot (-8) is approximately Right (-2).
-}
cubeRoot :: forall a. RealFloat a => a -> Either MathError a
cubeRoot a
x = Int -> a -> Either MathError a
forall a. RealFloat a => Int -> a -> Either MathError a
nthRoot Int
3 a
x

{- | Return 1/sqrt(x) for finite strictly positive x. Zero/negative/non-finite inputs fail,
as does a non-finite result. Example: reciprocalSqrt 4 == Right 0.5.
-}
reciprocalSqrt :: forall a. RealFloat a => a -> Either MathError a
reciprocalSqrt a
x =
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
x]
        Either MathError () -> Either MathError a -> Either MathError a
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> if a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0
            then MathError -> Either MathError a
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"reciprocal square root needs positive input")
            else a -> Either MathError a
forall a. RealFloat a => a -> Either MathError a
result (a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a
forall a. Floating a => a -> a
sqrt a
x)

{- | Given positive integer degree n then finite x, return its signed real n-th root. Even
roots of negative x and nonpositive degrees fail. Example: nthRoot 3 (-8) is
approximately Right (-2).
-}
nthRoot :: (RealFloat a) => Int -> a -> Either MathError a
nthRoot :: forall a. RealFloat a => Int -> a -> Either MathError a
nthRoot Int
n a
x = do
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
x]
    if Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0 Bool -> Bool -> Bool
|| (a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 Bool -> Bool -> Bool
&& Int -> Bool
forall a. Integral a => a -> Bool
even Int
n)
        then MathError -> Either MathError a
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"invalid real root")
        else if a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 then a -> Either MathError a
forall a b. b -> Either a b
Right a
x else a -> Either MathError a
forall a. RealFloat a => a -> Either MathError a
result (a -> a
forall a. Num a => a -> a
signum a
x a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
exp (a -> a
forall a. Floating a => a -> a
log (a -> a
forall a. Num a => a -> a
abs a
x) a -> a -> a
forall a. Fractional a => a -> a -> a
/ Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n))

expChecked, logChecked :: (RealFloat a) => a -> Either MathError a

{- | Compute exp of a finite input, reporting non-finite output as Left NonFiniteResult.
Finite underflow to zero is allowed. Example: expChecked 0 == Right 1.
-}
expChecked :: forall a. RealFloat a => a -> Either MathError a
expChecked a
x = [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
x] Either MathError () -> Either MathError a -> Either MathError a
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> a -> Either MathError a
forall a. RealFloat a => a -> Either MathError a
result (a -> a
forall a. Floating a => a -> a
exp a
x)

{- | Compute the natural logarithm of finite positive x. Zero/negative/non-finite input
fails. Example: logChecked 1 == Right 0.
-}
logChecked :: forall a. RealFloat a => a -> Either MathError a
logChecked a
x =
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
x] Either MathError () -> Either MathError a -> Either MathError a
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> if a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0 then MathError -> Either MathError a
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"logarithm needs positive input") else a -> Either MathError a
forall a. RealFloat a => a -> Either MathError a
result (a -> a
forall a. Floating a => a -> a
log a
x)

expm1, log1p :: (RealFloat a) => a -> Either MathError a

{- | Compute exp(x)-1, using a convergent series near zero to reduce cancellation. Reject
non-finite input/output; preserve signed zero. Example: expm1 0 == Right 0. This is not
a correct-rounding guarantee.
-}
expm1 :: forall a. RealFloat a => a -> Either MathError a
expm1 a
x = do
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
x]
    if a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 then a -> Either MathError a
forall a b. b -> Either a b
Right a
x else a -> Either MathError a
forall a. RealFloat a => a -> Either MathError a
result (if a -> a
forall a. Num a => a -> a
abs a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0.5 then Int -> a -> a -> a
expSeries Int
1 a
x a
x else a -> a
forall a. Floating a => a -> a
exp a
x a -> a -> a
forall a. Num a => a -> a -> a
- a
1)
  where
    expSeries :: Int -> a -> a -> a
expSeries Int
k a
term a
acc
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
256 Bool -> Bool -> Bool
|| a
next a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
acc = a
acc
        | Bool
otherwise = Int -> a -> a -> a
expSeries (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) a
term' a
next
      where
        term' :: a
term' = a
term a -> a -> a
forall a. Num a => a -> a -> a
* a
x a -> a -> a
forall a. Fractional a => a -> a -> a
/ Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 :: Int)
        next :: a
next = a
acc a -> a -> a
forall a. Num a => a -> a -> a
+ a
term'

{- | Compute log(1+x), using a series near zero. Require finite x greater than -1; preserve
signed zero and reject non-finite output. Example: log1p 0 == Right 0.
-}
log1p :: forall a. RealFloat a => a -> Either MathError a
log1p a
x = do
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
x]
    if a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= -a
1
        then MathError -> Either MathError a
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"log1p needs input greater than -1")
        else
            if a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
                then a -> Either MathError a
forall a b. b -> Either a b
Right a
x
                else a -> Either MathError a
forall a. RealFloat a => a -> Either MathError a
result (if a -> a
forall a. Num a => a -> a
abs a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0.5 then Int -> a -> a -> a
series Int
1 a
y a
y else a -> a
forall a. Floating a => a -> a
log (a
1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
x))
  where
    y :: a
y = a
x a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
2 a -> a -> a
forall a. Num a => a -> a -> a
+ a
x)
    series :: Int -> a -> a -> a
series Int
k a
term a
acc
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
256 Bool -> Bool -> Bool
|| a
next a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
acc = a
2 a -> a -> a
forall a. Num a => a -> a -> a
* a
acc
        | Bool
otherwise = Int -> a -> a -> a
series (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) a
term' a
next
      where
        term' :: a
term' = a
term a -> a -> a
forall a. Num a => a -> a -> a
* a
y a -> a -> a
forall a. Num a => a -> a -> a
* a
y
        next :: a
next = a
acc a -> a -> a
forall a. Num a => a -> a -> a
+ a
term' a -> a -> a
forall a. Fractional a => a -> a -> a
/ Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 :: Int)

{- | Return sqrt(a*a+b*b) with scaling to avoid avoidable intermediate overflow. Require
finite inputs and result. Example: hypot 3 4 == Right 5.
-}
hypot :: (RealFloat a) => a -> a -> Either MathError a
hypot :: forall a. RealFloat a => a -> a -> Either MathError a
hypot a
a a
b = do
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
a, a
b]
    let hi :: a
hi = a -> a -> a
forall a. Ord a => a -> a -> a
max (a -> a
forall a. Num a => a -> a
abs a
a) (a -> a
forall a. Num a => a -> a
abs a
b); lo :: a
lo = a -> a -> a
forall a. Ord a => a -> a -> a
min (a -> a
forall a. Num a => a -> a
abs a
a) (a -> a
forall a. Num a => a -> a
abs a
b)
    if a
hi a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 then a -> Either MathError a
forall a b. b -> Either a b
Right a
0 else a -> Either MathError a
forall a. RealFloat a => a -> Either MathError a
result (a
hi a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
sqrt (a
1 a -> a -> a
forall a. Num a => a -> a -> a
+ (a
lo a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
hi) a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)))

floorChecked
    , ceilChecked
    , roundChecked
    , truncateChecked
    , bankersRound ::
        (RealFloat a) => a -> Either MathError Integer

{- | Convert a finite real to the greatest Integer no larger than it. Non-finite input
fails. Example: floorChecked (-1.2) == Right (-2).
-}
floorChecked :: forall a. RealFloat a => a -> Either MathError Integer
floorChecked a
x = [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
x] Either MathError ()
-> Either MathError Integer -> Either MathError Integer
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Integer -> Either MathError Integer
forall a b. b -> Either a b
Right (a -> Integer
forall b. Integral b => a -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor a
x)

{- | Convert a finite real to the least Integer no smaller than it. Non-finite input fails.
Example: ceilChecked 1.2 == Right 2.
-}
ceilChecked :: forall a. RealFloat a => a -> Either MathError Integer
ceilChecked a
x = [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
x] Either MathError ()
-> Either MathError Integer -> Either MathError Integer
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Integer -> Either MathError Integer
forall a b. b -> Either a b
Right (a -> Integer
forall b. Integral b => a -> b
forall a b. (RealFrac a, Integral b) => a -> b
ceiling a
x)

{- | Round a finite real to Integer with ties to even, as bankersRound. Non-finite input
fails. Example: roundChecked 2.5 == Right 2.
-}
roundChecked :: forall a. RealFloat a => a -> Either MathError Integer
roundChecked = a -> Either MathError Integer
forall a. RealFloat a => a -> Either MathError Integer
bankersRound

{- | Convert finite real to Integer toward zero. Non-finite input fails. Example:
truncateChecked (-1.8) == Right (-1).
-}
truncateChecked :: forall a. RealFloat a => a -> Either MathError Integer
truncateChecked a
x = [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
x] Either MathError ()
-> Either MathError Integer -> Either MathError Integer
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Integer -> Either MathError Integer
forall a b. b -> Either a b
Right (a -> Integer
forall b. Integral b => a -> b
forall a b. (RealFrac a, Integral b) => a -> b
truncate a
x)

{- | Round a finite real to nearest Integer, resolving exact half ties to even. Non-finite
input fails. Example: bankersRound 3.5 == Right 4.
-}
bankersRound :: forall a. RealFloat a => a -> Either MathError Integer
bankersRound a
x = [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
x] Either MathError ()
-> Either MathError Integer -> Either MathError Integer
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Integer -> Either MathError Integer
forall a b. b -> Either a b
Right (a -> Integer
forall b. Integral b => a -> b
forall a b. (RealFrac a, Integral b) => a -> b
round a
x)

{- | Decimal-place rounding with exact rational scaling and ties to even. Negative
digits round to tens, hundreds, etc. Digits are bounded to avoid huge allocation.
-}
roundDigits :: (RealFloat a) => Int -> a -> Either MathError a
roundDigits :: forall a. RealFloat a => Int -> a -> Either MathError a
roundDigits Int
n a
x = do
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
x]
    if Integer -> Integer
forall a. Num a => a -> a
abs (Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
n) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
10000
        then MathError -> Either MathError a
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"rounding digit count exceeds 10000")
        else
            let scale :: Rational
scale = if Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 then Integer -> Rational
forall a. Num a => Integer -> a
fromInteger (Integer
10 Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
n) else Rational
1 Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Integer -> Rational
forall a. Num a => Integer -> a
fromInteger (Integer
10 Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ Int -> Int
forall a. Num a => a -> a
negate Int
n)
             in a -> Either MathError a
forall a. RealFloat a => a -> Either MathError a
result (Rational -> a
forall a. Fractional a => Rational -> a
fromRational (Integer -> Rational
forall a. Num a => Integer -> a
fromInteger (Rational -> Integer
forall b. Integral b => Rational -> b
forall a b. (RealFrac a, Integral b) => a -> b
round (a -> Rational
forall a. Real a => a -> Rational
toRational a
x Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
scale)) Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
scale))

{- | Split finite x into its truncated Integer and signed fractional remainder. Non-finite
input fails. Example: integerFraction (-1.5) == Right (-1,-0.5).
-}
integerFraction :: (RealFloat a) => a -> Either MathError (Integer, a)
integerFraction :: forall a. RealFloat a => a -> Either MathError (Integer, a)
integerFraction a
x = do
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
x]
    let n :: Integer
n = a -> Integer
forall b. Integral b => a -> b
forall a b. (RealFrac a, Integral b) => a -> b
truncate a
x
    (Integer, a) -> Either MathError (Integer, a)
forall a b. b -> Either a b
Right (Integer
n, if a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 then a
x else a
x a -> a -> a
forall a. Num a => a -> a -> a
- Integer -> a
forall a. Num a => Integer -> a
fromInteger Integer
n)

{- | Given integral value then strictly positive modulus, return the nonnegative remainder.
Nonpositive modulus fails. Example: positiveModulo (-1) 5 == Right 4.
-}
positiveModulo :: (Integral a) => a -> a -> Either MathError a
positiveModulo :: forall a. Integral a => a -> a -> Either MathError a
positiveModulo a
_ a
m | a
m a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0 = MathError -> Either MathError a
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"modulus must be positive")
positiveModulo a
a a
m = a -> Either MathError a
forall a b. b -> Either a b
Right (a
a a -> a -> a
forall a. Integral a => a -> a -> a
`mod` a
m)

absoluteError, relativeError :: (RealFloat a) => a -> a -> Either MathError a

{- | Given actual then expected finite values, return abs(actual-expected); reject
non-finite input/result. Example: absoluteError 3 2 == Right 1.
-}
absoluteError :: forall a. RealFloat a => a -> a -> Either MathError a
absoluteError a
actual a
expected = [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
actual, a
expected] Either MathError () -> Either MathError a -> Either MathError a
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> a -> Either MathError a
forall a. RealFloat a => a -> Either MathError a
result (a -> a
forall a. Num a => a -> a
abs (a
actual a -> a -> a
forall a. Num a => a -> a -> a
- a
expected))

{- | Given actual then expected, return absolute error divided by abs(expected). Expected
must be nonzero and inputs/results finite. Example: relativeError 3 2 == Right 0.5.
-}
relativeError :: forall a. RealFloat a => a -> a -> Either MathError a
relativeError a
actual a
expected =
    if a
expected a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
        then MathError -> Either MathError a
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"zero reference for relative error")
        else a -> a -> Either MathError a
forall a. RealFloat a => a -> a -> Either MathError a
absoluteError a
actual a
expected Either MathError a
-> (a -> Either MathError a) -> Either MathError a
forall a b.
Either MathError a
-> (a -> Either MathError b) -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= \a
e -> a -> Either MathError a
forall a. RealFloat a => a -> Either MathError a
result (a
e a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a
forall a. Num a => a -> a
abs a
expected)

nonnegative :: Integer -> Either MathError ()
nonnegative :: Integer -> Either MathError ()
nonnegative Integer
n
    | Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 = MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"negative integer argument")
    | Bool
otherwise = () -> Either MathError ()
forall a b. b -> Either a b
Right ()
factorialChecked
    , doubleFactorial
    , fibonacciChecked
    , lucas
    , catalan ::
        Integer -> Either MathError Integer

{- | Return exact n! for nonnegative Integer n; negatives fail with DomainError. Example:
factorialChecked 5 == Right 120. No fixed-size arithmetic overflow occurs.
-}
factorialChecked :: Integer -> Either MathError Integer
factorialChecked Integer
n = Integer -> Either MathError ()
nonnegative Integer
n Either MathError ()
-> Either MathError Integer -> Either MathError Integer
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Integer -> Either MathError Integer
forall a b. b -> Either a b
Right ([Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer
1 .. Integer
n])

{- | Return n*(n-2)*... exactly for n >= 0; (-1)!! is also 1. Other negative n fail.
Example: doubleFactorial 5 == Right 15.
-}
doubleFactorial :: Integer -> Either MathError Integer
doubleFactorial Integer
n
    | Integer
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== -Integer
1 = Integer -> Either MathError Integer
forall a b. b -> Either a b
Right Integer
1
    | Bool
otherwise = Integer -> Either MathError ()
nonnegative Integer
n Either MathError ()
-> Either MathError Integer -> Either MathError Integer
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Integer -> Either MathError Integer
forall a b. b -> Either a b
Right ([Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer
n, Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
2 .. Integer
1])

risingFactorial, fallingFactorial :: Integer -> Integer -> Either MathError Integer

{- | Given integer x then nonnegative length n, return x*(x+1)*...*(x+n-1). Negative length
fails; length zero gives 1. Example: risingFactorial 3 2 == Right 12.
-}
risingFactorial :: Integer -> Integer -> Either MathError Integer
risingFactorial Integer
x Integer
n = Integer -> Either MathError ()
nonnegative Integer
n Either MathError ()
-> Either MathError Integer -> Either MathError Integer
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Integer -> Either MathError Integer
forall a b. b -> Either a b
Right ([Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer
x .. Integer
x Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
1])

{- | Given integer x then nonnegative length n, return x*(x-1)*...*(x-n+1). Negative length
fails; length zero gives 1. Example: fallingFactorial 3 2 == Right 6.
-}
fallingFactorial :: Integer -> Integer -> Either MathError Integer
fallingFactorial Integer
x Integer
n = Integer -> Either MathError ()
nonnegative Integer
n Either MathError ()
-> Either MathError Integer -> Either MathError Integer
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Integer -> Either MathError Integer
forall a b. b -> Either a b
Right ([Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer
x Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1 .. Integer
x])

{- | Return exact zero-based Fibonacci number using fast doubling (F0=0,F1=1). Negative
index fails. Example: fibonacciChecked 10 == Right 55.
-}
fibonacciChecked :: Integer -> Either MathError Integer
fibonacciChecked Integer
n = Integer -> Either MathError ()
nonnegative Integer
n Either MathError ()
-> Either MathError Integer -> Either MathError Integer
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Integer -> Either MathError Integer
forall a b. b -> Either a b
Right ((Integer, Integer) -> Integer
forall a b. (a, b) -> a
fst (Integer -> (Integer, Integer)
fibPair Integer
n))

fibPair :: Integer -> (Integer, Integer)
fibPair :: Integer -> (Integer, Integer)
fibPair Integer
0 = (Integer
0, Integer
1)
fibPair Integer
n =
    let (Integer
a, Integer
b) = Integer -> (Integer, Integer)
fibPair (Integer
n Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`quot` Integer
2); c :: Integer
c = Integer
a Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* (Integer
2 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
b Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
a); d :: Integer
d = Integer
a Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
a Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
b Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
b
     in if Integer -> Bool
forall a. Integral a => a -> Bool
even Integer
n then (Integer
c, Integer
d) else (Integer
d, Integer
c Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
d)

{- | Return exact zero-based Lucas number (L0=2,L1=1). Negative index fails. Example: lucas
5 == Right 11.
-}
lucas :: Integer -> Either MathError Integer
lucas Integer
n = Integer -> Either MathError ()
nonnegative Integer
n Either MathError ()
-> Either MathError Integer -> Either MathError Integer
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> let (Integer
a, Integer
b) = Integer -> (Integer, Integer)
fibPair Integer
n in Integer -> Either MathError Integer
forall a b. b -> Either a b
Right (Integer
2 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
b Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
a)

permutations, combinations :: Integer -> Integer -> Either MathError Integer

{- | Given nonnegative n,k, return n!/(n-k)!; k>n gives zero. Negative inputs fail. Example:
permutations 4 2 == Right 12.
-}
permutations :: Integer -> Integer -> Either MathError Integer
permutations Integer
n Integer
k = do
    Integer -> Either MathError ()
nonnegative Integer
n
    Integer -> Either MathError ()
nonnegative Integer
k
    if Integer
k Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
n then Integer -> Either MathError Integer
forall a b. b -> Either a b
Right Integer
0 else Integer -> Either MathError Integer
forall a b. b -> Either a b
Right ([Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
k Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1 .. Integer
n])

{- | Given nonnegative n,k, return the exact binomial coefficient; k>n gives zero. Negative
inputs fail. Example: combinations 4 2 == Right 6.
-}
combinations :: Integer -> Integer -> Either MathError Integer
combinations Integer
n Integer
k = do
    Integer -> Either MathError ()
nonnegative Integer
n
    Integer -> Either MathError ()
nonnegative Integer
k
    if Integer
k Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
n
        then Integer -> Either MathError Integer
forall a b. b -> Either a b
Right Integer
0
        else let j :: Integer
j = Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
min Integer
k (Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
k) in Integer -> Either MathError Integer
forall a b. b -> Either a b
Right ([Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
j Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1 .. Integer
n] Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`quot` [Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer
1 .. Integer
j])

{- | Given nonnegative group sizes, return factorial(sum sizes)/product(factorial sizes).
Negative entries fail. Example: multinomial [2,1] == Right 3.
-}
multinomial :: [Integer] -> Either MathError Integer
multinomial :: [Integer] -> Either MathError Integer
multinomial [Integer]
ks = do
    (Integer -> Either MathError ())
-> [Integer] -> Either MathError ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ Integer -> Either MathError ()
nonnegative [Integer]
ks
    f <- Integer -> Either MathError Integer
factorialChecked ([Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Integer]
ks)
    fs <- mapM factorialChecked ks
    Right (f `quot` product fs)

{- | Return exact n-th Catalan number for nonnegative Integer n; negative input fails.
Example: catalan 3 == Right 5.
-}
catalan :: Integer -> Either MathError Integer
catalan Integer
n = Integer -> Either MathError ()
nonnegative Integer
n Either MathError ()
-> Either MathError Integer -> Either MathError Integer
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> (Integer -> Integer)
-> Either MathError Integer -> Either MathError Integer
forall a b. (a -> b) -> Either MathError a -> Either MathError b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`quot` (Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1)) (Integer -> Integer -> Either MathError Integer
combinations (Integer
2 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
n) Integer
n)

-- | Bernoulli convention B1 = -1/2. Akiyama-Tanigawa exact rational recurrence.
bernoulli :: Int -> Either MathError Rational
bernoulli :: Int -> Either MathError Rational
bernoulli Int
n
    | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = MathError -> Either MathError Rational
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"negative Bernoulli index")
    | Bool
otherwise = Rational -> Either MathError Rational
forall a b. b -> Either a b
Right (if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 then -Rational
1 Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
2 else [Rational] -> Rational
forall a. Num a => [a] -> a
first (([Rational] -> Int -> [Rational])
-> [Rational] -> [Int] -> [Rational]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
L.foldl' [Rational] -> Int -> [Rational]
forall {a}. Integral a => [Ratio a] -> Int -> [Ratio a]
advance [] [Int
0 .. Int
n]))
  where
    first :: [a] -> a
first (a
x : [a]
_) = a
x
    first [] = a
0
    advance :: [Ratio a] -> Int -> [Ratio a]
advance [Ratio a]
acc Int
m = [Ratio a] -> Int -> [Ratio a]
forall {a}. Num a => [a] -> Int -> [a]
update ([Ratio a]
acc [Ratio a] -> [Ratio a] -> [Ratio a]
forall a. [a] -> [a] -> [a]
++ [a
1 a -> a -> Ratio a
forall a. Integral a => a -> a -> Ratio a
% Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)]) Int
m
    update :: [a] -> Int -> [a]
update [a]
xs Int
0 = [a]
xs
    update [a]
xs Int
j = [a] -> Int -> [a]
update (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) [a]
xs [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
j a -> a -> a
forall a. Num a => a -> a -> a
* ([a]
xs [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) a -> a -> a
forall a. Num a => a -> a -> a
- [a]
xs [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
j)] [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
j [a]
xs) (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)

{- | Return exact Rational sum 1+1/2+...+1/n; n=0 gives 0 and negatives fail. Example:
harmonic 2 == Right (3/2).
-}
harmonic :: Integer -> Either MathError Rational
harmonic :: Integer -> Either MathError Rational
harmonic Integer
n = Integer -> Either MathError ()
nonnegative Integer
n Either MathError ()
-> Either MathError Rational -> Either MathError Rational
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Rational -> Either MathError Rational
forall a b. b -> Either a b
Right ([Rational] -> Rational
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Integer
1 Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% Integer
k | Integer
k <- [Integer
1 .. Integer
n]])

cumulativeSum, prefixSum, cumulativeProduct :: (Num a) => [a] -> [a]

{- | Return inclusive running sums without a leading zero; output has input length.
Arithmetic follows Num. Example: cumulativeSum [1,2,3] == [1,3,6].
-}
cumulativeSum :: forall a. Num a => [a] -> [a]
cumulativeSum = Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
1 ([a] -> [a]) -> ([a] -> [a]) -> [a] -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a -> a) -> a -> [a] -> [a]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl a -> a -> a
forall a. Num a => a -> a -> a
(+) a
0

{- | Return running sums including a leading zero; output length is input length plus one.
Example: prefixSum [1,2] == [0,1,3].
-}
prefixSum :: forall a. Num a => [a] -> [a]
prefixSum = (a -> a -> a) -> a -> [a] -> [a]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl a -> a -> a
forall a. Num a => a -> a -> a
(+) a
0

{- | Return inclusive running products without a leading one. Arithmetic follows Num.
Example: cumulativeProduct [2,3,4] == [2,6,24].
-}
cumulativeProduct :: forall a. Num a => [a] -> [a]
cumulativeProduct = Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
1 ([a] -> [a]) -> ([a] -> [a]) -> [a] -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a -> a) -> a -> [a] -> [a]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl a -> a -> a
forall a. Num a => a -> a -> a
(*) a
1

{- | Software single-rounding FMA for finite inputs using exact 'Rational'. This
avoids intermediate overflow and cancellation; it does not change CPU flags.
-}
fusedMultiplyAdd :: (RealFloat a) => a -> a -> a -> Either MathError a
fusedMultiplyAdd :: forall a. RealFloat a => a -> a -> a -> Either MathError a
fusedMultiplyAdd a
a a
b a
c = do
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
a, a
b, a
c]
    let r :: Rational
r = a -> Rational
forall a. Real a => a -> Rational
toRational a
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* a -> Rational
forall a. Real a => a -> Rational
toRational a
b Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ a -> Rational
forall a. Real a => a -> Rational
toRational a
c
    a -> Either MathError a
forall a. RealFloat a => a -> Either MathError a
result (if Rational
r Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational
0 Bool -> Bool -> Bool
&& a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
b a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 Bool -> Bool -> Bool
&& a
c a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 then a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
b a -> a -> a
forall a. Num a => a -> a -> a
+ a
c else Rational -> a
forall a. Fractional a => Rational -> a
fromRational Rational
r)

saturatingAdd, saturatingSubtract, saturatingMultiply :: (Integral a, Bounded a) => a -> a -> a

{- | Add bounded integers through an exact Integer intermediate and clamp to the destination
minBound/maxBound. Example: saturatingAdd (127 :: Int8) 1 == 127.
-}
saturatingAdd :: forall a. (Integral a, Bounded a) => a -> a -> a
saturatingAdd a
a a
b = Integer -> a
forall a. (Integral a, Bounded a) => Integer -> a
saturate (a -> Integer
forall a. Integral a => a -> Integer
toInteger a
a Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ a -> Integer
forall a. Integral a => a -> Integer
toInteger a
b)

{- | Subtract bounded integers through Integer and clamp to type bounds instead of wrapping.
Example: saturatingSubtract (0 :: Word8) 1 == 0.
-}
saturatingSubtract :: forall a. (Integral a, Bounded a) => a -> a -> a
saturatingSubtract a
a a
b = Integer -> a
forall a. (Integral a, Bounded a) => Integer -> a
saturate (a -> Integer
forall a. Integral a => a -> Integer
toInteger a
a Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- a -> Integer
forall a. Integral a => a -> Integer
toInteger a
b)

{- | Multiply bounded integers through Integer and clamp to type bounds instead of wrapping.
Example: saturatingMultiply (100 :: Int8) 2 == 127.
-}
saturatingMultiply :: forall a. (Integral a, Bounded a) => a -> a -> a
saturatingMultiply a
a a
b = Integer -> a
forall a. (Integral a, Bounded a) => Integer -> a
saturate (a -> Integer
forall a. Integral a => a -> Integer
toInteger a
a Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* a -> Integer
forall a. Integral a => a -> Integer
toInteger a
b)

saturate :: (Integral a, Bounded a) => Integer -> a
saturate :: forall a. (Integral a, Bounded a) => Integer -> a
saturate Integer
n = a
x
  where
    x :: a
x = Integer -> a
forall a. Num a => Integer -> a
fromInteger (Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
max (a -> Integer
forall a. Integral a => a -> Integer
toInteger (a
forall a. Bounded a => a
minBound a -> a -> a
forall a. a -> a -> a
`asTypeOf` a
x)) (Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
min (a -> Integer
forall a. Integral a => a -> Integer
toInteger (a
forall a. Bounded a => a
maxBound a -> a -> a
forall a. a -> a -> a
`asTypeOf` a
x)) Integer
n))
checkedAdd
    , checkedSubtract
    , checkedMultiply ::
        (Integral a, Bounded a) => a -> a -> Either MathError a

{- | Add bounded integers using an exact intermediate. Out-of-range results return Left
DomainError, never wrap. Example: checkedAdd (2 :: Int8) 3 == Right 5.
-}
checkedAdd :: forall a. (Integral a, Bounded a) => a -> a -> Either MathError a
checkedAdd a
a a
b = Integer -> Either MathError a
forall a. (Integral a, Bounded a) => Integer -> Either MathError a
boundedResult (a -> Integer
forall a. Integral a => a -> Integer
toInteger a
a Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ a -> Integer
forall a. Integral a => a -> Integer
toInteger a
b)

{- | Subtract bounded integers using an exact intermediate, returning Left DomainError on
underflow/overflow. Example: checkedSubtract (0 :: Word8) 1 fails.
-}
checkedSubtract :: forall a. (Integral a, Bounded a) => a -> a -> Either MathError a
checkedSubtract a
a a
b = Integer -> Either MathError a
forall a. (Integral a, Bounded a) => Integer -> Either MathError a
boundedResult (a -> Integer
forall a. Integral a => a -> Integer
toInteger a
a Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- a -> Integer
forall a. Integral a => a -> Integer
toInteger a
b)

{- | Multiply bounded integers using an exact intermediate and reject out-of-range results.
Example: checkedMultiply (100 :: Int8) 2 fails.
-}
checkedMultiply :: forall a. (Integral a, Bounded a) => a -> a -> Either MathError a
checkedMultiply a
a a
b = Integer -> Either MathError a
forall a. (Integral a, Bounded a) => Integer -> Either MathError a
boundedResult (a -> Integer
forall a. Integral a => a -> Integer
toInteger a
a Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* a -> Integer
forall a. Integral a => a -> Integer
toInteger a
b)

boundedResult :: (Integral a, Bounded a) => Integer -> Either MathError a
boundedResult :: forall a. (Integral a, Bounded a) => Integer -> Either MathError a
boundedResult Integer
n = Either MathError a
r
  where
    r :: Either MathError a
r =
        if Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< a -> Integer
forall a. Integral a => a -> Integer
toInteger (a
forall a. Bounded a => a
minBound a -> a -> a
forall a. a -> a -> a
`asTypeOf` a
x) Bool -> Bool -> Bool
|| Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> a -> Integer
forall a. Integral a => a -> Integer
toInteger (a
forall a. Bounded a => a
maxBound a -> a -> a
forall a. a -> a -> a
`asTypeOf` a
x)
            then MathError -> Either MathError a
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"integer overflow")
            else a -> Either MathError a
forall a b. b -> Either a b
Right a
x
    x :: a
x = Integer -> a
forall a. Num a => Integer -> a
fromInteger Integer
n

-- | Report finite overflow and any nonzero exact product rounded to zero.
checkedFloatMultiply :: (RealFloat a) => a -> a -> Either MathError a
checkedFloatMultiply :: forall a. RealFloat a => a -> a -> Either MathError a
checkedFloatMultiply a
a a
b = do
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
inputs [a
a, a
b]
    v <- a -> Either MathError a
forall a. RealFloat a => a -> Either MathError a
result (a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
b)
    if v == 0 && a /= 0 && b /= 0 then Left (DomainError "floating-point underflow") else Right v