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