{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ScopedTypeVariables #-}
module Numeric.Rounded.Hardware.Interval.ElementaryFunctions where
import Control.Exception (assert)
import Numeric.Rounded.Hardware.Internal
import Numeric.Rounded.Hardware.Interval.Class
sqrtI :: (IsInterval i, RoundedSqrt (EndPoint i)) => i -> i
sqrtI = withEndPoints $ \x y -> case intervalSqrt x y of (u, v) -> makeInterval u v
{-# INLINE sqrtI #-}
expP :: forall i. (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RealFloatConstants (EndPoint i), RoundedFractional (EndPoint i)) => EndPoint i -> i
expP x | isInfinite x = if x > 0
then makeInterval (Rounded maxFinite) (Rounded positiveInfinity)
else makeInterval (Rounded 0) (Rounded minPositive)
expP x = let a = round x
b = x - fromInteger a
b' = singleton b
series :: Int -> i -> i
series n acc | n == 0 = acc
| otherwise = series (n-1) $ 1 + acc * b' / fromIntegral n
in assert (fromInteger a + b == x && abs b <= 0.5) $
(makeInterval exp1_down exp1_up)^^a * series 15 (makeInterval expm1_2_down exp1_2_up)
{-# INLINABLE expP #-}
expI :: (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RealFloatConstants (EndPoint i), RoundedFractional (EndPoint i)) => i -> i
expI = withEndPoints (\(Rounded x) (Rounded y) -> hull (expP x) (expP y))
{-# INLINABLE expI #-}
expm1P :: forall i. (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RealFloatConstants (EndPoint i), RoundedFractional (EndPoint i)) => EndPoint i -> i
expm1P x | -0.5 <= x && x <= 0.5 = let b' = singleton x
series :: Int -> i -> i
series n acc | n == 1 = acc * b'
| otherwise = series (n-1) $ 1 + acc * b' / fromIntegral n
in series 15 (makeInterval expm1_2_down exp1_2_up)
| otherwise = expP x - 1
{-# INLINABLE expm1P #-}
expm1I :: (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RealFloatConstants (EndPoint i), RoundedFractional (EndPoint i)) => i -> i
expm1I = withEndPoints (\(Rounded x) (Rounded y) -> hull (expm1P x) (expm1P y))
{-# INLINABLE expm1I #-}
logP :: forall i. (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RealFloatConstants (EndPoint i)) => EndPoint i -> i
logP x
| floatRadix (undefined :: EndPoint i) /= 2 = error "Unsupported float radix"
| x < 0 = error "Negative logarithm"
| x == 0 = makeInterval (Rounded negativeInfinity) (Rounded (-maxFinite))
| isInfinite x = makeInterval (Rounded maxFinite) (Rounded positiveInfinity)
| isNaN x = error "NaN"
| otherwise = let m :: Integer
n :: Int
(m,n) = decodeFloat x
a0, b, bm1 :: i
a0 = fromIntegral (n + d)
x' :: EndPoint i
x' = encodeFloat m (- d)
(a,b) | 0.5 <= x' && x' <= getRounded two_minus_sqrt2_down = (a0 - 1, singleton x' * 2)
| getRounded two_minus_sqrt2_up <= x' && x' <= 2 * getRounded sqrt2m1_down = (a0 - 0.5, singleton x' * sqrt2_iv)
| 2 * getRounded sqrt2m1_up <= x' && x' < 1 = (a0, singleton x')
| otherwise = error "interval log: internal error"
bm1 = b - 1
series :: Int -> i -> i
series k acc | k == 0 = bm1 * acc
| otherwise = series (k-1) $ recip (fromIntegral k) - bm1 * acc
in a * log2_iv + series 21 (hull 1 b ^^ (-22 :: Int) * bm1 / fromInteger 22)
where
d = floatDigits (undefined :: EndPoint i)
log2_iv :: i
log2_iv = makeInterval log2_down log2_up
sqrt2_iv :: i
sqrt2_iv = makeInterval sqrt2_down sqrt2_up
{-# INLINABLE logP #-}
logI :: (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RealFloatConstants (EndPoint i)) => i -> i
logI = withEndPoints (\(Rounded x) (Rounded y) -> hull (logP x) (logP y))
{-# INLINABLE logI #-}
log1pP :: forall i. (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RealFloatConstants (EndPoint i)) => EndPoint i -> i
log1pP x | - getRounded three_minus_2sqrt2_down <= x && x <= getRounded three_minus_2sqrt2_down =
let x' :: i
x' = singleton x
series :: Int -> i -> i
series k acc | k == 0 = x' * acc
| otherwise = series (k-1) $ recip (fromIntegral k) - x' * acc
in series 21 (hull 1 (x' + 1) ^^ (-22 :: Int) * x' / fromInteger 22)
| otherwise = logI (singleton x + 1)
{-# INLINABLE log1pP #-}
log1pI :: (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RealFloatConstants (EndPoint i)) => i -> i
log1pI = withEndPoints (\(Rounded x) (Rounded y) -> hull (log1pP x) (log1pP y))
{-# INLINABLE log1pI #-}
sin_small :: forall i. (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RoundedRing (EndPoint i), RealFloatConstants (EndPoint i)) => i -> i
sin_small x = let xx = x * x
series :: Int -> i -> i
series k acc | k == 0 = x * acc
| otherwise = series (k-2) $ 1 - xx * acc / fromIntegral (k * (k+1))
in series 18 (makeInterval (-1) 1)
{-# INLINABLE sin_small #-}
cos_small :: forall i. (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RoundedRing (EndPoint i), RealFloatConstants (EndPoint i)) => i -> i
cos_small x = let xx = x * x
series :: Int -> i -> i
series k acc | k == 1 = acc
| otherwise = series (k-2) $ 1 - xx * acc / fromIntegral ((k-1) * (k-2))
in series 17 (makeInterval (-1) 1)
{-# INLINABLE cos_small #-}
sinP :: forall i. (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RoundedRing (EndPoint i), RealFloatConstants (EndPoint i)) => i -> i
sinP x | x `weaklyLess` - three_pi_iv / 4 = - sin_small (x + pi_iv)
| x `weaklyLess` - pi_iv / 2 = - cos_small (- pi_iv / 2 - x)
| x `weaklyLess` - pi_iv / 4 = - cos_small (x + pi_iv / 2)
| x `weaklyLess` 0 = - sin_small (- x)
| x `weaklyLess` pi_iv / 4 = sin_small x
| x `weaklyLess` pi_iv / 2 = cos_small (pi_iv / 2 - x)
| x `weaklyLess` three_pi_iv / 4 = cos_small (x - pi_iv / 2)
| otherwise = sin_small (pi_iv - x)
where
pi_iv :: i
pi_iv = makeInterval pi_down pi_up
three_pi_iv :: i
three_pi_iv = makeInterval three_pi_down three_pi_up
{-# INLINABLE sinP #-}
cosP :: forall i. (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RoundedRing (EndPoint i), RealFloatConstants (EndPoint i)) => i -> i
cosP x | x `weaklyLess` - three_pi_iv / 4 = - cos_small (x + pi_iv)
| x `weaklyLess` - pi_iv / 2 = - sin_small (- pi_iv / 2 - x)
| x `weaklyLess` - pi_iv / 4 = sin_small (x + pi_iv / 2)
| x `weaklyLess` 0 = cos_small (- x)
| x `weaklyLess` pi_iv / 4 = cos_small x
| x `weaklyLess` pi_iv / 2 = sin_small (pi_iv / 2 - x)
| x `weaklyLess` three_pi_iv / 4 = - sin_small (x - pi_iv / 2)
| otherwise = - cos_small (pi_iv - x)
where
pi_iv :: i
pi_iv = makeInterval pi_down pi_up
three_pi_iv :: i
three_pi_iv = makeInterval three_pi_down three_pi_up
{-# INLINABLE cosP #-}
sinI :: forall i. (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RoundedRing (EndPoint i), RealFloatConstants (EndPoint i)) => i -> i
sinI t = flip withEndPoints t $ \(Rounded x) (Rounded y) ->
if isInfinite x || isInfinite y
then wholeRange
else flip withEndPoints (singleton x / (2 * pi_iv)) $ \(Rounded x0) (Rounded _) ->
let n = round x0
t' = t - 2 * pi_iv * fromInteger n
in flip withEndPoints t' $ \(Rounded x') (Rounded y') ->
if y' - x' >= getRounded (2 * pi_up)
then wholeRange
else
let include_minus_1 = minus_half_pi_iv `subset` t' || three_pi_2_iv `subset` t'
include_plus_1 = pi_iv / 2 `subset` t' || five_pi_2_iv `subset` t'
u = hull (sinP $ singleton x') $ sinP (if y <= getRounded pi_down then singleton y' else singleton y' - 2 * pi_iv)
v | include_minus_1 = hull (-1) u
| otherwise = u
w | include_plus_1 = hull 1 v
| otherwise = v
in intersection wholeRange w
where
pi_iv :: i
pi_iv = makeInterval pi_down pi_up
minus_half_pi_iv :: i
minus_half_pi_iv = - pi_iv / 2
three_pi_2_iv :: i
three_pi_2_iv = makeInterval three_pi_down three_pi_up / 2
five_pi_2_iv :: i
five_pi_2_iv = makeInterval five_pi_down five_pi_up / 2
wholeRange :: i
wholeRange = makeInterval (-1) 1
{-# INLINABLE sinI #-}
cosI :: forall i. (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RoundedRing (EndPoint i), RealFloatConstants (EndPoint i)) => i -> i
cosI t = flip withEndPoints t $ \(Rounded x) (Rounded y) ->
if isInfinite x || isInfinite y
then wholeRange
else flip withEndPoints (singleton x / (2 * pi_iv)) $ \(Rounded x0) (Rounded _) ->
let n = round x0
t' = t - 2 * pi_iv * fromInteger n
in flip withEndPoints t' $ \(Rounded x') (Rounded y') ->
if y' - x' >= getRounded (2 * pi_up)
then wholeRange
else
let include_minus_1 = -pi_iv `subset` t' || pi_iv `subset` t' || three_pi_iv `subset` t'
include_plus_1 = 0 `subset` t' || 2 * pi_iv `subset` t'
u = hull (cosP $ singleton x') $ cosP (if y <= getRounded pi_down then singleton y' else singleton y' - 2 * pi_iv)
v | include_minus_1 = hull (-1) u
| otherwise = u
w | include_plus_1 = hull 1 v
| otherwise = v
in intersection wholeRange w
where
pi_iv :: i
pi_iv = makeInterval pi_down pi_up
three_pi_iv :: i
three_pi_iv = makeInterval three_pi_down three_pi_up
wholeRange :: i
wholeRange = makeInterval (-1) 1
{-# INLINABLE cosI #-}
tanI :: forall i. (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RoundedRing (EndPoint i), RealFloatConstants (EndPoint i)) => i -> i
tanI t = flip withEndPoints t $ \(Rounded x) (Rounded y) ->
if isInfinite x || isInfinite y
then wholeRange
else flip withEndPoints (t / pi_iv) $ \(Rounded x0) (Rounded _) ->
let n = round x0
t' = t - pi_iv * fromInteger n
in flip withEndPoints t' $ \(Rounded x') (Rounded y') ->
if y' >= getRounded pi_up / 2
then wholeRange
else let lb = sinP (singleton x') / cosP (singleton x')
ub = sinP (singleton y') / cosP (singleton y')
in hull lb ub
where
pi_iv :: i
pi_iv = makeInterval pi_down pi_up
wholeRange :: i
wholeRange = makeInterval (Rounded negativeInfinity) (Rounded positiveInfinity)
{-# INLINABLE tanI #-}
atan_small :: forall i. (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RoundedRing (EndPoint i), RealFloatConstants (EndPoint i)) => i -> i
atan_small x = let n = 39
in series n (makeInterval (-1) 1 / fromIntegral n)
where
xx = x * x
series :: Int -> i -> i
series k acc | k == 1 = x * acc
| otherwise = series (k-2) $ recip (fromIntegral (k-2)) - xx * acc
{-# INLINABLE atan_small #-}
atanP :: forall i. (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RoundedRing (EndPoint i), RealFloatConstants (EndPoint i)) => EndPoint i -> i
atanP x | getRounded (1 + sqrt2_up) <= x = pi_iv / 2 - atan_small (recip x')
| getRounded sqrt2m1_up <= x = pi_iv / 4 + atan_small ((x' - 1) / (x' + 1))
| - getRounded sqrt2m1_down <= x = atan_small x'
| - getRounded (sqrt2_down + 1) <= x = - pi_iv / 4 + atan_small ((1 + x') / (1 - x'))
| otherwise = - pi_iv / 2 - atan_small (recip x')
where
x' = singleton x
pi_iv :: i
pi_iv = makeInterval pi_down pi_up
{-# INLINABLE atanP #-}
atanI :: forall i. (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RoundedRing (EndPoint i), RealFloatConstants (EndPoint i)) => i -> i
atanI = withEndPoints (\(Rounded x) (Rounded y) -> hull (atanP x) (atanP y))
{-# INLINABLE atanI #-}
asinP :: forall i. (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RoundedRing (EndPoint i), RoundedSqrt (EndPoint i), RealFloatConstants (EndPoint i)) => EndPoint i -> i
asinP x | x < -1 || 1 < x = error "asin"
| x == -1 = - pi_iv / 2
| x == 1 = pi_iv / 2
| otherwise = atanI (x' / sqrtI (1 - x'*x'))
where
x' = singleton x
pi_iv :: i
pi_iv = makeInterval pi_down pi_up
{-# INLINABLE asinP #-}
asinI :: forall i. (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RoundedRing (EndPoint i), RoundedSqrt (EndPoint i), RealFloatConstants (EndPoint i)) => i -> i
asinI = withEndPoints (\(Rounded x) (Rounded y) -> hull (asinP x) (asinP y))
{-# INLINABLE asinI #-}
acosP :: forall i. (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RoundedRing (EndPoint i), RoundedSqrt (EndPoint i), RealFloatConstants (EndPoint i)) => EndPoint i -> i
acosP x | x < -1 || 1 < x = error "asin"
| x == -1 = pi_iv
| x == 1 = 0
| otherwise = case x' / sqrtI (1 - x'*x') of
y' | one_plus_sqrt2 `weaklyLess` y' -> atan_small (recip y')
| sqrt2_minus_one `weaklyLess` y' -> pi_iv / 4 - atan_small ((y' - 1) / (y' + 1))
| - sqrt2_minus_one `weaklyLess` y' -> pi_iv / 2 - atan_small y'
| - one_plus_sqrt2 `weaklyLess` y' -> three_pi_iv / 4 - atan_small ((1 + y') / (1 - y'))
| otherwise -> pi_iv + atan_small (recip y')
where
x' = singleton x
pi_iv :: i
pi_iv = makeInterval pi_down pi_up
three_pi_iv :: i
three_pi_iv = makeInterval three_pi_down three_pi_up
one_plus_sqrt2 :: i
one_plus_sqrt2 = 1 + makeInterval sqrt2_down sqrt2_up
sqrt2_minus_one :: i
sqrt2_minus_one = makeInterval sqrt2m1_down sqrt2m1_up
{-# INLINABLE acosP #-}
acosI :: forall i. (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RoundedRing (EndPoint i), RoundedSqrt (EndPoint i), RealFloatConstants (EndPoint i)) => i -> i
acosI = withEndPoints (\(Rounded x) (Rounded y) -> hull (acosP x) (acosP y))
{-# INLINABLE acosI #-}
sinhP :: (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RealFloatConstants (EndPoint i), RoundedFractional (EndPoint i)) => EndPoint i -> i
sinhP x | x >= 0 = let y = expP x
in (y - recip y) / 2
| otherwise = let y = expP (- x)
in (recip y - y) / 2
{-# INLINABLE sinhP #-}
sinhI :: (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RealFloatConstants (EndPoint i), RoundedFractional (EndPoint i)) => i -> i
sinhI = withEndPoints (\(Rounded x) (Rounded y) -> hull (sinhP x) (sinhP y))
{-# INLINABLE sinhI #-}
coshP :: (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RealFloatConstants (EndPoint i), RoundedFractional (EndPoint i)) => EndPoint i -> i
coshP x | x >= 0 = let y = expP x
in (y + recip y) / 2
| otherwise = let y = expP (- x)
in (recip y + y) / 2
{-# INLINABLE coshP #-}
coshI :: (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RealFloatConstants (EndPoint i), RoundedFractional (EndPoint i)) => i -> i
coshI = withEndPoints $ \(Rounded x) (Rounded y) ->
let z = hull (coshP x) (coshP y)
in if x <= 0 && 0 <= y
then hull 0 z
else z
{-# INLINABLE coshI #-}
tanhP :: (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RealFloatConstants (EndPoint i), RoundedFractional (EndPoint i)) => EndPoint i -> i
tanhP x | -0.5 <= x && x <= 0.5 = sinhP x / coshP x
| 0 < x = 1 - 2 / (1 + expP (2 * x))
| otherwise = 2 / (1 + expP (- 2 * x)) - 1
{-# INLINABLE tanhP #-}
tanhI :: (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RealFloatConstants (EndPoint i), RoundedFractional (EndPoint i)) => i -> i
tanhI = withEndPoints $ \(Rounded x) (Rounded y) -> hull (tanhP x) (tanhP y)
{-# INLINABLE tanhI #-}
asinhP :: (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RealFloatConstants (EndPoint i), RoundedSqrt (EndPoint i)) => EndPoint i -> i
asinhP x = let x' = singleton x
in logI (x' + sqrtI (1 + x' ^ (2 :: Int)))
{-# INLINABLE asinhP #-}
asinhI :: (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RealFloatConstants (EndPoint i), RoundedSqrt (EndPoint i)) => i -> i
asinhI = withEndPoints $ \(Rounded x) (Rounded y) -> hull (asinhP x) (asinhP y)
{-# INLINABLE asinhI #-}
acoshP :: (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RealFloatConstants (EndPoint i), RoundedSqrt (EndPoint i)) => EndPoint i -> i
acoshP x | x < 1 = error "acosh: domain"
| otherwise = let x' = singleton x
in logI (x' + sqrtI (x' ^ (2 :: Int) - 1))
{-# INLINABLE acoshP #-}
acoshI :: (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RealFloatConstants (EndPoint i), RoundedSqrt (EndPoint i)) => i -> i
acoshI = withEndPoints $ \(Rounded x) (Rounded y) -> hull (acoshP x) (acoshP y)
{-# INLINABLE acoshI #-}
atanhP :: (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RealFloatConstants (EndPoint i)) => EndPoint i -> i
atanhP x | x < -1 || 1 < x = error "atanh: domain"
| x == -1 = - makeInterval (Rounded maxFinite) (Rounded positiveInfinity)
| x == 1 = makeInterval (Rounded maxFinite) (Rounded positiveInfinity)
| otherwise = let x' = singleton x
in logI ((1 + x') / (1 - x')) / 2
{-# INLINABLE atanhP #-}
atanhI :: (IsInterval i, Fractional i, Eq (EndPoint i), RealFloat (EndPoint i), RealFloatConstants (EndPoint i)) => i -> i
atanhI = withEndPoints $ \(Rounded x) (Rounded y) -> hull (atanhP x) (atanhP y)
{-# INLINABLE atanhI #-}