diff --git a/Data/Number/MPFR.hs b/Data/Number/MPFR.hs new file mode 100644 --- /dev/null +++ b/Data/Number/MPFR.hs @@ -0,0 +1,117 @@ +{-| + Module : Data.Number.MPFR + Description : top level + Copyright : (c) Aleš Bizjak + License : BSD3 + + Maintainer : ales.bizjak0@gmail.com + Stability : experimental + Portability : non-portable + + /Naming/ + + - functions ending with _ (underscore) usually return a pair (MPFR, Int), where + Int is a return value of a corresponding mpfr_ function. See the MPFR manual for + a description of return values. + + - the same functions without the _ return just the MPFR. + + - mpfr_ prefix in functions is removed + + - _ui and ui_ in function becomes w (stands for Word). For example mpfr_add_ui becomes addw. + + - si_ and _si in functions becomes i (stands for Int). + + - comparison functions which have _p appended loose it. For example mpfr_less_p becomes less. + + /Instances/ + + Eq + + - NaN \/= Nan, + + - Infinity = Infinity, + + - \-Infinity = -Infinity + + - otherwise normal comparison + + + + Ord + + - compare NaN _ = GT + + - compare _ NaN = GT + + - infinity < _ = False + + - \-infinity > _ = False + + - NaN [\<,\>,\>=,<=] _ = False + + This mimics the behaviour of built in Haskell Float and Double. + + + + Num + + Operations defined in Num class will be computed exactly (no precision is lost). + This isn't particularly useful as precision grows fairly quickly and everything becomes + slow so it preferably shouldn't be used. + + + + /This module should always be imported qualified./ + +-} + +{-# INCLUDE #-} +{-# INCLUDE #-} + + +module Data.Number.MPFR ( + module Data.Number.MPFR.Base +) where + +import Data.Number.MPFR.Base + +import Data.Number.MPFR.Internal + +import Data.Maybe + +import Data.Ratio + +instance Num MPFR 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 + abs d = absD Zero (getPrec d) d + signum d = fromInt Zero minPrec (fromMaybe (-1) (sgn d)) + fromInteger i = fromIntegerA Zero (checkPrec $ binprec i) i + -- TODO works only partially + +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 +-} + + +instance Real MPFR where + toRational d = n % 2 ^ e + where (n', e') = decompose d + (n, e) = if e' >= 0 then ((n' * 2 ^ e'), 0) + else (n', - e') + diff --git a/Data/Number/MPFR/Arithmetic.hs b/Data/Number/MPFR/Arithmetic.hs new file mode 100644 --- /dev/null +++ b/Data/Number/MPFR/Arithmetic.hs @@ -0,0 +1,224 @@ +{-| + Module : Data.Number.MPFR.Arithmetic + Description : wrappers for basic arithmetic functions + Copyright : (c) Aleš Bizjak + License : BSD3 + + Maintainer : ales.bizjak0@gmail.com + Stability : experimental + Portability : non-portable + + Basic arithmetic functions. See MPFR manual for detailed documentation. +-} +{-# INCLUDE #-} +{-# INCLUDE #-} + +module Data.Number.MPFR.Arithmetic + where + +import Data.Number.MPFR.Internal + +add :: RoundMode -> Precision -> MPFR -> MPFR -> MPFR +add r p d1 d2 = fst $ add_ r p d1 d2 + +addw :: RoundMode -> Precision -> MPFR -> Word -> MPFR +addw r p d1 d = fst $ addw_ r p d1 d + +addi :: RoundMode -> Precision -> MPFR -> Int -> MPFR +addi r p d1 d = fst $ addi_ r p d1 d + +sub :: RoundMode -> Precision -> MPFR -> MPFR -> MPFR +sub r p d1 d2 = fst $ sub_ r p d1 d2 + +subw :: RoundMode -> Precision -> MPFR -> Word -> MPFR +subw r p d1 d = fst $ subw_ r p d1 d + +subi :: RoundMode -> Precision -> MPFR -> Int -> MPFR +subi r p d1 d = fst $ subi_ r p d1 d + +wsub :: RoundMode -> Precision -> Word -> MPFR -> MPFR +wsub r p d d1 = fst $ wsub_ r p d d1 + +isub :: RoundMode -> Precision -> Int -> MPFR -> MPFR +isub r p d d1 = fst $ isub_ r p d d1 + +mul :: RoundMode -> Precision -> MPFR -> MPFR -> MPFR +mul r p d1 d2 = fst $ mul_ r p d1 d2 + +mulw :: RoundMode -> Precision -> MPFR -> Word -> MPFR +mulw r p d1 d = fst $ mulw_ r p d1 d + +muli :: RoundMode -> Precision -> MPFR -> Int -> MPFR +muli r p d1 d = fst $ muli_ r p d1 d + +sqr :: RoundMode -> Precision -> MPFR -> MPFR +sqr r p d = fst $ sqr_ r p d + +div :: RoundMode -> Precision -> MPFR -> MPFR -> MPFR +div r p d1 d2 = fst $ div_ r p d1 d2 + +divw :: RoundMode -> Precision -> MPFR -> Word -> MPFR +divw r p d1 d = fst $ divw_ r p d1 d + +divi :: RoundMode -> Precision -> MPFR -> Int -> MPFR +divi r p d1 d = fst $ divi_ r p d1 d + +wdiv :: RoundMode -> Precision -> Word -> MPFR -> MPFR +wdiv r p d d1 = fst $ wdiv_ r p d d1 + +idiv :: RoundMode -> Precision -> Int -> MPFR -> MPFR +idiv r p d d1 = fst $ idiv_ r p d d1 + +sqrt :: RoundMode -> Precision -> MPFR -> MPFR +sqrt r p d = fst $ sqrt_ r p d + +sqrtw :: RoundMode -> Precision -> Word -> MPFR +sqrtw r p d = fst $ sqrtw_ r p d + +cbrt :: RoundMode -> Precision -> MPFR -> MPFR +cbrt r p d = fst $ cbrt_ r p d + +root :: RoundMode -> Precision -> MPFR -> Word -> MPFR +root r p d n = fst $ root_ r p d n + +pow :: RoundMode -> Precision -> MPFR -> MPFR -> MPFR +pow r p d1 d2 = fst $ pow_ r p d1 d2 + +poww :: RoundMode -> Precision -> MPFR -> Word -> MPFR +poww r p d1 d2 = fst $ poww_ r p d1 d2 + +powi :: RoundMode -> Precision -> MPFR -> Int -> MPFR +powi r p d1 d2 = fst $ powi_ r p d1 d2 + +wpoww :: RoundMode -> Precision -> Word -> Word -> MPFR +wpoww r p d1 d2 = fst $ wpoww_ r p d1 d2 + +wpow :: RoundMode -> Precision -> Word -> MPFR -> MPFR +wpow r p d1 d2 = fst $ wpow_ r p d1 d2 + +neg :: RoundMode -> Precision -> MPFR -> MPFR +neg r p d = fst $ neg_ r p d + +absD :: RoundMode -> Precision -> MPFR -> MPFR +absD r p d = fst $ absD_ r p d + +dim :: RoundMode -> Precision -> MPFR -> MPFR -> MPFR +dim r p d1 d2 = fst $ dim_ r p d1 d2 + +mul2w :: RoundMode -> Precision -> MPFR -> Word -> MPFR +mul2w r p d1 d2 = fst $ mul2w_ r p d1 d2 + +mul2i :: RoundMode -> Precision -> MPFR -> Int -> MPFR +mul2i r p d1 d2 = fst $ mul2i_ r p d1 d2 + +div2w :: RoundMode -> Precision -> MPFR -> Word -> MPFR +div2w r p d1 d2 = fst $ div2w_ r p d1 d2 + +div2i :: RoundMode -> Precision -> MPFR -> Int -> MPFR +div2i r p d1 d2 = fst $ div2i_ r p d1 d2 + +add_ :: RoundMode -> Precision -> MPFR -> MPFR -> (MPFR,Int) +add_ r p d1 d2 = withMPFRsBA r p d1 d2 mpfr_add + +addw_ :: RoundMode -> Precision -> MPFR -> Word -> (MPFR, Int) +addw_ r p d1 d = withMPFRBAui r p d1 (fromIntegral d) mpfr_add_ui + +addi_ :: RoundMode -> Precision -> MPFR -> Int -> (MPFR, Int) +addi_ r p d1 d = withMPFRBAsi r p d1 (fromIntegral d) mpfr_add_si + +sub_ :: RoundMode -> Precision -> MPFR -> MPFR -> (MPFR,Int) +sub_ r p d1 d2 = withMPFRsBA r p d1 d2 mpfr_sub + +subw_ :: RoundMode -> Precision -> MPFR -> Word -> (MPFR, Int) +subw_ r p d1 d = withMPFRBAui r p d1 (fromIntegral d) mpfr_sub_ui + +subi_ :: RoundMode -> Precision -> MPFR -> Int -> (MPFR, Int) +subi_ r p d1 d = withMPFRBAsi r p d1 (fromIntegral d) mpfr_sub_si + +wsub_ :: RoundMode -> Precision -> Word -> MPFR -> (MPFR, Int) +wsub_ r p d d1 = withMPFRBAiu r p (fromIntegral d) d1 mpfr_ui_sub + +isub_ :: RoundMode -> Precision -> Int -> MPFR -> (MPFR, Int) +isub_ r p d d1 = withMPFRBAis r p (fromIntegral d) d1 mpfr_si_sub + +mul_ :: RoundMode -> Precision -> MPFR -> MPFR -> (MPFR,Int) +mul_ r p d1 d2 = withMPFRsBA r p d1 d2 mpfr_mul + +mulw_ :: RoundMode -> Precision -> MPFR -> Word -> (MPFR, Int) +mulw_ r p d1 d = withMPFRBAui r p d1 (fromIntegral d) mpfr_mul_ui + +muli_ :: RoundMode -> Precision -> MPFR -> Int -> (MPFR, Int) +muli_ r p d1 d = withMPFRBAsi r p d1 (fromIntegral d) mpfr_mul_si + +sqr_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +sqr_ r p d = withMPFR r p d mpfr_sqr + +div_ :: RoundMode -> Precision -> MPFR -> MPFR -> (MPFR, Int) +div_ r p d1 d2 = withMPFRsBA r p d1 d2 mpfr_div + +divw_ :: RoundMode -> Precision -> MPFR -> Word -> (MPFR, Int) +divw_ r p d1 d = withMPFRBAui r p d1 (fromIntegral d) mpfr_div_ui + +divi_ :: RoundMode -> Precision -> MPFR -> Int -> (MPFR, Int) +divi_ r p d1 d = withMPFRBAsi r p d1 (fromIntegral d) mpfr_div_si + +wdiv_ :: RoundMode -> Precision -> Word -> MPFR -> (MPFR, Int) +wdiv_ r p d d1 = withMPFRBAiu r p (fromIntegral d) d1 mpfr_ui_div + +idiv_ :: RoundMode -> Precision -> Int -> MPFR -> (MPFR, Int) +idiv_ r p d d1 = withMPFRBAis r p (fromIntegral d) d1 mpfr_si_div + +sqrt_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +sqrt_ r p d = withMPFR r p d mpfr_sqrt + +sqrtw_ :: RoundMode -> Precision -> Word -> (MPFR, Int) +sqrtw_ r p d = withMPFRUI r p d mpfr_sqrt_ui + +cbrt_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +cbrt_ r p d = withMPFR r p d mpfr_cbrt + +root_ :: RoundMode -> Precision -> MPFR -> Word -> (MPFR, Int) +root_ r p d n = withMPFRBAui r p d (fromIntegral n) mpfr_root + +pow_ :: RoundMode -> Precision -> MPFR -> MPFR -> (MPFR, Int) +pow_ r p d1 d2 = withMPFRsBA r p d1 d2 mpfr_pow + +poww_ :: RoundMode -> Precision -> MPFR -> Word -> (MPFR , Int) +poww_ r p d1 d2 = withMPFRBAui r p d1 (fromIntegral d2) mpfr_pow_ui + +powi_ :: RoundMode -> Precision -> MPFR -> Int -> (MPFR , Int) +powi_ r p d1 d2 = withMPFRBAsi r p d1 (fromIntegral d2) mpfr_pow_si + +wpoww_ :: RoundMode -> Precision -> Word -> Word -> (MPFR , 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 -> MPFR -> (MPFR , Int) +wpow_ r p d1 d2 = withMPFRBAiu r p (fromIntegral d1) d2 mpfr_ui_pow + +neg_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +neg_ r p d = withMPFR r p d mpfr_neg + +absD_ :: RoundMode -> Precision -> MPFR -> (MPFR , Int) +absD_ r p d = withMPFR r p d mpfr_abs + +dim_ :: RoundMode -> Precision -> MPFR -> MPFR -> (MPFR, Int) +dim_ r p d1 d2 = withMPFRsBA r p d1 d2 mpfr_dim + +mul2w_ :: RoundMode -> Precision -> MPFR -> Word -> (MPFR, Int) +mul2w_ r p d1 d2 = withMPFRBAui r p d1 (fromIntegral d2) mpfr_mul_2ui + +mul2i_ :: RoundMode -> Precision -> MPFR -> Int -> (MPFR, Int) +mul2i_ r p d1 d2 = withMPFRBAsi r p d1 (fromIntegral d2) mpfr_mul_2si + +div2w_ :: RoundMode -> Precision -> MPFR -> Word -> (MPFR, Int) +div2w_ r p d1 d2 = withMPFRBAui r p d1 (fromIntegral d2) mpfr_div_2ui + +div2i_ :: RoundMode -> Precision -> MPFR -> Int -> (MPFR, Int) +div2i_ r p d1 d2 = withMPFRBAsi r p d1 (fromIntegral d2) mpfr_div_2si diff --git a/Data/Number/MPFR/Assignment.hs b/Data/Number/MPFR/Assignment.hs new file mode 100644 --- /dev/null +++ b/Data/Number/MPFR/Assignment.hs @@ -0,0 +1,172 @@ +{-| + Module : Data.Number.MPFR.Assignment + Description : wrappers for assignment functions + Copyright : (c) Aleš Bizjak + License : BSD3 + + Maintainer : ales.bizjak0@gmail.com + Stability : experimental + Portability : non-portable + + Conversion from basic Haskell types to MPFR. See MPFR manual for detailed documentation. +-} + +{-# INCLUDE #-} +{-# INCLUDE #-} + + +module Data.Number.MPFR.Assignment where + +import Data.Number.MPFR.Internal + +import Data.Number.MPFR.Arithmetic + +set :: RoundMode -> Precision -> MPFR -> MPFR +set r p d = fst $ set_ r p d + +set_ :: RoundMode -> Precision -> MPFR -> (MPFR, 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) + +fromWord :: RoundMode -> Precision -> Word -> MPFR +fromWord r p d = fst $ fromWord_ r p d + +fromInt :: RoundMode -> Precision -> Int -> MPFR +fromInt r p d = fst $ fromInt_ r p d + +fromDouble :: RoundMode -> Precision -> Double -> MPFR +fromDouble r p d = fst $ fromDouble_ r p d + +fromWord_ :: RoundMode -> Precision -> Word -> (MPFR, 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) + +fromInt_ :: RoundMode -> Precision -> Int -> (MPFR, 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) + +fromDouble_ :: RoundMode -> Precision -> Double -> (MPFR, 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) + +-- x * 2 ^ y +int2w :: RoundMode -> Precision -> Word -> Int -> MPFR +int2w r p i e = fst $ int2w_ r p i e + +int2i :: RoundMode -> Precision -> Int -> Int -> MPFR +int2i r p i e = fst $ int2i_ r p i e + +int2w_ :: RoundMode -> Precision -> Word -> Int -> (MPFR, 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) + +int2i_ :: RoundMode -> Precision -> Int -> Int -> (MPFR, 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) + + +stringToMPFR :: RoundMode -> Precision + -> Word -- ^ Base + -> String -> MPFR +stringToMPFR r p b d = fst $ stringToMPFR_ r p b d + +stringToMPFR_ :: RoundMode -> Precision + -> Word -- ^ Base + -> String -> (MPFR, Int) +stringToMPFR_ 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 + r2 <- mpfr_set_str p1 p2 (fromIntegral b) ((fromIntegral . fromEnum) r) + r1 <- peekP p1 fp + return (r1, fromIntegral r2) + +strtofr :: RoundMode -> Precision + -> Word -- ^ base + -> String -> (MPFR, String) +strtofr r p b d = case strtofr_ r p b d of + (a, b', _) -> (a,b') + +strtofr_ :: RoundMode -> Precision + -> Word -- ^ base + -> String -> (MPFR, String, Int) +strtofr_ 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 + alloca $ \p3 -> do + r3 <- mpfr_strtofr p1 p2 p3 (fromIntegral b) ((fromIntegral . fromEnum) r) + p3' <- peek p3 + r2 <- peekCString p3' + r1 <- peekP p1 fp + return (r1, r2, fromIntegral r3) + + +setInf :: Precision -> Int -> MPFR +setInf p i = 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 + mpfr_set_inf p1 (fromIntegral i) + peekP p1 fp + +setNaN :: Precision -> MPFR +setNaN p = 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 + mpfr_set_nan p1 + peekP p1 fp + +fromIntegerA :: RoundMode -> Precision -> Integer -> MPFR +fromIntegerA r p d = stringToMPFR r p 10 (show d) + +compose :: RoundMode -> Precision -> (Integer, Int) -> MPFR +compose r p (i, ii) = div2i r p (fromIntegerA r p i) ii + +-- | @stringToMPFR@ with default rounding to Near. +fromString :: String -> Precision -> Word -> MPFR +fromString s p b = stringToMPFR Near p b s diff --git a/Data/Number/MPFR/Base.hs b/Data/Number/MPFR/Base.hs new file mode 100644 --- /dev/null +++ b/Data/Number/MPFR/Base.hs @@ -0,0 +1,29 @@ +{-# INCLUDE #-} +{-# INCLUDE #-} + +module Data.Number.MPFR.Base ( + module Data.Number.MPFR.Assignment, + module Data.Number.MPFR.Conversion, + module Data.Number.MPFR.Arithmetic, + module Data.Number.MPFR.Comparison, + module Data.Number.MPFR.Special, + module Data.Number.MPFR.Integer, + module Data.Number.MPFR.Misc, + RoundMode (Near, Up, Down, Zero), + MPFR, Precision, Exp, Dyadic + ) +where + +import Data.Number.MPFR.Assignment +import Data.Number.MPFR.Conversion +import Data.Number.MPFR.Arithmetic +import Data.Number.MPFR.Comparison +import Data.Number.MPFR.Special +import Data.Number.MPFR.Integer +import Data.Number.MPFR.Misc + + +import Data.Number.MPFR.Internal + +type Dyadic = MPFR + diff --git a/Data/Number/MPFR/Comparison.hs b/Data/Number/MPFR/Comparison.hs new file mode 100644 --- /dev/null +++ b/Data/Number/MPFR/Comparison.hs @@ -0,0 +1,127 @@ +{-| + Module : Data.Number.MPFR.Comparison + Description : wrappers for comparison functions + Copyright : (c) Aleš Bizjak + License : BSD3 + + Maintainer : ales.bizjak0@gmail.com + Stability : experimental + Portability : non-portable + + Comparison functions. All the functions that return Maybe Ordering return Nothing + when one of the operands is NaN and Just _ otherwise. + +-} + +{-# INCLUDE #-} +{-# INCLUDE #-} + +module Data.Number.MPFR.Comparison where + +import Data.Number.MPFR.Misc + +import Data.Number.MPFR.Internal + +import Prelude hiding (isNaN) + +import Data.Maybe + +cmp :: MPFR -> MPFR -> Maybe Ordering +cmp mp1 mp2 = if isNaN mp1 || isNaN mp2 then Nothing + else Just (compare (withMPFRBB mp1 mp2 mpfr_cmp) 0) + +cmpw :: MPFR -> Word -> Maybe Ordering +cmpw mp1 w = if isNaN mp1 then Nothing else Just (compare (unsafePerformIO go) 0) + where go = with mp1 $ \p -> mpfr_cmp_ui p (fromIntegral w) + +cmpi :: MPFR -> Int -> Maybe Ordering +cmpi mp1 i = if isNaN mp1 then Nothing else Just (compare (unsafePerformIO go) 0) + where go = with mp1 $ \p -> mpfr_cmp_si p (fromIntegral i) + +cmpd :: MPFR -> Double -> Maybe Ordering +cmpd mp1 d = unsafePerformIO go + where go = do mpfr_clear_erangeflag + with mp1 $ \p -> do + r1 <- mpfr_cmp_d p (realToFrac d) + r2 <- mpfr_erangeflag_p + if r2 == 0 then return (Just (compare r1 0)) + else do mpfr_clear_erangeflag + return Nothing + +cmp2w :: MPFR -> Word -> Exp -> Maybe Ordering +cmp2w d w e = unsafePerformIO go + where go = do mpfr_clear_erangeflag + with d $ \p -> do + r1 <- mpfr_cmp_ui_2exp p (fromIntegral w) e + r2 <- mpfr_erangeflag_p + if r2 == 0 then return (Just (compare r1 0)) + else do mpfr_clear_erangeflag + return Nothing + +cmp2i :: MPFR -> Int -> Exp -> Maybe Ordering +cmp2i d w e = unsafePerformIO go + where go = do mpfr_clear_erangeflag + with d $ \p -> do + r1 <- mpfr_cmp_si_2exp p (fromIntegral w) e + r2 <- mpfr_erangeflag_p + if r2 == 0 then return (Just (compare r1 0)) + else do mpfr_clear_erangeflag + return Nothing +cmpabs :: MPFR -> MPFR -> Maybe Ordering +cmpabs mp1 mp2 = if isNaN mp1 || isNaN mp2 then Nothing + else Just (compare (withMPFRBB mp1 mp2 mpfr_cmpabs) 0) + +isNaN :: MPFR -> Bool +isNaN d = withMPFRB d mpfr_nan_p /= 0 + +isInfinite :: MPFR -> Bool +isInfinite d = withMPFRB d mpfr_inf_p /= 0 + +isNumber :: MPFR -> Bool +isNumber d = withMPFRB d mpfr_number_p /= 0 + +isZero :: MPFR -> Bool +isZero d = withMPFRB d mpfr_zero_p /= 0 + +sgn :: MPFR -> Maybe Int +sgn mp1 = case (cmpw mp1 0) of + Nothing -> Nothing + Just x -> Just (pred . fromEnum $ x) + +-- TODO Maybe Bool???? +greater :: MPFR -> MPFR -> Bool +greater d1 d2 = withMPFRBB d1 d2 mpfr_greater_p /= 0 + +greatereq :: MPFR -> MPFR -> Bool +greatereq d1 d2 = withMPFRBB d1 d2 mpfr_greaterequal_p /= 0 + +less :: MPFR -> MPFR -> Bool +less d1 d2 = withMPFRBB d1 d2 mpfr_less_p /= 0 + +lesseq :: MPFR -> MPFR -> Bool +lesseq d1 d2 = withMPFRBB d1 d2 mpfr_lessequal_p /= 0 + +lessgreater :: MPFR -> MPFR -> Maybe Bool +lessgreater d1 d2 = if isNaN d1 || isNaN d2 then Nothing + else Just (withMPFRBB d1 d2 mpfr_lessgreater_p /= 0) + +equal :: MPFR -> MPFR -> Bool +equal d1 d2 = withMPFRBB d1 d2 mpfr_equal_p /= 0 + +unordered :: MPFR -> MPFR -> Maybe Bool +unordered d1 d2 = if isNaN d1 || isNaN d2 then Nothing + else Just (withMPFRBB d1 d2 mpfr_unordered_p /= 0) + + +instance Eq MPFR where + (==) = equal + +instance Ord MPFR where + compare d d' = fromMaybe GT (cmp d d') + (<) = less + (<=) = lesseq + (>) = greater + (>=) = greatereq + max d d' = maxD Zero (maxPrec d d') d d' + min d d' = minD Zero (maxPrec d d') d d' + diff --git a/Data/Number/MPFR/Conversion.hs b/Data/Number/MPFR/Conversion.hs new file mode 100644 --- /dev/null +++ b/Data/Number/MPFR/Conversion.hs @@ -0,0 +1,126 @@ +{-| + Module : Data.Number.MPFR.Conversion + Description : wrappers for conversion functions + Copyright : (c) Aleš Bizjak + License : BSD3 + + Maintainer : ales.bizjak0@gmail.com + Stability : experimental + Portability : non-portable + + Conversion from MPFR to basic Haskell types. See MPFR manual for detailed documentation. +-} + +{-# INCLUDE #-} +{-# INCLUDE #-} + + +module Data.Number.MPFR.Conversion where + +import Data.Number.MPFR.Internal +import Data.Number.MPFR.Misc + +import Data.List (isInfixOf) + +toDouble :: RoundMode -> MPFR -> Double +toDouble r mp1 = (realToFrac . unsafePerformIO) go + where go = with mp1 $ \p -> mpfr_get_d p ((fromIntegral . fromEnum) r) + +toDouble2exp :: RoundMode -> MPFR -> (Double, Int) +toDouble2exp r mp1 = unsafePerformIO go + where go = do with mp1 $ \p1 -> do + alloca $ \p2 -> do + r1 <- mpfr_get_d_2exp p2 p1 ((fromIntegral . fromEnum) r) + r2 <- peek p2 + return (realToFrac r1, fromIntegral r2) + +toInt :: RoundMode -> MPFR -> Int +toInt r mp1 = (fromIntegral . unsafePerformIO) go + where go = with mp1 $ \p -> mpfr_get_si p ((fromIntegral . fromEnum) r) + +toWord :: RoundMode -> MPFR -> Word +toWord r mp1 = (fromIntegral . unsafePerformIO) go + where go = with mp1 $ \p -> mpfr_get_ui p ((fromIntegral . fromEnum) r) + + +mpfrToString :: RoundMode + -> Word -- ^ number of decimals + -> Word -- ^ base + -> MPFR -> (String, Exp) +mpfrToString r n b mp1 = unsafePerformIO go + where go = with mp1 $ \p1 -> do + alloca $ \p2 -> do + p3 <- mpfr_get_str nullPtr p2 (fromIntegral b) (fromIntegral n) p1 ((fromIntegral . fromEnum) r) + r1 <- peekCString p3 + r2 <- peek p2 + mpfr_free_str p3 + return (r1, r2) + +fitsULong :: RoundMode -> MPFR -> Bool +fitsULong r d = withMPFRF d r mpfr_fits_ulong_p /= 0 + +fitsSLong :: RoundMode -> MPFR -> Bool +fitsSLong r d = withMPFRF d r mpfr_fits_slong_p /= 0 + +fitsUInt :: RoundMode -> MPFR -> Bool +fitsUInt r d = withMPFRF d r mpfr_fits_uint_p /= 0 + +fitsSInt :: RoundMode -> MPFR -> Bool +fitsSInt r d = withMPFRF d r mpfr_fits_sint_p /= 0 + +fitsUShort :: RoundMode -> MPFR -> Bool +fitsUShort r d = withMPFRF d r mpfr_fits_ushort_p /= 0 + +fitsSShort :: RoundMode -> MPFR -> Bool +fitsSShort r d = withMPFRF d r mpfr_fits_sshort_p /= 0 + +-- TODO +decompose :: MPFR -> (Integer, Exp) +decompose d = (getMantissa d, getExp d - fromIntegral (Prelude.ceiling (fromIntegral (getPrec d) / fromIntegral bitsPerMPLimb :: Double) * bitsPerMPLimb)) + +-- | Output a string in base 10 rounded to Near in exponential form. +toStringExp :: Word -- ^ number of digits + -> MPFR -> String +toStringExp dec d = + if isInfixOf "NaN" ss then "NaN" + else if isInfixOf "Inf" ss then s ++ "Infinity" + else s ++ case e > 0 of + True -> case Prelude.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') = mpfrToString Near n 10 d + e = fromIntegral e' + n = max dec 5 + (s, ss) = case head str of + '-' -> ("-", tail str) + _ -> ("" , str) + backtrim = reverse . dropWhile (== '0') . reverse + +-- | Output a string in base 10 rounded to Near. The difference from @toStringExp@ is that +-- it won't output in exponential form if it is sensible to do so. +toString :: Word -> MPFR -> String +toString dec d = + if isInfixOf "NaN" ss then "NaN" + else if isInfixOf "Inf" ss then s ++ "Infinity" + else 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') = mpfrToString Near n 10 d + n = max dec 5 + e = fromIntegral e' + (s, ss) = case head str of + '-' -> ("-", tail str) + _ -> ("" , str) + backtrim = reverse . dropWhile (== '0') . reverse + +instance Show MPFR where + show = toStringExp 16 diff --git a/Data/Number/MPFR/Down.hs b/Data/Number/MPFR/Down.hs new file mode 100644 --- /dev/null +++ b/Data/Number/MPFR/Down.hs @@ -0,0 +1,81 @@ +{-| + Module : Data.Number.MPFR.Down + Description : top level + Copyright : (c) Aleš Bizjak + License : BSD3 + + Maintainer : ales.bizjak0@gmail.com + Stability : experimental + Portability : non-portable + + This module defines instances Num, Real, Fractional, Floating and RealFrac of MPFR. + Operations are rounded with RoundMode Down and computed with max precision of two + operands or with the precision of the operand. Otherwise it is equivalent to + Data.Number.MPFR +-} +{-# INCLUDE #-} +{-# INCLUDE #-} + +module Data.Number.MPFR.Down ( + module Data.Number.MPFR.Base +) +where + +import Data.Number.MPFR.Base + +import Data.Number.MPFR.Internal + +import Data.Maybe + +import Data.Ratio + +instance Num MPFR where + d + d' = add Down (maxPrec d d') d d' + d - d' = sub Down (maxPrec d d') d d' + d * d' = mul Down (maxPrec d d') d d' + negate d = neg Down (getPrec d) d + abs d = absD Down (getPrec d) d + signum d = fromInt Down minPrec (fromMaybe (-1) (sgn d)) + fromInteger i = fromIntegerA Zero (checkPrec $ binprec i) i + -- TODO works only partially + +instance Real MPFR where + toRational d = n % 2 ^ e + where (n', e') = decompose d + (n, e) = if e' >= 0 then ((n' * 2 ^ e'), 0) + else (n', - e') + +instance Fractional MPFR where + d / d' = Data.Number.MPFR.Base.div Down (maxPrec d d') d d' + fromRational r = (fromInteger n) / (fromInteger d) + where n = numerator r + d = denominator r + recip d = one / d + +instance Floating MPFR where + pi = Data.Number.MPFR.Base.pi Down 53 + exp d = Data.Number.MPFR.Base.exp Down (getPrec d) d + log d = Data.Number.MPFR.Base.log Down (getPrec d) d + sqrt d = Data.Number.MPFR.Base.sqrt Down (getPrec d) d + (**) d d' = Data.Number.MPFR.Base.pow Down (maxPrec d d') d d' + logBase d d' = Prelude.log d' / Prelude.log d + sin d = Data.Number.MPFR.Base.sin Down (getPrec d) d + cos d = Data.Number.MPFR.Base.cos Down (getPrec d) d + tan d = Data.Number.MPFR.Base.tan Down (getPrec d) d + asin d = Data.Number.MPFR.Base.asin Down (getPrec d) d + acos d = Data.Number.MPFR.Base.acos Down (getPrec d) d + atan d = Data.Number.MPFR.Base.atan Down (getPrec d) d + sinh d = Data.Number.MPFR.Base.sinh Down (getPrec d) d + cosh d = Data.Number.MPFR.Base.cosh Down (getPrec d) d + tanh d = Data.Number.MPFR.Base.tanh Down (getPrec d) d + asinh d = Data.Number.MPFR.Base.asinh Down (getPrec d) d + acosh d = Data.Number.MPFR.Base.acosh Down (getPrec d) d + atanh d = Data.Number.MPFR.Base.atanh Down (getPrec d) d + +instance RealFrac MPFR where + properFraction d = (fromIntegral n, f) + where r = toRational d + m = numerator r + e = denominator r + n = quot m e + f = frac Down (getPrec d) d diff --git a/Data/Number/MPFR/FFIhelper.hsc b/Data/Number/MPFR/FFIhelper.hsc new file mode 100644 --- /dev/null +++ b/Data/Number/MPFR/FFIhelper.hsc @@ -0,0 +1,727 @@ +{-# LANGUAGE ForeignFunctionInterface, DeriveDataTypeable #-} +#include +#include + +module Data.Number.MPFR.FFIhelper where + +import Data.Word + +import Data.Int + +import Foreign.C.String(CString) +import Foreign.C.Types(CULong, CLong, CInt, CUInt, CDouble, CChar) +import Foreign.Ptr(FunPtr, Ptr) + +import Foreign.Storable +import Foreign.ForeignPtr (ForeignPtr, withForeignPtr) + +import Data.Typeable(Typeable) + +data RoundMode = Near | Zero | Up | Down | GMP_RND_MAX | GMP_RNDNA + +instance Enum RoundMode where + fromEnum Near = #{const GMP_RNDN} + fromEnum Zero = #{const GMP_RNDZ} + fromEnum Up = #{const GMP_RNDU} + fromEnum Down = #{const GMP_RNDD} + fromEnum GMP_RND_MAX = #{const GMP_RND_MAX} + fromEnum GMP_RNDNA = #{const GMP_RNDNA} + + toEnum #{const GMP_RNDN} = Near + toEnum #{const GMP_RNDZ} = Zero + toEnum #{const GMP_RNDU} = Up + toEnum #{const GMP_RNDD} = Down + toEnum #{const GMP_RND_MAX} = GMP_RND_MAX + toEnum (#{const GMP_RNDNA}) = GMP_RNDNA + toEnum i = error $ "RoundMode.toEnum called with illegal argument :" ++ show i + + +data MPFR = MP { precision :: !CPrecision, + sign :: !Sign, + exponent :: !Exp, + limbs :: !(ForeignPtr Limb) +} deriving (Typeable) + +instance Storable MPFR where + sizeOf _ = #size __mpfr_struct + alignment _ = (undefined :: Int) + peek = error "MPFR.peek: Not needed and not applicable" + poke p (MP prec s e fp) = do #{poke __mpfr_struct, _mpfr_prec} p prec + #{poke __mpfr_struct, _mpfr_sign} p s + #{poke __mpfr_struct, _mpfr_exp} p e + withForeignPtr fp $ \p1 -> #{poke __mpfr_struct, _mpfr_d} p p1 + +peekP :: Ptr MPFR -> ForeignPtr Limb -> IO MPFR +peekP p fp = do r11 <- #{peek __mpfr_struct, _mpfr_prec} p + r21 <- #{peek __mpfr_struct, _mpfr_sign} p + r22 <- #{peek __mpfr_struct, _mpfr_exp} p + return (MP r11 r21 r22 fp) + +bitsPerMPLimb :: Int +bitsPerMPLimb = 8 * #size mp_limb_t + +type CRoundMode = CInt + +type Limb = #type mp_limb_t + +type Sign = #type mpfr_sign_t + +type CPrecision = #type mpfr_prec_t + +type Exp = #type mp_exp_t + +type MpSize = #type mp_size_t + +-- utility functions from chsmpfr.h +foreign import ccall unsafe "initS" + initS :: CPrecision -> IO (Ptr MPFR) + +foreign import ccall unsafe "&clear" + clear :: FunPtr (Ptr MPFR -> IO ()) + +-------------------- +foreign import ccall unsafe "mpfr_get_prec_wrap" + mpfr_get_prec :: Ptr MPFR -> IO CPrecision + +---------------------------------------------------------------- + +-- assignment functions +foreign import ccall unsafe "mpfr_set_wrap" + mpfr_set :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_set_ui_wrap" + mpfr_set_ui :: Ptr MPFR -> CULong -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_set_si_wrap" + mpfr_set_si :: Ptr MPFR -> CLong -> CRoundMode -> IO CInt + +--TODO set_uj, set_sj + +foreign import ccall unsafe "mpfr_set_d" + mpfr_set_d :: Ptr MPFR -> CDouble -> CRoundMode -> IO CInt + +--foreign import ccall unsafe "mfpr_set_ld" + -- mpfr_set_ld :: Ptr MPFR -> #{type long double} -> CRoundMode -> IO CInt +--long double does not seem to be supported + +--TODO set_decimal64, set_z, set_q, set_f + +foreign import ccall unsafe "mpfr_set_ui_2exp" + mpfr_set_ui_2exp :: Ptr MPFR -> CULong -> Exp -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_set_si_2exp" + mpfr_set_si_2exp :: Ptr MPFR -> CLong -> Exp -> CRoundMode -> IO CInt + +--TODO set_uj_2exp, set_sj_2exp + +foreign import ccall unsafe "mpfr_set_str" + mpfr_set_str :: Ptr MPFR -> CString -> CInt -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_strtofr" + mpfr_strtofr :: Ptr MPFR -> CString -> Ptr (Ptr CChar) -> CInt -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_set_inf" + mpfr_set_inf :: Ptr MPFR -> CInt -> IO () + +foreign import ccall unsafe "mpfr_set_nan" + mpfr_set_nan :: Ptr MPFR -> IO () + +foreign import ccall unsafe "mpfr_swap" + mpfr_swap :: Ptr MPFR -> Ptr MPFR -> IO () + +-- THINK combined initialization and assignment functions are non-applicable +-- with custom interface? + +-------------------------------------------------------------------------------- + + +-- conversion functions +foreign import ccall unsafe "mpfr_get_d" + mpfr_get_d :: Ptr MPFR -> CRoundMode -> IO CDouble + +-- TODO get_decimal64 + +foreign import ccall unsafe "mpfr_get_d_2exp" + mpfr_get_d_2exp :: Ptr CLong -> Ptr MPFR -> CRoundMode -> IO CDouble + +-- TODO get_ld_2exp +-- !!!!!!! next 4 set erange flags +foreign import ccall unsafe "mpfr_get_si" + mpfr_get_si :: Ptr MPFR -> CRoundMode -> IO CLong + +foreign import ccall unsafe "mpfr_get_ui" + mpfr_get_ui :: Ptr MPFR -> CRoundMode -> IO CULong + +{- +foreign import ccall unsafe "mpfr_get_sj_wrap" + mpfr_get_sj :: Ptr MPFR -> CRoundMode -> IO #type intmax_t + +foreign import ccall unsafe "mpfr_get_uj_wrap" + mpft_get_uj :: Ptr MPFR -> CRoundMode -> IO #type uintmax_t +-} +--TODO get_z_exp, get_z, get_f, + +foreign import ccall unsafe "mpfr_get_str" + mpfr_get_str :: CString -> Ptr Exp -> CInt -> CUInt -> Ptr MPFR -> CRoundMode -> IO CString + +foreign import ccall unsafe "mpfr_free_str" + mpfr_free_str :: CString -> IO () + +foreign import ccall unsafe "mpfr_fits_ulong_p" + mpfr_fits_ulong_p :: Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_fits_slong_p" + mpfr_fits_slong_p :: Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_fits_uint_p" + mpfr_fits_uint_p :: Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_fits_sint_p" + mpfr_fits_sint_p :: Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_fits_ushort_p" + mpfr_fits_ushort_p :: Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_fits_sshort_p" + mpfr_fits_sshort_p :: Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_fits_intmax_p" + mpfr_fits_intmax_p :: Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_fits_uintmax_p" + mpfr_fits_uintmax_p :: Ptr MPFR -> CRoundMode -> IO CInt + + +------------------------------------------------------------------------------- + +-- basic arithmetic functions + +foreign import ccall unsafe "mpfr_add" + mpfr_add :: Ptr MPFR -> Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_add_ui" + mpfr_add_ui :: Ptr MPFR -> Ptr MPFR -> CULong -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_add_si" + mpfr_add_si :: Ptr MPFR -> Ptr MPFR -> CLong -> CRoundMode -> IO CInt + +-- TODO add_z, add_q + +foreign import ccall unsafe "mpfr_sub" + mpfr_sub :: Ptr MPFR -> Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_ui_sub" + mpfr_ui_sub :: Ptr MPFR -> CULong -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_sub_ui" + mpfr_sub_ui :: Ptr MPFR -> Ptr MPFR -> CULong -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_si_sub" + mpfr_si_sub :: Ptr MPFR -> CLong -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_sub_si" + mpfr_sub_si :: Ptr MPFR -> Ptr MPFR -> CLong -> CRoundMode -> IO CInt + +--TODO sub_z, sub_q + +foreign import ccall unsafe "mpfr_mul" + mpfr_mul :: Ptr MPFR -> Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_mul_ui" + mpfr_mul_ui :: Ptr MPFR -> Ptr MPFR -> CULong -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_mul_si" + mpfr_mul_si :: Ptr MPFR -> Ptr MPFR -> CLong -> CRoundMode -> IO CInt + +--TODO mul_z, mul_q + +foreign import ccall unsafe "mpfr_sqr" + mpfr_sqr :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_div" + mpfr_div :: Ptr MPFR -> Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_ui_div" + mpfr_ui_div :: Ptr MPFR -> CULong -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_div_ui" + mpfr_div_ui :: Ptr MPFR -> Ptr MPFR -> CULong -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_si_div" + mpfr_si_div :: Ptr MPFR -> CLong -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_div_si" + mpfr_div_si :: Ptr MPFR -> Ptr MPFR -> CLong -> CRoundMode -> IO CInt + +-- TODO div_z, div_q + +foreign import ccall unsafe "mpfr_sqrt" + mpfr_sqrt :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_sqrt_ui" + mpfr_sqrt_ui :: Ptr MPFR -> CULong -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_cbrt" + mpfr_cbrt :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_root" + mpfr_root :: Ptr MPFR -> Ptr MPFR -> CULong -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_pow" + mpfr_pow :: Ptr MPFR -> Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_pow_ui" + mpfr_pow_ui :: Ptr MPFR -> Ptr MPFR -> CULong -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_pow_si" + mpfr_pow_si :: Ptr MPFR -> Ptr MPFR -> CLong -> CRoundMode -> IO CInt + +--TODO pow_z + +foreign import ccall unsafe "mpfr_ui_pow_ui" + mpfr_ui_pow_ui :: Ptr MPFR -> CULong -> CULong -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_ui_pow" + mpfr_ui_pow :: Ptr MPFR -> CULong -> Ptr MPFR -> CRoundMode -> IO CInt + + +foreign import ccall unsafe "mpfr_neg" + mpfr_neg :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_abs_wrap" + mpfr_abs :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_dim" + mpfr_dim :: Ptr MPFR -> Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_mul_2ui" + mpfr_mul_2ui :: Ptr MPFR -> Ptr MPFR -> CULong -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_mul_2si" + mpfr_mul_2si :: Ptr MPFR -> Ptr MPFR -> CLong -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_div_2ui" + mpfr_div_2ui :: Ptr MPFR -> Ptr MPFR -> CULong -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_div_2si" + mpfr_div_2si :: Ptr MPFR -> Ptr MPFR -> CLong -> CRoundMode -> IO CInt + + + +-------------------------------------------------------------------------------- +-- comparison functions +-- !!!!!!!! these set erange flags +foreign import ccall unsafe "mpfr_cmp_wrap" + mpfr_cmp :: Ptr MPFR -> Ptr MPFR -> IO CInt + +foreign import ccall unsafe "mpfr_cmp_ui_wrap" + mpfr_cmp_ui :: Ptr MPFR -> CULong -> IO CInt + +foreign import ccall unsafe "mpfr_cmp_si_wrap" + mpfr_cmp_si :: Ptr MPFR -> CLong -> IO CInt + +foreign import ccall unsafe "mpfr_cmp_d" + mpfr_cmp_d :: Ptr MPFR -> CDouble -> IO CInt + +--TODO cmp_ld, cmp_z, cmp_q, cmp_f + +foreign import ccall unsafe "mpfr_cmp_ui_2exp" + mpfr_cmp_ui_2exp :: Ptr MPFR -> CULong -> Exp -> IO CInt + +foreign import ccall unsafe "mpfr_cmp_si_2exp" + mpfr_cmp_si_2exp :: Ptr MPFR -> CLong -> Exp -> IO CInt + +foreign import ccall unsafe "mpfr_cmpabs" + mpfr_cmpabs :: Ptr MPFR -> Ptr MPFR -> IO CInt + +foreign import ccall unsafe "mpfr_nan_p_wrap" + mpfr_nan_p :: Ptr MPFR -> IO CInt + +foreign import ccall unsafe "mpfr_inf_p_wrap" + mpfr_inf_p :: Ptr MPFR -> IO CInt + +foreign import ccall unsafe "mpfr_number_p" + mpfr_number_p :: Ptr MPFR -> IO CInt + +foreign import ccall unsafe "mpfr_zero_p_wrap" + mpfr_zero_p :: Ptr MPFR -> IO CInt + +foreign import ccall unsafe "mpfr_sgn_wrap" + mpfr_sgn :: Ptr MPFR -> IO CInt + +foreign import ccall unsafe "mpfr_greater_p" + mpfr_greater_p :: Ptr MPFR -> Ptr MPFR -> IO CInt + +foreign import ccall unsafe "mpfr_greaterequal_p" + mpfr_greaterequal_p :: Ptr MPFR -> Ptr MPFR -> IO CInt + +foreign import ccall unsafe "mpfr_less_p" + mpfr_less_p :: Ptr MPFR -> Ptr MPFR -> IO CInt + +foreign import ccall unsafe "mpfr_lessequal_p" + mpfr_lessequal_p :: Ptr MPFR -> Ptr MPFR -> IO CInt + +foreign import ccall unsafe "mpfr_lessgreater_p" + mpfr_lessgreater_p :: Ptr MPFR -> Ptr MPFR -> IO CInt + +foreign import ccall unsafe "mpfr_equal_p" + mpfr_equal_p :: Ptr MPFR -> Ptr MPFR -> IO CInt + +foreign import ccall unsafe "mpfr_unordered_p" + mpfr_unordered_p :: Ptr MPFR -> Ptr MPFR -> IO CInt + +-- special functions + +foreign import ccall unsafe "mpfr_log" + mpfr_log :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_log2" + mpfr_log2 :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_log10" + mpfr_log10 :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_exp" + mpfr_exp :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_exp2" + mpfr_exp2 :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_exp10" + mpfr_exp10 :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_sin" + mpfr_sin :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_cos" + mpfr_cos :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_tan" + mpfr_tan :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_sec" + mpfr_sec :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_csc" + mpfr_csc :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_cot" + mpfr_cot :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_sin_cos" + mpfr_sin_cos :: Ptr MPFR -> Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + + +foreign import ccall unsafe "mpfr_asin" + mpfr_asin :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_acos" + mpfr_acos :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_atan" + mpfr_atan :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_atan2" + mpfr_atan2 :: Ptr MPFR -> Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + + +foreign import ccall unsafe "mpfr_cosh" + mpfr_cosh :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_sinh" + mpfr_sinh :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_tanh" + mpfr_tanh :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + + +foreign import ccall unsafe "mpfr_sech" + mpfr_sech :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_csch" + mpfr_csch :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_coth" + mpfr_coth :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_asinh" + mpfr_asinh :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_acosh" + mpfr_acosh :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_atanh" + mpfr_atanh :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_fac_ui" + mpfr_fac_ui :: Ptr MPFR -> CULong -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_log1p" + mpfr_log1p :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_expm1" + mpfr_expm1 :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_eint" + mpfr_eint :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_gamma" + mpfr_gamma :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_lngamma" + mpfr_lngamma :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_lgamma" + mpfr_lgamma :: Ptr MPFR -> Ptr CInt -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_zeta" + mpfr_zeta :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_zeta_ui" + mpfr_zeta_ui :: Ptr MPFR -> CULong -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_erf" + mpfr_erf :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_erfc" + mpfr_erfc :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_j0" + mpfr_j0 :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_j1" + mpfr_j1 :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_jn" + mpfr_jn :: Ptr MPFR -> CLong -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_y0" + mpfr_y0 :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_y1" + mpfr_y1 :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_yn" + mpfr_yn :: Ptr MPFR -> CLong -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_fma" + mpfr_fma :: Ptr MPFR -> Ptr MPFR -> Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_fms" + mpfr_fms :: Ptr MPFR -> Ptr MPFR -> Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + + +foreign import ccall unsafe "mpfr_agm" + mpfr_agm :: Ptr MPFR -> Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_hypot" + mpfr_hypot :: Ptr MPFR -> Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +-- constants +foreign import ccall unsafe "mpfr_const_log2" + mpfr_const_log2 :: Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_const_pi" + mpfr_const_pi :: Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_const_euler" + mpfr_const_euler :: Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_const_catalan" + mpfr_const_catalan :: Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_free_cache" + mpfr_free_cache :: IO () + +foreign import ccall unsafe "mpfr_sum" + mpfr_sum :: Ptr MPFR -> Ptr (Ptr MPFR) -> CULong -> CRoundMode -> IO CInt + +-- TODO input and output functions + +-- integer related functions + +foreign import ccall unsafe "mpfr_rint" + mpfr_rint :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_ceil_wrap" + mpfr_ceil :: Ptr MPFR -> Ptr MPFR -> IO CInt + +foreign import ccall unsafe "mpfr_floor_wrap" + mpfr_floor :: Ptr MPFR -> Ptr MPFR -> IO CInt + +foreign import ccall unsafe "mpfr_round_wrap" + mpfr_round :: Ptr MPFR -> Ptr MPFR -> IO CInt + +foreign import ccall unsafe "mpfr_trunc_wrap" + mpfr_trunc :: Ptr MPFR -> Ptr MPFR -> IO CInt + +foreign import ccall unsafe "mpfr_rint_ceil" + mpfr_rint_ceil :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_rint_floor" + mpfr_rint_floor :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_rint_round" + mpfr_rint_round :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_rint_trunc" + mpfr_rint_trunc :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_frac" + mpfr_frac :: Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_remainder" + mpfr_remainder :: Ptr MPFR -> Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_remquo" + mpfr_remquo :: Ptr MPFR -> Ptr CLong -> Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_integer_p" + mpfr_integer_p :: Ptr MPFR -> IO CInt + +-------------------- +-- miscellaneus functions + +foreign import ccall unsafe "mpfr_nexttoward" + mpfr_nexttoward :: Ptr MPFR -> Ptr MPFR -> IO () + +foreign import ccall unsafe "mpfr_nextabove" + mpfr_nextabove :: Ptr MPFR -> IO () + +foreign import ccall unsafe "mpfr_nextbelow" + mpfr_nextbelow :: Ptr MPFR -> IO () + +foreign import ccall unsafe "mpfr_min" + mpfr_min :: Ptr MPFR -> Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_max" + mpfr_max :: Ptr MPFR -> Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +-- TODO urandomb + +foreign import ccall unsafe "mpfr_random2" + mpfr_random2 :: Ptr MPFR -> MpSize -> Exp -> IO () + + +foreign import ccall unsafe "mpfr_get_exp_wrap" + mpfr_get_exp :: Ptr MPFR -> IO Exp + +foreign import ccall unsafe "mpfr_set_exp" + mpfr_set_exp :: Ptr MPFR -> Exp -> IO CInt + +foreign import ccall unsafe "mpfr_signbit_wrap" + mpfr_signbit :: Ptr MPFR -> IO CInt + +foreign import ccall unsafe "mpfr_setsign_wrap" + mpfr_setsign :: Ptr MPFR -> Ptr MPFR -> CInt -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_copysign_wrap" + mpfr_copysign :: Ptr MPFR -> Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt + +--------------------------------------------------------------- +-- rounding mode related functions + +foreign import ccall unsafe "mpfr_get_emin" + mpfr_get_emin :: IO Exp + +foreign import ccall unsafe "mpfr_get_emax" + mpfr_get_emax :: IO Exp + +foreign import ccall unsafe "mpfr_set_emin" + mpfr_set_emin :: Exp -> IO CInt + +foreign import ccall unsafe "mpfr_set_emax" + mpfr_set_emax :: Exp -> IO CInt + +foreign import ccall unsafe "mpfr_get_emin_min" + mpfr_get_emin_min :: IO Exp + +foreign import ccall unsafe "mpfr_get_emin_max" + mpfr_get_emin_max :: IO Exp + +foreign import ccall unsafe "mpfr_get_emax_min" + mpfr_get_emax_min :: IO Exp + +foreign import ccall unsafe "mpfr_get_emax_max" + mpfr_get_emax_max :: IO Exp + +foreign import ccall unsafe "mpfr_check_range" + mpfr_check_range :: Ptr MPFR -> CInt -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_subnormalize" + mpfr_subnormalize :: Ptr MPFR -> CInt -> CRoundMode -> IO CInt + +foreign import ccall unsafe "mpfr_clear_underflow" + mpfr_clear_underflow :: IO () + +foreign import ccall unsafe "mpfr_clear_overflow" + mpfr_clear_overflow :: IO () + +foreign import ccall unsafe "mpfr_clear_nanflag" + mpfr_clear_nanflag :: IO () + +foreign import ccall unsafe "mpfr_clear_inexflag" + mpfr_clear_inexflag :: IO () + +foreign import ccall unsafe "mpfr_clear_erangeflag" + mpfr_clear_erangeflag :: IO () + +foreign import ccall unsafe "mpfr_set_underflow" + mpfr_set_underflow :: IO () + +foreign import ccall unsafe "mpfr_set_overflow" + mpfr_set_overflow :: IO () + +foreign import ccall unsafe "mpfr_set_nanflag" + mpfr_set_nanflag :: IO () + +foreign import ccall unsafe "mpfr_set_inexflag" + mpfr_set_inexflag :: IO () + +foreign import ccall unsafe "mpfr_set_erangeflag" + mpfr_set_erangeflag :: IO () + +foreign import ccall unsafe "mpfr_clear_flags" + mpfr_clear_flags :: IO () + + +foreign import ccall unsafe "mpfr_underflow_p" + mpfr_underflow_p :: IO CInt + +foreign import ccall unsafe "mpfr_overflow_p" + mpfr_overflow_p :: IO CInt + +foreign import ccall unsafe "mpfr_nanflag_p" + mpfr_nanflag_p :: IO CInt + +foreign import ccall unsafe "mpfr_inexflag_p" + mpfr_inexflag_p :: IO CInt + +foreign import ccall unsafe "mpfr_erangeflag_p" + mpfr_erangeflag_p :: IO CInt + +--------------------------------------------------------------- +-- custom interface +foreign import ccall unsafe "mpfr_custom_get_size_wrap" + mpfr_custom_get_size :: CPrecision -> IO #{type size_t} + +foreign import ccall unsafe "mpfr_custom_init_wrap" + mpfr_custom_init :: Ptr #{type mp_limb_t} -> CPrecision -> IO () + +foreign import ccall unsafe "mpfr_custom_init_set_wrap" + mpfr_custom_init_set :: Ptr MPFR -> CInt -> Exp -> CPrecision -> Ptr Limb -> IO () + +foreign import ccall unsafe "mpfr_custom_get_kind_wrap" + mpfr_custom_get_kind :: Ptr MPFR -> IO CInt + +foreign import ccall unsafe "mpfr_custom_get_mantissa_wrap" + mpfr_custom_get_mantissa :: Ptr MPFR -> IO (Ptr Limb) + +foreign import ccall unsafe "mpfr_custom_get_exp_wrap" + mpfr_custom_get_exp :: Ptr MPFR -> IO Exp + +foreign import ccall unsafe "mpfr_custom_move_wrap" + mpfr_custom_move :: Ptr MPFR -> Ptr #{type mp_limb_t} -> IO () + +------------------------------------------------- + diff --git a/Data/Number/MPFR/Integer.hs b/Data/Number/MPFR/Integer.hs new file mode 100644 --- /dev/null +++ b/Data/Number/MPFR/Integer.hs @@ -0,0 +1,106 @@ +{-| + Module : Data.Number.MPFR.Integer + Description : wrappers for integer related functions + Copyright : (c) Aleš Bizjak + License : BSD3 + + Maintainer : ales.bizjak0@gmail.com + Stability : experimental + Portability : non-portable + + Integer related functions. See MPFR manual for detailed documentation. +-} +{-# INCLUDE #-} +{-# INCLUDE #-} + + +module Data.Number.MPFR.Integer where + +import Data.Number.MPFR.Internal + +rint :: RoundMode -> Precision -> MPFR -> MPFR +rint r p d = fst $ rint_ r p d + +ceil :: Precision -> MPFR -> MPFR +ceil p d = fst $ ceil_ p d + +floor :: Precision -> MPFR -> MPFR +floor p d = fst $ floor_ p d + +round :: Precision -> MPFR -> MPFR +round p d = fst $ round_ p d + +trunc :: Precision -> MPFR -> MPFR +trunc p d = fst $ trunc_ p d + +rintCeil :: RoundMode -> Precision -> MPFR -> MPFR +rintCeil r p d = fst $ rintCeil_ r p d + +rintFloor :: RoundMode -> Precision -> MPFR -> MPFR +rintFloor r p d = fst $ rintFloor_ r p d + +rintRound :: RoundMode -> Precision -> MPFR -> MPFR +rintRound r p d = fst $ rintRound_ r p d + +rintTrunc :: RoundMode -> Precision -> MPFR -> MPFR +rintTrunc r p d = fst $ rintTrunc_ r p d + +frac :: RoundMode -> Precision -> MPFR -> MPFR +frac r p d = fst $ frac_ r p d + +remainder :: RoundMode -> Precision -> MPFR -> MPFR -> MPFR +remainder r p d d' = fst $ remainder_ r p d d' + +remquo :: RoundMode -> Precision -> MPFR -> MPFR -> (MPFR, Int) +remquo r p d d' = case remquo_ r p d d' of + (a, b, _) -> (a, b) + +rint_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +rint_ r p d = withMPFR r p d mpfr_rint + +ceil_ :: Precision -> MPFR -> (MPFR, Int) +ceil_ p d = withMPFRR p d mpfr_ceil + +floor_ :: Precision -> MPFR -> (MPFR, Int) +floor_ p d = withMPFRR p d mpfr_floor + +round_ :: Precision -> MPFR -> (MPFR, Int) +round_ p d = withMPFRR p d mpfr_round + +trunc_ :: Precision -> MPFR -> (MPFR, Int) +trunc_ p d = withMPFRR p d mpfr_trunc + +rintCeil_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +rintCeil_ r p d = withMPFR r p d mpfr_rint_ceil + +rintFloor_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +rintFloor_ r p d = withMPFR r p d mpfr_rint_floor + +rintRound_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +rintRound_ r p d = withMPFR r p d mpfr_rint_round + +rintTrunc_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +rintTrunc_ r p d = withMPFR r p d mpfr_rint_trunc + +frac_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +frac_ r p d = withMPFR r p d mpfr_frac + +remainder_ :: RoundMode -> Precision -> MPFR -> MPFR -> (MPFR,Int) +remainder_ r p d d' = withMPFRsBA r p d d' mpfr_remainder + +remquo_ :: RoundMode -> Precision -> MPFR -> MPFR -> (MPFR, Int, Int) +remquo_ r p d 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 + with d $ \p2 -> do + with d' $ \p3 -> do + alloca $ \p4 -> do + r3 <- mpfr_remquo p1 p4 p2 p3 ((fromIntegral . fromEnum) r) + r1 <- peekP p1 fp + r2 <- peek p4 + return (r1, fromIntegral r2, fromIntegral r3) + +isInteger :: MPFR -> Bool +isInteger d = withMPFRB d mpfr_integer_p /= 0 diff --git a/Data/Number/MPFR/Internal.hs b/Data/Number/MPFR/Internal.hs new file mode 100644 --- /dev/null +++ b/Data/Number/MPFR/Internal.hs @@ -0,0 +1,204 @@ +{-# INCLUDE #-} +{-# INCLUDE #-} +{-# LANGUAGE BangPatterns #-} + +module Data.Number.MPFR.Internal ( + module FFIhelper, + withMPFRsBA, withMPFRBAui, withMPFRBAiu, withMPFRBAsi, withMPFRBAis, + withMPFRB, withMPFRP, withMPFR, withMPFRBB, withMPFRC, + withMPFRF, withMPFRUI, withMPFRR, checkPrec, getMantissa', binprec, + unsafePerformIO, peek, Ptr, nullPtr, mallocForeignPtrBytes, with, + withForeignPtr, CInt, CLong, CULong, withCString, peekCString, alloca, + peekArray, shiftL, Word, minPrec, + + Precision +) +where + +import Data.Number.MPFR.FFIhelper as FFIhelper + +import Foreign.C(CInt, CLong, CULong, withCString, peekCString) +import Foreign.Marshal(alloca, peekArray) +import Foreign(unsafePerformIO, peek, Ptr, nullPtr, mallocForeignPtrBytes, with, withForeignPtr) + +import Data.Bits(shiftL) + +import Data.Word(Word) + +type Precision = Word + + +-- these are helper functions, only for internal use +{-# INLINE withMPFRsBA #-} +withMPFRsBA :: RoundMode -> Precision -> MPFR -> MPFR + -> (Ptr MPFR -> Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt) + -> (MPFR, Int) +withMPFRsBA 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 withMPFRBAui #-} +withMPFRBAui :: RoundMode -> Precision -> MPFR -> CULong + -> (Ptr MPFR -> Ptr MPFR -> CULong -> CRoundMode -> IO CInt) + -> (MPFR, Int) +withMPFRBAui 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 withMPFRBAsi #-} +withMPFRBAsi :: RoundMode -> Precision -> MPFR -> CLong + -> (Ptr MPFR -> Ptr MPFR -> CLong -> CRoundMode -> IO CInt) + -> (MPFR, Int) +withMPFRBAsi 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 withMPFRBAiu #-} +withMPFRBAiu :: RoundMode -> Precision -> CULong -> MPFR + -> (Ptr MPFR -> CULong -> Ptr MPFR -> CRoundMode -> IO CInt) + -> (MPFR, Int) +withMPFRBAiu 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 withMPFRBAis #-} +withMPFRBAis :: RoundMode -> Precision -> CLong -> MPFR + -> (Ptr MPFR -> CLong -> Ptr MPFR -> CRoundMode -> IO CInt) + -> (MPFR, Int) +withMPFRBAis 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 withMPFRB #-} +withMPFRB :: MPFR -> (Ptr MPFR -> IO CInt) -> CInt +withMPFRB mp1 !f = unsafePerformIO go + where go = with mp1 $! \p1 -> f p1 + +withMPFRP :: MPFR -> (Ptr MPFR -> IO CPrecision) -> CPrecision +withMPFRP mp1 !f = unsafePerformIO go + where go = with mp1 $! \p1 -> f p1 + +{-# INLINE withMPFR #-} +withMPFR :: RoundMode -> Precision -> MPFR + -> (Ptr MPFR -> Ptr MPFR -> CRoundMode -> IO CInt) + -> (MPFR, Int) +withMPFR 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 withMPFRBB #-} +withMPFRBB :: MPFR -> MPFR + -> (Ptr MPFR -> Ptr MPFR -> IO CInt) + -> CInt +withMPFRBB mp1 mp2 f = unsafePerformIO go + where go = do with mp1 $! \p1 -> do + with mp2 $! \p2 -> do + f p1 p2 + +{-# INLINE withMPFRC #-} +withMPFRC :: RoundMode -> Precision -> + (Ptr MPFR -> CRoundMode -> IO CInt) -> (MPFR, Int) +withMPFRC 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) + +withMPFRF :: MPFR -> RoundMode + -> (Ptr MPFR -> CRoundMode -> IO CInt) + -> Int +withMPFRF !mp1 r f = (fromIntegral . unsafePerformIO) go + where go = do with mp1 $! \p1 -> f p1 ((fromIntegral . fromEnum) r) + +{-# INLINE withMPFRUI #-} +withMPFRUI :: RoundMode -> Precision -> Word + -> (Ptr MPFR -> CULong -> CRoundMode -> IO CInt) + -> (MPFR, Int) +withMPFRUI r p 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 + r2 <- f p1 (fromIntegral d) ((fromIntegral . fromEnum) r) + r1 <- peekP p1 fp + return (r1, fromIntegral r2) +{-# INLINE withMPFRR #-} +withMPFRR :: Precision -> MPFR + -> (Ptr MPFR -> Ptr MPFR -> IO CInt) + -> (MPFR, Int) +withMPFRR p 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 d $! \p2 -> do + r2 <- f p1 p2 + r1 <- peekP p1 fp + return (r1, fromIntegral r2) + + +checkPrec :: Precision -> Precision +checkPrec = max minPrec + +getMantissa' :: MPFR -> [Limb] +getMantissa' (MP p _ _ p1) = unsafePerformIO go + where go = do withForeignPtr p1 $! \pt -> do + arr <- peekArray (Prelude.ceiling ((fromIntegral p ::Double) / fromIntegral bitsPerMPLimb)) pt ; + return arr + +{- TODO: this is inefficient +binprec :: Integer -> Precision +binprec i = length (takeWhile (/= 0) (iterate (flip shiftR 1) i) +-} +-- TODO +binprec :: Integer -> Precision +binprec d = Prelude.floor (logBase 2 (fromInteger (if d >= 0 then d else -d)) :: Double) + 1 + +--one :: MPFR +--one = fromWord Near minPrec 1 + +--zero :: MPFR +--ggzero = fromWord Near minPrec 0 + +minPrec :: Precision +minPrec = 32 diff --git a/Data/Number/MPFR/Misc.hs b/Data/Number/MPFR/Misc.hs new file mode 100644 --- /dev/null +++ b/Data/Number/MPFR/Misc.hs @@ -0,0 +1,117 @@ +{-| + Module : Data.Number.MPFR.Misc + Description : wrappers for miscellaneous functions + Copyright : (c) Aleš Bizjak + License : BSD3 + + Maintainer : ales.bizjak0@gmail.com + Stability : experimental + Portability : non-portable + + Functions that don't belong anywhere else. See MPFR manual for detailed documentation. +-} +{-# INCLUDE #-} +{-# INCLUDE #-} + +module Data.Number.MPFR.Misc where + +import Data.Number.MPFR.Internal + +import Data.Number.MPFR.Assignment + + +nextToward :: MPFR -> MPFR -> MPFR +nextToward mp1 mp2 = 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 + with mp2 $ \p3 -> do + _ <- mpfr_set p1 p2 ((fromIntegral . fromEnum) Near) + mpfr_nexttoward p1 p3 + peekP p1 fp + + +nextAbove :: MPFR -> MPFR +nextAbove 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_nextabove p1 + peekP p1 fp + +nextBelow :: MPFR -> MPFR +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 + +maxD :: RoundMode -> Precision -> MPFR -> MPFR -> MPFR +maxD r p d1 d2 = fst $ maxD_ r p d1 d2 + +minD :: RoundMode -> Precision -> MPFR -> MPFR -> MPFR +minD r p d1 d2 = fst $ minD_ r p d1 d2 + +random2 :: Precision -> MpSize -> Exp -> IO MPFR +random2 p m e = do ls <- mpfr_custom_get_size (fromIntegral p) + fp <- mallocForeignPtrBytes (fromIntegral ls) + let dummy = MP (fromIntegral p) 0 0 fp + with dummy $ \p1 -> do + mpfr_random2 p1 m e + peekP p1 fp + +getExp :: MPFR -> Exp +getExp d = (fromIntegral . unsafePerformIO) go + where go = do with d $ \p1 -> mpfr_custom_get_exp p1 + +setExp :: MPFR -> Exp -> MPFR +setExp d e = unsafePerformIO go + where go = do let p = fromIntegral (getPrec d) + ls <- mpfr_custom_get_size p + fp <- mallocForeignPtrBytes (fromIntegral ls) + let dummy = MP p 0 0 fp + with dummy $ \p1 -> do + with d $ \p2 -> do + mpfr_set p1 p2 ((fromIntegral . fromEnum) Near) + mpfr_set_exp p1 e + peekP p1 fp + +signbit :: MPFR -> Bool +signbit d = withMPFRB d mpfr_signbit /= 0 + +maxD_ :: RoundMode -> Precision -> MPFR -> MPFR -> (MPFR, Int) +maxD_ r p d1 d2 = withMPFRsBA r p d1 d2 mpfr_max + +minD_ :: RoundMode -> Precision -> MPFR -> MPFR -> (MPFR, Int) +minD_ r p d1 d2 = withMPFRsBA r p d1 d2 mpfr_min + +getPrec :: MPFR -> Precision +getPrec d = fromIntegral (withMPFRP d mpfr_get_prec) + +-- | getMantissa and getExp return values such that +-- +-- > d = getMantissa d * 2^(getExp d - Prelude.ceiling ((getPrec d) / bitsPerMPLimb)* bitsPerMPLimb ) +getMantissa :: MPFR -> Integer +getMantissa d = toInteger (sign d) * h + where (h, _) = foldl (\(a,b) c -> (a + (toInteger c) `shiftL` b, b + bitsPerMPLimb)) (0,0) (getMantissa' d) + +one :: MPFR +one = fromWord Near minPrec 1 + +zero :: MPFR +zero = fromWord Near minPrec 0 + +maxPrec :: MPFR -> MPFR -> Precision +maxPrec d d' = max (getPrec d) (getPrec d') diff --git a/Data/Number/MPFR/Near.hs b/Data/Number/MPFR/Near.hs new file mode 100644 --- /dev/null +++ b/Data/Number/MPFR/Near.hs @@ -0,0 +1,82 @@ +{-| + Module : Data.Number.MPFR.Near + Description : top level + Copyright : (c) Aleš Bizjak + License : BSD3 + + Maintainer : ales.bizjak0@gmail.com + Stability : experimental + Portability : non-portable + + This module defines instances Num, Real, Fractional, Floating and RealFrac of MPFR. + Operations are rounded with RoundMode Near and computed with max precision of two + operands or with the precision of the operand. Otherwise it is equivalent to + Data.Number.MPFR + +-} +{-# INCLUDE #-} +{-# INCLUDE #-} + +module Data.Number.MPFR.Near ( + module Data.Number.MPFR.Base +) +where + +import Data.Number.MPFR.Base + +import Data.Number.MPFR.Internal + +import Data.Maybe + +import Data.Ratio + +instance Num MPFR where + d + d' = add Near (maxPrec d d') d d' + d - d' = sub Near (maxPrec d d') d d' + d * d' = mul Near (maxPrec d d') d d' + negate d = neg Near (getPrec d) d + abs d = absD Near (getPrec d) d + signum d = fromInt Near minPrec (fromMaybe (-1) (sgn d)) + fromInteger i = fromIntegerA Zero (checkPrec $ binprec i) i + -- TODO works only partially + +instance Real MPFR where + toRational d = n % 2 ^ e + where (n', e') = decompose d + (n, e) = if e' >= 0 then ((n' * 2 ^ e'), 0) + else (n', - e') + +instance Fractional MPFR where + d / d' = Data.Number.MPFR.Base.div Up (maxPrec d d') d d' + fromRational r = (fromInteger n) / (fromInteger d) + where n = numerator r + d = denominator r + recip d = one / d + +instance Floating MPFR where + pi = Data.Number.MPFR.Base.pi Near 53 + exp d = Data.Number.MPFR.Base.exp Near (getPrec d) d + log d = Data.Number.MPFR.Base.log Near (getPrec d) d + sqrt d = Data.Number.MPFR.Base.sqrt Near (getPrec d) d + (**) d d' = Data.Number.MPFR.Base.pow Near (maxPrec d d') d d' + logBase d d' = Prelude.log d' / Prelude.log d + sin d = Data.Number.MPFR.Base.sin Near (getPrec d) d + cos d = Data.Number.MPFR.Base.cos Near (getPrec d) d + tan d = Data.Number.MPFR.Base.tan Near (getPrec d) d + asin d = Data.Number.MPFR.Base.asin Near (getPrec d) d + acos d = Data.Number.MPFR.Base.acos Near (getPrec d) d + atan d = Data.Number.MPFR.Base.atan Near (getPrec d) d + sinh d = Data.Number.MPFR.Base.sinh Near (getPrec d) d + cosh d = Data.Number.MPFR.Base.cosh Near (getPrec d) d + tanh d = Data.Number.MPFR.Base.tanh Near (getPrec d) d + asinh d = Data.Number.MPFR.Base.asinh Near (getPrec d) d + acosh d = Data.Number.MPFR.Base.acosh Near (getPrec d) d + atanh d = Data.Number.MPFR.Base.atanh Near (getPrec d) d + +instance RealFrac MPFR where + properFraction d = (fromIntegral n, f) + where r = toRational d + m = numerator r + e = denominator r + n = quot m e + f = frac Near (getPrec d) d diff --git a/Data/Number/MPFR/Special.hs b/Data/Number/MPFR/Special.hs new file mode 100644 --- /dev/null +++ b/Data/Number/MPFR/Special.hs @@ -0,0 +1,382 @@ +{-| + Module : Data.Number.MPFR.Special + Description : wrappers for special functions + Copyright : (c) Aleš Bizjak + License : BSD3 + + Maintainer : ales.bizjak0@gmail.com + Stability : experimental + Portability : non-portable + + Special functions. See MPFR manual for detailed documentation. +-} +{-# INCLUDE #-} +{-# INCLUDE #-} + +module Data.Number.MPFR.Special where + +import Data.Number.MPFR.Internal + +log :: RoundMode -> Precision -> MPFR -> MPFR +log r p d = fst $ log_ r p d + +log2 :: RoundMode -> Precision -> MPFR -> MPFR +log2 r p d = fst $ log2_ r p d + +log10 :: RoundMode -> Precision -> MPFR -> MPFR +log10 r p d = fst $ log10_ r p d + +exp :: RoundMode -> Precision -> MPFR -> MPFR +exp r p d = fst $ exp_ r p d + +exp2 :: RoundMode -> Precision -> MPFR -> MPFR +exp2 r p d = fst $ exp2_ r p d + +exp10 :: RoundMode -> Precision -> MPFR -> MPFR +exp10 r p d = fst $ exp10_ r p d + +sin :: RoundMode -> Precision -> MPFR -> MPFR +sin r p d = fst $ sin_ r p d + +cos :: RoundMode -> Precision -> MPFR -> MPFR +cos r p d = fst $ cos_ r p d + +tan :: RoundMode -> Precision -> MPFR -> MPFR +tan r p d = fst $ tan_ r p d + +sec :: RoundMode -> Precision -> MPFR -> MPFR +sec r p d = fst $ sec_ r p d + +csc :: RoundMode -> Precision -> MPFR -> MPFR +csc r p d = fst $ csc_ r p d + +cot :: RoundMode -> Precision -> MPFR -> MPFR +cot r p d = fst $ cot_ r p d + +sincos :: RoundMode + -> Precision -- ^ precision to compute sin + -> Precision -- ^ precision to compute cos + -> MPFR + -> (MPFR, MPFR) -- ^ return (sin x, cos x) +sincos r p p' d = case sincos_ r p p' d of + (a, b, _) -> (a, b) + +asin :: RoundMode -> Precision -> MPFR -> MPFR +asin r p d = fst $ asin_ r p d + +acos :: RoundMode -> Precision -> MPFR -> MPFR +acos r p d = fst $ acos_ r p d + +atan :: RoundMode -> Precision -> MPFR -> MPFR +atan r p d = fst $ atan_ r p d + +atan2 :: RoundMode -> Precision -> MPFR -> MPFR -> MPFR +atan2 r p d d' = fst $ atan2_ r p d d' + +sinh :: RoundMode -> Precision -> MPFR -> MPFR +sinh r p d = fst $ sinh_ r p d + +cosh :: RoundMode -> Precision -> MPFR -> MPFR +cosh r p d = fst $ cosh_ r p d + +tanh :: RoundMode -> Precision -> MPFR -> MPFR +tanh r p d = fst $ tanh_ r p d + +sech :: RoundMode -> Precision -> MPFR -> MPFR +sech r p d = fst $ sech_ r p d + +csch :: RoundMode -> Precision -> MPFR -> MPFR +csch r p d = fst $ csch_ r p d + +coth :: RoundMode -> Precision -> MPFR -> MPFR +coth r p d = fst $ coth_ r p d + +acosh :: RoundMode -> Precision -> MPFR -> MPFR +acosh r p d = fst $ acosh_ r p d + +asinh :: RoundMode -> Precision -> MPFR -> MPFR +asinh r p d = fst $ asinh_ r p d + +atanh :: RoundMode -> Precision -> MPFR -> MPFR +atanh r p d = fst $ atanh_ r p d + +facw :: RoundMode -> Precision -> Word -> MPFR +facw r p d = fst $ facw_ r p d + +log1p :: RoundMode -> Precision -> MPFR -> MPFR +log1p r p d = fst $ log1p_ r p d + +expm1 :: RoundMode -> Precision -> MPFR -> MPFR +expm1 r p d = fst $ expm1_ r p d + +eint :: RoundMode -> Precision -> MPFR -> MPFR +eint r p d = fst $ eint_ r p d + +gamma :: RoundMode -> Precision -> MPFR -> MPFR +gamma r p d = fst $ gamma_ r p d + +lngamma :: RoundMode -> Precision -> MPFR -> MPFR +lngamma r p d = fst $ lngamma_ r p d + +lgamma :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +lgamma r p d = case lgamma_ r p d of + (a, b, _) -> (a,b) + +zeta :: RoundMode -> Precision -> MPFR -> MPFR +zeta r p d = fst $ zeta_ r p d + +zetaw :: RoundMode -> Precision -> Word -> MPFR +zetaw r p d = fst $ zetaw_ r p d + +erf :: RoundMode -> Precision -> MPFR -> MPFR +erf r p d = fst $ erf_ r p d + +erfc :: RoundMode -> Precision -> MPFR -> MPFR +erfc r p d = fst $ erfc_ r p d + +j0 :: RoundMode -> Precision -> MPFR -> MPFR +j0 r p d = fst $ j0_ r p d + +j1 :: RoundMode -> Precision -> MPFR -> MPFR +j1 r p d = fst $ j1_ r p d + +jn :: RoundMode -> Precision -> Int -> MPFR -> MPFR +jn r p w d = fst $ jn_ r p w d + +y0 :: RoundMode -> Precision -> MPFR -> MPFR +y0 r p d = fst $ y0_ r p d + +y1 :: RoundMode -> Precision -> MPFR -> MPFR +y1 r p d = fst $ y1_ r p d + +yn :: RoundMode -> Precision -> Int -> MPFR -> MPFR +yn r p w d = fst $ yn_ r p w d + +fma :: RoundMode -> Precision -> MPFR -> MPFR -> MPFR -> MPFR +fma r p d1 d2 d3 = fst $ fma_ r p d1 d2 d3 + +fms :: RoundMode -> Precision -> MPFR -> MPFR -> MPFR -> MPFR +fms r p d1 d2 d3 = fst $ fms_ r p d1 d2 d3 + +agm :: RoundMode -> Precision -> MPFR -> MPFR -> MPFR +agm r p d1 d2 = fst $ agm_ r p d1 d2 + +hypot :: RoundMode -> Precision -> MPFR -> MPFR -> MPFR +hypot r p d1 d2 = fst $ hypot_ r p d1 d2 + +pi :: RoundMode -> Precision -> MPFR +pi r p = fst $ pi_ r p + +log2c :: RoundMode -> Precision -> MPFR +log2c r p = fst $ pi_ r p + +euler :: RoundMode -> Precision -> MPFR +euler r p = fst $ pi_ r p + +catalan :: RoundMode -> Precision -> MPFR +catalan r p = fst $ pi_ r p + + +log_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +log_ r p d = withMPFR r p d mpfr_log + +log2_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +log2_ r p d = withMPFR r p d mpfr_log2 + +log10_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +log10_ r p d = withMPFR r p d mpfr_log10 + +exp_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +exp_ r p d = withMPFR r p d mpfr_exp + +exp2_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +exp2_ r p d = withMPFR r p d mpfr_exp2 + +exp10_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +exp10_ r p d = withMPFR r p d mpfr_exp10 + +sin_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +sin_ r p d = withMPFR r p d mpfr_sin + +cos_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +cos_ r p d = withMPFR r p d mpfr_cos + +tan_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +tan_ r p d = withMPFR r p d mpfr_tan + +sec_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +sec_ r p d = withMPFR r p d mpfr_sec + +csc_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +csc_ r p d = withMPFR r p d mpfr_csc + +cot_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +cot_ r p d = withMPFR r p d mpfr_cot + + +sincos_ :: RoundMode + -> Precision -- ^ precision to compute sin + -> Precision -- ^ precision to compute cos + -> MPFR + -> (MPFR, MPFR, Int) +sincos_ r p 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 + ls' <- mpfr_custom_get_size (fromIntegral p') + fp' <- mallocForeignPtrBytes (fromIntegral ls') + let dummy' = MP (fromIntegral p') 0 0 fp' + with dummy $ \p1 -> do + with dummy' $ \p2 -> do + with d $ \p3 -> do + r3 <- mpfr_sin_cos p1 p2 p3 ((fromIntegral . fromEnum) r) + r1 <- peekP p1 fp + r2 <- peekP p2 fp' + return (r1, r2, fromIntegral r3) + +asin_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +asin_ r p d = withMPFR r p d mpfr_asin + +acos_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +acos_ r p d = withMPFR r p d mpfr_acos + +atan_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +atan_ r p d = withMPFR r p d mpfr_atan + +atan2_ :: RoundMode -> Precision -> MPFR -> MPFR -> (MPFR, Int) +atan2_ r p d d' = withMPFRsBA r p d d' mpfr_atan2 + +sinh_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +sinh_ r p d = withMPFR r p d mpfr_sinh + +cosh_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +cosh_ r p d = withMPFR r p d mpfr_cosh + +tanh_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +tanh_ r p d = withMPFR r p d mpfr_tanh + +sech_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +sech_ r p d = withMPFR r p d mpfr_sech + +csch_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +csch_ r p d = withMPFR r p d mpfr_csch + +coth_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +coth_ r p d = withMPFR r p d mpfr_coth + +acosh_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +acosh_ r p d = withMPFR r p d mpfr_acosh + +asinh_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +asinh_ r p d = withMPFR r p d mpfr_asinh + +atanh_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +atanh_ r p d = withMPFR r p d mpfr_atanh + +facw_ :: RoundMode -> Precision -> Word -> (MPFR, Int) +facw_ r p w = withMPFRUI r p w mpfr_fac_ui + +log1p_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +log1p_ r p d = withMPFR r p d mpfr_log1p + +expm1_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +expm1_ r p d = withMPFR r p d mpfr_expm1 + +eint_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +eint_ r p d = withMPFR r p d mpfr_eint + +gamma_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +gamma_ r p d = withMPFR r p d mpfr_gamma + +lngamma_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +lngamma_ r p d = withMPFR r p d mpfr_lngamma + +lgamma_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int, Int) +lgamma_ 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 + with d $ \p2 -> do + alloca $ \p3 -> do + r3 <- mpfr_lgamma p1 p3 p2 ((fromIntegral . fromEnum) r) + r2 <- peek p3 + r1 <- peekP p1 fp + return (r1, fromIntegral r2, fromIntegral r3) + +zeta_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +zeta_ r p d = withMPFR r p d mpfr_zeta + +zetaw_ :: RoundMode -> Precision -> Word -> (MPFR, Int) +zetaw_ r p d = withMPFRUI r p d mpfr_zeta_ui + +erf_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +erf_ r p d = withMPFR r p d mpfr_erf + +erfc_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +erfc_ r p d = withMPFR r p d mpfr_erfc + +j0_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +j0_ r p d = withMPFR r p d mpfr_j0 + +j1_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +j1_ r p d = withMPFR r p d mpfr_j1 + +jn_ :: RoundMode -> Precision -> Int -> MPFR -> (MPFR, Int) +jn_ r p i d = withMPFRBAis r p (fromIntegral i) d mpfr_jn + +y0_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +y0_ r p d = withMPFR r p d mpfr_y0 + +y1_ :: RoundMode -> Precision -> MPFR -> (MPFR, Int) +y1_ r p d = withMPFR r p d mpfr_y1 + +yn_ :: RoundMode -> Precision -> Int -> MPFR -> (MPFR, Int) +yn_ r p i d = withMPFRBAis r p (fromIntegral i) d mpfr_yn + +fma_ :: RoundMode -> Precision -> MPFR -> MPFR -> MPFR -> (MPFR, 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 -> MPFR -> MPFR -> MPFR -> (MPFR, 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) + +agm_ :: RoundMode -> Precision -> MPFR -> MPFR -> (MPFR,Int) +agm_ r p d1 d2 = withMPFRsBA r p d1 d2 mpfr_agm + +hypot_ :: RoundMode -> Precision -> MPFR -> MPFR -> (MPFR,Int) +hypot_ r p d1 d2 = withMPFRsBA r p d1 d2 mpfr_hypot + +pi_ :: RoundMode -> Precision -> (MPFR, Int) +pi_ r p = withMPFRC r p mpfr_const_pi + +log2c_ :: RoundMode -> Precision -> (MPFR, Int) +log2c_ r p = withMPFRC r p mpfr_const_log2 + +euler_ :: RoundMode -> Precision -> (MPFR, Int) +euler_ r p = withMPFRC r p mpfr_const_euler + +catalan_ :: RoundMode -> Precision -> (MPFR, Int) +catalan_ r p = withMPFRC r p mpfr_const_catalan + +freeCache :: IO () +freeCache = mpfr_free_cache diff --git a/Data/Number/MPFR/Up.hs b/Data/Number/MPFR/Up.hs new file mode 100644 --- /dev/null +++ b/Data/Number/MPFR/Up.hs @@ -0,0 +1,83 @@ +{-| + Module : Data.Number.MPFR.Up + Description : top level + Copyright : (c) Aleš Bizjak + License : BSD3 + + Maintainer : ales.bizjak0@gmail.com + Stability : experimental + Portability : non-portable + + This module defines instances Num, Real, Fractional, Floating and RealFrac of MPFR. + Operations are rounded with RoundMode Up and computed with max precision of two + operands or with the precision of the operand. Otherwise it is equivalent to + Data.Number.MPFR +-} + +{-# INCLUDE #-} +{-# INCLUDE #-} + + +module Data.Number.MPFR.Up ( + module Data.Number.MPFR.Base +) +where + +import Data.Number.MPFR.Base + +import Data.Number.MPFR.Internal + +import Data.Maybe + +import Data.Ratio + +instance Num MPFR where + d + d' = add Up (maxPrec d d') d d' + d - d' = sub Up (maxPrec d d') d d' + d * d' = mul Up (maxPrec d d') d d' + negate d = neg Up (getPrec d) d + abs d = absD Up (getPrec d) d + signum d = fromInt Up minPrec (fromMaybe (-1) (sgn d)) + fromInteger i = fromIntegerA Zero (checkPrec $ binprec i) i + -- TODO works only partially + +instance Real MPFR where + toRational d = n % 2 ^ e + where (n', e') = decompose d + (n, e) = if e' >= 0 then ((n' * 2 ^ e'), 0) + else (n', - e') + +instance Fractional MPFR where + d / d' = Data.Number.MPFR.Base.div Up (maxPrec d d') d d' + fromRational r = (fromInteger n) / (fromInteger d) + where n = numerator r + d = denominator r + recip d = one / d + +instance Floating MPFR where + pi = Data.Number.MPFR.Base.pi Up 53 + exp d = Data.Number.MPFR.Base.exp Up (getPrec d) d + log d = Data.Number.MPFR.Base.log Up (getPrec d) d + sqrt d = Data.Number.MPFR.Base.sqrt Up (getPrec d) d + (**) d d' = Data.Number.MPFR.Base.pow Up (maxPrec d d') d d' + logBase d d' = Prelude.log d' / Prelude.log d + sin d = Data.Number.MPFR.Base.sin Up (getPrec d) d + cos d = Data.Number.MPFR.Base.cos Up (getPrec d) d + tan d = Data.Number.MPFR.Base.tan Up (getPrec d) d + asin d = Data.Number.MPFR.Base.asin Up (getPrec d) d + acos d = Data.Number.MPFR.Base.acos Up (getPrec d) d + atan d = Data.Number.MPFR.Base.atan Up (getPrec d) d + sinh d = Data.Number.MPFR.Base.sinh Up (getPrec d) d + cosh d = Data.Number.MPFR.Base.cosh Up (getPrec d) d + tanh d = Data.Number.MPFR.Base.tanh Up (getPrec d) d + asinh d = Data.Number.MPFR.Base.asinh Up (getPrec d) d + acosh d = Data.Number.MPFR.Base.acosh Up (getPrec d) d + atanh d = Data.Number.MPFR.Base.atanh Up (getPrec d) d + +instance RealFrac MPFR where + properFraction d = (fromIntegral n, f) + where r = toRational d + m = numerator r + e = denominator r + n = quot m e + f = frac Up (getPrec d) d diff --git a/Data/Number/MPFR/Zero.hs b/Data/Number/MPFR/Zero.hs new file mode 100644 --- /dev/null +++ b/Data/Number/MPFR/Zero.hs @@ -0,0 +1,83 @@ +{-| + Module : Data.Number.MPFR.Zero + Description : top level + Copyright : (c) Aleš Bizjak + License : BSD3 + + Maintainer : ales.bizjak0@gmail.com + Stability : experimental + Portability : non-portable + + This module defines instances Num, Real, Fractional, Floating and RealFrac of MPFR. + Operations are rounded with RoundMode Zero and computed with max precision of two + operands or with the precision of the operand. Otherwise it is equivalent to + Data.Number.MPFR + +-} + +{-# INCLUDE #-} +{-# INCLUDE #-} + +module Data.Number.MPFR.Zero ( + module Data.Number.MPFR.Base +) +where + +import Data.Number.MPFR.Base + +import Data.Number.MPFR.Internal + +import Data.Maybe + +import Data.Ratio + +instance Num MPFR where + d + d' = add Zero (maxPrec d d') d d' + d - d' = sub Zero (maxPrec d d') d d' + d * d' = mul Zero (maxPrec d d') d d' + negate d = neg Zero (getPrec d) d + abs d = absD Zero (getPrec d) d + signum d = fromInt Zero minPrec (fromMaybe (-1) (sgn d)) + fromInteger i = fromIntegerA Zero (checkPrec $ binprec i) i + -- TODO works only partially + +instance Real MPFR where + toRational d = n % 2 ^ e + where (n', e') = decompose d + (n, e) = if e' >= 0 then ((n' * 2 ^ e'), 0) + else (n', - e') + +instance Fractional MPFR where + d / d' = Data.Number.MPFR.Base.div Zero (maxPrec d d') d d' + fromRational r = (fromInteger n) / (fromInteger d) + where n = numerator r + d = denominator r + recip d = one / d + +instance Floating MPFR where + pi = Data.Number.MPFR.Base.pi Zero 53 + exp d = Data.Number.MPFR.Base.exp Zero (getPrec d) d + log d = Data.Number.MPFR.Base.log Zero (getPrec d) d + sqrt d = Data.Number.MPFR.Base.sqrt Zero (getPrec d) d + (**) d d' = Data.Number.MPFR.Base.pow Zero (maxPrec d d') d d' + logBase d d' = Prelude.log d' / Prelude.log d + sin d = Data.Number.MPFR.Base.sin Zero (getPrec d) d + cos d = Data.Number.MPFR.Base.cos Zero (getPrec d) d + tan d = Data.Number.MPFR.Base.tan Zero (getPrec d) d + asin d = Data.Number.MPFR.Base.asin Zero (getPrec d) d + acos d = Data.Number.MPFR.Base.acos Zero (getPrec d) d + atan d = Data.Number.MPFR.Base.atan Zero (getPrec d) d + sinh d = Data.Number.MPFR.Base.sinh Zero (getPrec d) d + cosh d = Data.Number.MPFR.Base.cosh Zero (getPrec d) d + tanh d = Data.Number.MPFR.Base.tanh Zero (getPrec d) d + asinh d = Data.Number.MPFR.Base.asinh Zero (getPrec d) d + acosh d = Data.Number.MPFR.Base.acosh Zero (getPrec d) d + atanh d = Data.Number.MPFR.Base.atanh Zero (getPrec d) d + +instance RealFrac MPFR where + properFraction d = (fromIntegral n, f) + where r = toRational d + m = numerator r + e = denominator r + n = quot m e + f = frac Zero (getPrec d) d diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,31 @@ +Copyright Aleš Bizjak 2008 + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Aleš Bizjak nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + diff --git a/README b/README new file mode 100644 --- /dev/null +++ b/README @@ -0,0 +1,7 @@ +Haskell MPFR binding. + +It was successfully tested with MPFR versions 2.3.1 and 2.3.2 and with +GHC versions 6.8.3 and 6.8.2 on (Ubuntu) Linux and on Windows (MPFR was built with +MSYS). + +If you find bugs or have any questions please write to . diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,3 @@ +#!/usr/bin/env runhaskell +> import Distribution.Simple +> main = defaultMain diff --git a/cbits/chsmpfr.c b/cbits/chsmpfr.c new file mode 100644 --- /dev/null +++ b/cbits/chsmpfr.c @@ -0,0 +1,131 @@ +#include "chsmpfr.h" + + +mpfr_ptr initS (const mp_prec_t prec) { + mpfr_ptr rVal = malloc (sizeof(__mpfr_struct)); + mp_limb_t *limb = (mp_limb_t*)malloc (mpfr_custom_get_size(prec)); + mpfr_custom_init(limb, prec); + mpfr_custom_init_set(rVal, MPFR_NAN_KIND, 0, prec, limb); + return rVal; +} + +void clear (const mpfr_ptr p) { + free (p->_mpfr_d); + free (p); +} + +int mpfr_nan_p_wrap(const mpfr_ptr p) { + return mpfr_nan_p(p); +} + +int mpfr_inf_p_wrap(const mpfr_ptr p) { + return mpfr_inf_p(p); +} +int mpfr_zero_p_wrap(const mpfr_ptr p) { + return mpfr_inf_p(p); +} + +int mpfr_set_wrap(const mpfr_ptr p1, const mpfr_ptr p2, mp_rnd_t r) { + return mpfr_set(p1, p2, r); +} + +int mpfr_abs_wrap(const mpfr_ptr p1, const mpfr_ptr p2, mp_rnd_t r) { + return mpfr_abs(p1, p2, r); +} + +int mpfr_cmp_wrap (const mpfr_ptr p1 , const mpfr_ptr p2) { + return mpfr_cmp(p1, p2); +} + +int mpfr_cmp_si_wrap (const mpfr_ptr p1, signed long int p2) { + return mpfr_cmp_si(p1, p2); +} + +int mpfr_cmp_ui_wrap (const mpfr_ptr p1, unsigned long int p2) { + return mpfr_cmp_ui (p1, p2); +} + +int mpfr_sgn_wrap (const mpfr_ptr p1) { + return mpfr_sgn (p1); +} + +int mpfr_set_si_wrap (const mpfr_ptr p, long int si, mp_rnd_t r) { + return mpfr_set_si(p, si, r); +} + +int mpfr_set_ui_wrap (const mpfr_ptr p, unsigned long int si, mp_rnd_t r) { + return mpfr_set_ui(p, si, r); +} + +int mpfr_ceil_wrap (const mpfr_ptr p, const mpfr_ptr p2) { + return mpfr_ceil(p, p2); +} + +int mpfr_floor_wrap (const mpfr_ptr p, const mpfr_ptr p2) { + return mpfr_floor(p, p2); +} + +int mpfr_round_wrap (const mpfr_ptr p, const mpfr_ptr p2) { + return mpfr_round(p, p2); +} + +int mpfr_trunc_wrap (const mpfr_ptr p, const mpfr_ptr p2) { + return mpfr_trunc(p, p2); +} + +mp_prec_t mpfr_get_prec_wrap (const mpfr_ptr e) { + return mpfr_get_prec(e); +} + +mp_exp_t mpfr_get_exp_wrap (const mpfr_ptr p1) { + return mpfr_get_exp (p1); +} + +int mpfr_signbit_wrap (const mpfr_ptr p1) { + return mpfr_signbit (p1); +} + +int mpfr_setsign_wrap (const mpfr_ptr p1, const mpfr_ptr p2, int p3, mp_rnd_t p4) { + return mpfr_setsign (p1, p2, p3, p4); +} + +int mpfr_copysign_wrap (const mpfr_ptr p1, const mpfr_ptr p2, const mpfr_ptr p3, mp_rnd_t p4) { + return mpfr_copysign (p1, p2, p3, p4); +} + +size_t mpfr_custom_get_size_wrap (mp_prec_t p1) { + return mpfr_custom_get_size (p1); +} + +void mpfr_custom_init_wrap (void *p1 , mp_prec_t p2) { + mpfr_custom_init (p1, p2); +} + +void mpfr_custom_init_set_wrap (mpfr_ptr p1, int p2, mp_exp_t p3, mp_prec_t p4, void *p5) { + mpfr_custom_init_set (p1, p2, p3, p4, p5); +} + +int mpfr_custom_get_kind_wrap (mpfr_ptr p1) { + return mpfr_custom_get_kind (p1); +} + +void * mpfr_custom_get_mantissa_wrap (const mpfr_ptr p) { + return mpfr_custom_get_mantissa(p); +} + +mp_exp_t mpfr_custom_get_exp_wrap(const mpfr_ptr p) { + return mpfr_custom_get_exp(p); +} + +void mpfr_custom_move_wrap (mpfr_ptr p1, void *p2 ) { + mpfr_custom_move(p1, p2); +} +/* +intmax_t mpfr_get_sj_wrap (mpfr_ptr p1, mp_rnd_t p2) { + return mpfr_get_sj(p1, p2); +} + +uintmax_t mpfr_get_uj_wrap (mpfr_ptr p1, mp_rnd_t p2) { + return mpfr_get_uj (p1, p2); +} +*/ diff --git a/cbits/chsmpfr.h b/cbits/chsmpfr.h new file mode 100644 --- /dev/null +++ b/cbits/chsmpfr.h @@ -0,0 +1,73 @@ +#include +#include +#include + +mpfr_ptr initS(const mp_prec_t ); + +void clear (const mpfr_ptr ) ; + +// these functions are defined as macros and so haskell ffi +// can't work with them directly + +int mpfr_nan_p_wrap(const mpfr_ptr) ; + +int mpfr_inf_p_wrap(const mpfr_ptr) ; + +int mpfr_zero_p_wrap(const mpfr_ptr) ; + +int mpfr_set_wrap(const mpfr_ptr p1, const mpfr_ptr p2, mp_rnd_t r) ; + +int mpfr_abs_wrap(const mpfr_ptr, const mpfr_ptr, mp_rnd_t) ; + +int mpfr_set_si_wrap (const mpfr_ptr, long int, mp_rnd_t) ; + +int mpfr_set_ui_wrap (const mpfr_ptr, unsigned long int, mp_rnd_t) ; + +int mpfr_cmp_wrap (const mpfr_ptr, const mpfr_ptr) ; + +int mpfr_cmp_si_wrap (const mpfr_ptr, signed long int ) ; + +int mpfr_cmp_ui_wrap (const mpfr_ptr, unsigned long int) ; + +int mpfr_sgn_wrap (const mpfr_ptr) ; + +int mpfr_ceil_wrap (const mpfr_ptr , const mpfr_ptr ) ; + +int mpfr_floor_wrap (const mpfr_ptr , const mpfr_ptr ) ; + +int mpfr_round_wrap (const mpfr_ptr , const mpfr_ptr ) ; + +int mpfr_trunc_wrap (const mpfr_ptr , const mpfr_ptr ) ; + +mp_prec_t mpfr_get_prec_wrap (const mpfr_ptr ) ; + + +mp_exp_t mpfr_get_exp_wrap (const mpfr_ptr ) ; + +int mpfr_sign_bit_wrap (const mpfr_ptr ) ; + +int mpfr_setsign_wrap (const mpfr_ptr , const mpfr_ptr, int , mp_rnd_t ) ; + +int mpfr_copysign_wrap (const mpfr_ptr , const mpfr_ptr , const mpfr_ptr , mp_rnd_t ) ; + +int mpfr_signbit_wrap (mpfr_ptr ) ; + +size_t mpfr_custom_get_size_wrap (mp_prec_t) ; + +void mpfr_custom_init_wrap (void * , mp_prec_t) ; + +void mpfr_custom_init_set_wrap (const mpfr_ptr , int , mp_exp_t , mp_prec_t , void *) ; +/* +intmax_t mpfr_get_sj_wrap (mpfr_ptr, mp_rnd_t ); + +uintmax_t mpfr_get_uj_wrap (mpfr_ptr, mp_rnd_t ); +*/ + + +int mpfr_custom_get_kind_wrap (const mpfr_ptr ) ; + +void * mpfr_custom_get_mantissa_wrap (const mpfr_ptr ) ; + +mp_exp_t mpfr_custom_get_exp_wrap(const mpfr_ptr ) ; + +void mpfr_custom_move_wrap (const mpfr_ptr , void * ) ; diff --git a/hmpfr.cabal b/hmpfr.cabal new file mode 100644 --- /dev/null +++ b/hmpfr.cabal @@ -0,0 +1,50 @@ +name: hmpfr +version: 0.1 +synopsis: Haskell binding to MPFR library +description: Haskell binding to MPFR library. Tested with MPFR 2.3.1 and 2.3.2. + Some simple examples of usage can be found in test/Demo.hs. +category: Data, Math +license: BSD3 +license-file: LICENSE +Stability: experimental +tested-with: GHC == 6.8.2 GHC == 6.8.3 +author: Aleš Bizjak +maintainer: Aleš Bizjak +Homepage: http://code.haskell.org/hmpfr/ +build-Depends: base + +Other-modules: Data.Number.MPFR.FFIhelper + Data.Number.MPFR.Internal + Data.Number.MPFR.Base + + +Exposed-modules: Data.Number.MPFR.Assignment + Data.Number.MPFR.Conversion + Data.Number.MPFR.Arithmetic + Data.Number.MPFR.Comparison + Data.Number.MPFR.Special + Data.Number.MPFR.Integer + Data.Number.MPFR.Misc + + Data.Number.MPFR.Near + Data.Number.MPFR.Up + Data.Number.MPFR.Down + Data.Number.MPFR.Zero + + Data.Number.MPFR + +Extra-source-files: test/Demo.hs + +Data-files: README + +build-type: Simple +build-tools: hsc2hs +GHC-options: -Wall -fno-warn-orphans -O2 +GHC-prof-options: -prof -auto-all +include-dirs: cbits +includes: chsmpfr.h mpfr.h +install-includes: chsmpfr.h +c-sources: cbits/chsmpfr.c + +extra-libraries: mpfr + diff --git a/test/Demo.hs b/test/Demo.hs new file mode 100644 --- /dev/null +++ b/test/Demo.hs @@ -0,0 +1,38 @@ +module Demo where + +import qualified Data.Number.MPFR as M + +-- compute the sum from 1 to n with precision of p bits rounded to Near +s :: M.Precision -> Int -> M.MPFR +s p n = s' 1 0 + where s' k acc | k <= n = s' (succ k) (M.add M.Near p acc (M.fromInt M.Near 32 k)) + | otherwise = acc + +s' :: M.Precision -> Int -> M.MPFR +s' p = foldl (M.addi M.Near p) 0 . enumFromTo 1 + +-- compute pi with precision of n bits +pi' n = M.pi M.Near n + +-- compute pi and get an indicator where the result is rounded +pi'' n = fst (M.mpfrToString M.Near 0 10 p) ++ case compare i 0 of + GT -> " result is rounded up" + EQ -> " result is exact" + LT -> " result is rounded down" + where (p, i) = M.pi_ M.Near n + + + +-- sum up first n terms of a Taylor series for e with precision p +e :: M.Precision -> Int -> M.MPFR +e p n = e' 1 1 1 + where e' k acc acc' | k == n = acc + | otherwise= e' (succ k) (M.add M.Down p acc (M.div M.Down p M.one acc'')) acc'' + where acc'' = M.muli M.Up p acc' k + +-- factorial of n with p bits rounded to Near +fac :: M.Precision -> Int -> M.MPFR +fac p n = s' 1 1 + where s' k acc | k <= n = s' (succ k) (M.muli M.Near p acc k) + | otherwise = acc +