{-# LINE 1 "src/Data/Number/Flint/Support/D/Interval/FFI.hsc" #-}
module Data.Number.Flint.Support.D.Interval.FFI (
Di (..)
, CDi (..)
, di_interval
, arb_get_di
, arb_set_di
, di_print
, di_randtest2
, di_randtest
, di_neg
, 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
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
1forall a. Fractional a => a -> a -> a
/CDouble
0 :: CDouble
data Di = Di {-# UNPACK #-} !(ForeignPtr CDi)
data CDi = CDi CDouble CDouble deriving Int -> CDi -> ShowS
[CDi] -> ShowS
CDi -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [CDi] -> ShowS
$cshowList :: [CDi] -> ShowS
show :: CDi -> String
$cshow :: CDi -> String
showsPrec :: Int -> CDi -> ShowS
$cshowsPrec :: Int -> 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" #-}
di_interval :: CDouble -> CDouble -> CDi
di_interval :: CDouble -> CDouble -> CDi
di_interval CDouble
a CDouble
b =
if CDouble
a forall a. Ord a => a -> a -> Bool
<= CDouble
b
then CDouble -> CDouble -> CDi
CDi CDouble
a CDouble
b
else forall a. HasCallStack => String -> a
error forall a b. (a -> b) -> a -> b
$ forall r. PrintfType r => String -> r
printf String
"di_interval endpoints %g, %g not ordered.\n"
(forall a b. (Real a, Fractional b) => a -> b
realToFrac CDouble
a :: 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 forall a. Ord a => a -> a -> Bool
<= CDouble
1e300 then
CDouble
x forall a. Num a => a -> a -> a
- (CDouble
1e-300 forall a. Num a => a -> a -> a
+ if CDouble
x forall a. Ord a => a -> a -> Bool
< CDouble
0 then -CDouble
x else CDouble
x) forall a. Num a => a -> a -> a
* CDouble
4.440892098500626e-16
else
if CDouble
x 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 forall a. Ord a => a -> a -> Bool
>= -CDouble
1e300 then
CDouble
x forall a. Num a => a -> a -> a
+ (CDouble
1e-300 forall a. Num a => a -> a -> a
+ if CDouble
x forall a. Ord a => a -> a -> Bool
< CDouble
0 then -CDouble
x else CDouble
x) forall a. Num a => a -> a -> a
* CDouble
4.440892098500626e-16
else
if CDouble
x forall a. Eq a => a -> a -> Bool
/= CDouble
x then CDouble
d_inf else -CDouble
1e300
arb_get_di :: Ptr CArb -> IO CDi
arb_get_di :: Ptr CArb -> IO CDi
arb_get_di Ptr CArb
x = do
(Arf
_, CDi
result) <- forall {a}. (Ptr CArf -> IO a) -> IO (Arf, a)
withNewArf 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
forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ CDouble -> CDouble -> CDi
CDi CDouble
a CDouble
b
forall (m :: * -> *) a. Monad m => a -> m a
return CDi
result
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
forall {a}. (Ptr CArf -> IO a) -> IO (Arf, a)
withNewArf forall a b. (a -> b) -> a -> b
$ \Ptr CArf
t -> do
forall {a}. (Ptr CArf -> IO a) -> IO (Arf, a)
withNewArf 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
forall (m :: * -> *) a. Monad m => a -> m a
return ()
di_print :: CDi -> IO ()
di_print :: CDi -> IO ()
di_print (CDi CDouble
a CDouble
b) = do
String -> IO ()
putStr forall a b. (a -> b) -> a -> b
$ forall r. PrintfType r => String -> r
printf String
"[%.17g, %.17g]" (forall a b. (Real a, Fractional b) => a -> b
realToFrac CDouble
a :: Double)
(forall a b. (Real a, Fractional b) => a -> b
realToFrac CDouble
b :: Double)
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
forall (m :: * -> *) a. Monad m => a -> m a
return CDouble
x
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
forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ if CDouble
a 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
di_neg :: CDi -> CDi
di_neg :: CDi -> CDi
di_neg (CDi CDouble
a CDouble
b) = CDouble -> CDouble -> CDi
CDi (-CDouble
b) CDouble
a
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
aforall a. Num a => a -> a -> a
+CDouble
a')) (CDouble -> CDouble
_di_above (CDouble
bforall a. Num a => a -> a -> a
+CDouble
b'))
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
aforall a. Num a => a -> a -> a
-CDouble
b')) (CDouble -> CDouble
_di_above (CDouble
bforall a. Num a => a -> a -> a
-CDouble
a'))
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 forall a. Ord a => a -> a -> Bool
> CDouble
0 Bool -> Bool -> Bool
&& CDouble
ya forall a. Ord a => a -> a -> Bool
> CDouble
0 = (CDouble
xaforall a. Num a => a -> a -> a
*CDouble
ya, CDouble
xbforall a. Num a => a -> a -> a
*CDouble
yb)
| CDouble
xa forall a. Ord a => a -> a -> Bool
> CDouble
0 Bool -> Bool -> Bool
&& CDouble
yb forall a. Ord a => a -> a -> Bool
< CDouble
0 = (CDouble
xbforall a. Num a => a -> a -> a
*CDouble
ya, CDouble
xaforall a. Num a => a -> a -> a
*CDouble
yb)
| CDouble
xb forall a. Ord a => a -> a -> Bool
< CDouble
0 Bool -> Bool -> Bool
&& CDouble
ya forall a. Ord a => a -> a -> Bool
> CDouble
0 = (CDouble
xaforall a. Num a => a -> a -> a
*CDouble
yb, CDouble
xbforall a. Num a => a -> a -> a
*CDouble
ya)
| CDouble
xb forall a. Ord a => a -> a -> Bool
< CDouble
0 Bool -> Bool -> Bool
&& CDouble
yb forall a. Ord a => a -> a -> Bool
< CDouble
0 = (CDouble
xbforall a. Num a => a -> a -> a
*CDouble
yb, CDouble
xaforall a. Num a => a -> a -> a
*CDouble
ya)
| CDouble
a forall a. Eq a => a -> a -> Bool
/= CDouble
a Bool -> Bool -> Bool
|| CDouble
b forall a. Eq a => a -> a -> Bool
/= CDouble
b Bool -> Bool -> Bool
|| CDouble
c forall a. Eq a => a -> a -> Bool
/= CDouble
c Bool -> Bool -> Bool
|| CDouble
d forall a. Eq a => a -> a -> Bool
/= CDouble
d = (-CDouble
d_inf, CDouble
d_inf)
| Bool
otherwise = (forall a. Ord a => a -> a -> a
min (forall a. Ord a => a -> a -> a
min CDouble
a CDouble
b) (forall a. Ord a => a -> a -> a
min CDouble
c CDouble
d), forall a. Ord a => a -> a -> a
max (forall a. Ord a => a -> a -> a
max CDouble
a CDouble
b) (forall a. Ord a => a -> a -> a
max CDouble
c CDouble
d))
where
a :: CDouble
a = CDouble
xa forall a. Num a => a -> a -> a
* CDouble
ya
b :: CDouble
b = CDouble
xa forall a. Num a => a -> a -> a
* CDouble
yb
c :: CDouble
c = CDouble
xb forall a. Num a => a -> a -> a
* CDouble
ya
d :: CDouble
d = CDouble
xb forall a. Num a => a -> a -> a
* CDouble
yb
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 forall a. Ord a => a -> a -> Bool
> CDouble
0 Bool -> Bool -> Bool
&& CDouble
xa forall a. Ord a => a -> a -> Bool
>= CDouble
0 = (CDouble
xaforall a. Fractional a => a -> a -> a
/CDouble
yb, CDouble
xbforall a. Fractional a => a -> a -> a
/CDouble
ya)
| CDouble
ya forall a. Ord a => a -> a -> Bool
> CDouble
0 Bool -> Bool -> Bool
&& CDouble
xb forall a. Ord a => a -> a -> Bool
<= CDouble
0 = (CDouble
xaforall a. Fractional a => a -> a -> a
/CDouble
ya, CDouble
xbforall a. Fractional a => a -> a -> a
/CDouble
yb)
| CDouble
ya forall a. Ord a => a -> a -> Bool
> CDouble
0 = (CDouble
xaforall a. Fractional a => a -> a -> a
/CDouble
ya, CDouble
xbforall a. Fractional a => a -> a -> a
/CDouble
ya)
| CDouble
yb forall a. Ord a => a -> a -> Bool
< CDouble
0 Bool -> Bool -> Bool
&& CDouble
xa forall a. Ord a => a -> a -> Bool
>= CDouble
0 = (CDouble
xbforall a. Fractional a => a -> a -> a
/CDouble
yb, CDouble
xaforall a. Fractional a => a -> a -> a
/CDouble
ya)
| CDouble
yb forall a. Ord a => a -> a -> Bool
< CDouble
0 Bool -> Bool -> Bool
&& CDouble
xb forall a. Ord a => a -> a -> Bool
<= CDouble
0 = (CDouble
xbforall a. Fractional a => a -> a -> a
/CDouble
ya, CDouble
xaforall a. Fractional a => a -> a -> a
/CDouble
yb)
| CDouble
yb forall a. Ord a => a -> a -> Bool
<CDouble
0 = (CDouble
xbforall a. Fractional a => a -> a -> a
/CDouble
yb, CDouble
xaforall a. Fractional a => a -> a -> a
/CDouble
yb)
| Bool
otherwise = (-CDouble
d_inf, CDouble
d_inf)
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 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 forall a. Ord a => a -> a -> Bool
>= CDouble
0 = (CDouble
aforall a. Num a => a -> a -> a
*CDouble
a, CDouble
bforall a. Num a => a -> a -> a
*CDouble
b)
| CDouble
b forall a. Ord a => a -> a -> Bool
<= CDouble
0 = (CDouble
bforall a. Num a => a -> a -> a
*CDouble
b, CDouble
aforall a. Num a => a -> a -> a
*CDouble
a)
| Bool
otherwise = (CDouble
0, forall a. Ord a => a -> a -> a
max (CDouble
aforall a. Num a => a -> a -> a
*CDouble
a) (CDouble
bforall a. Num a => a -> a -> a
*CDouble
b))
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 :: 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 :: 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 :: 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 :: 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 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 :: CDi -> CDi
di_fast_mid :: CDi -> CDi
di_fast_mid (CDi CDouble
a CDouble
b)
| CDouble
a forall a. Eq a => a -> a -> Bool
== -CDouble
d_inf Bool -> Bool -> Bool
|| CDouble
b 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 :: CDi -> CDouble
di_fast_ubound_radius :: CDi -> CDouble
di_fast_ubound_radius (CDi CDouble
a CDouble
b) = CDouble -> CDouble
_di_above (CDouble
0.5 forall a. Num a => a -> a -> a
* (CDouble
b forall a. Num a => a -> a -> a
-CDouble
a))