bugfix: make fromInteger and fromRational work with full precision
[hsqd:hsqd.git] / Numeric / QD / QuadDouble.hs
1 {-# LANGUAGE BangPatterns, DeriveDataTypeable #-}
2 module Numeric.QD.QuadDouble
3   ( QuadDouble(QuadDouble)
4   , toDouble
5   , fromDouble
6   , toDoubleDouble
7   , fromDoubleDouble
8   , sqr
9   ) where
10
11 import Foreign (Ptr, alloca, allocaArray, castPtr, Storable(..), unsafePerformIO, with)
12 import Foreign.C (CDouble, CInt, peekCString)
13 import Data.Ratio ((%))
14 import Data.Typeable (Typeable(..))
15
16 import Numeric.QD.DoubleDouble (DoubleDouble(DoubleDouble))
17 import Numeric.QD.QuadDouble.Raw
18   ( c_qd_add
19   , c_qd_sub
20   , c_qd_mul
21   , c_qd_div
22   , c_qd_pi
23   , c_qd_exp
24   , c_qd_sqrt
25   , c_qd_log
26   , c_qd_sin
27   , c_qd_cos
28   , c_qd_tan
29   , c_qd_asin
30   , c_qd_acos
31   , c_qd_atan
32   , c_qd_sinh
33   , c_qd_cosh
34   , c_qd_tanh
35   , c_qd_asinh
36   , c_qd_acosh
37   , c_qd_atanh
38   , c_qd_comp
39   , c_qd_swrite
40   , c_qd_neg
41   , c_qd_abs
42   , c_qd_ceil
43   , c_qd_floor
44   , c_qd_atan2
45   , c_qd_sqr
46   )
47
48 data QuadDouble = QuadDouble {-# UNPACK #-} !CDouble {-# UNPACK #-} !CDouble {-# UNPACK #-} !CDouble {-# UNPACK #-} !CDouble deriving Typeable
49
50 toDouble :: QuadDouble -> Double
51 toDouble !(QuadDouble a _ _ _) = realToFrac a
52
53 fromDouble :: Double -> QuadDouble
54 fromDouble !a = QuadDouble (realToFrac a) 0 0 0
55
56 toDoubleDouble :: QuadDouble -> DoubleDouble
57 toDoubleDouble !(QuadDouble a b _ _) = DoubleDouble a b
58
59 fromDoubleDouble :: DoubleDouble -> QuadDouble
60 fromDoubleDouble !(DoubleDouble a b) = QuadDouble a b 0 0
61
62 instance Eq QuadDouble where
63   (!a) == (!b) = a `compare` b == EQ
64   (!a) /= (!b) = a `compare` b /= EQ
65
66 instance Ord QuadDouble where
67   (!a) `compare` (!b) = unsafePerformIO $ with a $ \p -> with b $ \q -> alloca $ \r -> do
68                           c_qd_comp (castPtr p) (castPtr q) (castPtr r)
69                           !i <- peek r
70                           return $ i `compare` (0 :: CInt)
71
72 instance Show QuadDouble where
73   show !a =  unsafePerformIO $ with a $ \p -> allocaArray (fromIntegral l) $ \r -> do
74               c_qd_swrite (castPtr p) (l`div`2) (castPtr r) l
75               peekCString r
76             where l = 128
77   
78 instance Num QuadDouble where
79   (+) = lift_qd_qd_qd c_qd_add
80   (*) = lift_qd_qd_qd c_qd_mul
81   (-) = lift_qd_qd_qd c_qd_sub
82   negate = lift_qd_qd c_qd_neg
83   abs = lift_qd_qd c_qd_abs
84   signum !a = case a `compare` 0 of { LT -> -1 ; EQ -> 0 ; GT -> 1 }
85   fromInteger !i = fromRational (i % 1)
86
87 sqr :: QuadDouble -> QuadDouble
88 sqr = lift_qd_qd c_qd_sqr
89
90 instance Fractional QuadDouble where
91   (/) = lift_qd_qd_qd c_qd_div
92   recip !b = 1 / b
93   fromRational !k = let a = fromRational k
94                         ka = k  - toRational a
95                         b = fromRational ka
96                         kb = ka - toRational b
97                         c = fromRational kb
98                         kc = kb - toRational c
99                         d = fromRational kc
100                     in  QuadDouble a b c d
101
102 instance Real QuadDouble where
103   toRational (QuadDouble a b c d) = toRational a + toRational b + toRational c + toRational d
104
105 instance RealFrac QuadDouble where
106   properFraction = error "Numeric.QD.QuadDouble.properFraction: not yet implemented" -- FIXME
107   truncate = error "Numeric.QD.QuadDouble.truncate: not yet implemented" -- FIXME
108   round = error "Numeric.QD.QuadDouble.round: not yet implemented" -- FIXME
109   ceiling = ceiling . toDouble . lift_qd_qd c_qd_ceil
110   floor = floor . toDouble . lift_qd_qd c_qd_floor
111
112 instance Floating QuadDouble where
113   pi = unsafePerformIO $ alloca $ \r -> c_qd_pi (castPtr r) >> peek r
114   exp = lift_qd_qd c_qd_exp
115   sqrt = lift_qd_qd c_qd_sqrt
116   log = lift_qd_qd c_qd_log
117   sin = lift_qd_qd c_qd_sin
118   cos = lift_qd_qd c_qd_cos
119   tan = lift_qd_qd c_qd_tan
120   asin = lift_qd_qd c_qd_asin
121   acos = lift_qd_qd c_qd_acos
122   atan = lift_qd_qd c_qd_atan
123   sinh = lift_qd_qd c_qd_sinh
124   cosh = lift_qd_qd c_qd_cosh
125   tanh = lift_qd_qd c_qd_tanh
126   asinh = lift_qd_qd c_qd_asinh
127   acosh = lift_qd_qd c_qd_acosh
128   atanh = lift_qd_qd c_qd_atanh
129
130 instance RealFloat QuadDouble where
131   floatRadix _ = 2
132   floatDigits _ = 4 * floatDigits (error "Numeric.QD.QuadDouble.floatDigits" :: CDouble)
133   floatRange _ = floatRange (error "Numeric.QD.QuadDouble.floatRange" :: CDouble)
134   decodeFloat = error "Numeric.QD.QuadDouble.decodeFloat: not yet implemented" -- FIXME
135   encodeFloat = error "Numeric.QD.QuadDouble.encodeFloat: not yet implemented" -- FIXME
136   exponent _ = error "Numeric.QD.QuadDouble.exponent: not yet implemented" -- FIXME
137   significand _ = error "Numeric.QD.QuadDouble.significand: not yet implemented" -- FIXME
138   scaleFloat !n !x = x * 2 ** fromIntegral n
139   isNaN !(QuadDouble a b c d) = isNaN a || isNaN b || isNaN c || isNaN d
140   isInfinite !(QuadDouble a b c d) = isInfinite a || isInfinite b || isInfinite c || isInfinite d
141   isDenormalized !(QuadDouble a b c d) = isDenormalized a || isDenormalized b || isDenormalized c || isDenormalized d
142   isNegativeZero !(QuadDouble a b c d) = isNegativeZero a || (a == 0 && (isNegativeZero b || (b == 0 && (isNegativeZero c || (c == 0 && isNegativeZero d)))))
143   isIEEE _ = False -- FIXME what does this imply?
144   atan2 = lift_qd_qd_qd c_qd_atan2
145
146 -- instance Enum QuadDouble -- FIXME
147 -- instance Read QuadDouble -- FIXME
148
149 instance Storable QuadDouble where
150   sizeOf _ = 4 * sizeOf (error "Numeric.QD.QuadDouble.sizeOf" :: CDouble)
151   alignment _ = alignment (error "Numeric.QD.QuadDouble.alignment" :: CDouble)
152   peek !p = do
153     let !q = castPtr p
154     a <- peekElemOff q 0
155     b <- peekElemOff q 1
156     c <- peekElemOff q 2
157     d <- peekElemOff q 3
158     return $ QuadDouble a b c d
159   poke !p !(QuadDouble a b c d) = do
160     let !q = castPtr p
161     pokeElemOff q 0 a
162     pokeElemOff q 1 b
163     pokeElemOff q 2 c
164     pokeElemOff q 3 d
165
166 lift_qd_qd :: (Ptr CDouble -> Ptr CDouble -> IO ()) -> QuadDouble -> QuadDouble
167 lift_qd_qd f !a = unsafePerformIO $ with a $ \p -> alloca $ \r -> f (castPtr p) (castPtr r) >> peek r
168
169 lift_qd_qd_qd :: (Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()) -> QuadDouble -> QuadDouble -> QuadDouble
170 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