qd 0.1 → 0.2
raw patch · 5 files changed
+135/−86 lines, 5 files
Files
- Numeric/QD.hs +2/−2
- Numeric/QD/DoubleDouble.hs +34/−27
- Numeric/QD/QuadDouble.hs +35/−28
- TODO +0/−24
- qd.cabal +64/−5
Numeric/QD.hs view
@@ -6,6 +6,6 @@ ) where import Numeric.QD.Bits-import Numeric.QD.DoubleDouble hiding (toDouble, fromDouble)+import Numeric.QD.DoubleDouble hiding (toDouble, fromDouble, sqr) import Numeric.QD.FPU-import Numeric.QD.QuadDouble hiding (toDouble, fromDouble, toDoubleDouble, fromDoubleDouble)+import Numeric.QD.QuadDouble hiding (toDouble, fromDouble, toDoubleDouble, fromDoubleDouble, sqr)
Numeric/QD/DoubleDouble.hs view
@@ -1,8 +1,9 @@-{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE BangPatterns, DeriveDataTypeable #-} module Numeric.QD.DoubleDouble ( DoubleDouble(DoubleDouble) , toDouble , fromDouble+ , sqr ) where import Foreign (Ptr, alloca, allocaArray, castPtr, Storable(..), unsafePerformIO, with)@@ -39,28 +40,29 @@ , c_dd_ceil , c_dd_floor , c_dd_atan2+ , c_dd_sqr ) -data DoubleDouble = DoubleDouble !CDouble !CDouble deriving Typeable+data DoubleDouble = DoubleDouble {-# UNPACK #-} !CDouble {-# UNPACK #-} !CDouble deriving Typeable toDouble :: DoubleDouble -> Double-toDouble (DoubleDouble a _) = realToFrac a+toDouble !(DoubleDouble a _) = realToFrac a fromDouble :: Double -> DoubleDouble-fromDouble a = DoubleDouble (realToFrac a) 0+fromDouble !a = DoubleDouble (realToFrac a) 0 instance Eq DoubleDouble where- a == b = a `compare` b == EQ- a /= b = a `compare` b /= EQ+ (!a) == (!b) = a `compare` b == EQ+ (!a) /= (!b) = a `compare` b /= EQ instance Ord DoubleDouble where- a `compare` b = unsafePerformIO $ with a $ \p -> with b $ \q -> alloca $ \r -> do- c_dd_comp (castPtr p) (castPtr q) (castPtr r)- i <- peek r- return $ i `compare` (0 :: CInt)+ (!a) `compare` (!b) = unsafePerformIO $ with a $ \p -> with b $ \q -> alloca $ \r -> do+ c_dd_comp (castPtr p) (castPtr q) (castPtr r)+ !i <- peek r+ return $ i `compare` (0 :: CInt) instance Show DoubleDouble where- show a = unsafePerformIO $ with a $ \p -> allocaArray (fromIntegral l) $ \r -> do+ 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@@ -71,16 +73,19 @@ (-) = lift_dd_dd_dd c_dd_sub 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+ 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 +sqr :: DoubleDouble -> DoubleDouble+sqr = lift_dd_dd c_dd_sqr+ 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+ recip !b = 1 / b+ fromRational !k = let { !n = numerator k ; !d = denominator k } in fromInteger n / fromInteger d instance Real DoubleDouble where- toRational (DoubleDouble a b) = toRational a + toRational b+ toRational !(DoubleDouble a b) = toRational a + toRational b instance RealFrac DoubleDouble where properFraction = error "Numeric.QD.DoubleDouble.properFraction: not yet implemented" -- FIXME@@ -115,11 +120,11 @@ 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- isNaN (DoubleDouble a b) = isNaN a || isNaN b- isInfinite (DoubleDouble a b) = isInfinite a || isInfinite b- isDenormalized (DoubleDouble a b) = isDenormalized a || isDenormalized b- isNegativeZero (DoubleDouble a b) = isNegativeZero a || (a == 0 && isNegativeZero b)+ scaleFloat !n !x = x * 2 ** fromIntegral n+ isNaN !(DoubleDouble a b) = isNaN a || isNaN b+ isInfinite !(DoubleDouble a b) = isInfinite a || isInfinite b+ isDenormalized !(DoubleDouble a b) = isDenormalized a || isDenormalized b+ isNegativeZero !(DoubleDouble a b) = isNegativeZero a || (a == 0 && isNegativeZero b) isIEEE _ = False -- FIXME what does this imply? atan2 = lift_dd_dd_dd c_dd_atan2 @@ -129,18 +134,20 @@ instance Storable DoubleDouble where sizeOf _ = 2 * sizeOf (error "Numeric.QD.DoubleDouble.sizeOf" :: CDouble) alignment _ = alignment (error "Numeric.QD.DoubleDouble.alignment" :: CDouble)- peek p = do- let q = castPtr p+ peek !p = do+ let !q = castPtr p a <- peekElemOff q 0 b <- peekElemOff q 1 return $ DoubleDouble a b- poke p (DoubleDouble a b) = do- let q = castPtr p+ poke !p !(DoubleDouble a b) = do+ let !q = castPtr p pokeElemOff q 0 a pokeElemOff q 1 b lift_dd_dd :: (Ptr CDouble -> Ptr CDouble -> IO ()) -> DoubleDouble -> DoubleDouble-lift_dd_dd f a = unsafePerformIO $ with a $ \p -> alloca $ \r -> f (castPtr p) (castPtr r) >> peek r+{-# 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-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+{-# 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
@@ -1,10 +1,11 @@-{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE BangPatterns, DeriveDataTypeable #-} module Numeric.QD.QuadDouble ( QuadDouble(QuadDouble) , toDouble , fromDouble , toDoubleDouble , fromDoubleDouble+ , sqr ) where import Foreign (Ptr, alloca, allocaArray, castPtr, Storable(..), unsafePerformIO, with)@@ -42,34 +43,35 @@ , c_qd_ceil , c_qd_floor , c_qd_atan2+ , c_qd_sqr ) -data QuadDouble = QuadDouble !CDouble !CDouble !CDouble !CDouble deriving Typeable+data QuadDouble = QuadDouble {-# UNPACK #-} !CDouble {-# UNPACK #-} !CDouble {-# UNPACK #-} !CDouble {-# UNPACK #-} !CDouble deriving Typeable toDouble :: QuadDouble -> Double-toDouble (QuadDouble a _ _ _) = realToFrac a+toDouble !(QuadDouble a _ _ _) = realToFrac a fromDouble :: Double -> QuadDouble-fromDouble a = QuadDouble (realToFrac a) 0 0 0+fromDouble !a = QuadDouble (realToFrac a) 0 0 0 toDoubleDouble :: QuadDouble -> DoubleDouble-toDoubleDouble (QuadDouble a b _ _) = DoubleDouble a b+toDoubleDouble !(QuadDouble a b _ _) = DoubleDouble a b fromDoubleDouble :: DoubleDouble -> QuadDouble-fromDoubleDouble (DoubleDouble a b) = QuadDouble a b 0 0+fromDoubleDouble !(DoubleDouble a b) = QuadDouble a b 0 0 instance Eq QuadDouble where- a == b = a `compare` b == EQ- a /= b = a `compare` b /= EQ+ (!a) == (!b) = a `compare` b == EQ+ (!a) /= (!b) = a `compare` b /= EQ instance Ord QuadDouble where- a `compare` b = unsafePerformIO $ with a $ \p -> with b $ \q -> alloca $ \r -> do- c_qd_comp (castPtr p) (castPtr q) (castPtr r)- i <- peek r- return $ i `compare` (0 :: CInt)+ (!a) `compare` (!b) = unsafePerformIO $ with a $ \p -> with b $ \q -> alloca $ \r -> do+ c_qd_comp (castPtr p) (castPtr q) (castPtr r)+ !i <- peek r+ return $ i `compare` (0 :: CInt) instance Show QuadDouble where- show a = unsafePerformIO $ with a $ \p -> allocaArray (fromIntegral l) $ \r -> do+ 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@@ -80,13 +82,16 @@ (-) = lift_qd_qd_qd c_qd_sub 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+ 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 +sqr :: QuadDouble -> QuadDouble+sqr = lift_qd_qd c_qd_sqr+ 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+ recip !b = 1 / b+ fromRational !k = let { !n = numerator k ; !d = denominator k } in fromInteger n / fromInteger d instance Real QuadDouble where toRational (QuadDouble a b c d) = toRational a + toRational b + toRational c + toRational d@@ -124,11 +129,11 @@ 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- 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- isNegativeZero (QuadDouble a b c d) = isNegativeZero a || (a == 0 && (isNegativeZero b || (b == 0 && (isNegativeZero c || (c == 0 && isNegativeZero d)))))+ scaleFloat !n !x = x * 2 ** fromIntegral n+ 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+ isNegativeZero !(QuadDouble a b c d) = isNegativeZero a || (a == 0 && (isNegativeZero b || (b == 0 && (isNegativeZero c || (c == 0 && isNegativeZero d))))) isIEEE _ = False -- FIXME what does this imply? atan2 = lift_qd_qd_qd c_qd_atan2 @@ -138,22 +143,24 @@ instance Storable QuadDouble where sizeOf _ = 4 * sizeOf (error "Numeric.QD.QuadDouble.sizeOf" :: CDouble) alignment _ = alignment (error "Numeric.QD.QuadDouble.alignment" :: CDouble)- peek p = do- let q = castPtr p+ peek !p = do+ let !q = castPtr p a <- peekElemOff q 0 b <- peekElemOff q 1 c <- peekElemOff q 2 d <- peekElemOff q 3 return $ QuadDouble a b c d- poke p (QuadDouble a b c d) = do- let q = castPtr p+ poke !p !(QuadDouble a b c d) = do+ let !q = castPtr p pokeElemOff q 0 a pokeElemOff q 1 b pokeElemOff q 2 c pokeElemOff q 3 d lift_qd_qd :: (Ptr CDouble -> Ptr CDouble -> IO ()) -> QuadDouble -> QuadDouble-lift_qd_qd f a = unsafePerformIO $ with a $ \p -> alloca $ \r -> f (castPtr p) (castPtr r) >> peek r+{-# 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-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+{-# 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
− TODO
@@ -1,24 +0,0 @@-Present but undefined:-- properFraction- truncate- round- decodeFloat- encodeFloat- exponent- significand--Not yet implemented:-- instance Enum- instance Read--Need to be verified:-- isIEEE--Cannot be implemented:-- instance Data -- no instance for CDouble--More documentation!
qd.cabal view
@@ -1,9 +1,67 @@ Name: qd-Version: 0.1-Synopsis: double-double and quad-double type via libqd-Description: This package supports both a double-double datatype (approx. 32 decimal digits)+Version: 0.2+Synopsis: double-double and quad-double number type via libqd+Description:+ This package supports both a double-double datatype (approx. 32 decimal digits) and a quad-double datatype (approx. 64 decimal digits), using libqd (which is- implemented in C++ with C and Fortran wrappers).+ implemented in C++ with C and Fortran wrappers). To compile this package you+ 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@.+ .+ Additional note: libqd depends on 64bit doubles, while some FPU architectures+ use 80bit. It is highly recommended to compile with -fno-excess-precision and+ set the FPU control words to avoid erroneous behaviour, perhaps by doing+ something like this at the start of your program:+ .+ > import Foreign (nullPtr)+ > import GHC.Conc (forkOnIO, numCapabilities)+ > import Numeric.QD.FPU.Raw (fpu_fix_start)+ > main :: IO ()+ > main = do+ > mapM_ (flip forkOnIO $ fpu_fix_start nullPtr) [ 0 .. numCapabilities - 1 ]+ > -- ...+ License: BSD3 License-file: LICENSE Author: Claude Heiland-Allen@@ -11,7 +69,6 @@ Category: Math Build-type: Simple Cabal-version: >=1.2-Extra-source-files: TODO Library Build-depends: base >= 4 && < 5@@ -26,3 +83,5 @@ Numeric.QD.QuadDouble Numeric.QD.QuadDouble.Raw Numeric.QD.Raw+ GHC-options: -O2 -Wall -fno-excess-precision+ GHC-prof-options: -O2 -Wall -fno-excess-precision -prof -auto-all -caf-all