{-# LANGUAGE DeriveDataTypeable #-} {-# LANGUAGE Trustworthy #-} {-# LANGUAGE StandaloneDeriving #-} {-# LANGUAGE DeriveGeneric #-} {-# LANGUAGE DeriveTraversable #-} ----------------------------------------------------------------------------- -- | -- Module : Data.Complex -- Copyright : (c) The University of Glasgow 2001 -- License : BSD-style (see the file libraries/base/LICENSE) -- -- Maintainer : libraries@haskell.org -- Stability : provisional -- Portability : portable -- -- Complex numbers. -- ----------------------------------------------------------------------------- module Data.Complex ( -- * Rectangular form Complex((:+)) , realPart , imagPart -- * Polar form , mkPolar , cis , polar , magnitude , phase -- * Conjugate , conjugate ) where import GHC.Base (Applicative (..)) import GHC.Generics (Generic, Generic1) import GHC.Float (Floating(..)) import Data.Data (Data) import Foreign (Storable, castPtr, peek, poke, pokeElemOff, peekElemOff, sizeOf, alignment) infix 6 :+ -- ----------------------------------------------------------------------------- -- The Complex type -- | Complex numbers are an algebraic type. -- -- For a complex number @z@, @'abs' z@ is a number with the magnitude of @z@, -- but oriented in the positive real direction, whereas @'signum' z@ -- has the phase of @z@, but unit magnitude. -- -- The 'Foldable' and 'Traversable' instances traverse the real part first. data Complex a = !a :+ !a -- ^ forms a complex number from its real and imaginary -- rectangular components. deriving (Eq, Show, Read, Data, Generic, Generic1 , Functor, Foldable, Traversable) -- ----------------------------------------------------------------------------- -- Functions over Complex -- | Extracts the real part of a complex number. realPart :: Complex a -> a realPart (x :+ _) = x -- | Extracts the imaginary part of a complex number. imagPart :: Complex a -> a imagPart (_ :+ y) = y -- | The conjugate of a complex number. {-# SPECIALISE conjugate :: Complex Double -> Complex Double #-} conjugate :: Num a => Complex a -> Complex a conjugate (x:+y) = x :+ (-y) -- | Form a complex number from polar components of magnitude and phase. {-# SPECIALISE mkPolar :: Double -> Double -> Complex Double #-} mkPolar :: Floating a => a -> a -> Complex a mkPolar r theta = r * cos theta :+ r * sin theta -- | @'cis' t@ is a complex value with magnitude @1@ -- and phase @t@ (modulo @2*'pi'@). {-# SPECIALISE cis :: Double -> Complex Double #-} cis :: Floating a => a -> Complex a cis theta = cos theta :+ sin theta -- | The function 'polar' takes a complex number and -- returns a (magnitude, phase) pair in canonical form: -- the magnitude is nonnegative, and the phase in the range @(-'pi', 'pi']@; -- if the magnitude is zero, then so is the phase. {-# SPECIALISE polar :: Complex Double -> (Double,Double) #-} polar :: (RealFloat a) => Complex a -> (a,a) polar z = (magnitude z, phase z) -- | The nonnegative magnitude of a complex number. {-# SPECIALISE magnitude :: Complex Double -> Double #-} magnitude :: (RealFloat a) => Complex a -> a magnitude (x:+y) = scaleFloat k (sqrt (sqr (scaleFloat mk x) + sqr (scaleFloat mk y))) where k = max (exponent x) (exponent y) mk = - k sqr z = z * z -- | The phase of a complex number, in the range @(-'pi', 'pi']@. -- If the magnitude is zero, then so is the phase. {-# SPECIALISE phase :: Complex Double -> Double #-} phase :: (RealFloat a) => Complex a -> a phase (0 :+ 0) = 0 -- SLPJ July 97 from John Peterson phase (x:+y) = atan2 y x -- ----------------------------------------------------------------------------- -- Instances of Complex -- | @since 2.01 instance (RealFloat a) => Num (Complex a) where {-# SPECIALISE instance Num (Complex Float) #-} {-# SPECIALISE instance Num (Complex Double) #-} (x:+y) + (x':+y') = (x+x') :+ (y+y') (x:+y) - (x':+y') = (x-x') :+ (y-y') (x:+y) * (x':+y') = (x*x'-y*y') :+ (x*y'+y*x') negate (x:+y) = negate x :+ negate y abs z = magnitude z :+ 0 signum (0:+0) = 0 signum z@(x:+y) = x/r :+ y/r where r = magnitude z fromInteger n = fromInteger n :+ 0 -- | @since 2.01 instance (RealFloat a) => Fractional (Complex a) where {-# SPECIALISE instance Fractional (Complex Float) #-} {-# SPECIALISE instance Fractional (Complex Double) #-} (x:+y) / (x':+y') = (x*x''+y*y'') / d :+ (y*x''-x*y'') / d where x'' = scaleFloat k x' y'' = scaleFloat k y' k = - max (exponent x') (exponent y') d = x'*x'' + y'*y'' fromRational a = fromRational a :+ 0 -- | @since 2.01 instance (RealFloat a) => Floating (Complex a) where {-# SPECIALISE instance Floating (Complex Float) #-} {-# SPECIALISE instance Floating (Complex Double) #-} pi = pi :+ 0 exp (x:+y) = expx * cos y :+ expx * sin y where expx = exp x log z = log (magnitude z) :+ phase z x ** y = case (x,y) of (_ , (0:+0)) -> 1 :+ 0 ((0:+0), (exp_re:+_)) -> case compare exp_re 0 of GT -> 0 :+ 0 LT -> inf :+ 0 EQ -> nan :+ nan ((re:+im), (exp_re:+_)) | (isInfinite re || isInfinite im) -> case compare exp_re 0 of GT -> inf :+ 0 LT -> 0 :+ 0 EQ -> nan :+ nan | otherwise -> exp (log x * y) where inf = 1/0 nan = 0/0 sqrt (0:+0) = 0 sqrt z@(x:+y) = u :+ (if y < 0 then -v else v) where (u,v) = if x < 0 then (v',u') else (u',v') v' = abs y / (u'*2) u' = sqrt ((magnitude z + abs x) / 2) sin (x:+y) = sin x * cosh y :+ cos x * sinh y cos (x:+y) = cos x * cosh y :+ (- sin x * sinh y) tan (x:+y) = (sinx*coshy:+cosx*sinhy)/(cosx*coshy:+(-sinx*sinhy)) where sinx = sin x cosx = cos x sinhy = sinh y coshy = cosh y sinh (x:+y) = cos y * sinh x :+ sin y * cosh x cosh (x:+y) = cos y * cosh x :+ sin y * sinh x tanh (x:+y) = (cosy*sinhx:+siny*coshx)/(cosy*coshx:+siny*sinhx) where siny = sin y cosy = cos y sinhx = sinh x coshx = cosh x asin z@(x:+y) = y':+(-x') where (x':+y') = log (((-y):+x) + sqrt (1 - z*z)) acos z = y'':+(-x'') where (x'':+y'') = log (z + ((-y'):+x')) (x':+y') = sqrt (1 - z*z) atan z@(x:+y) = y':+(-x') where (x':+y') = log (((1-y):+x) / sqrt (1+z*z)) asinh z = log (z + sqrt (1+z*z)) -- Take care to allow (-1)::Complex, fixing #8532 acosh z = log (z + (sqrt \$ z+1) * (sqrt \$ z-1)) atanh z = 0.5 * log ((1.0+z) / (1.0-z)) log1p x@(a :+ b) | abs a < 0.5 && abs b < 0.5 , u <- 2*a + a*a + b*b = log1p (u/(1 + sqrt(u+1))) :+ atan2 (1 + a) b | otherwise = log (1 + x) {-# INLINE log1p #-} expm1 x@(a :+ b) | a*a + b*b < 1 , u <- expm1 a , v <- sin (b/2) , w <- -2*v*v = (u*w + u + w) :+ (u+1)*sin b | otherwise = exp x - 1 {-# INLINE expm1 #-} -- | @since 4.8.0.0 instance Storable a => Storable (Complex a) where sizeOf a = 2 * sizeOf (realPart a) alignment a = alignment (realPart a) peek p = do q <- return \$ castPtr p r <- peek q i <- peekElemOff q 1 return (r :+ i) poke p (r :+ i) = do q <-return \$ (castPtr p) poke q r pokeElemOff q 1 i -- | @since 4.9.0.0 instance Applicative Complex where pure a = a :+ a f :+ g <*> a :+ b = f a :+ g b liftA2 f (x :+ y) (a :+ b) = f x a :+ f y b -- | @since 4.9.0.0 instance Monad Complex where a :+ b >>= f = realPart (f a) :+ imagPart (f b) -- ----------------------------------------------------------------------------- -- Rules on Complex {-# RULES "realToFrac/a->Complex Double" realToFrac = \x -> realToFrac x :+ (0 :: Double) "realToFrac/a->Complex Float" realToFrac = \x -> realToFrac x :+ (0 :: Float) #-}