{-# LINE 1 "src/Data/Number/Flint/Support/D/Interval/FFI.hsc" #-}
{-|
module      :  Data.Number.Flint.Support.D.Interval.FFI
copyright   :  (c) 2022 Hartmut Monien
license     :  GNU GPL, version 2 or above (see LICENSE)
maintainer  :  hmonien@uni-bonn.de
-}
module Data.Number.Flint.Support.D.Interval.FFI (
  -- * Double-precision interval arithmetic and helpers
    Di (..)
  , CDi (..)
  -- * Basic manipulation
  , di_interval
  , arb_get_di
  , arb_set_di
  , di_print
  , di_randtest2
  , di_randtest
  -- * Arithmetic
  , di_neg
  -- * Fast arithmetic
  , di_fast_add
  , di_fast_sub
  , di_fast_mul
  , di_fast_div
  , di_fast_sqr
  , di_fast_add_d
  , di_fast_sub_d
  , di_fast_mul_d
  , di_fast_div_d
  , di_fast_log_nonnegative
  , di_fast_mid
  , di_fast_ubound_radius
) where

-- Double-precision interval arithmetic and helpers ----------------------------

import System.IO.Unsafe

import Foreign.Ptr
import Foreign.ForeignPtr
import Foreign.C.Types
import Foreign.C.String
import Foreign.Storable

import Text.Printf

import Data.Number.Flint.Flint
import Data.Number.Flint.Arb
import Data.Number.Flint.Arb.Types
import Data.Number.Flint.Arb.Arf
import Data.Number.Flint.Arb.Mag
import Data.Number.Flint.Support.D.Extras




d_inf :: CDouble
d_inf = CDouble
1CDouble -> CDouble -> CDouble
forall a. Fractional a => a -> a -> a
/CDouble
0 :: CDouble 

-- di_t ------------------------------------------------------------------------

data Di = Di {-# UNPACK #-} !(ForeignPtr CDi)
data CDi = CDi CDouble CDouble deriving Int -> CDi -> ShowS
[CDi] -> ShowS
CDi -> String
(Int -> CDi -> ShowS)
-> (CDi -> String) -> ([CDi] -> ShowS) -> Show CDi
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> CDi -> ShowS
showsPrec :: Int -> CDi -> ShowS
$cshow :: CDi -> String
show :: CDi -> String
$cshowList :: [CDi] -> ShowS
showList :: [CDi] -> ShowS
Show

instance Storable CDi where
  sizeOf :: CDi -> Int
sizeOf    CDi
_ = (Int
16)
{-# LINE 66 "src/Data/Number/Flint/Support/D/Interval/FFI.hsc" #-}
  alignment _ = 8
{-# LINE 67 "src/Data/Number/Flint/Support/D/Interval/FFI.hsc" #-}
  peek ptr = CDi
    <$> (\hsc_ptr -> peekByteOff hsc_ptr 0) ptr
{-# LINE 69 "src/Data/Number/Flint/Support/D/Interval/FFI.hsc" #-}
    <*> (\hsc_ptr -> peekByteOff hsc_ptr 8) ptr
{-# LINE 70 "src/Data/Number/Flint/Support/D/Interval/FFI.hsc" #-}
  poke ptr (CDi a b) = do
    (\hsc_ptr -> pokeByteOff hsc_ptr 0) ptr a
{-# LINE 72 "src/Data/Number/Flint/Support/D/Interval/FFI.hsc" #-}
    (\hsc_ptr -> pokeByteOff hsc_ptr 8) ptr b
{-# LINE 73 "src/Data/Number/Flint/Support/D/Interval/FFI.hsc" #-}
    

-- Basic manipulation ----------------------------------------------------------

-- | /di_interval/ /a/ /b/ 
--
-- Returns the interval \([a, b]\). We require that the endpoints are
-- ordered and not NaN.
di_interval :: CDouble -> CDouble -> CDi
di_interval :: CDouble -> CDouble -> CDi
di_interval CDouble
a CDouble
b =
  if CDouble
a CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
<= CDouble
b
    then CDouble -> CDouble -> CDi
CDi CDouble
a CDouble
b
    else String -> CDi
forall a. HasCallStack => String -> a
error (String -> CDi) -> String -> CDi
forall a b. (a -> b) -> a -> b
$ String -> Double -> Double -> String
forall r. PrintfType r => String -> r
printf String
"di_interval endpoints %g, %g not ordered.\n"
                        (CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac CDouble
a :: Double)
                        (CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac CDouble
b :: Double)

_di_below :: CDouble -> CDouble
_di_below CDouble
x =
  if CDouble
x CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
<= CDouble
1e300 then
    CDouble
x CDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
- (CDouble
1e-300 CDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
+ if CDouble
x CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
< CDouble
0 then -CDouble
x else CDouble
x) CDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
* CDouble
4.440892098500626e-16
  else
    if CDouble
x CDouble -> CDouble -> Bool
forall a. Eq a => a -> a -> Bool
/= CDouble
x then -CDouble
d_inf else CDouble
1e300

_di_above :: CDouble -> CDouble
_di_above CDouble
x =
  if CDouble
x CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
>= -CDouble
1e300 then
    CDouble
x CDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
+ (CDouble
1e-300 CDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
+ if CDouble
x CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
< CDouble
0 then -CDouble
x else CDouble
x) CDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
* CDouble
4.440892098500626e-16
  else
    if CDouble
x CDouble -> CDouble -> Bool
forall a. Eq a => a -> a -> Bool
/= CDouble
x then CDouble
d_inf else -CDouble
1e300
    
-- | /arb_get_di/ /x/ 
--
-- Returns the ball /x/ converted to a double-precision interval.
arb_get_di :: Ptr CArb -> IO CDi
arb_get_di :: Ptr CArb -> IO CDi
arb_get_di Ptr CArb
x = do
  (Arf
_, CDi
result) <- (Ptr CArf -> IO CDi) -> IO (Arf, CDi)
forall {a}. (Ptr CArf -> IO a) -> IO (Arf, a)
withNewArf ((Ptr CArf -> IO CDi) -> IO (Arf, CDi))
-> (Ptr CArf -> IO CDi) -> IO (Arf, CDi)
forall a b. (a -> b) -> a -> b
$ \Ptr CArf
t -> do
    Ptr CArf -> Ptr CArb -> CLong -> IO ()
arb_get_lbound_arf Ptr CArf
t Ptr CArb
x CLong
53
    CDouble
a <- Ptr CArf -> ArfRnd -> IO CDouble
arf_get_d Ptr CArf
t ArfRnd
arf_rnd_floor
    Ptr CArf -> Ptr CArb -> CLong -> IO ()
arb_get_ubound_arf Ptr CArf
t Ptr CArb
x CLong
53
    CDouble
b <- Ptr CArf -> ArfRnd -> IO CDouble
arf_get_d Ptr CArf
t ArfRnd
arf_rnd_ceil
    CDi -> IO CDi
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (CDi -> IO CDi) -> CDi -> IO CDi
forall a b. (a -> b) -> a -> b
$ CDouble -> CDouble -> CDi
CDi CDouble
a CDouble
b
  CDi -> IO CDi
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return CDi
result

-- | /arb_set_di/ /res/ /x/ /prec/ 
--
-- Sets the ball /res/ to the double-precision interval /x/, rounded to
-- /prec/ bits.
arb_set_di :: Ptr CArb -> CDi -> CLong -> IO ()
arb_set_di :: Ptr CArb -> CDi -> CLong -> IO ()
arb_set_di Ptr CArb
res (CDi CDouble
a CDouble
b) CLong
prec = do
  (Ptr CArf -> IO (Arf, ())) -> IO (Arf, (Arf, ()))
forall {a}. (Ptr CArf -> IO a) -> IO (Arf, a)
withNewArf ((Ptr CArf -> IO (Arf, ())) -> IO (Arf, (Arf, ())))
-> (Ptr CArf -> IO (Arf, ())) -> IO (Arf, (Arf, ()))
forall a b. (a -> b) -> a -> b
$ \Ptr CArf
t -> do
    (Ptr CArf -> IO ()) -> IO (Arf, ())
forall {a}. (Ptr CArf -> IO a) -> IO (Arf, a)
withNewArf ((Ptr CArf -> IO ()) -> IO (Arf, ()))
-> (Ptr CArf -> IO ()) -> IO (Arf, ())
forall a b. (a -> b) -> a -> b
$ \Ptr CArf
u -> do
      Ptr CArf -> CDouble -> IO ()
arf_set_d Ptr CArf
t CDouble
a
      Ptr CArf -> CDouble -> IO ()
arf_set_d Ptr CArf
u CDouble
b
      Ptr CArb -> Ptr CArf -> Ptr CArf -> CLong -> IO ()
arb_set_interval_arf Ptr CArb
res Ptr CArf
t Ptr CArf
u CLong
prec
  () -> IO ()
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return ()

-- | /di_print/ /x/ 
--
-- Prints /x/ to standard output. This simply prints decimal
-- representations of the floating-point endpoints; the decimals are not
-- guaranteed to be rounded outward.
di_print :: CDi -> IO ()
di_print :: CDi -> IO ()
di_print (CDi CDouble
a CDouble
b) = do
  String -> IO ()
putStr (String -> IO ()) -> String -> IO ()
forall a b. (a -> b) -> a -> b
$ String -> Double -> Double -> String
forall r. PrintfType r => String -> r
printf String
"[%.17g, %.17g]" (CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac CDouble
a :: Double)
                                   (CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac CDouble
b :: Double)
  
-- | /d_randtest2/ /state/ 
--
-- Returns a random non-NaN @double@ with any exponent. The value can be
-- infinite or subnormal.
di_randtest2 :: Ptr CFRandState -> IO CDouble
di_randtest2 :: Ptr CFRandState -> IO CDouble
di_randtest2 Ptr CFRandState
state = do
  CDouble
x <- Ptr CFRandState -> IO CDouble
d_randtest Ptr CFRandState
state
  CDouble -> IO CDouble
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return CDouble
x
  
-- | /di_randtest/ /state/ 
--
-- Returns an interval with random endpoints.
di_randtest :: Ptr CFRandState -> IO CDi
di_randtest :: Ptr CFRandState -> IO CDi
di_randtest Ptr CFRandState
state = do
  CDouble
a <- Ptr CFRandState -> IO CDouble
d_randtest Ptr CFRandState
state
  CDouble
b <- Ptr CFRandState -> IO CDouble
d_randtest Ptr CFRandState
state
  CDi -> IO CDi
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (CDi -> IO CDi) -> CDi -> IO CDi
forall a b. (a -> b) -> a -> b
$ if CDouble
a CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
> CDouble
b then CDouble -> CDouble -> CDi
CDi CDouble
b CDouble
a else CDouble -> CDouble -> CDi
CDi CDouble
a CDouble
b
    

-- Arithmetic ------------------------------------------------------------------

-- | /di_neg/ /x/ 
--
-- Returns the exact negation of /x/.
di_neg :: CDi -> CDi
di_neg :: CDi -> CDi
di_neg (CDi CDouble
a CDouble
b) = CDouble -> CDouble -> CDi
CDi (-CDouble
b) CDouble
a

-- Fast arithmetic -------------------------------------------------------------

-- The following methods perform fast but sloppy interval arithmetic: we
-- manipulate the endpoints with default rounding and then add or subtract
-- generic perturbations regardless of whether the operations were exact.
-- It is currently assumed that the CPU rounding mode is to nearest.
--
-- | /di_fast_add/ /x/ /y/ 
di_fast_add :: CDi -> CDi -> CDi
di_fast_add :: CDi -> CDi -> CDi
di_fast_add (CDi CDouble
a CDouble
b) (CDi CDouble
a' CDouble
b') = CDouble -> CDouble -> CDi
CDi (CDouble -> CDouble
_di_below (CDouble
aCDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
+CDouble
a')) (CDouble -> CDouble
_di_above (CDouble
bCDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
+CDouble
b'))
  
-- | /di_fast_sub/ /x/ /y/ 
di_fast_sub :: CDi -> CDi -> CDi
di_fast_sub :: CDi -> CDi -> CDi
di_fast_sub (CDi CDouble
a CDouble
b) (CDi CDouble
a' CDouble
b') = CDouble -> CDouble -> CDi
CDi (CDouble -> CDouble
_di_below (CDouble
aCDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
-CDouble
b')) (CDouble -> CDouble
_di_above (CDouble
bCDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
-CDouble
a'))

-- | /di_fast_mul/ /x/ /y/ 
di_fast_mul :: CDi -> CDi -> CDi
di_fast_mul :: CDi -> CDi -> CDi
di_fast_mul (CDi CDouble
xa CDouble
xb) (CDi CDouble
ya CDouble
yb) = CDouble -> CDouble -> CDi
CDi (CDouble -> CDouble
_di_below CDouble
u) (CDouble -> CDouble
_di_above CDouble
v) where
  (CDouble
u, CDouble
v) 
    | CDouble
xa CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
> CDouble
0 Bool -> Bool -> Bool
&& CDouble
ya CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
> CDouble
0 = (CDouble
xaCDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
*CDouble
ya, CDouble
xbCDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
*CDouble
yb)
    | CDouble
xa CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
> CDouble
0 Bool -> Bool -> Bool
&& CDouble
yb CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
< CDouble
0 = (CDouble
xbCDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
*CDouble
ya, CDouble
xaCDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
*CDouble
yb)
    | CDouble
xb CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
< CDouble
0 Bool -> Bool -> Bool
&& CDouble
ya CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
> CDouble
0 = (CDouble
xaCDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
*CDouble
yb, CDouble
xbCDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
*CDouble
ya)
    | CDouble
xb CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
< CDouble
0 Bool -> Bool -> Bool
&& CDouble
yb CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
< CDouble
0 = (CDouble
xbCDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
*CDouble
yb, CDouble
xaCDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
*CDouble
ya)
    | CDouble
a CDouble -> CDouble -> Bool
forall a. Eq a => a -> a -> Bool
/= CDouble
a Bool -> Bool -> Bool
|| CDouble
b CDouble -> CDouble -> Bool
forall a. Eq a => a -> a -> Bool
/= CDouble
b Bool -> Bool -> Bool
|| CDouble
c CDouble -> CDouble -> Bool
forall a. Eq a => a -> a -> Bool
/= CDouble
c Bool -> Bool -> Bool
|| CDouble
d CDouble -> CDouble -> Bool
forall a. Eq a => a -> a -> Bool
/= CDouble
d = (-CDouble
d_inf, CDouble
d_inf)
    | Bool
otherwise = (CDouble -> CDouble -> CDouble
forall a. Ord a => a -> a -> a
min (CDouble -> CDouble -> CDouble
forall a. Ord a => a -> a -> a
min CDouble
a CDouble
b) (CDouble -> CDouble -> CDouble
forall a. Ord a => a -> a -> a
min CDouble
c CDouble
d), CDouble -> CDouble -> CDouble
forall a. Ord a => a -> a -> a
max (CDouble -> CDouble -> CDouble
forall a. Ord a => a -> a -> a
max CDouble
a CDouble
b) (CDouble -> CDouble -> CDouble
forall a. Ord a => a -> a -> a
max CDouble
c CDouble
d))
    where
      a :: CDouble
a = CDouble
xa CDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
* CDouble
ya
      b :: CDouble
b = CDouble
xa CDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
* CDouble
yb
      c :: CDouble
c = CDouble
xb CDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
* CDouble
ya
      d :: CDouble
d = CDouble
xb CDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
* CDouble
yb
  
-- | /di_fast_div/ /x/ /y/ 
--
-- Returns the sum, difference, product or quotient of /x/ and /y/.
-- Division by zero is currently defined to return \([-\infty, +\infty]\).
di_fast_div :: CDi -> CDi -> CDi
di_fast_div :: CDi -> CDi -> CDi
di_fast_div (CDi CDouble
xa CDouble
xb) (CDi CDouble
ya CDouble
yb) = CDouble -> CDouble -> CDi
CDi (CDouble -> CDouble
_di_below CDouble
u) (CDouble -> CDouble
_di_above CDouble
v) where
  (CDouble
u, CDouble
v)
    | CDouble
ya CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
> CDouble
0 Bool -> Bool -> Bool
&& CDouble
xa CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
>= CDouble
0 = (CDouble
xaCDouble -> CDouble -> CDouble
forall a. Fractional a => a -> a -> a
/CDouble
yb, CDouble
xbCDouble -> CDouble -> CDouble
forall a. Fractional a => a -> a -> a
/CDouble
ya)
    | CDouble
ya CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
> CDouble
0 Bool -> Bool -> Bool
&& CDouble
xb CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
<= CDouble
0 = (CDouble
xaCDouble -> CDouble -> CDouble
forall a. Fractional a => a -> a -> a
/CDouble
ya, CDouble
xbCDouble -> CDouble -> CDouble
forall a. Fractional a => a -> a -> a
/CDouble
yb)
    | CDouble
ya CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
> CDouble
0            = (CDouble
xaCDouble -> CDouble -> CDouble
forall a. Fractional a => a -> a -> a
/CDouble
ya, CDouble
xbCDouble -> CDouble -> CDouble
forall a. Fractional a => a -> a -> a
/CDouble
ya)
    | CDouble
yb CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
< CDouble
0 Bool -> Bool -> Bool
&& CDouble
xa CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
>= CDouble
0 = (CDouble
xbCDouble -> CDouble -> CDouble
forall a. Fractional a => a -> a -> a
/CDouble
yb, CDouble
xaCDouble -> CDouble -> CDouble
forall a. Fractional a => a -> a -> a
/CDouble
ya)
    | CDouble
yb CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
< CDouble
0 Bool -> Bool -> Bool
&& CDouble
xb CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
<= CDouble
0 = (CDouble
xbCDouble -> CDouble -> CDouble
forall a. Fractional a => a -> a -> a
/CDouble
ya, CDouble
xaCDouble -> CDouble -> CDouble
forall a. Fractional a => a -> a -> a
/CDouble
yb)
    | CDouble
yb CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
<CDouble
0             = (CDouble
xbCDouble -> CDouble -> CDouble
forall a. Fractional a => a -> a -> a
/CDouble
yb, CDouble
xaCDouble -> CDouble -> CDouble
forall a. Fractional a => a -> a -> a
/CDouble
yb)
    | Bool
otherwise = (-CDouble
d_inf, CDouble
d_inf)

-- | /di_fast_sqr/ /x/ 
--
-- Returns the square of /x/. The output is clamped to be nonnegative.
di_fast_sqr ::  CDi -> CDi
di_fast_sqr :: CDi -> CDi
di_fast_sqr (CDi CDouble
a CDouble
b) =
  CDouble -> CDouble -> CDi
CDi (if CDouble
a CDouble -> CDouble -> Bool
forall a. Eq a => a -> a -> Bool
/= CDouble
0 then CDouble -> CDouble
_di_below CDouble
u else CDouble
u) (CDouble -> CDouble
_di_above CDouble
b) where
  (CDouble
u, CDouble
v)
    | CDouble
a CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
>= CDouble
0 = (CDouble
aCDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
*CDouble
a, CDouble
bCDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
*CDouble
b)
    | CDouble
b CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
<= CDouble
0 = (CDouble
bCDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
*CDouble
b, CDouble
aCDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
*CDouble
a)
    | Bool
otherwise = (CDouble
0, CDouble -> CDouble -> CDouble
forall a. Ord a => a -> a -> a
max (CDouble
aCDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
*CDouble
a) (CDouble
bCDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
*CDouble
b))

-- | /di_fast_add_d/ /x/ /y/ 
di_fast_add_d :: CDi -> CDouble -> CDi
di_fast_add_d :: CDi -> CDouble -> CDi
di_fast_add_d CDi
x CDouble
y = CDi -> CDi -> CDi
di_fast_add CDi
x (CDouble -> CDouble -> CDi
di_interval CDouble
y CDouble
y)
-- -- | /di_fast_sub_d/ /x/ /y/ 
di_fast_sub_d :: CDi -> CDouble -> CDi
di_fast_sub_d :: CDi -> CDouble -> CDi
di_fast_sub_d CDi
x CDouble
y = CDi -> CDi -> CDi
di_fast_sub CDi
x (CDouble -> CDouble -> CDi
di_interval CDouble
y CDouble
y)
-- | /di_fast_mul_d/ /x/ /y/
di_fast_mul_d :: CDi -> CDouble -> CDi
di_fast_mul_d :: CDi -> CDouble -> CDi
di_fast_mul_d CDi
x CDouble
y = CDi -> CDi -> CDi
di_fast_mul CDi
x (CDouble -> CDouble -> CDi
di_interval CDouble
y CDouble
y)
-- | /di_fast_div_d/ /x/ /y/
-- Arithmetic with an exact @double@ operand.
di_fast_div_d :: CDi -> CDouble -> CDi
di_fast_div_d :: CDi -> CDouble -> CDi
di_fast_div_d CDi
x CDouble
y = CDi -> CDi -> CDi
di_fast_div CDi
x (CDouble -> CDouble -> CDi
di_interval CDouble
y CDouble
y)

-- | /di_fast_log_nonnegative/ /x/ 
--
-- Returns an enclosure of \(\log(x)\). The lower endpoint of /x/ is
-- rounded up to 0 if it is negative.
di_fast_log_nonnegative :: CDi -> CDi
di_fast_log_nonnegative :: CDi -> CDi
di_fast_log_nonnegative (CDi CDouble
a CDouble
b) = CDouble -> CDouble -> CDi
CDi CDouble
a' CDouble
b' where
  a' :: CDouble
a' = if CDouble
a CDouble -> CDouble -> Bool
forall a. Ord a => a -> a -> Bool
<= CDouble
0 then (-CDouble
d_inf) else CDouble -> CDouble
mag_d_log_lower_bound CDouble
a
  b' :: CDouble
b' = CDouble -> CDouble
mag_d_log_upper_bound CDouble
b

-- | /di_fast_mid/ /x/ 
--
-- Returns an enclosure of the midpoint of /x/.
di_fast_mid :: CDi -> CDi
di_fast_mid :: CDi -> CDi
di_fast_mid (CDi CDouble
a CDouble
b)
  | CDouble
a CDouble -> CDouble -> Bool
forall a. Eq a => a -> a -> Bool
== -CDouble
d_inf Bool -> Bool -> Bool
|| CDouble
b CDouble -> CDouble -> Bool
forall a. Eq a => a -> a -> Bool
== CDouble
d_inf = CDouble -> CDouble -> CDi
di_interval (-CDouble
d_inf) CDouble
d_inf
  | Bool
otherwise = CDi -> CDouble -> CDi
di_fast_mul_d (CDi -> CDi -> CDi
di_fast_add (CDouble -> CDouble -> CDi
di_interval CDouble
a CDouble
a)
                                           (CDouble -> CDouble -> CDi
di_interval CDouble
b CDouble
b)) CDouble
0.5
                                           
-- | /di_fast_ubound_radius/ /x/ 
--
-- Returns an upper bound for the radius of /x/.
di_fast_ubound_radius :: CDi -> CDouble
di_fast_ubound_radius :: CDi -> CDouble
di_fast_ubound_radius (CDi CDouble
a CDouble
b) = CDouble -> CDouble
_di_above (CDouble
0.5 CDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
* (CDouble
b CDouble -> CDouble -> CDouble
forall a. Num a => a -> a -> a
-CDouble
a))