Remove Hugs98 specific code
[packages/base.git] / Data / Complex.hs
1 {-# LANGUAGE Trustworthy #-}
2 {-# LANGUAGE CPP, DeriveDataTypeable #-}
3 #ifdef __GLASGOW_HASKELL__
4 {-# LANGUAGE StandaloneDeriving #-}
5 #endif
6
7 -----------------------------------------------------------------------------
8 -- |
9 -- Module : Data.Complex
10 -- Copyright : (c) The University of Glasgow 2001
11 -- License : BSD-style (see the file libraries/base/LICENSE)
12 --
13 -- Maintainer : libraries@haskell.org
14 -- Stability : provisional
15 -- Portability : portable
16 --
17 -- Complex numbers.
18 --
19 -----------------------------------------------------------------------------
20
21 module Data.Complex
22 (
23 -- * Rectangular form
24 Complex((:+))
25
26 , realPart
27 , imagPart
28 -- * Polar form
29 , mkPolar
30 , cis
31 , polar
32 , magnitude
33 , phase
34 -- * Conjugate
35 , conjugate
36
37 ) where
38
39 import Prelude
40
41 import Data.Typeable
42 #ifdef __GLASGOW_HASKELL__
43 import Data.Data (Data)
44 #endif
45
46 infix 6 :+
47
48 -- -----------------------------------------------------------------------------
49 -- The Complex type
50
51 -- | Complex numbers are an algebraic type.
52 --
53 -- For a complex number @z@, @'abs' z@ is a number with the magnitude of @z@,
54 -- but oriented in the positive real direction, whereas @'signum' z@
55 -- has the phase of @z@, but unit magnitude.
56 data Complex a
57 = !a :+ !a -- ^ forms a complex number from its real and imaginary
58 -- rectangular components.
59 # if __GLASGOW_HASKELL__
60 deriving (Eq, Show, Read, Data)
61 # else
62 deriving (Eq, Show, Read)
63 # endif
64
65 -- -----------------------------------------------------------------------------
66 -- Functions over Complex
67
68 -- | Extracts the real part of a complex number.
69 realPart :: (RealFloat a) => Complex a -> a
70 realPart (x :+ _) = x
71
72 -- | Extracts the imaginary part of a complex number.
73 imagPart :: (RealFloat a) => Complex a -> a
74 imagPart (_ :+ y) = y
75
76 -- | The conjugate of a complex number.
77 {-# SPECIALISE conjugate :: Complex Double -> Complex Double #-}
78 conjugate :: (RealFloat a) => Complex a -> Complex a
79 conjugate (x:+y) = x :+ (-y)
80
81 -- | Form a complex number from polar components of magnitude and phase.
82 {-# SPECIALISE mkPolar :: Double -> Double -> Complex Double #-}
83 mkPolar :: (RealFloat a) => a -> a -> Complex a
84 mkPolar r theta = r * cos theta :+ r * sin theta
85
86 -- | @'cis' t@ is a complex value with magnitude @1@
87 -- and phase @t@ (modulo @2*'pi'@).
88 {-# SPECIALISE cis :: Double -> Complex Double #-}
89 cis :: (RealFloat a) => a -> Complex a
90 cis theta = cos theta :+ sin theta
91
92 -- | The function 'polar' takes a complex number and
93 -- returns a (magnitude, phase) pair in canonical form:
94 -- the magnitude is nonnegative, and the phase in the range @(-'pi', 'pi']@;
95 -- if the magnitude is zero, then so is the phase.
96 {-# SPECIALISE polar :: Complex Double -> (Double,Double) #-}
97 polar :: (RealFloat a) => Complex a -> (a,a)
98 polar z = (magnitude z, phase z)
99
100 -- | The nonnegative magnitude of a complex number.
101 {-# SPECIALISE magnitude :: Complex Double -> Double #-}
102 magnitude :: (RealFloat a) => Complex a -> a
103 magnitude (x:+y) = scaleFloat k
104 (sqrt (sqr (scaleFloat mk x) + sqr (scaleFloat mk y)))
105 where k = max (exponent x) (exponent y)
106 mk = - k
107 sqr z = z * z
108
109 -- | The phase of a complex number, in the range @(-'pi', 'pi']@.
110 -- If the magnitude is zero, then so is the phase.
111 {-# SPECIALISE phase :: Complex Double -> Double #-}
112 phase :: (RealFloat a) => Complex a -> a
113 phase (0 :+ 0) = 0 -- SLPJ July 97 from John Peterson
114 phase (x:+y) = atan2 y x
115
116
117 -- -----------------------------------------------------------------------------
118 -- Instances of Complex
119
120 #include "Typeable.h"
121 INSTANCE_TYPEABLE1(Complex,complexTc,"Complex")
122
123 instance (RealFloat a) => Num (Complex a) where
124 {-# SPECIALISE instance Num (Complex Float) #-}
125 {-# SPECIALISE instance Num (Complex Double) #-}
126 (x:+y) + (x':+y') = (x+x') :+ (y+y')
127 (x:+y) - (x':+y') = (x-x') :+ (y-y')
128 (x:+y) * (x':+y') = (x*x'-y*y') :+ (x*y'+y*x')
129 negate (x:+y) = negate x :+ negate y
130 abs z = magnitude z :+ 0
131 signum (0:+0) = 0
132 signum z@(x:+y) = x/r :+ y/r where r = magnitude z
133 fromInteger n = fromInteger n :+ 0
134
135 instance (RealFloat a) => Fractional (Complex a) where
136 {-# SPECIALISE instance Fractional (Complex Float) #-}
137 {-# SPECIALISE instance Fractional (Complex Double) #-}
138 (x:+y) / (x':+y') = (x*x''+y*y'') / d :+ (y*x''-x*y'') / d
139 where x'' = scaleFloat k x'
140 y'' = scaleFloat k y'
141 k = - max (exponent x') (exponent y')
142 d = x'*x'' + y'*y''
143
144 fromRational a = fromRational a :+ 0
145
146 instance (RealFloat a) => Floating (Complex a) where
147 {-# SPECIALISE instance Floating (Complex Float) #-}
148 {-# SPECIALISE instance Floating (Complex Double) #-}
149 pi = pi :+ 0
150 exp (x:+y) = expx * cos y :+ expx * sin y
151 where expx = exp x
152 log z = log (magnitude z) :+ phase z
153
154 sqrt (0:+0) = 0
155 sqrt z@(x:+y) = u :+ (if y < 0 then -v else v)
156 where (u,v) = if x < 0 then (v',u') else (u',v')
157 v' = abs y / (u'*2)
158 u' = sqrt ((magnitude z + abs x) / 2)
159
160 sin (x:+y) = sin x * cosh y :+ cos x * sinh y
161 cos (x:+y) = cos x * cosh y :+ (- sin x * sinh y)
162 tan (x:+y) = (sinx*coshy:+cosx*sinhy)/(cosx*coshy:+(-sinx*sinhy))
163 where sinx = sin x
164 cosx = cos x
165 sinhy = sinh y
166 coshy = cosh y
167
168 sinh (x:+y) = cos y * sinh x :+ sin y * cosh x
169 cosh (x:+y) = cos y * cosh x :+ sin y * sinh x
170 tanh (x:+y) = (cosy*sinhx:+siny*coshx)/(cosy*coshx:+siny*sinhx)
171 where siny = sin y
172 cosy = cos y
173 sinhx = sinh x
174 coshx = cosh x
175
176 asin z@(x:+y) = y':+(-x')
177 where (x':+y') = log (((-y):+x) + sqrt (1 - z*z))
178 acos z = y'':+(-x'')
179 where (x'':+y'') = log (z + ((-y'):+x'))
180 (x':+y') = sqrt (1 - z*z)
181 atan z@(x:+y) = y':+(-x')
182 where (x':+y') = log (((1-y):+x) / sqrt (1+z*z))
183
184 asinh z = log (z + sqrt (1+z*z))
185 acosh z = log (z + (z+1) * sqrt ((z-1)/(z+1)))
186 atanh z = 0.5 * log ((1.0+z) / (1.0-z))
187