qd 0.2 → 0.4
raw patch · 8 files changed
+98/−114 lines, 8 files
Files
- LICENSE +1/−1
- Numeric/QD.hs +1/−3
- Numeric/QD/Bits.hs +0/−5
- Numeric/QD/Bits/Raw.hs +0/−8
- Numeric/QD/DoubleDouble.hs +42/−25
- Numeric/QD/QuadDouble.hs +47/−26
- Numeric/QD/Raw.hs +1/−3
- qd.cabal +6/−43
LICENSE view
@@ -1,4 +1,4 @@-Copyright (c)2010, Claude Heiland-Allen+Copyright (c)2010,2011, Claude Heiland-Allen All rights reserved.
Numeric/QD.hs view
@@ -1,11 +1,9 @@ module Numeric.QD- ( module Numeric.QD.Bits- , module Numeric.QD.DoubleDouble+ ( module Numeric.QD.DoubleDouble , module Numeric.QD.FPU , module Numeric.QD.QuadDouble ) where -import Numeric.QD.Bits import Numeric.QD.DoubleDouble hiding (toDouble, fromDouble, sqr) import Numeric.QD.FPU import Numeric.QD.QuadDouble hiding (toDouble, fromDouble, toDoubleDouble, fromDoubleDouble, sqr)
− Numeric/QD/Bits.hs
@@ -1,5 +0,0 @@-module Numeric.QD.Bits- (- ) where---- nothing here yet
− Numeric/QD/Bits/Raw.hs
@@ -1,8 +0,0 @@-{-# LANGUAGE ForeignFunctionInterface #-}-module Numeric.QD.Bits.Raw- ( get_double_expn- ) where--import Foreign.C (CDouble, CInt)--foreign import ccall unsafe "qd/bits.h get_double_expn" get_double_expn :: CDouble -> IO CInt
Numeric/QD/DoubleDouble.hs view
@@ -6,10 +6,12 @@ , sqr ) where -import Foreign (Ptr, alloca, allocaArray, castPtr, Storable(..), unsafePerformIO, with)-import Foreign.C (CDouble, CInt, peekCString)-import Data.Ratio (denominator, numerator)+import Foreign (Ptr, alloca, castPtr, Storable(..), unsafePerformIO, with)+import Foreign.C (CDouble, CInt)+import Data.Ratio ((%))+import Data.Bits (shiftL, shiftR) import Data.Typeable (Typeable(..))+import Numeric (showFloat, readSigned, readFloat) import Numeric.QD.DoubleDouble.Raw ( c_dd_add@@ -33,10 +35,10 @@ , c_dd_acosh , c_dd_atanh , c_dd_comp- , c_dd_swrite , c_dd_neg , c_dd_abs- , c_dd_copy_d+ , c_dd_aint+ , c_dd_nint , c_dd_ceil , c_dd_floor , c_dd_atan2@@ -62,11 +64,11 @@ return $ i `compare` (0 :: CInt) instance Show DoubleDouble where- show !a = unsafePerformIO $ with a $ \p -> allocaArray (fromIntegral l) $ \r -> do- c_dd_swrite (castPtr p) (l`div`2) (castPtr r) l- peekCString r- where l = 64- + show = flip showFloat ""++instance Read DoubleDouble where+ readsPrec _ = readSigned readFloat+ instance Num DoubleDouble where (+) = lift_dd_dd_dd c_dd_add (*) = lift_dd_dd_dd c_dd_mul@@ -74,7 +76,7 @@ negate = lift_dd_dd c_dd_neg abs = lift_dd_dd c_dd_abs signum !a = case a `compare` 0 of { LT -> -1 ; EQ -> 0 ; GT -> 1 }- fromInteger !i = unsafePerformIO $ alloca $ \r -> c_dd_copy_d (fromInteger i) (castPtr r) >> peek r+ fromInteger !i = fromRational (i % 1) sqr :: DoubleDouble -> DoubleDouble sqr = lift_dd_dd c_dd_sqr@@ -82,17 +84,21 @@ instance Fractional DoubleDouble where (/) = lift_dd_dd_dd c_dd_div recip !b = 1 / b- fromRational !k = let { !n = numerator k ; !d = denominator k } in fromInteger n / fromInteger d+ fromRational !k = let a = fromRational k+ k' = k - toRational a+ b = fromRational k'+ in DoubleDouble a b instance Real DoubleDouble where toRational !(DoubleDouble a b) = toRational a + toRational b instance RealFrac DoubleDouble where- properFraction = error "Numeric.QD.DoubleDouble.properFraction: not yet implemented" -- FIXME- truncate = error "Numeric.QD.DoubleDouble.truncate: not yet implemented" -- FIXME- round = error "Numeric.QD.DoubleDouble.round: not yet implemented" -- FIXME- ceiling = ceiling . toDouble . lift_dd_dd c_dd_ceil- floor = floor . toDouble . lift_dd_dd c_dd_floor+ properFraction k = let (n, r) = properFraction (toRational k)+ in (n, fromRational r)+ truncate = truncate . toRational . lift_dd_dd c_dd_aint+ round = round . toRational . lift_dd_dd c_dd_nint+ ceiling = ceiling . toRational . lift_dd_dd c_dd_ceil+ floor = floor . toRational . lift_dd_dd c_dd_floor instance Floating DoubleDouble where pi = unsafePerformIO $ alloca $ \r -> c_dd_pi (castPtr r) >> peek r@@ -116,11 +122,25 @@ floatRadix _ = 2 floatDigits _ = 2 * floatDigits (error "Numeric.QD.DoubleDouble.floatDigits" :: CDouble) floatRange _ = floatRange (error "Numeric.QD.DoubleDouble.floatRange" :: CDouble)- decodeFloat = error "Numeric.QD.DoubleDouble.decodeFloat: not yet implemented" -- FIXME- encodeFloat = error "Numeric.QD.DoubleDouble.encodeFloat: not yet implemented" -- FIXME- exponent _ = error "Numeric.QD.DoubleDouble.exponent: not yet implemented" -- FIXME- significand _ = error "Numeric.QD.DoubleDouble.significand: not yet implemented" -- FIXME- scaleFloat !n !x = x * 2 ** fromIntegral n+ decodeFloat !x = case toRational x of+ 0 -> (0, 0)+ r -> let k = floor $ fromIntegral ff - logBase (fromIntegral $ floatRadix x) (abs x)+ i = round $ (fromIntegral $ floatRadix x) ^ k * r+ fixup m e =+ if abs m < mMin+ then fixup (m `shiftL` 1) (e - 1)+ else if abs m >= mMax+ then fixup (m `shiftR` 1) (e + 1)+ else (m, e)+ mMin = 1 `shiftL` (ff - 1)+ mMax = 1 `shiftL` ff+ ff = floatDigits x+ g = -k+ in fixup i g+ encodeFloat m e = scaleFloat e (fromInteger m)+ -- exponent _ = -- use default implementation+ -- significand _ = -- use default implementation+ scaleFloat !n !(DoubleDouble a b) = DoubleDouble (scaleFloat n a) (scaleFloat n b) isNaN !(DoubleDouble a b) = isNaN a || isNaN b isInfinite !(DoubleDouble a b) = isInfinite a || isInfinite b isDenormalized !(DoubleDouble a b) = isDenormalized a || isDenormalized b@@ -129,7 +149,6 @@ atan2 = lift_dd_dd_dd c_dd_atan2 -- instance Enum DoubleDouble -- FIXME--- instance Read DoubleDouble -- FIXME instance Storable DoubleDouble where sizeOf _ = 2 * sizeOf (error "Numeric.QD.DoubleDouble.sizeOf" :: CDouble)@@ -145,9 +164,7 @@ pokeElemOff q 1 b lift_dd_dd :: (Ptr CDouble -> Ptr CDouble -> IO ()) -> DoubleDouble -> DoubleDouble-{-# INLINE lift_dd_dd #-} lift_dd_dd f !a = unsafePerformIO $ with a $ \p -> alloca $ \r -> f (castPtr p) (castPtr r) >> peek r lift_dd_dd_dd :: (Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()) -> DoubleDouble -> DoubleDouble -> DoubleDouble-{-# INLINE lift_dd_dd_dd #-} lift_dd_dd_dd f !a !b = unsafePerformIO $ with a $ \p -> with b $ \q -> alloca $ \r -> f (castPtr p) (castPtr q) (castPtr r) >> peek r
Numeric/QD/QuadDouble.hs view
@@ -8,10 +8,12 @@ , sqr ) where -import Foreign (Ptr, alloca, allocaArray, castPtr, Storable(..), unsafePerformIO, with)-import Foreign.C (CDouble, CInt, peekCString)-import Data.Ratio (denominator, numerator)+import Foreign (Ptr, alloca, castPtr, Storable(..), unsafePerformIO, with)+import Foreign.C (CDouble, CInt)+import Data.Bits (shiftL, shiftR)+import Data.Ratio ((%)) import Data.Typeable (Typeable(..))+import Numeric (showFloat, readSigned, readFloat) import Numeric.QD.DoubleDouble (DoubleDouble(DoubleDouble)) import Numeric.QD.QuadDouble.Raw@@ -36,10 +38,10 @@ , c_qd_acosh , c_qd_atanh , c_qd_comp- , c_qd_swrite , c_qd_neg , c_qd_abs- , c_qd_copy_d+ , c_qd_aint+ , c_qd_nint , c_qd_ceil , c_qd_floor , c_qd_atan2@@ -70,12 +72,6 @@ !i <- peek r return $ i `compare` (0 :: CInt) -instance Show QuadDouble where- show !a = unsafePerformIO $ with a $ \p -> allocaArray (fromIntegral l) $ \r -> do- c_qd_swrite (castPtr p) (l`div`2) (castPtr r) l- peekCString r- where l = 128- instance Num QuadDouble where (+) = lift_qd_qd_qd c_qd_add (*) = lift_qd_qd_qd c_qd_mul@@ -83,7 +79,7 @@ negate = lift_qd_qd c_qd_neg abs = lift_qd_qd c_qd_abs signum !a = case a `compare` 0 of { LT -> -1 ; EQ -> 0 ; GT -> 1 }- fromInteger !i = unsafePerformIO $ alloca $ \r -> c_qd_copy_d (fromInteger i) (castPtr r) >> peek r+ fromInteger !i = fromRational (i % 1) sqr :: QuadDouble -> QuadDouble sqr = lift_qd_qd c_qd_sqr@@ -91,17 +87,25 @@ instance Fractional QuadDouble where (/) = lift_qd_qd_qd c_qd_div recip !b = 1 / b- fromRational !k = let { !n = numerator k ; !d = denominator k } in fromInteger n / fromInteger d+ fromRational !k = let a = fromRational k+ ka = k - toRational a+ b = fromRational ka+ kb = ka - toRational b+ c = fromRational kb+ kc = kb - toRational c+ d = fromRational kc+ in QuadDouble a b c d instance Real QuadDouble where toRational (QuadDouble a b c d) = toRational a + toRational b + toRational c + toRational d instance RealFrac QuadDouble where- properFraction = error "Numeric.QD.QuadDouble.properFraction: not yet implemented" -- FIXME- truncate = error "Numeric.QD.QuadDouble.truncate: not yet implemented" -- FIXME- round = error "Numeric.QD.QuadDouble.round: not yet implemented" -- FIXME- ceiling = ceiling . toDouble . lift_qd_qd c_qd_ceil- floor = floor . toDouble . lift_qd_qd c_qd_floor+ properFraction k = let (n, r) = properFraction (toRational k)+ in (n, fromRational r)+ truncate = truncate . toRational . lift_qd_qd c_qd_aint+ round = round . toRational . lift_qd_qd c_qd_nint+ ceiling = ceiling . toRational . lift_qd_qd c_qd_ceil+ floor = floor . toRational . lift_qd_qd c_qd_floor instance Floating QuadDouble where pi = unsafePerformIO $ alloca $ \r -> c_qd_pi (castPtr r) >> peek r@@ -125,11 +129,25 @@ floatRadix _ = 2 floatDigits _ = 4 * floatDigits (error "Numeric.QD.QuadDouble.floatDigits" :: CDouble) floatRange _ = floatRange (error "Numeric.QD.QuadDouble.floatRange" :: CDouble)- decodeFloat = error "Numeric.QD.QuadDouble.decodeFloat: not yet implemented" -- FIXME- encodeFloat = error "Numeric.QD.QuadDouble.encodeFloat: not yet implemented" -- FIXME- exponent _ = error "Numeric.QD.QuadDouble.exponent: not yet implemented" -- FIXME- significand _ = error "Numeric.QD.QuadDouble.significand: not yet implemented" -- FIXME- scaleFloat !n !x = x * 2 ** fromIntegral n+ decodeFloat !x = case toRational x of+ 0 -> (0, 0)+ r -> let k = floor $ fromIntegral ff - logBase (fromIntegral $ floatRadix x) (abs x)+ i = round $ (fromIntegral $ floatRadix x) ^ k * r+ fixup m e =+ if abs m < mMin+ then fixup (m `shiftL` 1) (e - 1)+ else if abs m >= mMax+ then fixup (m `shiftR` 1) (e + 1)+ else (m, e)+ mMin = 1 `shiftL` (ff - 1)+ mMax = 1 `shiftL` ff+ ff = floatDigits x+ g = -k+ in fixup i g+ encodeFloat m e = scaleFloat e (fromInteger m)+ -- exponent _ = -- use default implementation+ -- significand _ = -- use default implementation+ scaleFloat !n !(QuadDouble a b c d) = QuadDouble (scaleFloat n a) (scaleFloat n b) (scaleFloat n c) (scaleFloat n d) isNaN !(QuadDouble a b c d) = isNaN a || isNaN b || isNaN c || isNaN d isInfinite !(QuadDouble a b c d) = isInfinite a || isInfinite b || isInfinite c || isInfinite d isDenormalized !(QuadDouble a b c d) = isDenormalized a || isDenormalized b || isDenormalized c || isDenormalized d@@ -138,8 +156,13 @@ atan2 = lift_qd_qd_qd c_qd_atan2 -- instance Enum QuadDouble -- FIXME--- instance Read QuadDouble -- FIXME +instance Show QuadDouble where+ show = flip showFloat ""++instance Read QuadDouble where+ readsPrec _ = readSigned readFloat+ instance Storable QuadDouble where sizeOf _ = 4 * sizeOf (error "Numeric.QD.QuadDouble.sizeOf" :: CDouble) alignment _ = alignment (error "Numeric.QD.QuadDouble.alignment" :: CDouble)@@ -158,9 +181,7 @@ pokeElemOff q 3 d lift_qd_qd :: (Ptr CDouble -> Ptr CDouble -> IO ()) -> QuadDouble -> QuadDouble-{-# INLINE lift_qd_qd #-} lift_qd_qd f !a = unsafePerformIO $ with a $ \p -> alloca $ \r -> f (castPtr p) (castPtr r) >> peek r lift_qd_qd_qd :: (Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()) -> QuadDouble -> QuadDouble -> QuadDouble-{-# INLINE lift_qd_qd_qd #-} lift_qd_qd_qd f !a !b = unsafePerformIO $ with a $ \p -> with b $ \q -> alloca $ \r -> f (castPtr p) (castPtr q) (castPtr r) >> peek r
Numeric/QD/Raw.hs view
@@ -1,11 +1,9 @@ module Numeric.QD.Raw- ( module Numeric.QD.Bits.Raw- , module Numeric.QD.DoubleDouble.Raw+ ( module Numeric.QD.DoubleDouble.Raw , module Numeric.QD.FPU.Raw , module Numeric.QD.QuadDouble.Raw ) where -import Numeric.QD.Bits.Raw import Numeric.QD.DoubleDouble.Raw import Numeric.QD.FPU.Raw import Numeric.QD.QuadDouble.Raw
qd.cabal view
@@ -1,5 +1,5 @@ Name: qd-Version: 0.2+Version: 0.4 Synopsis: double-double and quad-double number type via libqd Description: This package supports both a double-double datatype (approx. 32 decimal digits)@@ -8,46 +8,9 @@ need libqd to be installed. . @'Numeric.QD.DoubleDouble.DoubleDouble'@ and @'Numeric.QD.QuadDouble.QuadDouble'@- are strict tuples of @CDouble@s, with the following instances:- .- * @'Eq'@- .- * @'Floating'@- .- * @'Fractional'@- .- * @'Num'@- .- * @'Ord'@- .- * @'Real'@- .- * @'RealFloat'@- .- * @'RealFrac'@- .- * @'Show'@- .- * @'Storable'@- .- But note that the following functions (while present) are @'undefined'@:- .- * @'properFraction'@- .- * @'truncate'@- .- * @'round'@- .- * @'decodeFloat'@- .- * @'encodeFloat'@- .- * @'exponent'@- .- * @'significand'@- .- Non-crashing implementations of these are planned in a future update to- this package, as well as instances of @Enum@ and @Read@.+ are strict tuples of @CDouble@s, with instances of:+ @'Eq'@, @'Floating'@, @'Fractional'@, @'Num'@, @'Ord'@, @'Read'@, @'Real'@,+ @'RealFloat'@, @'RealFrac'@, @'Show'@, @'Storable'@, @'Typeable'@. . Additional note: libqd depends on 64bit doubles, while some FPU architectures use 80bit. It is highly recommended to compile with -fno-excess-precision and@@ -73,9 +36,9 @@ Library Build-depends: base >= 4 && < 5 Extra-Libraries: qd+ if os(linux)+ Extra-libraries: stdc++ Exposed-modules: Numeric.QD- Numeric.QD.Bits- Numeric.QD.Bits.Raw Numeric.QD.DoubleDouble Numeric.QD.DoubleDouble.Raw Numeric.QD.FPU