HERA (empty) → 0.2
raw patch · 12 files changed
+2691/−0 lines, 12 filesdep +basesetup-changed
Dependencies added: base
Files
- C/hsmpfr.c +131/−0
- C/hsmpfr.h +73/−0
- Data/Number/Ball.hs +354/−0
- Data/Number/Dyadic.hs +9/−0
- Data/Number/DyadicInterval.hs +181/−0
- Data/Number/FFIhelper.hsc +725/−0
- Data/Number/MPFR.hs +877/−0
- Data/Number/Real.hs +266/−0
- Data/Order.hs +10/−0
- HERA.cabal +29/−0
- LICENSE +30/−0
- Setup.hs +6/−0
+ C/hsmpfr.c view
@@ -0,0 +1,131 @@+#include "hsmpfr.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);+}+*/
+ C/hsmpfr.h view
@@ -0,0 +1,73 @@+#include <mpfr.h>+#include <malloc.h>+#include <inttypes.h>++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 * ) ;
+ Data/Number/Ball.hs view
@@ -0,0 +1,354 @@+{-# INCLUDE <mpfr.h> #-}+{-# INCLUDE <hsmpfr.h> #-}++module Data.Number.Ball (Ball(..), makeA, make, + normalizeBall,+ lower, upper, lower_, upper_,+ sgnLower, sgnUpper,+ width, compareB,+ below, contains, + intersectA, intersect,+ add, sub, neg, absB, mul, div, sqrt, exp, log,+ maxB, minB,+ fromDyadic, fromString, fromInt, fromWord )+where++import qualified Data.Number.Dyadic as D+import Data.Order++import Prelude hiding (div, sqrt, exp, log)+import Data.Word(Word)++-- | Ball represents a closed interval @[center-radius, center+radius] @+data Ball = Ball { center :: !D.Dyadic, -- ^ center of the ball+ radius :: !D.Dyadic -- ^ radius of the ball + }+{-+instance Show Ball where+ show b@(Ball c r) = "\ncenter = " ++ D.toString dc c ++ "\n" ++ "radius = " ++ + D.toString dr r+ where dc = min 60 $ (decimalPrec . correctDigits) b+ dr = D.getPrec r+ -}+instance Show Ball where+ show b@(Ball c r) = s ++ "[" ++ show go ++ "]"+ where go' = decimalPrec . correctDigits $ b+ go = let r' = D.getExp r in + if go' == 0 && r' < 0 then decimalPrec . fromIntegral . negate . succ $ r' else go'+ s = D.toString go c++-- | Precision of ball\'s radius. +radPrec :: D.Precision+radPrec = 32++-- | Create epsilon neighbourhood of d according to the number of accurate digits of d.+-- Specifically return m * 2 ^ (e - p - 1) +createEpsilon :: Int -- ^ m+ -> D.Dyadic -- ^ dyadic with magnitude e and precision p+ -> D.Dyadic+createEpsilon i d = D.int2i D.Zero radPrec i (if d == 0 then 0 else D.getExp d - (fromIntegral $ D.getPrec d) - 1)++-- | If first arugment \/= 0 then add to second argument the epsilon of the third.+addEpsilon :: Int -- ^ indicates whether correction is necessary + -> D.Dyadic -- ^ dyadic to be corrected+ -> D.Dyadic -- ^ dyadic which indicates the magnitude of correction+ -> D.Dyadic+addEpsilon e d x = if e /= 0 then D.add D.Up (D.getPrec d) d (createEpsilon 1 x)+ else d++-- | Make a ball from endpoints+makeA :: D.Precision -- ^ desired precision of the center+ -> D.Dyadic -- ^ left endpoint+ -> D.Dyadic -- ^ right endpoint+ -> Ball +makeA p d1 d2 = Ball cen rad+ where (c,e) = D.add_ D.Near p d1 d2+ cen = D.div2w D.Near p c 1+ r = D.sub D.Up radPrec d2 d1+ r' = D.div2w D.Up radPrec r 1+ rad = addEpsilon e r' c+++-- | Make a ball from endpoints so that no precision is lost. +make :: D.Dyadic -- ^ left endpoint+ -> D.Dyadic -- ^ right endpoint+ -> Ball +make d1 d2 = makeA (D.addPrec d1 d2) d1 d2+ +-- | Normalize the given ball's center to the specified precision.+-- Resulting ball might be larger.+normalizeBall :: D.Precision -> Ball -> Ball+normalizeBall p (Ball c r) = Ball c' r'+ where (c',e) = D.set_ D.Near p c+ r'' = D.set D.Up radPrec r+ r' = addEpsilon e r'' c'++-- | MakeA a ball from dyadic. Radius is 0 if desired precision is not smaller+-- than precision of dyadic.+fromDyadic :: D.Precision -> D.Dyadic -> Ball+fromDyadic p d = Ball c r+ where (c, e) = D.set_ D.Near p d+ r' = D.fromWord D.Up radPrec 0+ r = addEpsilon e r' c++-- | Similar to fromDyadic.+fromInt :: D.Precision -> Int -> Ball+fromInt p d = Ball c r+ where (c, e) = D.fromInt_ D.Near p d+ r' = D.fromWord D.Up radPrec 0+ r = addEpsilon e r' c++-- | Similar to fromInt.+fromWord :: D.Precision -> Word -> Ball+fromWord p = fromInt p . fromIntegral++-- | Lower endpoint of the ball rounded down to specified precision.+lower :: D.Precision -> Ball -> D.Dyadic +lower p (Ball c r) = D.sub D.Down p c r++-- | Upper endpoint of the ball rounded up to specified precision.+upper :: D.Precision -> Ball -> D.Dyadic +upper p (Ball c r) = D.add D.Up p c r++-- | Lower endpoint with precision of the center+lower_ :: Ball -> D.Dyadic+lower_ b@(Ball c _) = lower (D.getPrec c) b++-- | Upper endpoint with precision of the center+upper_ :: Ball -> D.Dyadic+upper_ b@(Ball c _) = upper (D.getPrec c) b++-- | Sign of lower endpoint of the ball. This should be faster than using @ signum (center b - radius b) @ +sgnLower :: Ball -> Int+sgnLower (Ball c r) = case compare c r of+ LT -> -1+ EQ -> 0+ _ -> 1++-- | Analogous to sgnLower.+sgnUpper :: Ball -> Int+sgnUpper (Ball c r) = case compare (D.neg D.Near (D.getPrec r) r) c of + LT -> 1+ EQ -> 0+ _ -> -1++-- | Upper bound on the width of the ball. @ 2 * radius b @ rounded up.+width :: Ball -> D.Dyadic+width (Ball _ r) = D.mul2w D.Up radPrec r 1++-- | Check if second ball is included in the first+below :: Ball -> Ball -> Bool+below a b = lower_ a <= lower_ b && upper_ a >= upper_ b++-- | Check if dyadic is element of the ball.+contains :: Ball -> D.Dyadic -> Bool+contains b d = lower_ b <= d && upper_ b >= d++-- | Returns an intersection of two balls. If balls are disjoint then computation fails with fail.+intersectA :: Monad m => D.Precision -- ^ precision of the resulting ball's center+ -> Ball -> Ball -> m Ball+intersectA p b1 b2 | l <= u = return $ makeA p l u+ | otherwise = fail "cannot intersect disjoint intervals"+ where l = D.maxD D.Down p (lower p b1) (lower p b2)+ u = D.minD D.Up p (upper p b1) (upper p b2)++-- | Intersection of two balls exactly (no precision is lost).+intersect :: Monad m => Ball -> Ball -> m Ball+intersect b1@(Ball c _) b2@(Ball c' _) = intersectA (D.addPrec c c') b1 b2++-- | Addition of two balls.+--+-- - @ center = center a + center b @+-- +-- - @ radius = radius a + radius b @+--+-- Rounding errors are added to the radius.+add :: D.Precision -> Ball -> Ball -> Ball+add p (Ball c r) (Ball c' r') = Ball cen rad+ where (cen, e) = D.add_ D.Near p c c'+ rad = D.add D.Up radPrec r' (addEpsilon e r cen)++-- | Subtraction of two balls.+--+-- - @ center = center a - center b @+-- +-- - @ radius = radius a + radius b @+--+-- Rounding errors are added to the radius.+sub :: D.Precision -> Ball -> Ball -> Ball+sub p (Ball c r) (Ball c' r') = Ball cen rad+ where (cen, e) = D.sub_ D.Near p c c'+ rad = D.add D.Up radPrec r' (addEpsilon e r cen)+-- | Negation of the ball. +--+-- - center = - center b rounded to specified precision.+-- +-- - radius is only modified for the rounding error.+neg :: D.Precision -> Ball -> Ball+neg p (Ball c r) = Ball c' r'+ where (c',e) = D.neg_ D.Near p c+ r' = addEpsilon e r c'++absB :: D.Precision -> Ball -> Ball+absB p b = if lb > 0 then normalizeBall p b+ else if ub < 0 then neg p b+ else makeA p 0 (max (negate lb) (ub))+ where lb = lower_ b+ ub = upper_ b+++-- | Multiplication of two balls. (centers of both balls are assumed positive)+--+-- - If none of the balls contains 0 then+--+-- @ center = center a * center b + radius a * radius b @+--+-- @ radius = center a * radius b + radius a * center b @+-- +-- - If one of the operands (left) contains 0+-- +-- @ center = center a * upper b @+--+-- @ radius = radius a * upper b @+-- +-- - If both of the balls contain 0+--+-- @ lower = min ((lower a) * (upper b)) ((lower b) * (upper a)) @+-- +-- @ upper = max ((lower a) * (lower b)) ((upper b) * (upper a)) @+--+-- Rounding errors are added to the radius.+mul :: D.Precision -> Ball -> Ball -> Ball+mul p b1 b2 = if D.sgn (center b1) * D.sgn (center b2) < 0 then neg p ret else ret+ where ret = mul' (absB p b1) (absB p b2)+ mul' b b' = case (sgnLower b, sgnLower b') of + (1, 1) -> nonzero b b'+ (1, _) -> leftzero b' b+ (_, 1) -> leftzero b b'+ _ -> bothzero b b'+ nonzero (Ball c r) (Ball c' r') = Ball cen rad + where r'' = D.fma D.Up radPrec c r' (D.mul D.Up radPrec c' r) + cr = D.mul D.Near (2 * radPrec) r r' + (cen, e) = D.fma_ D.Near p c c' cr+ rad = addEpsilon e r'' cen+ leftzero (Ball c r) b = Ball cen rad+ where (cen, e) = D.mul_ D.Near p c up+ rad = addEpsilon e (D.mul D.Up radPrec r up) cen+ up = upper p b+ bothzero b b' = makeA p l u+ where l = D.minD D.Down p l1 l2+ u = D.maxD D.Up p u1 u2+ l1 = D.mul D.Down p (lower p b) (upper p b')+ l2 = D.mul D.Down p (upper p b) (lower p b')+ u1 = D.mul D.Up p (lower p b) (lower p b')+ u2 = D.mul D.Up p (upper p b) (upper p b')++-- | Division of two balls+-- +-- - If radius is \"large\" then divide endpoints and makeA a ball from them.+--+-- - If radius is \"small\" then division can be optimized+--+-- - @ center = center a / center b @+-- +-- - @ (radius = radius a * center b + center a * radius b) / (center b * center b) + 2 * 2 ^ (e1 - e2 - p)@ +-- where @ p @ is precision of the result, @ e1 = getExp c1, e2 = getExp c2 @. This way the resulting interval is +-- guaranteed to be correct.+--+-- Rounding errors are added to the radius.+--+-- If divisor ball contains zero compuatation fails with fail.+div :: (Monad m) => D.Precision -> Ball -> Ball -> m Ball+div p b1 b2 = if sgnLower b2 > 0 then return (div' b1 b2)+ else if sgnUpper b2 < 0 then return (neg p (div' b1 (neg p b2)))+ else fail "Division by interval containing zero"+ where div' b b' = if smallRad b && smallRad b' then divSmall b b'+ else divLarge b b'+ -- radius is small if (a) it is 0 or if number of correct digits is + -- more than half of the end precision+ smallRad (Ball c r) = D.sgn r == 0 || 2 * (D.getExp c - D.getExp r) > fromIntegral p + 2+ divSmall (Ball c r) (Ball c' r') = Ball cen rad+ where cen = D.div D.Near p c c'+ r'' = D.fma D.Up radPrec c r' (D.mul D.Up radPrec c' r)+ (bsq, e) = D.sqr_ D.Down radPrec c'+ bsq' = if e == 0 then D.nextBelow bsq else bsq+ r''' = D.div D.Up radPrec r'' bsq'+ -- now r''' is at most 2 * 2 ^ (e1 - e2 - p) too small+ rad = D.add D.Up radPrec r''' (createEpsilon 3 cen)+ divLarge b b' = makeA p l u+ where l = D.div D.Down p l1 (if D.sgn l1 < 0 then l2 else u2)+ u = D.div D.Up p u1 (if D.sgn u1 < 0 then u2 else l2)+ l1 = lower p b+ l2 = lower p b'+ u1 = upper p b+ u2 = upper p b'++-- | Square root of a ball. If interval contains 0 then computation fails.+sqrt :: Monad m => D.Precision -> Ball -> m Ball+sqrt p b@(Ball c r) = if lower_ b < 0 then fail "Sqrt of a interval containing negative numbers"+ else if radSmall then return sqrt_small+ else return sqrt_large+ where sqrt_large = makeA p l u+ where l = D.sqrt D.Down p (D.sub D.Down p c r)+ u = D.sqrt D.Up p (D.add D.Up p c r)+ radSmall = D.sgn c /= 0 && (D.sgn r == 0 || D.getExp c `quot` 2 - D.getExp r > fromIntegral p)+ sqrt_small = Ball cen rad+ where (cen,e) = D.sqrt_ D.Near p c+ rad' = D.div D.Up radPrec r cen+ rad = addEpsilon e rad' cen++-- | @ e ^ b @+exp :: D.Precision -> Ball -> Ball+exp p (Ball c r) = makeA p l u+ where l = D.exp D.Down p (D.add D.Down p c r)+ u = D.exp D.Up p (D.add D.Up p c r)++-- | Natural logarithm of a ball. If interval contains 0 then computation fails.+log :: Monad m => D.Precision -> Ball -> m Ball+log p b@(Ball c r) = if lower_ b <= 0 then fail "Domain of log is R+"+ else return (makeA p l u)+ where l = D.log D.Down p (D.add D.Down p c r)+ u = D.log D.Up p (D.add D.Up p c r)++-- | Compare two balls.+--+-- - if upper a < lower b then Less+--+-- - if upper b < lower a then Greater +--+-- - otherwise balls are incomparable.+compareB :: Ball -> Ball -> POrdering+compareB b b' = if upper_ b < lower_ b' then Less+ else if lower_ b > upper_ b' then Greater+ else Incomparable++-- | Maximum of two balls, meaning:+--+-- - lower = max (lower a) (lower b) rounded down+--+-- - upper = max (upper a) (upper b) rounded up+maxB :: D.Precision -> Ball -> Ball -> Ball+maxB p b b' = makeA p l u+ where l = D.maxD D.Down p (lower p b) (lower p b')+ u = D.maxD D.Up p (upper p b) (upper p b')++-- | Analogous to maxB.+minB :: D.Precision -> Ball -> Ball -> Ball+minB p b b' = makeA p l u+ where l = D.minD D.Down p (lower p b) (lower p b')+ u = D.minD D.Up p (upper p b) (upper p b')++-- | Similar to fromDyadic.+fromString :: D.Precision -> String -> Ball+fromString p s = Ball cen rad+ where cen = D.fromString s p 10+ rad = createEpsilon 1 cen++decimalPrec :: Word -> Word+decimalPrec d = floor (fromIntegral d * logBase 10 2 :: Double)++correctDigits :: Ball -> Word+correctDigits (Ball c r) = case compare D.zero r of+ EQ -> (fromIntegral . D.getPrec) c+ LT -> let cd = D.getExp c - D.getExp r in fromIntegral (max 0 cd)+ _ -> error "Ball.correctDigits: radius should be a nonnegative, non-degenerate number"
+ Data/Number/Dyadic.hs view
@@ -0,0 +1,9 @@+module Data.Number.Dyadic (+ module Data.Number.MPFR,+ pow2+) where++import Data.Number.MPFR ++pow2 :: Int -> Dyadic+pow2 = int2i Near minPrec 1
+ Data/Number/DyadicInterval.hs view
@@ -0,0 +1,181 @@+{-# LANGUAGE BangPatterns #-}+module Data.Number.DyadicInterval(Interval,+ fromBallA, fromBall, make, makeA,+ below, contains, includes, + intersectA, intersect, + neg, add, mul, sub, div, sqrt, exp, log,+ compareI, maxI, minI,+ center, radius, lower, upper, width,+ fromDyadic, fromString, fromInt, fromWord,+ toString)+ where++import qualified Data.Number.Ball as B+import qualified Data.Number.Dyadic as D++import Data.Order++import Prelude hiding (div, sqrt, exp, log)+import Data.Word(Word)++import Control.Monad+ +-- | A wrapper around Ball allowing the results of operations like division+-- by interval containing zero to be represented and do not cause errors.+--+-- Nothing represents undefined interval.+type Interval = Maybe B.Ball++data Inclusion = Above | Below | NoInclusion++-- | Make an interval from a ball and normalize it to specified precision.+fromBallA :: D.Precision -> B.Ball -> Interval+fromBallA p b = Just (B.normalizeBall p b)++-- | Just make an interval from a ball.+fromBall :: B.Ball -> Interval+fromBall b = Just b++-- | Make an interval from two endpoints.+makeA :: D.Precision -- ^ precision of the interval's center+ -> D.Dyadic -- ^ left endpoint+ -> D.Dyadic -- ^ right endpoint+ -> Interval+makeA p l u = Just (B.makeA p l u)++-- | Make an interval from two endpoints so that no precision is lost.+make :: D.Dyadic -> D.Dyadic -> Interval+make l u = Just (B.make l u)++-- | Checks if second interval is inside the first. _|_ is above all.+below :: Interval -> Interval -> Bool+below (Just b) (Just b') = B.below b b'+below Nothing _ = True+below _ Nothing = False++-- | Checks if interval contains dyadic. _|_ contains everything.+contains :: Interval -> D.Dyadic -> Bool+contains Nothing _ = True+contains (Just b) d = B.contains b d++-- | Returns Below if second interval is inside first, Above if converse, NoInclusion otherwise.+includes :: Interval -> Interval -> Inclusion+includes i i' | below i i' = Below+ | below i' i = Above+ | otherwise = NoInclusion++-- | Return the intersection of two intervals. The resulting interval's center has specified precision.+-- +-- If one of the intervals is _|_ then just return the other (even if it is _|_).+intersectA :: D.Precision -> Interval -> Interval -> Interval+intersectA p (Just b) (Just b') = B.intersectA p b b'+intersectA _ Nothing i = i+intersectA _ i Nothing = i++-- | Return the intersection of two intervals so that no precision is lost.+intersect :: Interval -> Interval -> Interval+intersect (Just b) (Just b') = B.intersect b b'+intersect Nothing x = x+intersect x Nothing = x++-- | Negate the interval. neg _|_ = _|_.+neg :: D.Precision -> Interval -> Interval+neg p i = do b <- i+ return $! B.neg p b++{-# INLINE wrap #-}+wrap :: (D.Precision -> B.Ball -> B.Ball -> B.Ball)+ -> D.Precision -> Interval -> Interval -> Interval+wrap f p i i' = do !b <- i+ !b' <- i'+ return $! f p b b'++-- | Addition. If one of the arguments is _|_, so is the result.+add :: D.Precision -> Interval -> Interval -> Interval+add = wrap B.add++-- | Multiplication. If one of the arguments is _|_, so is the result+mul :: D.Precision -> Interval -> Interval -> Interval+mul = wrap B.mul++-- | Subtraction. If one of the arguments is _|_, so is the result+sub :: D.Precision -> Interval -> Interval -> Interval+sub = wrap B.sub++-- | Division. If one of the arguments is _|_ or divisor contains 0 then result is _|_.+div :: D.Precision -> Interval -> Interval -> Interval+div p i i' = do !b <- i+ !b' <- i'+ B.div p b b'+-- | Square root. If one argument is _|_ or interval contains 0 then result is _|_.+sqrt :: D.Precision -> Interval -> Interval+sqrt p i = do !b <- i+ B.sqrt p b++-- | Natural logarithm. If one argument is _|_ or interval contains 0 then result is _|_.+log :: D.Precision -> Interval -> Interval+log p i = do !b <- i+ B.log p b++-- | @ e ^ i @ If argument is _|_ so is the result.+exp :: D.Precision -> Interval -> Interval+exp p i = do !b <- i+ return $! (B.exp p b)++-- | Compare two intervals. If one of them is _|_ the result is incomparable, +-- otherwise result is comparison of balls.+compareI :: Interval -> Interval -> POrdering+compareI (Just b) (Just b') = B.compareB b b'+compareI _ _ = Incomparable++-- | Maximum of intervals. If one interval is _|_ so is the result.+maxI :: D.Precision -> Interval -> Interval -> Interval+maxI = wrap B.maxB++-- | Similar to maxI.+minI :: D.Precision -> Interval -> Interval -> Interval+minI = wrap B.minB++-- | Center of interval. Center on _|_ will result in fail.+center :: (Monad m) => Interval -> m D.Dyadic+center Nothing = fail "center of _|_ is not defined"+center (Just (B.Ball c _)) = return c++-- | Radius of interval. Radius on _|_ will result in fail.+radius :: (Monad m) => Interval -> m D.Dyadic+radius Nothing = fail "radius of _|_ is infinity"+radius (Just (B.Ball _ r)) = return r++-- | Lower endpoint of interval with precision of the center. +-- Lower on _|_ will result in fail.+lower :: (Monad m) => Interval -> m D.Dyadic+lower Nothing = fail "lower bound of _|_ is -infinity"+lower (Just b) = return (B.lower_ b)++-- | Upper endpoint of interval with precision of the center. +-- Upper on _|_ will result in fail.+upper :: (Monad m) => Interval -> m D.Dyadic+upper Nothing = fail "upper bound of _|_ is +infinity"+upper (Just b) = return (B.upper_ b)+++-- | Width of the interval. Widht on _|_ will result in fail.+width :: (Monad m) => Interval -> m D.Dyadic+width Nothing = fail "width of _|_ is infinity"+width (Just b) = return (B.width b)++fromDyadic :: D.Precision -> D.Dyadic -> Interval+fromDyadic p d = Just (B.fromDyadic p d)++fromString :: D.Precision -> String -> Interval+fromString p s = Just (B.fromString p s)++fromInt :: D.Precision -> Int -> Interval+fromInt p d = Just (B.fromInt p d)++fromWord :: D.Precision -> Word -> Interval+fromWord p d = Just (B.fromWord p d)++toString :: Interval -> String+toString Nothing = "_|_"+toString (Just b) = show b
+ Data/Number/FFIhelper.hsc view
@@ -0,0 +1,725 @@+{-# LANGUAGE ForeignFunctionInterface #-}+#include <hsmpfr.h>+#include <mpfr.h>++module Data.Number.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)++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_T = MP CPrecision Sign Exp !(ForeignPtr Limb)++instance Storable MPFR_T where+ sizeOf _ = #size __mpfr_struct+ alignment _ = (undefined :: Int)+ peek = error "Not needed and not applicable"+ poke p (MP prec s e fp) = do withForeignPtr fp $ \p1 -> do + #{poke __mpfr_struct, _mpfr_prec} p prec+ #{poke __mpfr_struct, _mpfr_sign} p s + #{poke __mpfr_struct, _mpfr_exp} p e+ #{poke __mpfr_struct, _mpfr_d} p p1++peekP :: Ptr MPFR_T -> ForeignPtr Limb -> IO MPFR_T+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++--data MPFR_T = MPFR_T++-- utility functions from hsmpfr.h+foreign import ccall unsafe "initS"+ initS :: CPrecision -> IO (Ptr MPFR_T)++foreign import ccall unsafe "&clear"+ clear :: FunPtr (Ptr MPFR_T -> IO ())++--------------------+foreign import ccall unsafe "mpfr_get_prec_wrap"+ mpfr_get_prec :: Ptr MPFR_T -> IO CPrecision ++----------------------------------------------------------------++-- assignment functions+foreign import ccall unsafe "mpfr_set_wrap"+ mpfr_set :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_set_ui_wrap"+ mpfr_set_ui :: Ptr MPFR_T -> CULong -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_set_si_wrap"+ mpfr_set_si :: Ptr MPFR_T -> CLong -> CRoundMode -> IO CInt++--TODO set_uj, set_sj++foreign import ccall unsafe "mpfr_set_d"+ mpfr_set_d :: Ptr MPFR_T -> CDouble -> CRoundMode -> IO CInt++--foreign import ccall unsafe "mfpr_set_ld"+ -- mpfr_set_ld :: Ptr MPFR_T -> #{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_T -> CULong -> CInt -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_set_si_2exp"+ mpfr_set_si_2exp :: Ptr MPFR_T -> CLong -> CInt -> CRoundMode -> IO CInt++--TODO set_uj_2exp, set_sj_2exp++foreign import ccall unsafe "mpfr_set_str"+ mpfr_set_str :: Ptr MPFR_T -> CString -> CInt -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_strtofr"+ mpfr_strtofr :: Ptr MPFR_T -> CString -> Ptr (Ptr CChar) -> CInt -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_set_inf"+ mpfr_set_inf :: Ptr MPFR_T -> CInt -> IO ()++foreign import ccall unsafe "mpfr_set_nan"+ mpfr_set_nan :: Ptr MPFR_T -> IO ()++foreign import ccall unsafe "mpfr_swap"+ mpfr_swap :: Ptr MPFR_T -> Ptr MPFR_T -> 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_T -> CRoundMode -> IO CDouble++-- TODO get_decimal64++foreign import ccall unsafe "mpfr_get_d_2exp"+ mpfr_get_d_2exp :: Ptr CLong -> Ptr MPFR_T -> CRoundMode -> IO CDouble++-- TODO get_ld_2exp+-- !!!!!!! next 4 set erange flags+foreign import ccall unsafe "mpfr_get_si" + mpfr_get_si :: Ptr MPFR_T -> CRoundMode -> IO CLong++foreign import ccall unsafe "mpfr_get_ui" + mpfr_get_ui :: Ptr MPFR_T -> CRoundMode -> IO CULong++{-+foreign import ccall unsafe "mpfr_get_sj_wrap"+ mpfr_get_sj :: Ptr MPFR_T -> CRoundMode -> IO #type intmax_t++foreign import ccall unsafe "mpfr_get_uj_wrap"+ mpft_get_uj :: Ptr MPFR_T -> 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 CInt -> CInt -> CUInt -> Ptr MPFR_T -> 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_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_fits_slong_p"+ mpfr_fits_slong_p :: Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_fits_uint_p"+ mpfr_fits_uint_p :: Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_fits_sint_p"+ mpfr_fits_sint_p :: Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_fits_ushort_p"+ mpfr_fits_ushort_p :: Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_fits_sshort_p"+ mpfr_fits_sshort_p :: Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_fits_intmax_p"+ mpfr_fits_intmax_p :: Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_fits_uintmax_p"+ mpfr_fits_uintmax_p :: Ptr MPFR_T -> CRoundMode -> IO CInt+++-------------------------------------------------------------------------------++-- basic arithmetic functions++foreign import ccall unsafe "mpfr_add"+ mpfr_add :: Ptr MPFR_T -> Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_add_ui"+ mpfr_add_ui :: Ptr MPFR_T -> Ptr MPFR_T -> CULong -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_add_si"+ mpfr_add_si :: Ptr MPFR_T -> Ptr MPFR_T -> CLong -> CRoundMode -> IO CInt++-- TODO add_z, add_q++foreign import ccall unsafe "mpfr_sub"+ mpfr_sub :: Ptr MPFR_T -> Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_ui_sub" + mpfr_ui_sub :: Ptr MPFR_T -> CULong -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_sub_ui"+ mpfr_sub_ui :: Ptr MPFR_T -> Ptr MPFR_T -> CULong -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_si_sub" + mpfr_si_sub :: Ptr MPFR_T -> CLong -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_sub_si"+ mpfr_sub_si :: Ptr MPFR_T -> Ptr MPFR_T -> CLong -> CRoundMode -> IO CInt++--TODO sub_z, sub_q++foreign import ccall unsafe "mpfr_mul"+ mpfr_mul :: Ptr MPFR_T -> Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt ++foreign import ccall unsafe "mpfr_mul_ui"+ mpfr_mul_ui :: Ptr MPFR_T -> Ptr MPFR_T -> CULong -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_mul_si"+ mpfr_mul_si :: Ptr MPFR_T -> Ptr MPFR_T -> CLong -> CRoundMode -> IO CInt++--TODO mul_z, mul_q++foreign import ccall unsafe "mpfr_sqr"+ mpfr_sqr :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_div"+ mpfr_div :: Ptr MPFR_T -> Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_ui_div"+ mpfr_ui_div :: Ptr MPFR_T -> CULong -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_div_ui"+ mpfr_div_ui :: Ptr MPFR_T -> Ptr MPFR_T -> CULong -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_si_div"+ mpfr_si_div :: Ptr MPFR_T -> CLong -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_div_si"+ mpfr_div_si :: Ptr MPFR_T -> Ptr MPFR_T -> CLong -> CRoundMode -> IO CInt++-- TODO div_z, div_q++foreign import ccall unsafe "mpfr_sqrt"+ mpfr_sqrt :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_sqrt_ui"+ mpfr_sqrt_ui :: Ptr MPFR_T -> CULong -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_cbrt"+ mpfr_cbrt :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_root"+ mpfr_root :: Ptr MPFR_T -> Ptr MPFR_T -> CULong -> CRoundMode -> IO CInt ++foreign import ccall unsafe "mpfr_pow"+ mpfr_pow :: Ptr MPFR_T -> Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt ++foreign import ccall unsafe "mpfr_pow_ui"+ mpfr_pow_ui :: Ptr MPFR_T -> Ptr MPFR_T -> CULong -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_pow_si"+ mpfr_pow_si :: Ptr MPFR_T -> Ptr MPFR_T -> CLong -> CRoundMode -> IO CInt++--TODO pow_z++foreign import ccall unsafe "mpfr_ui_pow_ui"+ mpfr_ui_pow_ui :: Ptr MPFR_T -> CULong -> CULong -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_ui_pow"+ mpfr_ui_pow :: Ptr MPFR_T -> CULong -> Ptr MPFR_T -> CRoundMode -> IO CInt+++foreign import ccall unsafe "mpfr_neg"+ mpfr_neg :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt ++foreign import ccall unsafe "mpfr_abs_wrap"+ mpfr_abs :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt ++foreign import ccall unsafe "mpfr_dim"+ mpfr_dim :: Ptr MPFR_T -> Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_mul_2ui"+ mpfr_mul_2ui :: Ptr MPFR_T -> Ptr MPFR_T -> CULong -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_mul_2si"+ mpfr_mul_2si :: Ptr MPFR_T -> Ptr MPFR_T -> CLong -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_div_2ui"+ mpfr_div_2ui :: Ptr MPFR_T -> Ptr MPFR_T -> CULong -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_div_2si"+ mpfr_div_2si :: Ptr MPFR_T -> Ptr MPFR_T -> CLong -> CRoundMode -> IO CInt++++--------------------------------------------------------------------------------+-- comparison functions+-- !!!!!!!! these set erange flags+foreign import ccall unsafe "mpfr_cmp_wrap"+ mpfr_cmp :: Ptr MPFR_T -> Ptr MPFR_T -> IO CInt++foreign import ccall unsafe "mpfr_cmp_ui_wrap"+ mpfr_cmp_ui :: Ptr MPFR_T -> CULong -> IO CInt++foreign import ccall unsafe "mpfr_cmp_si_wrap"+ mpfr_cmp_si :: Ptr MPFR_T -> CLong -> IO CInt++foreign import ccall unsafe "mpfr_cmp_d"+ mpfr_cmp_d :: Ptr MPFR_T -> 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_T -> CULong -> CInt -> IO CInt++foreign import ccall unsafe "mpfr_cmp_si_2exp"+ mpfr_cmp_si_2exp :: Ptr MPFR_T -> CLong -> CInt -> IO CInt++foreign import ccall unsafe "mpfr_cmpabs"+ mpfr_cmpabs :: Ptr MPFR_T -> Ptr MPFR_T -> IO CInt++foreign import ccall unsafe "mpfr_nan_p_wrap"+ mpfr_nan_p :: Ptr MPFR_T -> IO CInt++foreign import ccall unsafe "mpfr_inf_p_wrap"+ mpfr_inf_p :: Ptr MPFR_T -> IO CInt++foreign import ccall unsafe "mpfr_number_p"+ mpfr_number_p :: Ptr MPFR_T -> IO CInt++foreign import ccall unsafe "mpfr_zero_p_wrap"+ mpfr_zero_p :: Ptr MPFR_T -> IO CInt++foreign import ccall unsafe "mpfr_sgn_wrap"+ mpfr_sgn :: Ptr MPFR_T -> IO CInt ++foreign import ccall unsafe "mpfr_greater_p"+ mpfr_greater_p :: Ptr MPFR_T -> Ptr MPFR_T -> IO CInt++foreign import ccall unsafe "mpfr_greaterequal_p"+ mpfr_greaterequal_p :: Ptr MPFR_T -> Ptr MPFR_T -> IO CInt ++foreign import ccall unsafe "mpfr_less_p"+ mpfr_less_p :: Ptr MPFR_T -> Ptr MPFR_T -> IO CInt ++foreign import ccall unsafe "mpfr_lessequal_p"+ mpfr_lessequal_p :: Ptr MPFR_T -> Ptr MPFR_T -> IO CInt ++foreign import ccall unsafe "mpfr_lessgreater_p"+ mpfr_lessgreater_p :: Ptr MPFR_T -> Ptr MPFR_T -> IO CInt ++foreign import ccall unsafe "mpfr_equal_p"+ mpfr_equal_p :: Ptr MPFR_T -> Ptr MPFR_T -> IO CInt ++foreign import ccall unsafe "mpfr_unordered_p"+ mpfr_unordered_p :: Ptr MPFR_T -> Ptr MPFR_T -> IO CInt ++-- special functions ++foreign import ccall unsafe "mpfr_log"+ mpfr_log :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_log2"+ mpfr_log2 :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_log10"+ mpfr_log10 :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_exp"+ mpfr_exp :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_exp2"+ mpfr_exp2 :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_exp10"+ mpfr_exp10 :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_sin"+ mpfr_sin :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_cos"+ mpfr_cos :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_tan"+ mpfr_tan :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_sec"+ mpfr_sec :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_csc"+ mpfr_csc :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_cot"+ mpfr_cot :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_sin_cos"+ mpfr_sin_cos :: Ptr MPFR_T -> Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt+++foreign import ccall unsafe "mpfr_asin"+ mpfr_asin :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_acos"+ mpfr_acos :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_atan"+ mpfr_atan :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_atan2"+ mpfr_atan2 :: Ptr MPFR_T -> Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt+++foreign import ccall unsafe "mpfr_cosh"+ mpfr_cosh :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_sinh"+ mpfr_sinh :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_tanh"+ mpfr_tanh :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt+++foreign import ccall unsafe "mpfr_sech"+ mpfr_sech :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_csch"+ mpfr_csch :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_coth"+ mpfr_coth :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_asinh"+ mpfr_asinh :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_acosh"+ mpfr_acosh :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_atanh"+ mpfr_atanh :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_fac_ui"+ mpfr_fac_ui :: Ptr MPFR_T -> CULong -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_log1p"+ mpfr_log1p :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_expm1"+ mpfr_expm1 :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_eint"+ mpfr_eint :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_gamma"+ mpfr_gamma :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_lngamma"+ mpfr_lngamma :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_lgamma"+ mpfr_lgamma :: Ptr MPFR_T -> Ptr CInt -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_zeta"+ mpfr_zeta :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_zeta_ui"+ mpfr_zeta_ui :: Ptr MPFR_T -> CULong -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_erf"+ mpfr_erf :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_erfc"+ mpfr_erfc :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_j0"+ mpfr_j0 :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_j1"+ mpfr_j1 :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_jn"+ mpfr_jn :: Ptr MPFR_T -> CLong -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_y0"+ mpfr_y0 :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_y1"+ mpfr_y1 :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_yn"+ mpfr_yn :: Ptr MPFR_T -> CLong -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_fma"+ mpfr_fma :: Ptr MPFR_T -> Ptr MPFR_T -> Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt ++foreign import ccall unsafe "mpfr_fms"+ mpfr_fms :: Ptr MPFR_T -> Ptr MPFR_T -> Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt+++foreign import ccall unsafe "mpfr_agm"+ mpfr_agm :: Ptr MPFR_T -> Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt ++foreign import ccall unsafe "mpfr_hypot"+ mpfr_hypot :: Ptr MPFR_T -> Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt ++-- constants+foreign import ccall unsafe "mpfr_const_log2"+ mpfr_const_log2 :: Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_const_pi"+ mpfr_const_pi :: Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_const_euler"+ mpfr_const_euler :: Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_const_catalan"+ mpfr_const_catalan :: Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_free_cache"+ mpfr_free_cache :: IO ()++foreign import ccall unsafe "mpfr_sum"+ mpfr_sum :: Ptr MPFR_T -> Ptr (Ptr MPFR_T) -> CULong -> CRoundMode -> IO CInt++-- TODO input and output functions++-- integer related functions++foreign import ccall unsafe "mpfr_rint"+ mpfr_rint :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_ceil_wrap"+ mpfr_ceil :: Ptr MPFR_T -> Ptr MPFR_T -> IO CInt++foreign import ccall unsafe "mpfr_floor_wrap"+ mpfr_floor :: Ptr MPFR_T -> Ptr MPFR_T -> IO CInt++foreign import ccall unsafe "mpfr_round_wrap"+ mpfr_round :: Ptr MPFR_T -> Ptr MPFR_T -> IO CInt++foreign import ccall unsafe "mpfr_trunc_wrap"+ mpfr_trunc :: Ptr MPFR_T -> Ptr MPFR_T -> IO CInt+ +foreign import ccall unsafe "mpfr_rint_ceil"+ mpfr_rint_ceil :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_rint_floor"+ mpfr_rint_floor :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_rint_round"+ mpfr_rint_round :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_rint_trunc"+ mpfr_rint_trunc :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_frac"+ mpfr_frac :: Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_remainder" + mpfr_remainder :: Ptr MPFR_T -> Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_remquo" + mpfr_remquo :: Ptr MPFR_T -> Ptr CLong -> Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_integer_p"+ mpfr_integer_p :: Ptr MPFR_T -> IO CInt++--------------------+-- miscellaneus functions++foreign import ccall unsafe "mpfr_nexttoward"+ mpfr_nexttoward :: Ptr MPFR_T -> Ptr MPFR_T -> IO ()++foreign import ccall unsafe "mpfr_nextabove"+ mpfr_nextabove :: Ptr MPFR_T -> IO ()++foreign import ccall unsafe "mpfr_nextbelow"+ mpfr_nextbelow :: Ptr MPFR_T -> IO ()++foreign import ccall unsafe "mpfr_min"+ mpfr_min :: Ptr MPFR_T -> Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_max"+ mpfr_max :: Ptr MPFR_T -> Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt++-- TODO urandomb++foreign import ccall unsafe "mpfr_random2"+ mpfr_random2 :: Ptr MPFR_T -> MpSize -> Exp -> IO ()+++foreign import ccall unsafe "mpfr_get_exp_wrap"+ mpfr_get_exp :: Ptr MPFR_T -> IO Exp++foreign import ccall unsafe "mpfr_set_exp"+ mpfr_set_exp :: Ptr MPFR_T -> Exp -> IO CInt++foreign import ccall unsafe "mpfr_signbit_wrap"+ mpfr_signbit :: Ptr MPFR_T -> IO CInt++foreign import ccall unsafe "mpfr_setsign_wrap"+ mpfr_setsign :: Ptr MPFR_T -> Ptr MPFR_T -> CInt -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_copysign_wrap"+ mpfr_copysign :: Ptr MPFR_T -> Ptr MPFR_T -> Ptr MPFR_T -> 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_T -> CInt -> CRoundMode -> IO CInt++foreign import ccall unsafe "mpfr_subnormalize"+ mpfr_subnormalize :: Ptr MPFR_T -> 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_T -> CInt -> Exp -> CPrecision -> Ptr Limb -> IO ()++foreign import ccall unsafe "mpfr_custom_get_kind_wrap"+ mpfr_custom_get_kind :: Ptr MPFR_T -> IO CInt++foreign import ccall unsafe "mpfr_custom_get_mantissa_wrap"+ mpfr_custom_get_mantissa :: Ptr MPFR_T -> IO (Ptr Limb)++foreign import ccall unsafe "mpfr_custom_get_exp_wrap"+ mpfr_custom_get_exp :: Ptr MPFR_T -> IO CInt++foreign import ccall unsafe "mpfr_custom_move_wrap"+ mpfr_custom_move :: Ptr MPFR_T -> Ptr #{type mp_limb_t} -> IO ()++-------------------------------------------------+
+ Data/Number/MPFR.hs view
@@ -0,0 +1,877 @@+{-# LANGUAGE TypeSynonymInstances #-}+{-# INCLUDE <mpfr.h> #-}+{-# INCLUDE <hsmpfr.h> #-}+module Data.Number.MPFR (+-- | This module should always be imported qualified.++-- *** Naming +-- | - functions ending with _ return a pair (value, rounding indicator). +-- Rounding indicator indicates whether the result is rounded and in which+-- directon as described in the MPFR manual.+--+-- - the same functions without the _ return just the value. +--+-- - functions with added \"w\" correspond to MPFR _ui functions+--+-- - functions with added \"i\" correspond to MPFR _si functions+++-- *** Equality testing+-- | Equality works as follows: +-- +-- - NaN \/= Nan, +--+-- - Infinity = Infinity, +--+-- - \-Infinity = -Infinity+--+-- - otherwise normal comparison ++-- *** Ordering +-- | Ordering works as follows:+-- +-- - compare NaN _ = GT+--+-- - compare _ NaN = GT+--+-- - infinity < _ = False+--+-- - \-infinity > _ = False+--+-- - NaN [\<,\>,\>=,<=] _ = False+--+-- This mimics the behaviour of built in haskell Float and Double.++-- *** Num instance+-- | Operations defined in Num will be computed so that no precision is lost.++ Dyadic,+ Precision,+ RoundMode(Near, Zero, Up, Down),+ add, sub, mul, div, inverse,+ add_, sub_, mul_, div_,+ addw, addi, mulw, muli, divw, divi, wdiv, idiv, subw, subi, wsub, isub,+ addw_, addi_, mulw_, muli_, divw_, divi_, wdiv_, idiv_, subw_, subi_, wsub_, isub_,+ mul2w, mul2i, div2w, div2i, mul2w_, mul2i_, div2w_, div2i_,+ int2i, int2w, int2i_, int2w_,+ fma, fms, fma_, fms_, nextBelow,+ sqr, sqrt, root, pow, poww, powi, wpoww, wpow, + sqr_, sqrt_, root_, pow_, poww_, powi_, wpoww_, wpow_,+ exp, exp2, exp10, log, log2, log10, sinh, cosh, tanh,+ exp_, exp2_, exp10_, log_, log2_, log10_, sinh_, cosh_, tanh_,+ neg, absD, dim, neg_, absD_, dim_, + isNaN, isInfinite, isNumber, isZero, greater, greatereq, less, lesseq,+ equal, maxD, minD, maxD_, minD_, sgn, + dyadicToDouble, dyadicToWord, dyadicToInt, dyadicToString, decompose, toStringExp, toString,+ pi, log2c, euler, catalan, pi_, log2c_, euler_, catalan_,+ set, set_,+ fromDouble, fromInt, fromWord, fromDouble_, fromInt_, fromWord_, fromIntegerA, compose, fromString,+ getPrec, getMantissa, getExp, + minPrec, one, zero, addPrec+) where++++import Data.Number.FFIhelper++import Foreign.C(CInt, CLong, CULong, withCString, peekCString)+import Foreign.Marshal(alloca, peekArray)+import Foreign(unsafePerformIO, peek, Ptr, mallocForeignPtrBytes, with)++import Data.Bits(shiftL)++import Data.Word(Word)+import Prelude hiding (div, sqrt, read, isNaN, isInfinite, exp, log, sinh, cosh, tanh, pi)++type Dyadic = MPFR_T++type Precision = Word+++-- these are helper functions, only for internal use+{-# INLINE withDyadicsBA #-}+withDyadicsBA :: RoundMode -> Precision -> Dyadic -> Dyadic+ -> (Ptr MPFR_T -> Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt)+ -> (Dyadic, Int)+withDyadicsBA r p mp1 mp2 f = unsafePerformIO go+ where go = do ls <- mpfr_custom_get_size (fromIntegral p)+ fp <- mallocForeignPtrBytes (fromIntegral ls)+ let dummy = MP (fromIntegral p) 0 0 fp+ with dummy $ \p1 -> do+ with mp1 $ \p2 -> do+ with mp2 $ \p3 -> do+ r2 <- f p1 p2 p3 ((fromIntegral . fromEnum) r)+ r1 <- peekP p1 fp+ return (r1, fromIntegral r2)++{-# INLINE withDyadicBAui #-}+withDyadicBAui :: RoundMode -> Precision -> Dyadic -> CULong+ -> (Ptr MPFR_T -> Ptr MPFR_T -> CULong -> CRoundMode -> IO CInt)+ -> (Dyadic, Int) +withDyadicBAui r p mp1 d f = unsafePerformIO go+ where go = do ls <- mpfr_custom_get_size (fromIntegral p)+ fp <- mallocForeignPtrBytes (fromIntegral ls)+ let dummy = MP (fromIntegral p) 0 0 fp+ with dummy $ \p1 -> do+ with mp1 $ \p2 -> do+ r2 <- f p1 p2 d ((fromIntegral . fromEnum) r)+ r1 <- peekP p1 fp+ return (r1, fromIntegral r2)+ +{-# INLINE withDyadicBAsi #-}+withDyadicBAsi :: RoundMode -> Precision -> Dyadic -> CLong+ -> (Ptr MPFR_T -> Ptr MPFR_T -> CLong -> CRoundMode -> IO CInt)+ -> (Dyadic, Int)+withDyadicBAsi r p mp1 d f = unsafePerformIO go + where go = do ls <- mpfr_custom_get_size (fromIntegral p)+ fp <- mallocForeignPtrBytes (fromIntegral ls)+ let dummy = MP (fromIntegral p) 0 0 fp+ with dummy $ \p1 -> do+ with mp1 $ \p2 -> do+ r2 <- f p1 p2 d ((fromIntegral . fromEnum) r)+ r1 <- peekP p1 fp+ return (r1, fromIntegral r2)+ +{-# INLINE withDyadicBAiu #-}+withDyadicBAiu :: RoundMode -> Precision -> CULong -> Dyadic+ -> (Ptr MPFR_T -> CULong -> Ptr MPFR_T -> CRoundMode -> IO CInt)+ -> (Dyadic, Int) +withDyadicBAiu r p d mp1 f = unsafePerformIO go + where go = do ls <- mpfr_custom_get_size (fromIntegral p)+ fp <- mallocForeignPtrBytes (fromIntegral ls)+ let dummy = MP (fromIntegral p) 0 0 fp+ with dummy $ \p1 -> do+ with mp1 $ \p2 -> do+ r2 <- f p1 d p2 ((fromIntegral . fromEnum) r)+ r1 <- peekP p1 fp+ return (r1, fromIntegral r2)++{-# INLINE withDyadicBAis #-}+withDyadicBAis :: RoundMode -> Precision -> CLong -> Dyadic+ -> (Ptr MPFR_T -> CLong -> Ptr MPFR_T -> CRoundMode -> IO CInt)+ -> (Dyadic, Int) +withDyadicBAis r p d mp1 f = unsafePerformIO go+ where go = do ls <- mpfr_custom_get_size (fromIntegral p)+ fp <- mallocForeignPtrBytes (fromIntegral ls)+ let dummy = MP (fromIntegral p) 0 0 fp+ with dummy $ \p1 -> do+ with mp1 $ \p2 -> do+ r2 <- f p1 d p2 ((fromIntegral . fromEnum) r)+ r1 <- peekP p1 fp+ return (r1, fromIntegral r2)++{-# INLINE withDyadicB #-}+withDyadicB :: Dyadic -> (Ptr MPFR_T -> IO CInt) -> CInt +withDyadicB mp1 f = unsafePerformIO go+ where go = with mp1 $ \p1 -> f p1++withDyadicP :: Dyadic -> (Ptr MPFR_T -> IO CPrecision) -> CPrecision +withDyadicP mp1 f = unsafePerformIO go+ where go = with mp1 $ \p1 -> f p1++{-# INLINE withDyadic #-}+withDyadic :: RoundMode -> Precision -> Dyadic + -> (Ptr MPFR_T -> Ptr MPFR_T -> CRoundMode -> IO CInt) + -> (Dyadic, Int)+withDyadic r p mp1 f = unsafePerformIO go + where go = do ls <- mpfr_custom_get_size (fromIntegral p)+ fp <- mallocForeignPtrBytes (fromIntegral ls)+ let dummy = MP (fromIntegral p) 0 0 fp+ with dummy $ \p1 -> do+ with mp1 $ \p2 -> do+ r2 <- f p1 p2 ((fromIntegral . fromEnum) r)+ r1 <- peekP p1 fp+ return (r1, fromIntegral r2)+ +{-# INLINE withDyadicBB #-}+withDyadicBB :: Dyadic -> Dyadic + -> (Ptr MPFR_T -> Ptr MPFR_T -> IO CInt) + -> CInt +withDyadicBB mp1 mp2 f = unsafePerformIO go+ where go = do with mp1 $ \p1 -> do + with mp2 $ \p2 -> do + f p1 p2+ +{-# INLINE withDyadicC #-}+withDyadicC :: RoundMode -> Precision ->+ (Ptr MPFR_T -> CRoundMode -> IO CInt) -> (Dyadic, Int)+withDyadicC r p f = unsafePerformIO go+ where go = do ls <- mpfr_custom_get_size (fromIntegral p)+ fp <- mallocForeignPtrBytes (fromIntegral ls)+ let dummy = MP (fromIntegral p) 0 0 fp+ with dummy $ \p1 -> do+ r2 <- f p1 ((fromIntegral . fromEnum) r)+ r1 <- peekP p1 fp+ return (r1, fromIntegral r2)+ +checkPrec :: Precision -> Precision+checkPrec = max minPrec++stringToDyadic :: RoundMode -> Precision -> Word -> String -> Dyadic+stringToDyadic r p b d = unsafePerformIO go+ where go = do ls <- mpfr_custom_get_size (fromIntegral p)+ fp <- mallocForeignPtrBytes (fromIntegral ls)+ let dummy = MP (fromIntegral p) 0 0 fp+ with dummy $ \p1 -> do + withCString d $ \p2 -> do + _ <- mpfr_set_str p1 p2 (fromIntegral b) ((fromIntegral . fromEnum) r) + peekP p1 fp+++getMantissa' :: Dyadic -> [Limb]+getMantissa' mp1 = unsafePerformIO go+ where go = do with mp1 $ \p1 -> do + pt <- mpfr_custom_get_mantissa p1 + arr <- peekArray (ceiling ((fromIntegral p ::Double) / fromIntegral bitsPerMPLimb)) pt ;+ return arr + p = getPrec mp1++{- TODO: this is inefficient +binprec :: Integer -> Precision+binprec i = length (takeWhile (/= 0) (iterate (flip shiftR 1) i)+-}++binprec :: Integer -> Precision+binprec d = floor (logBase 2 (fromInteger (if d >= 0 then d else -d)) :: Double) + 1+++--------------------------------------------------------------------++-- pure wrappers for basic arithmetic operations++add :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic+add r p d1 d2 = fst $ add_ r p d1 d2 ++sub :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic+sub r p d1 d2 = fst $ sub_ r p d1 d2++mul :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic+mul r p d1 d2 = fst $ mul_ r p d1 d2++div :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic+div r p d1 d2 = fst $ div_ r p d1 d2++add_ :: RoundMode -> Precision -> Dyadic -> Dyadic -> (Dyadic,Int)+add_ r p d1 d2 = withDyadicsBA r p d1 d2 mpfr_add++sub_ :: RoundMode -> Precision -> Dyadic -> Dyadic -> (Dyadic,Int)+sub_ r p d1 d2 = withDyadicsBA r p d1 d2 mpfr_sub++mul_ :: RoundMode -> Precision -> Dyadic -> Dyadic -> (Dyadic,Int)+mul_ r p d1 d2 = withDyadicsBA r p d1 d2 mpfr_mul++div_ :: RoundMode -> Precision -> Dyadic -> Dyadic -> (Dyadic, Int)+div_ r p d1 d2 = withDyadicsBA r p d1 d2 mpfr_div+++inverse :: Dyadic -> Dyadic+inverse d = div Near (getPrec d) one d ++----------------------------------------------------------------++-- basic arithmetic operations with mixed operands++addw :: RoundMode -> Precision -> Dyadic -> Word -> Dyadic+addw r p d1 d = fst $ addw_ r p d1 d ++addi :: RoundMode -> Precision -> Dyadic -> Int -> Dyadic+addi r p d1 d = fst $ addi_ r p d1 d ++mulw :: RoundMode -> Precision -> Dyadic -> Word -> Dyadic+mulw r p d1 d = fst $ mulw_ r p d1 d ++muli :: RoundMode -> Precision -> Dyadic -> Int -> Dyadic+muli r p d1 d = fst $ muli_ r p d1 d ++divw :: RoundMode -> Precision -> Dyadic -> Word -> Dyadic+divw r p d1 d = fst $ divw_ r p d1 d ++divi :: RoundMode -> Precision -> Dyadic -> Int -> Dyadic+divi r p d1 d = fst $ divi_ r p d1 d ++wdiv :: RoundMode -> Precision -> Word -> Dyadic -> Dyadic+wdiv r p d d1 = fst $ wdiv_ r p d d1 ++idiv :: RoundMode -> Precision -> Int -> Dyadic -> Dyadic+idiv r p d d1 = fst $ idiv_ r p d d1 ++subw :: RoundMode -> Precision -> Dyadic -> Word -> Dyadic+subw r p d1 d = fst $ subw_ r p d1 d ++subi :: RoundMode -> Precision -> Dyadic -> Int -> Dyadic+subi r p d1 d = fst $ subi_ r p d1 d ++wsub :: RoundMode -> Precision -> Word -> Dyadic -> Dyadic+wsub r p d d1 = fst $ wsub_ r p d d1 ++isub :: RoundMode -> Precision -> Int -> Dyadic -> Dyadic+isub r p d d1 = fst $ isub_ r p d d1 ++addw_ :: RoundMode -> Precision -> Dyadic -> Word -> (Dyadic, Int)+addw_ r p d1 d = withDyadicBAui r p d1 (fromIntegral d) mpfr_add_ui++addi_ :: RoundMode -> Precision -> Dyadic -> Int -> (Dyadic, Int)+addi_ r p d1 d = withDyadicBAsi r p d1 (fromIntegral d) mpfr_add_si++mulw_ :: RoundMode -> Precision -> Dyadic -> Word -> (Dyadic, Int)+mulw_ r p d1 d = withDyadicBAui r p d1 (fromIntegral d) mpfr_mul_ui++muli_ :: RoundMode -> Precision -> Dyadic -> Int -> (Dyadic, Int)+muli_ r p d1 d = withDyadicBAsi r p d1 (fromIntegral d) mpfr_mul_si++divw_ :: RoundMode -> Precision -> Dyadic -> Word -> (Dyadic, Int)+divw_ r p d1 d = withDyadicBAui r p d1 (fromIntegral d) mpfr_div_ui++divi_ :: RoundMode -> Precision -> Dyadic -> Int -> (Dyadic, Int)+divi_ r p d1 d = withDyadicBAsi r p d1 (fromIntegral d) mpfr_div_si++wdiv_ :: RoundMode -> Precision -> Word -> Dyadic -> (Dyadic, Int)+wdiv_ r p d d1 = withDyadicBAiu r p (fromIntegral d) d1 mpfr_ui_div++idiv_ :: RoundMode -> Precision -> Int -> Dyadic -> (Dyadic, Int)+idiv_ r p d d1 = withDyadicBAis r p (fromIntegral d) d1 mpfr_si_div++subw_ :: RoundMode -> Precision -> Dyadic -> Word -> (Dyadic, Int)+subw_ r p d1 d = withDyadicBAui r p d1 (fromIntegral d) mpfr_sub_ui++subi_ :: RoundMode -> Precision -> Dyadic -> Int -> (Dyadic, Int)+subi_ r p d1 d = withDyadicBAsi r p d1 (fromIntegral d) mpfr_sub_si++wsub_ :: RoundMode -> Precision -> Word -> Dyadic -> (Dyadic, Int)+wsub_ r p d d1 = withDyadicBAiu r p (fromIntegral d) d1 mpfr_ui_sub++isub_ :: RoundMode -> Precision -> Int -> Dyadic -> (Dyadic, Int)+isub_ r p d d1 = withDyadicBAis r p (fromIntegral d) d1 mpfr_si_sub++----------------------------------------------------------++-- multiplication and division with 2 ^ x++mul2w :: RoundMode -> Precision -> Dyadic -> Word -> Dyadic+mul2w r p d1 d2 = fst $ mul2w_ r p d1 d2++mul2i :: RoundMode -> Precision -> Dyadic -> Int -> Dyadic+mul2i r p d1 d2 = fst $ mul2i_ r p d1 d2++div2w :: RoundMode -> Precision -> Dyadic -> Word -> Dyadic+div2w r p d1 d2 = fst $ div2w_ r p d1 d2++div2i :: RoundMode -> Precision -> Dyadic -> Int -> Dyadic+div2i r p d1 d2 = fst $ div2i_ r p d1 d2++mul2w_ :: RoundMode -> Precision -> Dyadic -> Word -> (Dyadic, Int)+mul2w_ r p d1 d2 = withDyadicBAui r p d1 (fromIntegral d2) mpfr_mul_2ui++mul2i_ :: RoundMode -> Precision -> Dyadic -> Int -> (Dyadic, Int)+mul2i_ r p d1 d2 = withDyadicBAsi r p d1 (fromIntegral d2) mpfr_mul_2si++div2w_ :: RoundMode -> Precision -> Dyadic -> Word -> (Dyadic, Int)+div2w_ r p d1 d2 = withDyadicBAui r p d1 (fromIntegral d2) mpfr_div_2ui++div2i_ :: RoundMode -> Precision -> Dyadic -> Int -> (Dyadic, Int)+div2i_ r p d1 d2 = withDyadicBAsi r p d1 (fromIntegral d2) mpfr_div_2si++----------------------------------------------------------++-- x * 2 ^ y+int2i :: RoundMode -> Precision -> Int -> Int -> Dyadic+int2i r p i e = fst $ int2i_ r p i e++int2w :: RoundMode -> Precision -> Word -> Int -> Dyadic+int2w r p i e = fst $ int2w_ r p i e++int2i_ :: RoundMode -> Precision -> Int -> Int -> (Dyadic, Int)+int2i_ r p i e = unsafePerformIO go+ where go = do ls <- mpfr_custom_get_size (fromIntegral p)+ fp <- mallocForeignPtrBytes (fromIntegral ls)+ let dummy = MP (fromIntegral p) 0 0 fp+ with dummy $ \p1 -> do+ r2 <- mpfr_set_si_2exp p1 (fromIntegral i) (fromIntegral e) ((fromIntegral . fromEnum) r)+ r1 <- peekP p1 fp+ return (r1, fromIntegral r2)++int2w_ :: RoundMode -> Precision -> Word -> Int -> (Dyadic, Int)+int2w_ r p i e = unsafePerformIO go+ where go = do ls <- mpfr_custom_get_size (fromIntegral p)+ fp <- mallocForeignPtrBytes (fromIntegral ls)+ let dummy = MP (fromIntegral p) 0 0 fp+ with dummy $ \p1 -> do+ r2 <- mpfr_set_ui_2exp p1 (fromIntegral i) (fromIntegral e) ((fromIntegral . fromEnum) r)+ r1 <- peekP p1 fp+ return (r1, fromIntegral r2)+ +----------------------------------------------------------+++-- Return d1 * d2 + d3+fma :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic -> Dyadic+fma r p d1 d2 d3 = fst $ fma_ r p d1 d2 d3++-- Return d1 * d2 - d3 +fms :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic -> Dyadic+fms r p d1 d2 d3 = fst $ fms_ r p d1 d2 d3+++fma_ :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic -> (Dyadic, Int)+fma_ r p mp1 mp2 mp3 = unsafePerformIO go+ where go = do ls <- mpfr_custom_get_size (fromIntegral p)+ fp <- mallocForeignPtrBytes (fromIntegral ls)+ let dummy = MP (fromIntegral p) 0 0 fp+ with dummy $ \p1 -> do+ with mp1 $ \p2 -> do + with mp2 $ \p3 -> do + with mp3 $ \p4 -> do + r2 <- mpfr_fma p1 p2 p3 p4 ((fromIntegral . fromEnum) r) + r1 <- peekP p1 fp + return (r1, fromIntegral r2)+++fms_ :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic -> (Dyadic, Int)+fms_ r p mp1 mp2 mp3 = unsafePerformIO go+ where go = do ls <- mpfr_custom_get_size (fromIntegral p)+ fp <- mallocForeignPtrBytes (fromIntegral ls)+ let dummy = MP (fromIntegral p) 0 0 fp+ with dummy $ \p1 -> do+ with mp1 $ \p2 -> do + with mp2 $ \p3 -> do + with mp3 $ \p4 -> do + r2 <- mpfr_fms p1 p2 p3 p4 ((fromIntegral . fromEnum) r) + r1 <- peekP p1 fp+ return (r1, fromIntegral r2)++nextBelow :: Dyadic -> Dyadic+nextBelow mp1 = unsafePerformIO go+ where go = do let p = fromIntegral (getPrec mp1)+ ls <- mpfr_custom_get_size p+ fp <- mallocForeignPtrBytes (fromIntegral ls)+ let dummy = MP p 0 0 fp+ with dummy $ \p1 -> do+ with mp1 $ \p2 -> do + _ <- mpfr_set p1 p2 ((fromIntegral . fromEnum) Near) + mpfr_nextbelow p1 + peekP p1 fp++----------------------------------------------------------+-- powers++sqr :: RoundMode -> Precision -> Dyadic -> Dyadic +sqr r p d = fst $ sqr_ r p d++sqrt :: RoundMode -> Precision -> Dyadic -> Dyadic+sqrt r p d = fst $ sqrt_ r p d++root :: RoundMode -> Precision -> Dyadic -> Word -> Dyadic+root r p d n = fst $ root_ r p d n++pow :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic+pow r p d1 d2 = fst $ pow_ r p d1 d2 ++poww :: RoundMode -> Precision -> Dyadic -> Word -> Dyadic +poww r p d1 d2 = fst $ poww_ r p d1 d2++powi :: RoundMode -> Precision -> Dyadic -> Int -> Dyadic +powi r p d1 d2 = fst $ powi_ r p d1 d2++wpoww :: RoundMode -> Precision -> Word -> Word -> Dyadic +wpoww r p d1 d2 = fst $ wpoww_ r p d1 d2++wpow :: RoundMode -> Precision -> Word -> Dyadic -> Dyadic +wpow r p d1 d2 = fst $ wpow_ r p d1 d2++sqr_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)+sqr_ r p d = withDyadic r p d mpfr_sqr++sqrt_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)+sqrt_ r p d = withDyadic r p d mpfr_sqrt+ +root_ :: RoundMode -> Precision -> Dyadic -> Word -> (Dyadic, Int)+root_ r p d n = withDyadicBAui r p d (fromIntegral n) mpfr_root++pow_ :: RoundMode -> Precision -> Dyadic -> Dyadic -> (Dyadic, Int)+pow_ r p d1 d2 = withDyadicsBA r p d1 d2 mpfr_pow ++poww_ :: RoundMode -> Precision -> Dyadic -> Word -> (Dyadic , Int)+poww_ r p d1 d2 = withDyadicBAui r p d1 (fromIntegral d2) mpfr_pow_ui++powi_ :: RoundMode -> Precision -> Dyadic -> Int -> (Dyadic , Int)+powi_ r p d1 d2 = withDyadicBAsi r p d1 (fromIntegral d2) mpfr_pow_si++wpoww_ :: RoundMode -> Precision -> Word -> Word -> (Dyadic , Int)+wpoww_ r p d1 d2 = unsafePerformIO go+ where go = do ls <- mpfr_custom_get_size (fromIntegral p)+ fp <- mallocForeignPtrBytes (fromIntegral ls)+ let dummy = MP (fromIntegral p) 0 0 fp+ with dummy $ \p1 -> do + r2 <- mpfr_ui_pow_ui p1 (fromIntegral d1) (fromIntegral d2) ((fromIntegral . fromEnum) r)+ r1 <- peekP p1 fp+ return (r1, fromIntegral r2)+ +wpow_ :: RoundMode -> Precision -> Word -> Dyadic -> (Dyadic , Int)+wpow_ r p d1 d2 = withDyadicBAiu r p (fromIntegral d1) d2 mpfr_ui_pow++-----------------------------------------------------------++-- transcendental functions++exp :: RoundMode -> Precision -> Dyadic -> Dyadic+exp r p d = fst $ exp_ r p d++exp2 :: RoundMode -> Precision -> Dyadic -> Dyadic+exp2 r p d = fst $ exp2_ r p d++exp10 :: RoundMode -> Precision -> Dyadic -> Dyadic+exp10 r p d = fst $ exp10_ r p d++log :: RoundMode -> Precision -> Dyadic -> Dyadic+log r p d = fst $ log_ r p d++log2 :: RoundMode -> Precision -> Dyadic -> Dyadic+log2 r p d = fst $ log2_ r p d++log10 :: RoundMode -> Precision -> Dyadic -> Dyadic+log10 r p d = fst $ log10_ r p d++sinh :: RoundMode -> Precision -> Dyadic -> Dyadic+sinh r p d = fst $ sinh_ r p d++cosh :: RoundMode -> Precision -> Dyadic -> Dyadic+cosh r p d = fst $ cosh_ r p d++tanh :: RoundMode -> Precision -> Dyadic -> Dyadic+tanh r p d = fst $ tanh_ r p d ++exp_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)+exp_ r p d = withDyadic r p d mpfr_exp++exp2_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)+exp2_ r p d = withDyadic r p d mpfr_exp2++exp10_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)+exp10_ r p d = withDyadic r p d mpfr_exp10++log_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)+log_ r p d = withDyadic r p d mpfr_log++log2_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)+log2_ r p d = withDyadic r p d mpfr_log2++log10_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)+log10_ r p d = withDyadic r p d mpfr_log10++sinh_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)+sinh_ r p d = withDyadic r p d mpfr_sinh++cosh_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)+cosh_ r p d = withDyadic r p d mpfr_cosh++tanh_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)+tanh_ r p d = withDyadic r p d mpfr_tanh++------------------------------------------------------------++neg :: RoundMode -> Precision -> Dyadic -> Dyadic+neg r p d = fst $ neg_ r p d++absD :: RoundMode -> Precision -> Dyadic -> Dyadic +absD r p d = fst $ absD_ r p d++dim :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic+dim r p d1 d2 = fst $ dim_ r p d1 d2 ++neg_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)+neg_ r p d = withDyadic r p d mpfr_neg++absD_ :: RoundMode -> Precision -> Dyadic -> (Dyadic , Int)+absD_ r p d = withDyadic r p d mpfr_abs++dim_ :: RoundMode -> Precision -> Dyadic -> Dyadic -> (Dyadic, Int)+dim_ r p d1 d2 = withDyadicsBA r p d1 d2 mpfr_dim++--------------------------------------------------------++-- comparison functions++isNaN :: Dyadic -> Bool+isNaN d = withDyadicB d mpfr_nan_p /= 0++isInfinite :: Dyadic -> Bool+isInfinite d = withDyadicB d mpfr_inf_p /= 0 ++isNumber :: Dyadic -> Bool+isNumber d = withDyadicB d mpfr_number_p /= 0 ++isZero :: Dyadic -> Bool+isZero d = withDyadicB d mpfr_zero_p /= 0++greater :: Dyadic -> Dyadic -> Bool+greater d1 d2 = withDyadicBB d1 d2 mpfr_greater_p /= 0++greatereq :: Dyadic -> Dyadic -> Bool+greatereq d1 d2 = withDyadicBB d1 d2 mpfr_greaterequal_p /= 0++less :: Dyadic -> Dyadic -> Bool+less d1 d2 = withDyadicBB d1 d2 mpfr_less_p /= 0++lesseq :: Dyadic -> Dyadic -> Bool+lesseq d1 d2 = withDyadicBB d1 d2 mpfr_lessequal_p /= 0++{- sets erange flag+lessgreater :: Dyadic -> Dyadic -> Bool+lessgreater d1 d2 = withDyadicBB d1 d2 mpfr_lessgreater_p /= 0+-}++equal :: Dyadic -> Dyadic -> Bool+equal d1 d2 = withDyadicBB d1 d2 mpfr_equal_p /= 0+++maxD :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic+maxD r p d1 d2 = fst $ maxD_ r p d1 d2++minD :: RoundMode -> Precision -> Dyadic -> Dyadic -> Dyadic+minD r p d1 d2 = fst $ minD_ r p d1 d2++maxD_ :: RoundMode -> Precision -> Dyadic -> Dyadic -> (Dyadic, Int)+maxD_ r p d1 d2 = withDyadicsBA r p d1 d2 mpfr_max++minD_ :: RoundMode -> Precision -> Dyadic -> Dyadic -> (Dyadic, Int)+minD_ r p d1 d2 = withDyadicsBA r p d1 d2 mpfr_min++sgn :: Dyadic -> Int +sgn d = case compare zero d of+ LT -> 1+ EQ -> 0+ _ -> -1++--conversion from dyadics to basic haskell types++dyadicToDouble :: RoundMode -> Dyadic -> Double+dyadicToDouble r mp1 = (realToFrac . unsafePerformIO) go+ where go = with mp1 $ \p -> mpfr_get_d p ((fromIntegral . fromEnum) r)++dyadicToWord :: RoundMode -> Dyadic -> Word+dyadicToWord r mp1 = (fromIntegral . unsafePerformIO) go+ where go = with mp1 $ \p -> mpfr_get_ui p ((fromIntegral . fromEnum) r)++dyadicToInt :: RoundMode -> Dyadic -> Int+dyadicToInt r mp1 = (fromIntegral . unsafePerformIO) go+ where go = with mp1 $ \p -> mpfr_get_si p ((fromIntegral . fromEnum) r)++dyadicToString :: RoundMode -> Word -- ^ number of significant digits + -> Word -- ^ base + -> Dyadic -> (String, Int)+dyadicToString r n b mp1 = unsafePerformIO go + where go = with mp1 $ \p1 -> do + alloca $ \p2 -> do + withCString (replicate (fromIntegral (n + 2)) '0') $ \p3 -> do + _ <- mpfr_get_str p3 p2 (fromIntegral b) (fromIntegral n) p1 ((fromIntegral . fromEnum) r)+ r1 <- peekCString p3 + r2 <- peek p2 + return (r1, fromIntegral r2) ++decompose :: Dyadic -> (Integer, Int)+decompose d = (getMantissa d, getExp d - ceiling (fromIntegral (getPrec d) / fromIntegral bitsPerMPLimb :: Double) * bitsPerMPLimb)+++toStringExp :: Word -> Dyadic -> String+toStringExp dec d = s ++ case e > 0 of+ True -> case floor (logBase 10 2 * fromIntegral (getExp d) :: Double) > dec of+ False -> take e ss ++ let bt = backtrim (drop e ss) in if null bt then "" else "." ++ bt+ True -> head ss : "." ++ let bt = (backtrim . tail) ss in if null bt then "0"+ else bt ++ "e" ++ show (pred e)+ False -> head ss : "." ++ (let bt = (backtrim . tail) ss in+ if null bt then "0" + else bt )+ ++ "e" ++ show (pred e)+ where (str, e) = dyadicToString Near n 10 d+ n = max dec 5+ (s, ss) = case head str of+ '-' -> ("-", tail str)+ _ -> ("" , str)+ backtrim = reverse . dropWhile (== '0') . reverse ++toString :: Word -> Dyadic -> String+toString dec d = s ++ case compare 0 e of+ LT -> take e ss ++ (let bt = all (== '0') (drop e ss) in if bt then "" else '.' : (drop e ss))+ ++ (if fromIntegral n - e < 0 then "e" ++ show (e - fromIntegral n) else "")+ GT -> let ee = fromIntegral dec + e in + if ee <= 0 then "0" else + head ss : "." ++ (backtrim . tail . take ee) ss ++ "e" ++ show (pred e)+ EQ -> "0." ++ let bt = all (== '0') ss in if bt then "0" else ss+ where (str, e) = dyadicToString Near n 10 d+ n = max dec 5+ (s, ss) = case head str of+ '-' -> ("-", tail str)+ _ -> ("" , str)+ backtrim = reverse . dropWhile (== '0') . reverse ++--asignment functions ++set :: RoundMode -> Precision -> Dyadic -> Dyadic+set r p d = fst $ set_ r p d+++set_ :: RoundMode -> Precision -> Dyadic -> (Dyadic, Int)+set_ r p mp1 = unsafePerformIO go+ where go = do ls <- mpfr_custom_get_size (fromIntegral p)+ fp <- mallocForeignPtrBytes (fromIntegral ls)+ let dummy = MP (fromIntegral p) 0 0 fp+ with dummy $ \p1 -> do+ with mp1 $ \p2 -> do + r2 <- mpfr_set p1 p2 ((fromIntegral . fromEnum) r) + r1 <- peekP p1 fp+ return (r1, fromIntegral r2)+ +------------------------------------+-- mpfr constants+pi :: RoundMode -> Precision -> Dyadic+pi r p = fst $ pi_ r p++log2c :: RoundMode -> Precision -> Dyadic+log2c r p = fst $ pi_ r p++euler :: RoundMode -> Precision -> Dyadic+euler r p = fst $ pi_ r p++catalan :: RoundMode -> Precision -> Dyadic+catalan r p = fst $ pi_ r p++pi_ :: RoundMode -> Precision -> (Dyadic, Int)+pi_ r p = withDyadicC r p mpfr_const_pi++log2c_ :: RoundMode -> Precision -> (Dyadic, Int)+log2c_ r p = withDyadicC r p mpfr_const_log2++euler_ :: RoundMode -> Precision -> (Dyadic, Int)+euler_ r p = withDyadicC r p mpfr_const_euler++catalan_ :: RoundMode -> Precision -> (Dyadic, Int)+catalan_ r p = withDyadicC r p mpfr_const_catalan++----------------------------------------------------------+-- conversion from basic haskell types to dyadics++fromDouble :: RoundMode -> Precision -> Double -> Dyadic+fromDouble r p d = fst $ fromDouble_ r p d++fromInt :: RoundMode -> Precision -> Int -> Dyadic+fromInt r p d = fst $ fromInt_ r p d++fromWord :: RoundMode -> Precision -> Word -> Dyadic+fromWord r p d = fst $ fromWord_ r p d++fromDouble_ :: RoundMode -> Precision -> Double -> (Dyadic, Int)+fromDouble_ r p d = unsafePerformIO go+ where go = do ls <- mpfr_custom_get_size (fromIntegral p)+ fp <- mallocForeignPtrBytes (fromIntegral ls)+ let dummy = MP (fromIntegral p) 0 0 fp+ with dummy $ \p1 -> do + r2 <- mpfr_set_d p1 (realToFrac d) ((fromIntegral . fromEnum) r)+ r1 <- peekP p1 fp+ return (r1, fromIntegral r2)++fromInt_ :: RoundMode -> Precision -> Int -> (Dyadic, Int)+fromInt_ r p d = unsafePerformIO go+ where go = do ls <- mpfr_custom_get_size (fromIntegral p)+ fp <- mallocForeignPtrBytes (fromIntegral ls)+ let dummy = MP (fromIntegral p) 0 0 fp+ with dummy $ \p1 -> do + r2 <- mpfr_set_si p1 (fromIntegral d) ((fromIntegral . fromEnum) r)+ r1 <- peekP p1 fp+ return (r1, fromIntegral r2)+++fromWord_ :: RoundMode -> Precision -> Word -> (Dyadic, Int)+fromWord_ r p d = unsafePerformIO go+ where go = do ls <- mpfr_custom_get_size (fromIntegral p)+ fp <- mallocForeignPtrBytes (fromIntegral ls)+ let dummy = MP (fromIntegral p) 0 0 fp+ with dummy $ \p1 -> do + r2 <- mpfr_set_ui p1 (fromIntegral d) ((fromIntegral . fromEnum) r)+ r1 <- peekP p1 fp+ return (r1, fromIntegral r2) +++fromIntegerA :: RoundMode -> Precision -> Integer -> Dyadic+fromIntegerA r p d = stringToDyadic r p 10 (show d)++compose :: RoundMode -> Precision -> (Integer, Int) -> Dyadic +compose r p (i, ii) = div2i r p (fromIntegerA r p i) ii++fromString :: String -> Precision -> Word -> Dyadic+fromString s p b = stringToDyadic Near p b s++---------------------------------------------------------++-- functions getting properties of dyadics++getPrec :: Dyadic -> Precision+getPrec d = fromIntegral (withDyadicP d mpfr_get_prec)++-- | getMantissa and getExp return values such that+--+-- > d = getMantissa d * 2^(getExp d - ceiling ((getPrec d) / bitsPerMPLimb)* bitsPerMPLimb )+getMantissa :: Dyadic -> Integer+getMantissa d = if d < zero then -h else h+ where (h, _) = foldl (\(a,b) c -> (a + (toInteger c) `shiftL` b, b + bitsPerMPLimb)) (0,0) (getMantissa' d) ++getExp :: Dyadic -> Int+getExp d = (fromIntegral . unsafePerformIO) go+ where go = do with d $ \p1 -> + mpfr_custom_get_exp p1++--------------------------------------------------------++-- some constants+minPrec :: Precision+minPrec = 32++one :: Dyadic +one = fromWord Near minPrec 1++zero :: Dyadic +zero = fromWord Near minPrec 0++-- instances++instance Eq Dyadic where+ (==) = equal++instance Ord Dyadic where+ (<) = less+ (<=) = lesseq+ (>) = greater+ (>=) = greatereq+ +instance Show Dyadic where+ show = toStringExp 16++-- these are exact operations, without rounding+instance Num Dyadic where+ d + d' = add Zero (addPrec d d') d d'+ d - d' = sub Zero (addPrec d d') d d'+ d * d' = mul Zero (getPrec d + getPrec d') d d'+ negate d = neg Zero (getPrec d) d + signum d = case compare d zero of + LT -> negate one+ EQ -> zero+ _ -> one+ abs d = absD Zero (getPrec d) d+ fromInteger i = fromIntegerA Zero (checkPrec $ binprec i) i++addPrec :: Dyadic -> Dyadic -> Precision+addPrec d1 d2 = fromIntegral (max (p1 + e1 - e3) (p2 + e2 - e3)) + 1+ where e1 = if d1 == 0 then 0 else getExp d1+ e2 = if d2 == 0 then 0 else getExp d2+ p1 = fromIntegral $ getPrec d1+ p2 = fromIntegral $ getPrec d2+ e3 = min e1 e2++{-+addPrec d1 d2 = max e1 e2 + 1 - min (e1 - p1) (p2 - e2)+ where e1 = if d1 == 0 then 0 else fromIntegral $ getExp d1+ e2 = if d2 == 0 then 0 else fromIntegral $ getExp d2+ p1 = getPrec d1+ p2 = getPrec d2+-}++
+ Data/Number/Real.hs view
@@ -0,0 +1,266 @@+module Data.Number.Real ( + -- | show x will output as much decimalas as+ -- a standard IEEE 754 double if possible.+++ -- | (==) and (/=) should not be used as x == y will diverge if+ -- two reals should be equal.++ CReal(), Nat, Chain,+ PBool (..),+ min, max, + lim, limRec, limRat, infSum, infSumRec,+ approx,+ pCompare, (<.), (>.), sqrt, exp, log,+ fromDyadic, fromInt, fromWord, fromString, toString, toStringDec+ ) where++import qualified Data.Number.DyadicInterval as DI+import qualified Data.Number.Ball as B+import qualified Data.Number.Dyadic as D++import Data.Order++import Data.Word(Word)+import Prelude hiding (min, max, sqrt, log, exp)++import Data.IORef(IORef, newIORef, writeIORef, readIORef)+import System.IO.Unsafe(unsafePerformIO)++import Data.Maybe(isNothing, fromMaybe)++import Data.Ratio(numerator, denominator)++type Nat = Word++type Chain = Nat -> DI.Interval++-- | Real number is represented as a chain of dyadic intervals which+-- are neither necessarily nested nor bounded away from 0.+--+-- On n-th stage computations are performed with precision of n bits.+data CReal = CReal { state :: IORef (Nat, DI.Interval),+ eval :: Nat -> CReal -> DI.Interval }++{-# INLINE make #-}+make :: Chain -> CReal+make c = CReal { state = unsafePerformIO $ newIORef (D.minPrec, c D.minPrec) ,+ eval = \n (CReal s _) -> unsafePerformIO $+ do (n', i) <- readIORef s+ if n' == n then return i+ else do let i' = c n+ writeIORef s (n, i')+ return i'+ } ++{-# INLINE represent #-}+represent :: (D.Precision -> DI.Interval -> DI.Interval) -> CReal -> CReal +represent f r = make (\n -> f n (eval r n r))++{-# INLINE represent2 #-}+represent2 :: (D.Precision -> DI.Interval -> DI.Interval -> DI.Interval)+ -> CReal -> CReal -> CReal+represent2 f r r' = make (\n -> f n (eval r n r) (eval r' n r'))+++max :: CReal -> CReal -> CReal+max = represent2 DI.maxI++min :: CReal -> CReal -> CReal+min = represent2 DI.minI++instance Eq CReal where+ r /= r' = or (map (isNothing . (\n -> DI.intersect (eval r n r) (eval r' n r'))) [1..]) ++instance Show CReal where+ show = toStringDec 16 ++instance Read CReal where+ readsPrec _ s = [(fromString s, "")]++instance Num CReal where+ (+) = represent2 DI.add+ (-) = represent2 DI.sub+ (*) = represent2 DI.mul+ negate = represent DI.neg+ abs r = max r (negate r)+ signum r = make $ \n -> case DI.compareI (eval r n r) (eval 0 n 0) of+ Less -> DI.fromInt D.minPrec (negate 1)+ Greater -> DI.fromInt D.minPrec 1+ _ -> DI.fromBall (B.Ball 0 1)+ fromInteger = fromDyadic . fromInteger++instance Fractional CReal where+ (/) = represent2 DI.div+ recip r = 1 / r+ fromRational r = fromIntegral (numerator r) / fromIntegral (denominator r)++sqrt :: CReal -> CReal+sqrt = represent DI.sqrt++exp :: CReal -> CReal+exp = represent DI.exp++log :: CReal -> CReal+log = represent DI.log++ +-- | A basic general limit which takes as arguments a sequence of reals and a sequence of +-- error bounds. +lim :: (Nat -> CReal) -- ^ Sequence+ -> (Nat -> CReal) -- ^ Error bounds+ -> CReal+lim am rm = make limStage+ where limStage n = foldl1 DI.intersect lst+ where lst = take (fromIntegral n) .+ map (\k -> let n' = if k < div n 2 then k else n+ (a, r) = (am k, rm k) -- get k-th element of the sequence + (an, rn) = (eval a n' a, eval r n' r) -- get the n-th approximation of the k-th element+ i = case (an, rn) of+ (Just b, Just b') -> DI.fromBall (B.Ball (B.center b) (B.radius b + B.upper_ b'))+ _ -> Nothing+ in i) $ [1..]++-- | Similar to lim, but can sometimes be more convenient for some sequences+limRec :: CReal -- ^ initial value+ -> (CReal -> Nat -> (CReal, CReal)) -- ^ a function which produces a pair, (next element, error estimate)+ -- from previous one and location+ -> CReal+limRec st f = make limStage+ where limStage n = limStage' 1 st (eval st n st)+ where limStage' k st' acc = + let (an, rn) = f st' k -- n-th element of the sequence+ (ak, rk) = (eval an n an, eval rn n rn) -- k-th approximation+ i = case (ak, rk) of+ (Just b, Just b') -> DI.fromBall (B.Ball (B.center b) (B.radius b + B.upper_ b'))+ _ -> Nothing+ in if k == n then DI.intersect acc i+ else limStage' (succ k) an (DI.intersect acc i)++-- | Limit of a sequence of rationals.+limRat :: (Nat -> D.Dyadic) -- ^ Sequence of dyadics+ -> (Nat -> D.Dyadic) -- ^ Sequence of (dyadic) error bounds+ -> CReal+limRat an rn = make (\n -> DI.fromBall (B.Ball (an n) (rn n)))+++-- | Computes an infinite sum of a series +infSum :: (Nat -> CReal) -- ^ Sequence of reals+ -> (Nat -> CReal) -- ^ Sequence of series remainders+ -> CReal+infSum am rm = make partialsum+ where partialsum k = psum 1 (eval a0 k a0) Nothing+ where psum n acc res = + let (an,rn) = (am n, rm n)+ err = eval rn k rn+ acc' = DI.add k acc (eval an k an)+ (res', p) = case (acc', err) of+ (Just b, Just b') -> + let (cac,rac) = (B.center b, B.radius b)+ (ler, uer) = (B.lower_ b', B.upper_ b')+ in (DI.intersect res (DI.fromBall (B.Ball cac (rac + uer))), rac <= ler)+ (Nothing, _) -> (Nothing, False)+ (_, Nothing) -> (res, True)+ in if p then psum (succ n) acc' res'+ else res'+ a0 = am 0+ ++-- | Similar to infSum but can sometimes be more convenient+-- Second argument is a_0+infSumRec :: CReal+ -> (CReal -> Nat -> (CReal, CReal)) + -> CReal+infSumRec st f = make partialsum+ where partialsum k = psum 1 (eval a0 k a0) Nothing a0+ where psum n acc res t = + let (an, rn) = f t n+ err = eval rn k rn+ acc' = DI.add k acc (eval an k an)+ (res', p) = case (acc', err) of+ (Just b, Just b') -> + let (cac,rac) = (B.center b, B.radius b)+ (ler, uer) = (B.lower_ b', B.upper_ b')+ in (DI.intersect res (DI.fromBall (B.Ball cac (rac + uer))), rac <= ler)+ (Nothing, _) -> (Nothing, False)+ (_, Nothing) -> (res, True)+ in if p then psum (succ n) acc' res' an+ else res'+ a0 = st++-- comparison functions++-- | @ pCompare x y @ returns a function @ Nat -> POrdering @ which+-- when applied to some @ n @ computes approximates with precision @ n @+-- and then compares the resulting intervals+pCompare :: CReal -> CReal -> Nat -> POrdering+pCompare r r' = \n -> DI.compareI (eval r n r) (eval r' n r')++-- | @ x \<. y @ is a function @ Nat -> PBool @ which, when+-- applied to some @ n @, computes the approximation with precision @ n @ +-- and then compares the intervals. If intervals are disjoint then result is +-- either PTrue or PFalse, otherwise result is Indeterminate.+infix 4 <.+(<.) :: CReal -> CReal -> Nat -> PBool+(<.) r r' = \n -> case pCompare r r' n of+ Less -> PTrue+ Greater -> PFalse+ _ -> Indeterminate++-- | Similar to (<.)+infix 4 >.+(>.) :: CReal -> CReal -> Nat -> PBool+(>.) r r' = \n -> case pCompare r r' n of+ Less -> PFalse+ Greater -> PTrue+ _ -> Indeterminate++-- | @ approx x n @ tries to compute a dyadic approximation to x so than @ |x - d| <= 10^(-n) @+-- If it succeeds it returns @ Right d @ where d is a dyadic rational, otherwise it returns+-- Left (d, n) where d is a dyadic rational and n is the number of accurate decimal places+--+-- Approx succeeds if result can be computed with precision less than the square of the number+-- of required bits of precision.+approx :: CReal -> Nat -> Either (D.Dyadic, Word) D.Dyadic+approx r k = approx' n+ where approx' :: Nat -> Either (D.Dyadic, Word) D.Dyadic+ approx' n' | cp >= fromIntegral n = Right c+ | threshold = Left (c, floor (logBase 10 2 * fromIntegral cp :: Double))+ | otherwise = approx' $ 2 * n'+ where cp = if r' == 0 then fromIntegral n + else let t = negate . D.getExp $ r'+ in if t >= 0 then t else 0+ B.Ball c r' = fromMaybe (B.Ball 0 (D.pow2 31)) (eval r n' r)+ threshold = n * n < n'+ n = ceiling ((logBase 2 10 :: Double) * fromIntegral k) + 1++fromDyadic :: D.Dyadic -> CReal+fromDyadic d = make $ \_ -> DI.fromBall (B.Ball d $ 0)++-- | fromInt should be preferred over fromIntegral where applicable+fromInt :: Int -> CReal+fromInt i = make $ \_ -> DI.fromBall $ B.Ball (D.fromInt D.Near 32 i) $ 0++-- | fromWord should be preferred over fromIntegral where applicable+fromWord :: Word -> CReal+fromWord i = make $ \_ -> DI.fromBall $ B.Ball (D.fromWord D.Near 32 i) $ 0+++fromString :: String -> CReal+fromString s = make (\_ -> let l = length s+ n = ceiling (logBase 2 10 * fromIntegral (if elem '.' s then pred l else l) :: Double)+ cen = D.fromString s n 10 in + DI.fromBall (B.Ball cen 0))+++-- | toStringDec tries to compute the result to the number of specified significand digits+toStringDec :: Nat -> CReal -> String +toStringDec n r = inf ++ s+ where (inf, s) = case approx r n of+ Right d -> ("",D.toString n d)+ Left (d,k) -> ("Could not compute to desired accuracy, only to " ++ show k ++ " significand digits : ",+ D.toString k d) ++-- | toString computes the result with specified precision.+toString :: Nat -> CReal -> String+toString n r = DI.toString (eval r n r)
+ Data/Order.hs view
@@ -0,0 +1,10 @@+module Data.Order where++-- | Partial ordering.+data POrdering = Less | Greater | Incomparable deriving (Eq, Show, Ord)++-- | Partial booleans+data PBool = PTrue -- ^ equivalent to True+ | PFalse -- ^ equivalent to False+ | Indeterminate -- ^ neither True nor False.+ deriving(Eq, Show, Ord)
+ HERA.cabal view
@@ -0,0 +1,29 @@+Name: HERA+Version: 0.2+Cabal-Version: >= 1.2+License: BSD3+License-file: LICENSE+Stability: experimental+Category: Math+Tested-With: GHC==6.8.2 GHC==6.8.3+Build-Depends: base+Exposed-modules: Data.Number.DyadicInterval+ Data.Number.Dyadic+ Data.Number.Ball+ Data.Number.Real+ Data.Order+ Data.Number.MPFR+ +Other-modules: Data.Number.FFIhelper+++build-type: Simple+build-tools: hsc2hs+GHC-options: -Wall -O2+include-dirs: C+includes: hsmpfr.h mpfr.h+install-includes: hsmpfr.h+c-sources: C/hsmpfr.c++extra-lib-dirs:+extra-libraries: mpfr
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) Ales Bizjak <ales.bizjak@student.fmf.uni-lj.si> 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:++1. Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++2. 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.++3. Neither the name of the author nor the names of his contributors+ may be used to endorse or promote products derived from this software+ without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``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 AUTHORS 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.
+ Setup.hs view
@@ -0,0 +1,6 @@+module Main where++import Distribution.Simple++main :: IO ()+main = defaultMain