ghc/lib/prelude/IComplex.hs

   1 -- Complex Numbers
   2
   3 module PreludeComplex where
   4
   5 import Cls
   6 import Core
   7
   8 import IDouble  -- instances
   9 import IChar
  10 import IFloat
  11 import IInt
  12 import IInteger
  13 import IList
  14 import List             ( (++), foldr )
  15 import Prel             ( (.), (&&), (||), (^), atan2, otherwise )
  16 import PS               ( _PackedString, _unpackPS )
  17 import Text
  18 import TyArray
  19 import TyComplex
  20
  21 -- infix  6  :+
  22
  23 -- data  (RealFloat a)     => Complex a = a :+ a  deriving (Eq,Binary,Text)
  24
  25 instance (Eq a) => Eq (Complex a) where
  26     (x :+ y) == (x2 :+ y2) = x == x2 && y == y2
  27     (x :+ y) /= (x2 :+ y2) = x /= x2 || y /= y2
  28
  29 instance (RealFloat a) => Num (Complex a) where
  30     (x:+y) + (x2:+y2)   =  (x+x2) :+ (y+y2)
  31     (x:+y) - (x2:+y2)   =  (x-x2) :+ (y-y2)
  32     (x:+y) * (x2:+y2)   =  (x*x2-y*y2) :+ (x*y2+y*x2)
  33     negate (x:+y)       =  negate x :+ negate y
  34     abs z               =  magnitude z :+ __i0
  35     signum 0            =  __i0
  36     signum z@(x:+y)     =  x/r :+ y/r  where { r = magnitude z }
  37     fromInteger n       =  fromInteger n :+ __i0
  38     fromInt n           =  fromInt n :+ __i0
  39
  40 instance (RealFloat a) => Fractional (Complex a) where
  41     (x:+y) / (x2:+y2)   =  (x*x3+y*y3) / d :+ (y*x3-x*y3) / d
  42                           where  x3 = scaleFloat k x2
  43                                  y3 = scaleFloat k y2
  44                                  k  = - max (exponent x2) (exponent y2)
  45                                  d  = x2*x3 + y2*y3
  46     fromRational a      =  fromRational a :+ __i0
  47     recip a             =  __i1 / a
  48
  49 instance (RealFloat a) => Floating (Complex a) where
  50     pi             =  pi :+ __i0
  51     exp (x:+y)     =  expx * cos y :+ expx * sin y
  52                       where expx = exp x
  53     log z          =  log (magnitude z) :+ phase z
  54
  55     (**) a b       =  exp (log a * b)
  56     logBase a b    =  log b / log a
  57
  58     sqrt z@(x:+y)  | z == __i0
  59                    =  __i0
  60                    | otherwise
  61                    =  u :+ (if y < __i0 then -v else v)
  62                       where (u,v) = if x < __i0 then (v2,u2) else (u2,v2)
  63                             v2    = abs y / (u2 * __i2)
  64                             u2    = sqrt ((magnitude z + abs x) / __i2)
  65
  66     sin (x:+y)     =  sin x * cosh y :+ cos x * sinh y
  67     cos (x:+y)     =  cos x * cosh y :+ (- sin x * sinh y)
  68     tan (x:+y)     =  (sinx*coshy:+cosx*sinhy)/(cosx*coshy:+(-sinx*sinhy))
  69                       where sinx  = sin x
  70                             cosx  = cos x
  71                             sinhy = sinh y
  72                             coshy = cosh y
  73
  74     sinh (x:+y)    =  cos y * sinh x :+ sin  y * cosh x
  75     cosh (x:+y)    =  cos y * cosh x :+ sin y * sinh x
  76     tanh (x:+y)    =  (cosy*sinhx:+siny*coshx)/(cosy*coshx:+siny*sinhx)
  77                       where siny  = sin y
  78                             cosy  = cos y
  79                             sinhx = sinh x
  80                             coshx = cosh x
  81
  82     asin z@(x:+y)  =  y2:+(-x2)
  83                       where (x2:+y2) = log (((-y):+x) + sqrt (__i1 - z*z))
  84     acos z@(x:+y)  =  y3:+(-x3)
  85                       where (x3:+y3) = log (z + ((-y2):+x2))
  86                             (x2:+y2)   = sqrt (__i1 - z*z)
  87     atan z@(x:+y)  =  y2:+(-x2)
  88                       where (x2:+y2) = log (((__i1 - y):+x) / sqrt (__i1 + z*z))
  89
  90     asinh z        =  log (z + sqrt (__i1 + z*z))
  91     acosh z        =  log (z + (z + __i1) * sqrt ((z - __i1)/(z + __i1)))
  92     atanh z        =  log ((__i1 + z) / sqrt (__i1 - z*z))
  93
  94
  95 instance (Text a) => Text (Complex a) where
  96
  97     -- magic fixity wired in: infix 6 :+
  98
  99     readsPrec p
 100       = readParen ( p > 6 )
 101           (\ r -> [ (x :+ y, s2) | (x,    s0) <- readsPrec 7 r,
 102                                    (":+", s1) <- lex s0,
 103                                    (y,    s2) <- readsPrec 7 s1 ])
 104     showsPrec d (a :+ b)
 105       = showParen (d > 6)
 106           (showsPrec 7 a . showString " :+ " . showsPrec 7 b)
 107
 108     readList    = _readList (readsPrec 0)
 109     showList    = _showList (showsPrec 0)
 110
 111 {-# SPECIALIZE instance Eq          (Complex Double) #-}
 112 {-# SPECIALIZE instance Num         (Complex Double) #-}
 113 {-# SPECIALIZE instance Fractional  (Complex Double) #-}
 114 {-# SPECIALIZE instance Floating    (Complex Double) #-}
 115 {-# SPECIALIZE instance Text        (Complex Double) #-}
 116
 117 #if defined(__UNBOXED_INSTANCES__)
 118
 119 {-# SPECIALIZE instance Eq          (Complex Double#) #-}
 120 {-# SPECIALIZE instance Num         (Complex Double#) #-}
 121 {-# SPECIALIZE instance Fractional  (Complex Double#) #-}
 122 {-# SPECIALIZE instance Floating    (Complex Double#) #-}
 123 {-# SPECIALIZE instance Text        (Complex Double#) #-}
 124
 125 #endif
 126
 127 -- ToDo: something for Binary
 128
 129 -- ToDo: Complex Double#  s/a{/a{Double#,?/
 130
 131 {-# GENERATE_SPECS realPart a{Double#} #-}
 132 realPart         :: Complex a -> a
 133 realPart (x:+y)  =  x
 134
 135 {-# GENERATE_SPECS imagPart a{Double#} #-}
 136 imagPart         :: Complex a -> a
 137 imagPart (x:+y)  =  y
 138
 139 {-# GENERATE_SPECS conjugate a{Double#,Double} #-}
 140 conjugate        :: (RealFloat a) => Complex a -> Complex a
 141 conjugate (x:+y) =  x :+ (-y)
 142
 143 {-# GENERATE_SPECS mkPolar a{Double#,Double} #-}
 144 mkPolar          :: (RealFloat a) => a -> a -> Complex a
 145 mkPolar r theta  =  r * cos theta :+ r * sin theta
 146
 147 {-# GENERATE_SPECS cis a{Double#,Double} #-}
 148 cis              :: (RealFloat a) => a -> Complex a
 149 cis theta        =  cos theta :+ sin theta
 150
 151 {-# GENERATE_SPECS polar a{Double#,Double} #-}
 152 polar            :: (RealFloat a) => Complex a -> (a,a)
 153 polar z          =  (magnitude z, phase z)
 154
 155 {-# GENERATE_SPECS magnitude a{Double#,Double} #-}
 156 magnitude :: (RealFloat a) => Complex a -> a
 157 magnitude (x:+y) =  scaleFloat k
 158                       (sqrt ((scaleFloat mk x)^2 + (scaleFloat mk y)^2))
 159                     where k  = max (exponent x) (exponent y)
 160                           mk = - k
 161
 162 {-# GENERATE_SPECS phase a{Double#,Double} #-}
 163 phase :: (RealFloat a) => Complex a -> a
 164 phase (x:+y) =  atan2 y x