qd (empty) → 0.1
raw patch · 14 files changed
+750/−0 lines, 14 filesdep +basesetup-changed
Dependencies added: base
Files
- LICENSE +30/−0
- Numeric/QD.hs +11/−0
- Numeric/QD/Bits.hs +5/−0
- Numeric/QD/Bits/Raw.hs +8/−0
- Numeric/QD/DoubleDouble.hs +146/−0
- Numeric/QD/DoubleDouble/Raw.hs +127/−0
- Numeric/QD/FPU.hs +22/−0
- Numeric/QD/FPU/Raw.hs +11/−0
- Numeric/QD/QuadDouble.hs +159/−0
- Numeric/QD/QuadDouble/Raw.hs +166/−0
- Numeric/QD/Raw.hs +11/−0
- Setup.hs +2/−0
- TODO +24/−0
- qd.cabal +28/−0
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c)2010, Claude Heiland-Allen++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Claude Heiland-Allen nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Numeric/QD.hs view
@@ -0,0 +1,11 @@+module Numeric.QD+ ( module Numeric.QD.Bits+ , module Numeric.QD.DoubleDouble+ , module Numeric.QD.FPU+ , module Numeric.QD.QuadDouble+ ) where++import Numeric.QD.Bits+import Numeric.QD.DoubleDouble hiding (toDouble, fromDouble)+import Numeric.QD.FPU+import Numeric.QD.QuadDouble hiding (toDouble, fromDouble, toDoubleDouble, fromDoubleDouble)
+ Numeric/QD/Bits.hs view
@@ -0,0 +1,5 @@+module Numeric.QD.Bits+ (+ ) where++-- nothing here yet
+ Numeric/QD/Bits/Raw.hs view
@@ -0,0 +1,8 @@+{-# 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
@@ -0,0 +1,146 @@+{-# LANGUAGE DeriveDataTypeable #-}+module Numeric.QD.DoubleDouble+ ( DoubleDouble(DoubleDouble)+ , toDouble+ , fromDouble+ ) where++import Foreign (Ptr, alloca, allocaArray, castPtr, Storable(..), unsafePerformIO, with)+import Foreign.C (CDouble, CInt, peekCString)+import Data.Ratio (denominator, numerator)+import Data.Typeable (Typeable(..))++import Numeric.QD.DoubleDouble.Raw+ ( c_dd_add+ , c_dd_sub+ , c_dd_mul+ , c_dd_div+ , c_dd_pi+ , c_dd_exp+ , c_dd_sqrt+ , c_dd_log+ , c_dd_sin+ , c_dd_cos+ , c_dd_tan+ , c_dd_asin+ , c_dd_acos+ , c_dd_atan+ , c_dd_sinh+ , c_dd_cosh+ , c_dd_tanh+ , c_dd_asinh+ , c_dd_acosh+ , c_dd_atanh+ , c_dd_comp+ , c_dd_swrite+ , c_dd_neg+ , c_dd_abs+ , c_dd_copy_d+ , c_dd_ceil+ , c_dd_floor+ , c_dd_atan2+ )++data DoubleDouble = DoubleDouble !CDouble !CDouble deriving Typeable++toDouble :: DoubleDouble -> Double+toDouble (DoubleDouble a _) = realToFrac a++fromDouble :: Double -> DoubleDouble+fromDouble a = DoubleDouble (realToFrac a) 0++instance Eq DoubleDouble where+ 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)++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+ +instance Num DoubleDouble where+ (+) = lift_dd_dd_dd c_dd_add+ (*) = lift_dd_dd_dd c_dd_mul+ (-) = 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++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++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++instance Floating DoubleDouble where+ pi = unsafePerformIO $ alloca $ \r -> c_dd_pi (castPtr r) >> peek r+ exp = lift_dd_dd c_dd_exp+ sqrt = lift_dd_dd c_dd_sqrt+ log = lift_dd_dd c_dd_log+ sin = lift_dd_dd c_dd_sin+ cos = lift_dd_dd c_dd_cos+ tan = lift_dd_dd c_dd_tan+ asin = lift_dd_dd c_dd_asin+ acos = lift_dd_dd c_dd_acos+ atan = lift_dd_dd c_dd_atan+ sinh = lift_dd_dd c_dd_sinh+ cosh = lift_dd_dd c_dd_cosh+ tanh = lift_dd_dd c_dd_tanh+ asinh = lift_dd_dd c_dd_asinh+ acosh = lift_dd_dd c_dd_acosh+ atanh = lift_dd_dd c_dd_atanh++instance RealFloat DoubleDouble where+ 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+ 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++-- instance Enum DoubleDouble -- FIXME+-- instance Read DoubleDouble -- FIXME++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+ a <- peekElemOff q 0+ b <- peekElemOff q 1+ return $ DoubleDouble a b+ 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++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
+ Numeric/QD/DoubleDouble/Raw.hs view
@@ -0,0 +1,127 @@+{-# LANGUAGE ForeignFunctionInterface #-}+module Numeric.QD.DoubleDouble.Raw+ ( c_dd_add+ , c_dd_add_d_dd+ , c_dd_add_dd_d+ , c_dd_sub+ , c_dd_sub_d_dd+ , c_dd_sub_dd_d+ , c_dd_mul+ , c_dd_mul_d_dd+ , c_dd_mul_dd_d+ , c_dd_div+ , c_dd_div_d_dd+ , c_dd_div_dd_d+ , c_dd_copy+ , c_dd_copy_d+ , c_dd_sqrt+ , c_dd_sqr+ , c_dd_abs+ , c_dd_npwr+ , c_dd_nroot+ , c_dd_nint+ , c_dd_aint+ , c_dd_floor+ , c_dd_ceil+ , c_dd_exp+ , c_dd_log+ , c_dd_log10+ , c_dd_sin+ , c_dd_cos+ , c_dd_tan+ , c_dd_asin+ , c_dd_acos+ , c_dd_atan+ , c_dd_atan2+ , c_dd_sinh+ , c_dd_cosh+ , c_dd_tanh+ , c_dd_asinh+ , c_dd_acosh+ , c_dd_atanh+ , c_dd_sincos+ , c_dd_sincosh+ , c_dd_read+ , c_dd_swrite+ , c_dd_neg+ , c_dd_comp+ , c_dd_comp_dd_d+ , c_dd_comp_d_dd+ , c_dd_pi+ , c_dd_write+ , c_dd_rand+ ) where++import Foreign (Ptr)+import Foreign.C (CChar, CDouble, CInt)++foreign import ccall unsafe "qd/c_dd.h c_dd_add" c_dd_add :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_add_d_dd" c_dd_add_d_dd :: CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_add_dd_d" c_dd_add_dd_d :: Ptr CDouble -> CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_dd.h c_dd_sub" c_dd_sub :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_sub_d_dd" c_dd_sub_d_dd :: CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_sub_dd_d" c_dd_sub_dd_d :: Ptr CDouble -> CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_dd.h c_dd_mul" c_dd_mul :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_mul_d_dd" c_dd_mul_d_dd :: CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_mul_dd_d" c_dd_mul_dd_d :: Ptr CDouble -> CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_dd.h c_dd_div" c_dd_div :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_div_d_dd" c_dd_div_d_dd :: CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_div_dd_d" c_dd_div_dd_d :: Ptr CDouble -> CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_dd.h c_dd_copy" c_dd_copy :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_copy_d" c_dd_copy_d :: CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_dd.h c_dd_sqrt" c_dd_sqrt :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_sqr" c_dd_sqr :: Ptr CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_dd.h c_dd_abs" c_dd_abs :: Ptr CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_dd.h c_dd_npwr" c_dd_npwr :: Ptr CDouble -> CInt -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_nroot" c_dd_nroot :: Ptr CDouble -> CInt -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_dd.h c_dd_nint" c_dd_nint :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_aint" c_dd_aint :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_floor" c_dd_floor :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_ceil" c_dd_ceil :: Ptr CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_dd.h c_dd_exp" c_dd_exp :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_log" c_dd_log :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_log10" c_dd_log10 :: Ptr CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_dd.h c_dd_sin" c_dd_sin :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_cos" c_dd_cos :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_tan" c_dd_tan :: Ptr CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_dd.h c_dd_asin" c_dd_asin :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_acos" c_dd_acos :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_atan" c_dd_atan :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_atan2" c_dd_atan2 :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_dd.h c_dd_sinh" c_dd_sinh :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_cosh" c_dd_cosh :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_tanh" c_dd_tanh :: Ptr CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_dd.h c_dd_asinh" c_dd_asinh :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_acosh" c_dd_acosh :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_atanh" c_dd_atanh :: Ptr CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_dd.h c_dd_sincos" c_dd_sincos :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_sincosh" c_dd_sincosh :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_dd.h c_dd_read" c_dd_read :: Ptr CChar -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_swrite" c_dd_swrite :: Ptr CDouble -> CInt -> Ptr CChar -> CInt -> IO ()++foreign import ccall unsafe "qd/c_dd.h c_dd_neg" c_dd_neg :: Ptr CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_dd.h c_dd_comp" c_dd_comp :: Ptr CDouble -> Ptr CDouble -> Ptr CInt -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_comp_dd_d" c_dd_comp_dd_d :: Ptr CDouble -> CDouble -> Ptr CInt -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_comp_d_dd" c_dd_comp_d_dd :: CDouble -> Ptr CDouble -> Ptr CInt -> IO ()++foreign import ccall unsafe "qd/c_dd.h c_dd_pi" c_dd_pi :: Ptr CDouble -> IO ()++-- these two are almost certainly not referentially transparent+foreign import ccall unsafe "qd/c_dd.h c_dd_write" c_dd_write :: Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_dd.h c_dd_rand" c_dd_rand :: Ptr CDouble -> IO ()
+ Numeric/QD/FPU.hs view
@@ -0,0 +1,22 @@+module Numeric.QD.FPU+ ( unsafePreservingFPU+ ) where++import Foreign.Marshal.Alloc (alloca)+import Control.Exception (bracket_)++import Numeric.QD.FPU.Raw (fpu_fix_start, fpu_fix_end)++-- |@'unsafePreservingFPU' f@ executes the computation @f@, ensuring+-- that the FPU control words are set to avoid problems from excess+-- precision. See the libqd documentation for further details.+-- +-- This function is unsafe in a threaded runtime as Haskell threads can+-- migrate between OS threads, moreover there is no checking for nested+-- calls which results in race conditions.+--+-- Some steps can be taken to mitigate some of this badness; perhaps+-- using (for example) @'GHC.Conc.forkOnIO'@ might help.++unsafePreservingFPU :: IO a -> IO a+unsafePreservingFPU f = alloca $ \p -> bracket_ (fpu_fix_start p) (fpu_fix_end p) f
+ Numeric/QD/FPU/Raw.hs view
@@ -0,0 +1,11 @@+{-# LANGUAGE ForeignFunctionInterface #-}+module Numeric.QD.FPU.Raw+ ( fpu_fix_start+ , fpu_fix_end+ ) where++import Foreign (Ptr)+import Foreign.C (CInt)++foreign import ccall unsafe "qd/fpu.h fpu_fix_start" fpu_fix_start :: Ptr CInt -> IO ()+foreign import ccall unsafe "qd/fpu.h fpu_fix_end" fpu_fix_end :: Ptr CInt -> IO ()
+ Numeric/QD/QuadDouble.hs view
@@ -0,0 +1,159 @@+{-# LANGUAGE DeriveDataTypeable #-}+module Numeric.QD.QuadDouble+ ( QuadDouble(QuadDouble)+ , toDouble+ , fromDouble+ , toDoubleDouble+ , fromDoubleDouble+ ) where++import Foreign (Ptr, alloca, allocaArray, castPtr, Storable(..), unsafePerformIO, with)+import Foreign.C (CDouble, CInt, peekCString)+import Data.Ratio (denominator, numerator)+import Data.Typeable (Typeable(..))++import Numeric.QD.DoubleDouble (DoubleDouble(DoubleDouble))+import Numeric.QD.QuadDouble.Raw+ ( c_qd_add+ , c_qd_sub+ , c_qd_mul+ , c_qd_div+ , c_qd_pi+ , c_qd_exp+ , c_qd_sqrt+ , c_qd_log+ , c_qd_sin+ , c_qd_cos+ , c_qd_tan+ , c_qd_asin+ , c_qd_acos+ , c_qd_atan+ , c_qd_sinh+ , c_qd_cosh+ , c_qd_tanh+ , c_qd_asinh+ , c_qd_acosh+ , c_qd_atanh+ , c_qd_comp+ , c_qd_swrite+ , c_qd_neg+ , c_qd_abs+ , c_qd_copy_d+ , c_qd_ceil+ , c_qd_floor+ , c_qd_atan2+ )++data QuadDouble = QuadDouble !CDouble !CDouble !CDouble !CDouble deriving Typeable++toDouble :: QuadDouble -> Double+toDouble (QuadDouble a _ _ _) = realToFrac a++fromDouble :: Double -> QuadDouble+fromDouble a = QuadDouble (realToFrac a) 0 0 0++toDoubleDouble :: QuadDouble -> DoubleDouble+toDoubleDouble (QuadDouble a b _ _) = DoubleDouble a b++fromDoubleDouble :: DoubleDouble -> QuadDouble+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++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)++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+ (-) = 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++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++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++instance Floating QuadDouble where+ pi = unsafePerformIO $ alloca $ \r -> c_qd_pi (castPtr r) >> peek r+ exp = lift_qd_qd c_qd_exp+ sqrt = lift_qd_qd c_qd_sqrt+ log = lift_qd_qd c_qd_log+ sin = lift_qd_qd c_qd_sin+ cos = lift_qd_qd c_qd_cos+ tan = lift_qd_qd c_qd_tan+ asin = lift_qd_qd c_qd_asin+ acos = lift_qd_qd c_qd_acos+ atan = lift_qd_qd c_qd_atan+ sinh = lift_qd_qd c_qd_sinh+ cosh = lift_qd_qd c_qd_cosh+ tanh = lift_qd_qd c_qd_tanh+ asinh = lift_qd_qd c_qd_asinh+ acosh = lift_qd_qd c_qd_acosh+ atanh = lift_qd_qd c_qd_atanh++instance RealFloat QuadDouble where+ 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+ 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++-- instance Enum QuadDouble -- FIXME+-- instance Read QuadDouble -- FIXME++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+ 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+ 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++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
+ Numeric/QD/QuadDouble/Raw.hs view
@@ -0,0 +1,166 @@+{-# LANGUAGE ForeignFunctionInterface #-}+module Numeric.QD.QuadDouble.Raw+ ( c_qd_add+ , c_qd_add_dd_qd+ , c_qd_add_qd_dd+ , c_qd_add_d_qd+ , c_qd_add_qd_d+ , c_qd_selfadd+ , c_qd_selfadd_dd+ , c_qd_selfadd_d+ , c_qd_sub+ , c_qd_sub_dd_qd+ , c_qd_sub_qd_dd+ , c_qd_sub_d_qd+ , c_qd_sub_qd_d+ , c_qd_selfsub+ , c_qd_selfsub_dd+ , c_qd_selfsub_d+ , c_qd_mul+ , c_qd_mul_dd_qd+ , c_qd_mul_qd_dd+ , c_qd_mul_d_qd+ , c_qd_mul_qd_d+ , c_qd_selfmul+ , c_qd_selfmul_dd+ , c_qd_selfmul_d+ , c_qd_div+ , c_qd_div_dd_qd+ , c_qd_div_qd_dd+ , c_qd_div_d_qd+ , c_qd_div_qd_d+ , c_qd_selfdiv+ , c_qd_selfdiv_dd+ , c_qd_selfdiv_d+ , c_qd_copy+ , c_qd_copy_dd+ , c_qd_copy_d+ , c_qd_sqrt+ , c_qd_sqr+ , c_qd_abs+ , c_qd_npwr+ , c_qd_nroot+ , c_qd_nint+ , c_qd_aint+ , c_qd_floor+ , c_qd_ceil+ , c_qd_exp+ , c_qd_log+ , c_qd_log10+ , c_qd_sin+ , c_qd_cos+ , c_qd_tan+ , c_qd_asin+ , c_qd_acos+ , c_qd_atan+ , c_qd_atan2+ , c_qd_sinh+ , c_qd_cosh+ , c_qd_tanh+ , c_qd_asinh+ , c_qd_acosh+ , c_qd_atanh+ , c_qd_sincos+ , c_qd_sincosh+ , c_qd_read+ , c_qd_swrite+ , c_qd_neg+ , c_qd_comp+ , c_qd_comp_qd_d+ , c_qd_comp_d_qd+ , c_qd_pi+ , c_qd_write+ , c_qd_rand+ ) where++import Foreign (Ptr)+import Foreign.C (CChar, CDouble, CInt)++foreign import ccall unsafe "qd/c_qd.h c_qd_add" c_qd_add :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ();+foreign import ccall unsafe "qd/c_qd.h c_qd_add_dd_qd" c_qd_add_dd_qd :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ();+foreign import ccall unsafe "qd/c_qd.h c_qd_add_qd_dd" c_qd_add_qd_dd :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ();+foreign import ccall unsafe "qd/c_qd.h c_qd_add_d_qd" c_qd_add_d_qd :: CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_add_qd_d" c_qd_add_qd_d :: Ptr CDouble -> CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_selfadd" c_qd_selfadd :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_selfadd_dd" c_qd_selfadd_dd :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_selfadd_d" c_qd_selfadd_d :: CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_qd.h c_qd_sub" c_qd_sub :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ();+foreign import ccall unsafe "qd/c_qd.h c_qd_sub_dd_qd" c_qd_sub_dd_qd :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ();+foreign import ccall unsafe "qd/c_qd.h c_qd_sub_qd_dd" c_qd_sub_qd_dd :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ();+foreign import ccall unsafe "qd/c_qd.h c_qd_sub_d_qd" c_qd_sub_d_qd :: CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_sub_qd_d" c_qd_sub_qd_d :: Ptr CDouble -> CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_selfsub" c_qd_selfsub :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_selfsub_dd" c_qd_selfsub_dd :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_selfsub_d" c_qd_selfsub_d :: CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_qd.h c_qd_mul" c_qd_mul :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ();+foreign import ccall unsafe "qd/c_qd.h c_qd_mul_dd_qd" c_qd_mul_dd_qd :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ();+foreign import ccall unsafe "qd/c_qd.h c_qd_mul_qd_dd" c_qd_mul_qd_dd :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ();+foreign import ccall unsafe "qd/c_qd.h c_qd_mul_d_qd" c_qd_mul_d_qd :: CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_mul_qd_d" c_qd_mul_qd_d :: Ptr CDouble -> CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_selfmul" c_qd_selfmul :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_selfmul_dd" c_qd_selfmul_dd :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_selfmul_d" c_qd_selfmul_d :: CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_qd.h c_qd_div" c_qd_div :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ();+foreign import ccall unsafe "qd/c_qd.h c_qd_div_dd_qd" c_qd_div_dd_qd :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ();+foreign import ccall unsafe "qd/c_qd.h c_qd_div_qd_dd" c_qd_div_qd_dd :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ();+foreign import ccall unsafe "qd/c_qd.h c_qd_div_d_qd" c_qd_div_d_qd :: CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_div_qd_d" c_qd_div_qd_d :: Ptr CDouble -> CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_selfdiv" c_qd_selfdiv :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_selfdiv_dd" c_qd_selfdiv_dd :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_selfdiv_d" c_qd_selfdiv_d :: CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_qd.h c_qd_copy" c_qd_copy :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_copy_dd" c_qd_copy_dd :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_copy_d" c_qd_copy_d :: CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_qd.h c_qd_sqrt" c_qd_sqrt :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_sqr" c_qd_sqr :: Ptr CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_qd.h c_qd_abs" c_qd_abs :: Ptr CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_qd.h c_qd_npwr" c_qd_npwr :: Ptr CDouble -> CInt -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_nroot" c_qd_nroot :: Ptr CDouble -> CInt -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_qd.h c_qd_nint" c_qd_nint :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_aint" c_qd_aint :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_floor" c_qd_floor :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_ceil" c_qd_ceil :: Ptr CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_qd.h c_qd_exp" c_qd_exp :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_log" c_qd_log :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_log10" c_qd_log10 :: Ptr CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_qd.h c_qd_sin" c_qd_sin :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_cos" c_qd_cos :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_tan" c_qd_tan :: Ptr CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_qd.h c_qd_asin" c_qd_asin :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_acos" c_qd_acos :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_atan" c_qd_atan :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_atan2" c_qd_atan2 :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ();++foreign import ccall unsafe "qd/c_qd.h c_qd_sinh" c_qd_sinh :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_cosh" c_qd_cosh :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_tanh" c_qd_tanh :: Ptr CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_qd.h c_qd_asinh" c_qd_asinh :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_acosh" c_qd_acosh :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_atanh" c_qd_atanh :: Ptr CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_qd.h c_qd_sincos" c_qd_sincos :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_sincosh" c_qd_sincosh :: Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()++foreign import ccall unsafe "qd/c_qd.h c_qd_read" c_qd_read :: Ptr CChar -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_swrite" c_qd_swrite :: Ptr CDouble -> CInt -> Ptr CChar -> CInt -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_neg" c_qd_neg :: Ptr CDouble -> Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_comp" c_qd_comp :: Ptr CDouble -> Ptr CDouble -> Ptr CInt -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_comp_qd_d" c_qd_comp_qd_d :: Ptr CDouble -> CDouble -> Ptr CInt -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_comp_d_qd" c_qd_comp_d_qd :: CDouble -> Ptr CDouble -> Ptr CInt -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_pi" c_qd_pi :: Ptr CDouble -> IO ()++-- these two are almost certainly not referentially transparent+foreign import ccall unsafe "qd/c_qd.h c_qd_write" c_qd_write :: Ptr CDouble -> IO ()+foreign import ccall unsafe "qd/c_qd.h c_qd_rand" c_qd_rand :: Ptr CDouble -> IO ()
+ Numeric/QD/Raw.hs view
@@ -0,0 +1,11 @@+module Numeric.QD.Raw+ ( module Numeric.QD.Bits.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
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ TODO view
@@ -0,0 +1,24 @@+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
@@ -0,0 +1,28 @@+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)+ and a quad-double datatype (approx. 64 decimal digits), using libqd (which is+ implemented in C++ with C and Fortran wrappers).+License: BSD3+License-file: LICENSE+Author: Claude Heiland-Allen+Maintainer: claudiusmaximus@goto10.org+Category: Math+Build-type: Simple+Cabal-version: >=1.2+Extra-source-files: TODO++Library+ Build-depends: base >= 4 && < 5+ Extra-Libraries: qd+ Exposed-modules: Numeric.QD+ Numeric.QD.Bits+ Numeric.QD.Bits.Raw+ Numeric.QD.DoubleDouble+ Numeric.QD.DoubleDouble.Raw+ Numeric.QD.FPU+ Numeric.QD.FPU.Raw+ Numeric.QD.QuadDouble+ Numeric.QD.QuadDouble.Raw+ Numeric.QD.Raw