diff --git a/C/hsmpfr.c b/C/hsmpfr.c new file mode 100644 --- /dev/null +++ b/C/hsmpfr.c @@ -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); +} +*/ diff --git a/C/hsmpfr.h b/C/hsmpfr.h new file mode 100644 --- /dev/null +++ b/C/hsmpfr.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/Data/Number/Ball.hs b/Data/Number/Ball.hs new file mode 100644 --- /dev/null +++ b/Data/Number/Ball.hs @@ -0,0 +1,354 @@ +{-# INCLUDE #-} +{-# INCLUDE #-} + +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" diff --git a/Data/Number/Dyadic.hs b/Data/Number/Dyadic.hs new file mode 100644 --- /dev/null +++ b/Data/Number/Dyadic.hs @@ -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 diff --git a/Data/Number/DyadicInterval.hs b/Data/Number/DyadicInterval.hs new file mode 100644 --- /dev/null +++ b/Data/Number/DyadicInterval.hs @@ -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 diff --git a/Data/Number/FFIhelper.hsc b/Data/Number/FFIhelper.hsc new file mode 100644 --- /dev/null +++ b/Data/Number/FFIhelper.hsc @@ -0,0 +1,725 @@ +{-# LANGUAGE ForeignFunctionInterface #-} +#include +#include + +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 () + +------------------------------------------------- + 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,877 @@ +{-# LANGUAGE TypeSynonymInstances #-} +{-# INCLUDE #-} +{-# INCLUDE #-} +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 +-} + + diff --git a/Data/Number/Real.hs b/Data/Number/Real.hs new file mode 100644 --- /dev/null +++ b/Data/Number/Real.hs @@ -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) diff --git a/Data/Order.hs b/Data/Order.hs new file mode 100644 --- /dev/null +++ b/Data/Order.hs @@ -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) diff --git a/HERA.cabal b/HERA.cabal new file mode 100644 --- /dev/null +++ b/HERA.cabal @@ -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 diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) Ales 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: + +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. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,6 @@ +module Main where + +import Distribution.Simple + +main :: IO () +main = defaultMain