HERA-0.2: Data/Number/MPFR.hs
{-# LANGUAGE TypeSynonymInstances #-}
{-# INCLUDE <mpfr.h> #-}
{-# INCLUDE <hsmpfr.h> #-}
module Data.Number.MPFR (
-- | This module should always be imported qualified.
-- *** Naming
-- | - functions ending with _ return a pair (value, rounding indicator).
-- Rounding indicator indicates whether the result is rounded and in which
-- directon as described in the MPFR manual.
--
-- - the same functions without the _ return just the value.
--
-- - functions with added \"w\" correspond to MPFR _ui functions
--
-- - functions with added \"i\" correspond to MPFR _si functions
-- *** Equality testing
-- | Equality works as follows:
--
-- - NaN \/= Nan,
--
-- - Infinity = Infinity,
--
-- - \-Infinity = -Infinity
--
-- - otherwise normal comparison
-- *** Ordering
-- | Ordering works as follows:
--
-- - compare NaN _ = GT
--
-- - compare _ NaN = GT
--
-- - infinity < _ = False
--
-- - \-infinity > _ = False
--
-- - NaN [\<,\>,\>=,<=] _ = False
--
-- This mimics the behaviour of built in haskell Float and Double.
-- *** Num instance
-- | Operations defined in Num will be computed so that no precision is lost.
Dyadic,
Precision,
RoundMode(Near, Zero, Up, Down),
add, sub, mul, div, inverse,
add_, sub_, mul_, div_,
addw, addi, mulw, muli, divw, divi, wdiv, idiv, subw, subi, wsub, isub,
addw_, addi_, mulw_, muli_, divw_, divi_, wdiv_, idiv_, subw_, subi_, wsub_, isub_,
mul2w, mul2i, div2w, div2i, mul2w_, mul2i_, div2w_, div2i_,
int2i, int2w, int2i_, int2w_,
fma, fms, fma_, fms_, nextBelow,
sqr, sqrt, root, pow, poww, powi, wpoww, wpow,
sqr_, sqrt_, root_, pow_, poww_, powi_, wpoww_, wpow_,
exp, exp2, exp10, log, log2, log10, sinh, cosh, tanh,
exp_, exp2_, exp10_, log_, log2_, log10_, sinh_, cosh_, tanh_,
neg, absD, dim, neg_, absD_, dim_,
isNaN, isInfinite, isNumber, isZero, greater, greatereq, less, lesseq,
equal, maxD, minD, maxD_, minD_, sgn,
dyadicToDouble, dyadicToWord, dyadicToInt, dyadicToString, decompose, toStringExp, toString,
pi, log2c, euler, catalan, pi_, log2c_, euler_, catalan_,
set, set_,
fromDouble, fromInt, fromWord, fromDouble_, fromInt_, fromWord_, fromIntegerA, compose, fromString,
getPrec, getMantissa, getExp,
minPrec, one, zero, addPrec
) where
import Data.Number.FFIhelper
import Foreign.C(CInt, CLong, CULong, withCString, peekCString)
import Foreign.Marshal(alloca, peekArray)
import Foreign(unsafePerformIO, peek, Ptr, mallocForeignPtrBytes, with)
import Data.Bits(shiftL)
import Data.Word(Word)
import Prelude hiding (div, sqrt, read, isNaN, isInfinite, exp, log, sinh, cosh, tanh, pi)
type Dyadic = MPFR_T
type Precision = Word
-- these are helper functions, only for internal use
{-# INLINE withDyadicsBA #-}
withDyadicsBA :: RoundMode -> Precision -> Dyadic -> Dyadic
-> (Ptr MPFR_T -> Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt)
-> (Dyadic, Int)
withDyadicsBA r p mp1 mp2 f = unsafePerformIO go
where go = do ls <- mpfr_custom_get_size (fromIntegral p)
fp <- mallocForeignPtrBytes (fromIntegral ls)
let dummy = MP (fromIntegral p) 0 0 fp
with dummy $ \p1 -> do
with mp1 $ \p2 -> do
with mp2 $ \p3 -> do
r2 <- f p1 p2 p3 ((fromIntegral . fromEnum) r)
r1 <- peekP p1 fp
return (r1, fromIntegral r2)
{-# INLINE withDyadicBAui #-}
withDyadicBAui :: RoundMode -> Precision -> Dyadic -> CULong
-> (Ptr MPFR_T -> Ptr MPFR_T -> CULong -> CRoundMode -> IO CInt)
-> (Dyadic, Int)
withDyadicBAui r p mp1 d f = unsafePerformIO go
where go = do ls <- mpfr_custom_get_size (fromIntegral p)
fp <- mallocForeignPtrBytes (fromIntegral ls)
let dummy = MP (fromIntegral p) 0 0 fp
with dummy $ \p1 -> do
with mp1 $ \p2 -> do
r2 <- f p1 p2 d ((fromIntegral . fromEnum) r)
r1 <- peekP p1 fp
return (r1, fromIntegral r2)
{-# INLINE withDyadicBAsi #-}
withDyadicBAsi :: RoundMode -> Precision -> Dyadic -> CLong
-> (Ptr MPFR_T -> Ptr MPFR_T -> CLong -> CRoundMode -> IO CInt)
-> (Dyadic, Int)
withDyadicBAsi r p mp1 d f = unsafePerformIO go
where go = do ls <- mpfr_custom_get_size (fromIntegral p)
fp <- mallocForeignPtrBytes (fromIntegral ls)
let dummy = MP (fromIntegral p) 0 0 fp
with dummy $ \p1 -> do
with mp1 $ \p2 -> do
r2 <- f p1 p2 d ((fromIntegral . fromEnum) r)
r1 <- peekP p1 fp
return (r1, fromIntegral r2)
{-# INLINE withDyadicBAiu #-}
withDyadicBAiu :: RoundMode -> Precision -> CULong -> Dyadic
-> (Ptr MPFR_T -> CULong -> Ptr MPFR_T -> CRoundMode -> IO CInt)
-> (Dyadic, Int)
withDyadicBAiu r p d mp1 f = unsafePerformIO go
where go = do ls <- mpfr_custom_get_size (fromIntegral p)
fp <- mallocForeignPtrBytes (fromIntegral ls)
let dummy = MP (fromIntegral p) 0 0 fp
with dummy $ \p1 -> do
with mp1 $ \p2 -> do
r2 <- f p1 d p2 ((fromIntegral . fromEnum) r)
r1 <- peekP p1 fp
return (r1, fromIntegral r2)
{-# INLINE withDyadicBAis #-}
withDyadicBAis :: RoundMode -> Precision -> CLong -> Dyadic
-> (Ptr MPFR_T -> CLong -> Ptr MPFR_T -> CRoundMode -> IO CInt)
-> (Dyadic, Int)
withDyadicBAis r p d mp1 f = unsafePerformIO go
where go = do ls <- mpfr_custom_get_size (fromIntegral p)
fp <- mallocForeignPtrBytes (fromIntegral ls)
let dummy = MP (fromIntegral p) 0 0 fp
with dummy $ \p1 -> do
with mp1 $ \p2 -> do
r2 <- f p1 d p2 ((fromIntegral . fromEnum) r)
r1 <- peekP p1 fp
return (r1, fromIntegral r2)
{-# INLINE withDyadicB #-}
withDyadicB :: Dyadic -> (Ptr MPFR_T -> IO CInt) -> CInt
withDyadicB mp1 f = unsafePerformIO go
where go = with mp1 $ \p1 -> f p1
withDyadicP :: Dyadic -> (Ptr MPFR_T -> IO CPrecision) -> CPrecision
withDyadicP mp1 f = unsafePerformIO go
where go = with mp1 $ \p1 -> f p1
{-# INLINE withDyadic #-}
withDyadic :: RoundMode -> Precision -> Dyadic
-> (Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt)
-> (Dyadic, Int)
withDyadic r p mp1 f = unsafePerformIO go
where go = do ls <- mpfr_custom_get_size (fromIntegral p)
fp <- mallocForeignPtrBytes (fromIntegral ls)
let dummy = MP (fromIntegral p) 0 0 fp
with dummy $ \p1 -> do
with mp1 $ \p2 -> do
r2 <- f p1 p2 ((fromIntegral . fromEnum) r)
r1 <- peekP p1 fp
return (r1, fromIntegral r2)
{-# INLINE withDyadicBB #-}
withDyadicBB :: Dyadic -> Dyadic
-> (Ptr MPFR_T -> Ptr MPFR_T -> IO CInt)
-> CInt
withDyadicBB mp1 mp2 f = unsafePerformIO go
where go = do with mp1 $ \p1 -> do
with mp2 $ \p2 -> do
f p1 p2
{-# INLINE withDyadicC #-}
withDyadicC :: RoundMode -> Precision ->
(Ptr MPFR_T -> CRoundMode -> IO CInt) -> (Dyadic, Int)
withDyadicC r p f = unsafePerformIO go
where go = do ls <- mpfr_custom_get_size (fromIntegral p)
fp <- mallocForeignPtrBytes (fromIntegral ls)
let dummy = MP (fromIntegral p) 0 0 fp
with dummy $ \p1 -> do
r2 <- f p1 ((fromIntegral . fromEnum) r)
r1 <- peekP p1 fp
return (r1, fromIntegral r2)
checkPrec :: Precision -> Precision
checkPrec = max minPrec
stringToDyadic :: RoundMode -> Precision -> Word -> String -> Dyadic
stringToDyadic r p b d = unsafePerformIO go
where go = do ls <- mpfr_custom_get_size (fromIntegral p)
fp <- mallocForeignPtrBytes (fromIntegral ls)
let dummy = MP (fromIntegral p) 0 0 fp
with dummy $ \p1 -> do
withCString d $ \p2 -> do
_ <- mpfr_set_str p1 p2 (fromIntegral b) ((fromIntegral . fromEnum) r)
peekP p1 fp
getMantissa' :: Dyadic -> [Limb]
getMantissa' mp1 = unsafePerformIO go
where go = do with mp1 $ \p1 -> do
pt <- mpfr_custom_get_mantissa p1
arr <- peekArray (ceiling ((fromIntegral p ::Double) / fromIntegral bitsPerMPLimb)) pt ;
return arr
p = getPrec mp1
{- TODO: this is inefficient
binprec :: Integer -> Precision
binprec i = length (takeWhile (/= 0) (iterate (flip shiftR 1) i)
-}
binprec :: Integer -> Precision
binprec d = floor (logBase 2 (fromInteger (if d >= 0 then d else -d)) :: Double) + 1
--------------------------------------------------------------------
-- pure wrappers for basic arithmetic operations
add :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic
add r p d1 d2 = fst $ add_ r p d1 d2
sub :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic
sub r p d1 d2 = fst $ sub_ r p d1 d2
mul :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic
mul r p d1 d2 = fst $ mul_ r p d1 d2
div :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic
div r p d1 d2 = fst $ div_ r p d1 d2
add_ :: RoundMode -> Precision -> Dyadic -> Dyadic -> (Dyadic,Int)
add_ r p d1 d2 = withDyadicsBA r p d1 d2 mpfr_add
sub_ :: RoundMode -> Precision -> Dyadic -> Dyadic -> (Dyadic,Int)
sub_ r p d1 d2 = withDyadicsBA r p d1 d2 mpfr_sub
mul_ :: RoundMode -> Precision -> Dyadic -> Dyadic -> (Dyadic,Int)
mul_ r p d1 d2 = withDyadicsBA r p d1 d2 mpfr_mul
div_ :: RoundMode -> Precision -> Dyadic -> Dyadic -> (Dyadic, Int)
div_ r p d1 d2 = withDyadicsBA r p d1 d2 mpfr_div
inverse :: Dyadic -> Dyadic
inverse d = div Near (getPrec d) one d
----------------------------------------------------------------
-- basic arithmetic operations with mixed operands
addw :: RoundMode -> Precision -> Dyadic -> Word -> Dyadic
addw r p d1 d = fst $ addw_ r p d1 d
addi :: RoundMode -> Precision -> Dyadic -> Int -> Dyadic
addi r p d1 d = fst $ addi_ r p d1 d
mulw :: RoundMode -> Precision -> Dyadic -> Word -> Dyadic
mulw r p d1 d = fst $ mulw_ r p d1 d
muli :: RoundMode -> Precision -> Dyadic -> Int -> Dyadic
muli r p d1 d = fst $ muli_ r p d1 d
divw :: RoundMode -> Precision -> Dyadic -> Word -> Dyadic
divw r p d1 d = fst $ divw_ r p d1 d
divi :: RoundMode -> Precision -> Dyadic -> Int -> Dyadic
divi r p d1 d = fst $ divi_ r p d1 d
wdiv :: RoundMode -> Precision -> Word -> Dyadic -> Dyadic
wdiv r p d d1 = fst $ wdiv_ r p d d1
idiv :: RoundMode -> Precision -> Int -> Dyadic -> Dyadic
idiv r p d d1 = fst $ idiv_ r p d d1
subw :: RoundMode -> Precision -> Dyadic -> Word -> Dyadic
subw r p d1 d = fst $ subw_ r p d1 d
subi :: RoundMode -> Precision -> Dyadic -> Int -> Dyadic
subi r p d1 d = fst $ subi_ r p d1 d
wsub :: RoundMode -> Precision -> Word -> Dyadic -> Dyadic
wsub r p d d1 = fst $ wsub_ r p d d1
isub :: RoundMode -> Precision -> Int -> Dyadic -> Dyadic
isub r p d d1 = fst $ isub_ r p d d1
addw_ :: RoundMode -> Precision -> Dyadic -> Word -> (Dyadic, Int)
addw_ r p d1 d = withDyadicBAui r p d1 (fromIntegral d) mpfr_add_ui
addi_ :: RoundMode -> Precision -> Dyadic -> Int -> (Dyadic, Int)
addi_ r p d1 d = withDyadicBAsi r p d1 (fromIntegral d) mpfr_add_si
mulw_ :: RoundMode -> Precision -> Dyadic -> Word -> (Dyadic, Int)
mulw_ r p d1 d = withDyadicBAui r p d1 (fromIntegral d) mpfr_mul_ui
muli_ :: RoundMode -> Precision -> Dyadic -> Int -> (Dyadic, Int)
muli_ r p d1 d = withDyadicBAsi r p d1 (fromIntegral d) mpfr_mul_si
divw_ :: RoundMode -> Precision -> Dyadic -> Word -> (Dyadic, Int)
divw_ r p d1 d = withDyadicBAui r p d1 (fromIntegral d) mpfr_div_ui
divi_ :: RoundMode -> Precision -> Dyadic -> Int -> (Dyadic, Int)
divi_ r p d1 d = withDyadicBAsi r p d1 (fromIntegral d) mpfr_div_si
wdiv_ :: RoundMode -> Precision -> Word -> Dyadic -> (Dyadic, Int)
wdiv_ r p d d1 = withDyadicBAiu r p (fromIntegral d) d1 mpfr_ui_div
idiv_ :: RoundMode -> Precision -> Int -> Dyadic -> (Dyadic, Int)
idiv_ r p d d1 = withDyadicBAis r p (fromIntegral d) d1 mpfr_si_div
subw_ :: RoundMode -> Precision -> Dyadic -> Word -> (Dyadic, Int)
subw_ r p d1 d = withDyadicBAui r p d1 (fromIntegral d) mpfr_sub_ui
subi_ :: RoundMode -> Precision -> Dyadic -> Int -> (Dyadic, Int)
subi_ r p d1 d = withDyadicBAsi r p d1 (fromIntegral d) mpfr_sub_si
wsub_ :: RoundMode -> Precision -> Word -> Dyadic -> (Dyadic, Int)
wsub_ r p d d1 = withDyadicBAiu r p (fromIntegral d) d1 mpfr_ui_sub
isub_ :: RoundMode -> Precision -> Int -> Dyadic -> (Dyadic, Int)
isub_ r p d d1 = withDyadicBAis r p (fromIntegral d) d1 mpfr_si_sub
----------------------------------------------------------
-- multiplication and division with 2 ^ x
mul2w :: RoundMode -> Precision -> Dyadic -> Word -> Dyadic
mul2w r p d1 d2 = fst $ mul2w_ r p d1 d2
mul2i :: RoundMode -> Precision -> Dyadic -> Int -> Dyadic
mul2i r p d1 d2 = fst $ mul2i_ r p d1 d2
div2w :: RoundMode -> Precision -> Dyadic -> Word -> Dyadic
div2w r p d1 d2 = fst $ div2w_ r p d1 d2
div2i :: RoundMode -> Precision -> Dyadic -> Int -> Dyadic
div2i r p d1 d2 = fst $ div2i_ r p d1 d2
mul2w_ :: RoundMode -> Precision -> Dyadic -> Word -> (Dyadic, Int)
mul2w_ r p d1 d2 = withDyadicBAui r p d1 (fromIntegral d2) mpfr_mul_2ui
mul2i_ :: RoundMode -> Precision -> Dyadic -> Int -> (Dyadic, Int)
mul2i_ r p d1 d2 = withDyadicBAsi r p d1 (fromIntegral d2) mpfr_mul_2si
div2w_ :: RoundMode -> Precision -> Dyadic -> Word -> (Dyadic, Int)
div2w_ r p d1 d2 = withDyadicBAui r p d1 (fromIntegral d2) mpfr_div_2ui
div2i_ :: RoundMode -> Precision -> Dyadic -> Int -> (Dyadic, Int)
div2i_ r p d1 d2 = withDyadicBAsi r p d1 (fromIntegral d2) mpfr_div_2si
----------------------------------------------------------
-- x * 2 ^ y
int2i :: RoundMode -> Precision -> Int -> Int -> Dyadic
int2i r p i e = fst $ int2i_ r p i e
int2w :: RoundMode -> Precision -> Word -> Int -> Dyadic
int2w r p i e = fst $ int2w_ r p i e
int2i_ :: RoundMode -> Precision -> Int -> Int -> (Dyadic, Int)
int2i_ r p i e = unsafePerformIO go
where go = do ls <- mpfr_custom_get_size (fromIntegral p)
fp <- mallocForeignPtrBytes (fromIntegral ls)
let dummy = MP (fromIntegral p) 0 0 fp
with dummy $ \p1 -> do
r2 <- mpfr_set_si_2exp p1 (fromIntegral i) (fromIntegral e) ((fromIntegral . fromEnum) r)
r1 <- peekP p1 fp
return (r1, fromIntegral r2)
int2w_ :: RoundMode -> Precision -> Word -> Int -> (Dyadic, Int)
int2w_ r p i e = unsafePerformIO go
where go = do ls <- mpfr_custom_get_size (fromIntegral p)
fp <- mallocForeignPtrBytes (fromIntegral ls)
let dummy = MP (fromIntegral p) 0 0 fp
with dummy $ \p1 -> do
r2 <- mpfr_set_ui_2exp p1 (fromIntegral i) (fromIntegral e) ((fromIntegral . fromEnum) r)
r1 <- peekP p1 fp
return (r1, fromIntegral r2)
----------------------------------------------------------
-- Return d1 * d2 + d3
fma :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic -> Dyadic
fma r p d1 d2 d3 = fst $ fma_ r p d1 d2 d3
-- Return d1 * d2 - d3
fms :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic -> Dyadic
fms r p d1 d2 d3 = fst $ fms_ r p d1 d2 d3
fma_ :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic -> (Dyadic, Int)
fma_ r p mp1 mp2 mp3 = unsafePerformIO go
where go = do ls <- mpfr_custom_get_size (fromIntegral p)
fp <- mallocForeignPtrBytes (fromIntegral ls)
let dummy = MP (fromIntegral p) 0 0 fp
with dummy $ \p1 -> do
with mp1 $ \p2 -> do
with mp2 $ \p3 -> do
with mp3 $ \p4 -> do
r2 <- mpfr_fma p1 p2 p3 p4 ((fromIntegral . fromEnum) r)
r1 <- peekP p1 fp
return (r1, fromIntegral r2)
fms_ :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic -> (Dyadic, Int)
fms_ r p mp1 mp2 mp3 = unsafePerformIO go
where go = do ls <- mpfr_custom_get_size (fromIntegral p)
fp <- mallocForeignPtrBytes (fromIntegral ls)
let dummy = MP (fromIntegral p) 0 0 fp
with dummy $ \p1 -> do
with mp1 $ \p2 -> do
with mp2 $ \p3 -> do
with mp3 $ \p4 -> do
r2 <- mpfr_fms p1 p2 p3 p4 ((fromIntegral . fromEnum) r)
r1 <- peekP p1 fp
return (r1, fromIntegral r2)
nextBelow :: Dyadic -> Dyadic
nextBelow mp1 = unsafePerformIO go
where go = do let p = fromIntegral (getPrec mp1)
ls <- mpfr_custom_get_size p
fp <- mallocForeignPtrBytes (fromIntegral ls)
let dummy = MP p 0 0 fp
with dummy $ \p1 -> do
with mp1 $ \p2 -> do
_ <- mpfr_set p1 p2 ((fromIntegral . fromEnum) Near)
mpfr_nextbelow p1
peekP p1 fp
----------------------------------------------------------
-- powers
sqr :: RoundMode -> Precision -> Dyadic -> Dyadic
sqr r p d = fst $ sqr_ r p d
sqrt :: RoundMode -> Precision -> Dyadic -> Dyadic
sqrt r p d = fst $ sqrt_ r p d
root :: RoundMode -> Precision -> Dyadic -> Word -> Dyadic
root r p d n = fst $ root_ r p d n
pow :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic
pow r p d1 d2 = fst $ pow_ r p d1 d2
poww :: RoundMode -> Precision -> Dyadic -> Word -> Dyadic
poww r p d1 d2 = fst $ poww_ r p d1 d2
powi :: RoundMode -> Precision -> Dyadic -> Int -> Dyadic
powi r p d1 d2 = fst $ powi_ r p d1 d2
wpoww :: RoundMode -> Precision -> Word -> Word -> Dyadic
wpoww r p d1 d2 = fst $ wpoww_ r p d1 d2
wpow :: RoundMode -> Precision -> Word -> Dyadic -> Dyadic
wpow r p d1 d2 = fst $ wpow_ r p d1 d2
sqr_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)
sqr_ r p d = withDyadic r p d mpfr_sqr
sqrt_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)
sqrt_ r p d = withDyadic r p d mpfr_sqrt
root_ :: RoundMode -> Precision -> Dyadic -> Word -> (Dyadic, Int)
root_ r p d n = withDyadicBAui r p d (fromIntegral n) mpfr_root
pow_ :: RoundMode -> Precision -> Dyadic -> Dyadic -> (Dyadic, Int)
pow_ r p d1 d2 = withDyadicsBA r p d1 d2 mpfr_pow
poww_ :: RoundMode -> Precision -> Dyadic -> Word -> (Dyadic , Int)
poww_ r p d1 d2 = withDyadicBAui r p d1 (fromIntegral d2) mpfr_pow_ui
powi_ :: RoundMode -> Precision -> Dyadic -> Int -> (Dyadic , Int)
powi_ r p d1 d2 = withDyadicBAsi r p d1 (fromIntegral d2) mpfr_pow_si
wpoww_ :: RoundMode -> Precision -> Word -> Word -> (Dyadic , Int)
wpoww_ r p d1 d2 = unsafePerformIO go
where go = do ls <- mpfr_custom_get_size (fromIntegral p)
fp <- mallocForeignPtrBytes (fromIntegral ls)
let dummy = MP (fromIntegral p) 0 0 fp
with dummy $ \p1 -> do
r2 <- mpfr_ui_pow_ui p1 (fromIntegral d1) (fromIntegral d2) ((fromIntegral . fromEnum) r)
r1 <- peekP p1 fp
return (r1, fromIntegral r2)
wpow_ :: RoundMode -> Precision -> Word -> Dyadic -> (Dyadic , Int)
wpow_ r p d1 d2 = withDyadicBAiu r p (fromIntegral d1) d2 mpfr_ui_pow
-----------------------------------------------------------
-- transcendental functions
exp :: RoundMode -> Precision -> Dyadic -> Dyadic
exp r p d = fst $ exp_ r p d
exp2 :: RoundMode -> Precision -> Dyadic -> Dyadic
exp2 r p d = fst $ exp2_ r p d
exp10 :: RoundMode -> Precision -> Dyadic -> Dyadic
exp10 r p d = fst $ exp10_ r p d
log :: RoundMode -> Precision -> Dyadic -> Dyadic
log r p d = fst $ log_ r p d
log2 :: RoundMode -> Precision -> Dyadic -> Dyadic
log2 r p d = fst $ log2_ r p d
log10 :: RoundMode -> Precision -> Dyadic -> Dyadic
log10 r p d = fst $ log10_ r p d
sinh :: RoundMode -> Precision -> Dyadic -> Dyadic
sinh r p d = fst $ sinh_ r p d
cosh :: RoundMode -> Precision -> Dyadic -> Dyadic
cosh r p d = fst $ cosh_ r p d
tanh :: RoundMode -> Precision -> Dyadic -> Dyadic
tanh r p d = fst $ tanh_ r p d
exp_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)
exp_ r p d = withDyadic r p d mpfr_exp
exp2_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)
exp2_ r p d = withDyadic r p d mpfr_exp2
exp10_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)
exp10_ r p d = withDyadic r p d mpfr_exp10
log_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)
log_ r p d = withDyadic r p d mpfr_log
log2_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)
log2_ r p d = withDyadic r p d mpfr_log2
log10_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)
log10_ r p d = withDyadic r p d mpfr_log10
sinh_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)
sinh_ r p d = withDyadic r p d mpfr_sinh
cosh_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)
cosh_ r p d = withDyadic r p d mpfr_cosh
tanh_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)
tanh_ r p d = withDyadic r p d mpfr_tanh
------------------------------------------------------------
neg :: RoundMode -> Precision -> Dyadic -> Dyadic
neg r p d = fst $ neg_ r p d
absD :: RoundMode -> Precision -> Dyadic -> Dyadic
absD r p d = fst $ absD_ r p d
dim :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic
dim r p d1 d2 = fst $ dim_ r p d1 d2
neg_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)
neg_ r p d = withDyadic r p d mpfr_neg
absD_ :: RoundMode -> Precision -> Dyadic -> (Dyadic , Int)
absD_ r p d = withDyadic r p d mpfr_abs
dim_ :: RoundMode -> Precision -> Dyadic -> Dyadic -> (Dyadic, Int)
dim_ r p d1 d2 = withDyadicsBA r p d1 d2 mpfr_dim
--------------------------------------------------------
-- comparison functions
isNaN :: Dyadic -> Bool
isNaN d = withDyadicB d mpfr_nan_p /= 0
isInfinite :: Dyadic -> Bool
isInfinite d = withDyadicB d mpfr_inf_p /= 0
isNumber :: Dyadic -> Bool
isNumber d = withDyadicB d mpfr_number_p /= 0
isZero :: Dyadic -> Bool
isZero d = withDyadicB d mpfr_zero_p /= 0
greater :: Dyadic -> Dyadic -> Bool
greater d1 d2 = withDyadicBB d1 d2 mpfr_greater_p /= 0
greatereq :: Dyadic -> Dyadic -> Bool
greatereq d1 d2 = withDyadicBB d1 d2 mpfr_greaterequal_p /= 0
less :: Dyadic -> Dyadic -> Bool
less d1 d2 = withDyadicBB d1 d2 mpfr_less_p /= 0
lesseq :: Dyadic -> Dyadic -> Bool
lesseq d1 d2 = withDyadicBB d1 d2 mpfr_lessequal_p /= 0
{- sets erange flag
lessgreater :: Dyadic -> Dyadic -> Bool
lessgreater d1 d2 = withDyadicBB d1 d2 mpfr_lessgreater_p /= 0
-}
equal :: Dyadic -> Dyadic -> Bool
equal d1 d2 = withDyadicBB d1 d2 mpfr_equal_p /= 0
maxD :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic
maxD r p d1 d2 = fst $ maxD_ r p d1 d2
minD :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic
minD r p d1 d2 = fst $ minD_ r p d1 d2
maxD_ :: RoundMode -> Precision -> Dyadic -> Dyadic -> (Dyadic, Int)
maxD_ r p d1 d2 = withDyadicsBA r p d1 d2 mpfr_max
minD_ :: RoundMode -> Precision -> Dyadic -> Dyadic -> (Dyadic, Int)
minD_ r p d1 d2 = withDyadicsBA r p d1 d2 mpfr_min
sgn :: Dyadic -> Int
sgn d = case compare zero d of
LT -> 1
EQ -> 0
_ -> -1
--conversion from dyadics to basic haskell types
dyadicToDouble :: RoundMode -> Dyadic -> Double
dyadicToDouble r mp1 = (realToFrac . unsafePerformIO) go
where go = with mp1 $ \p -> mpfr_get_d p ((fromIntegral . fromEnum) r)
dyadicToWord :: RoundMode -> Dyadic -> Word
dyadicToWord r mp1 = (fromIntegral . unsafePerformIO) go
where go = with mp1 $ \p -> mpfr_get_ui p ((fromIntegral . fromEnum) r)
dyadicToInt :: RoundMode -> Dyadic -> Int
dyadicToInt r mp1 = (fromIntegral . unsafePerformIO) go
where go = with mp1 $ \p -> mpfr_get_si p ((fromIntegral . fromEnum) r)
dyadicToString :: RoundMode -> Word -- ^ number of significant digits
-> Word -- ^ base
-> Dyadic -> (String, Int)
dyadicToString r n b mp1 = unsafePerformIO go
where go = with mp1 $ \p1 -> do
alloca $ \p2 -> do
withCString (replicate (fromIntegral (n + 2)) '0') $ \p3 -> do
_ <- mpfr_get_str p3 p2 (fromIntegral b) (fromIntegral n) p1 ((fromIntegral . fromEnum) r)
r1 <- peekCString p3
r2 <- peek p2
return (r1, fromIntegral r2)
decompose :: Dyadic -> (Integer, Int)
decompose d = (getMantissa d, getExp d - ceiling (fromIntegral (getPrec d) / fromIntegral bitsPerMPLimb :: Double) * bitsPerMPLimb)
toStringExp :: Word -> Dyadic -> String
toStringExp dec d = s ++ case e > 0 of
True -> case floor (logBase 10 2 * fromIntegral (getExp d) :: Double) > dec of
False -> take e ss ++ let bt = backtrim (drop e ss) in if null bt then "" else "." ++ bt
True -> head ss : "." ++ let bt = (backtrim . tail) ss in if null bt then "0"
else bt ++ "e" ++ show (pred e)
False -> head ss : "." ++ (let bt = (backtrim . tail) ss in
if null bt then "0"
else bt )
++ "e" ++ show (pred e)
where (str, e) = dyadicToString Near n 10 d
n = max dec 5
(s, ss) = case head str of
'-' -> ("-", tail str)
_ -> ("" , str)
backtrim = reverse . dropWhile (== '0') . reverse
toString :: Word -> Dyadic -> String
toString dec d = s ++ case compare 0 e of
LT -> take e ss ++ (let bt = all (== '0') (drop e ss) in if bt then "" else '.' : (drop e ss))
++ (if fromIntegral n - e < 0 then "e" ++ show (e - fromIntegral n) else "")
GT -> let ee = fromIntegral dec + e in
if ee <= 0 then "0" else
head ss : "." ++ (backtrim . tail . take ee) ss ++ "e" ++ show (pred e)
EQ -> "0." ++ let bt = all (== '0') ss in if bt then "0" else ss
where (str, e) = dyadicToString Near n 10 d
n = max dec 5
(s, ss) = case head str of
'-' -> ("-", tail str)
_ -> ("" , str)
backtrim = reverse . dropWhile (== '0') . reverse
--asignment functions
set :: RoundMode -> Precision -> Dyadic -> Dyadic
set r p d = fst $ set_ r p d
set_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)
set_ r p mp1 = unsafePerformIO go
where go = do ls <- mpfr_custom_get_size (fromIntegral p)
fp <- mallocForeignPtrBytes (fromIntegral ls)
let dummy = MP (fromIntegral p) 0 0 fp
with dummy $ \p1 -> do
with mp1 $ \p2 -> do
r2 <- mpfr_set p1 p2 ((fromIntegral . fromEnum) r)
r1 <- peekP p1 fp
return (r1, fromIntegral r2)
------------------------------------
-- mpfr constants
pi :: RoundMode -> Precision -> Dyadic
pi r p = fst $ pi_ r p
log2c :: RoundMode -> Precision -> Dyadic
log2c r p = fst $ pi_ r p
euler :: RoundMode -> Precision -> Dyadic
euler r p = fst $ pi_ r p
catalan :: RoundMode -> Precision -> Dyadic
catalan r p = fst $ pi_ r p
pi_ :: RoundMode -> Precision -> (Dyadic, Int)
pi_ r p = withDyadicC r p mpfr_const_pi
log2c_ :: RoundMode -> Precision -> (Dyadic, Int)
log2c_ r p = withDyadicC r p mpfr_const_log2
euler_ :: RoundMode -> Precision -> (Dyadic, Int)
euler_ r p = withDyadicC r p mpfr_const_euler
catalan_ :: RoundMode -> Precision -> (Dyadic, Int)
catalan_ r p = withDyadicC r p mpfr_const_catalan
----------------------------------------------------------
-- conversion from basic haskell types to dyadics
fromDouble :: RoundMode -> Precision -> Double -> Dyadic
fromDouble r p d = fst $ fromDouble_ r p d
fromInt :: RoundMode -> Precision -> Int -> Dyadic
fromInt r p d = fst $ fromInt_ r p d
fromWord :: RoundMode -> Precision -> Word -> Dyadic
fromWord r p d = fst $ fromWord_ r p d
fromDouble_ :: RoundMode -> Precision -> Double -> (Dyadic, Int)
fromDouble_ r p d = unsafePerformIO go
where go = do ls <- mpfr_custom_get_size (fromIntegral p)
fp <- mallocForeignPtrBytes (fromIntegral ls)
let dummy = MP (fromIntegral p) 0 0 fp
with dummy $ \p1 -> do
r2 <- mpfr_set_d p1 (realToFrac d) ((fromIntegral . fromEnum) r)
r1 <- peekP p1 fp
return (r1, fromIntegral r2)
fromInt_ :: RoundMode -> Precision -> Int -> (Dyadic, Int)
fromInt_ r p d = unsafePerformIO go
where go = do ls <- mpfr_custom_get_size (fromIntegral p)
fp <- mallocForeignPtrBytes (fromIntegral ls)
let dummy = MP (fromIntegral p) 0 0 fp
with dummy $ \p1 -> do
r2 <- mpfr_set_si p1 (fromIntegral d) ((fromIntegral . fromEnum) r)
r1 <- peekP p1 fp
return (r1, fromIntegral r2)
fromWord_ :: RoundMode -> Precision -> Word -> (Dyadic, Int)
fromWord_ r p d = unsafePerformIO go
where go = do ls <- mpfr_custom_get_size (fromIntegral p)
fp <- mallocForeignPtrBytes (fromIntegral ls)
let dummy = MP (fromIntegral p) 0 0 fp
with dummy $ \p1 -> do
r2 <- mpfr_set_ui p1 (fromIntegral d) ((fromIntegral . fromEnum) r)
r1 <- peekP p1 fp
return (r1, fromIntegral r2)
fromIntegerA :: RoundMode -> Precision -> Integer -> Dyadic
fromIntegerA r p d = stringToDyadic r p 10 (show d)
compose :: RoundMode -> Precision -> (Integer, Int) -> Dyadic
compose r p (i, ii) = div2i r p (fromIntegerA r p i) ii
fromString :: String -> Precision -> Word -> Dyadic
fromString s p b = stringToDyadic Near p b s
---------------------------------------------------------
-- functions getting properties of dyadics
getPrec :: Dyadic -> Precision
getPrec d = fromIntegral (withDyadicP d mpfr_get_prec)
-- | getMantissa and getExp return values such that
--
-- > d = getMantissa d * 2^(getExp d - ceiling ((getPrec d) / bitsPerMPLimb)* bitsPerMPLimb )
getMantissa :: Dyadic -> Integer
getMantissa d = if d < zero then -h else h
where (h, _) = foldl (\(a,b) c -> (a + (toInteger c) `shiftL` b, b + bitsPerMPLimb)) (0,0) (getMantissa' d)
getExp :: Dyadic -> Int
getExp d = (fromIntegral . unsafePerformIO) go
where go = do with d $ \p1 ->
mpfr_custom_get_exp p1
--------------------------------------------------------
-- some constants
minPrec :: Precision
minPrec = 32
one :: Dyadic
one = fromWord Near minPrec 1
zero :: Dyadic
zero = fromWord Near minPrec 0
-- instances
instance Eq Dyadic where
(==) = equal
instance Ord Dyadic where
(<) = less
(<=) = lesseq
(>) = greater
(>=) = greatereq
instance Show Dyadic where
show = toStringExp 16
-- these are exact operations, without rounding
instance Num Dyadic where
d + d' = add Zero (addPrec d d') d d'
d - d' = sub Zero (addPrec d d') d d'
d * d' = mul Zero (getPrec d + getPrec d') d d'
negate d = neg Zero (getPrec d) d
signum d = case compare d zero of
LT -> negate one
EQ -> zero
_ -> one
abs d = absD Zero (getPrec d) d
fromInteger i = fromIntegerA Zero (checkPrec $ binprec i) i
addPrec :: Dyadic -> Dyadic -> Precision
addPrec d1 d2 = fromIntegral (max (p1 + e1 - e3) (p2 + e2 - e3)) + 1
where e1 = if d1 == 0 then 0 else getExp d1
e2 = if d2 == 0 then 0 else getExp d2
p1 = fromIntegral $ getPrec d1
p2 = fromIntegral $ getPrec d2
e3 = min e1 e2
{-
addPrec d1 d2 = max e1 e2 + 1 - min (e1 - p1) (p2 - e2)
where e1 = if d1 == 0 then 0 else fromIntegral $ getExp d1
e2 = if d2 == 0 then 0 else fromIntegral $ getExp d2
p1 = getPrec d1
p2 = getPrec d2
-}