{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-}
{- |
Module      :  Number.Complex
Copyright   :  (c) The University of Glasgow 2001
License     :  BSD-style (see the file libraries/base/LICENSE)

Maintainer  :  numericprelude@henning-thielemann.de
Stability   :  provisional
Portability :  portable (?)

Complex numbers.
-}

module Number.Complex
	(
	-- * Cartesian form
	  T(real,imag)
        , fromReal

	, (+:)
	, (-:)
	-- * Polar form
	, fromPolar
	, cis
        , signum
	, toPolar
	, magnitude
	, phase
        , Polar
	, defltMagnitude
	, defltPhase
	-- * Conjugate
	, conjugate

        -- * Properties
        , propPolar

        -- * Auxiliary classes
        , Divisible(divide)
        , defltDiv
        , Power(power)
        , defltPow
        )  where

import qualified Number.Ratio as Ratio

import qualified Algebra.NormedSpace.Euclidean as NormedEuc
import qualified Algebra.NormedSpace.Sum       as NormedSum
import qualified Algebra.NormedSpace.Maximum   as NormedMax

import qualified Algebra.VectorSpace        as VectorSpace
import qualified Algebra.Module             as Module
import qualified Algebra.Vector             as Vector
import qualified Algebra.RealTranscendental as RealTrans
import qualified Algebra.Transcendental     as Trans
import qualified Algebra.Algebraic          as Algebraic
import qualified Algebra.Field              as Field
import qualified Algebra.Units              as Units
import qualified Algebra.PrincipalIdealDomain as PID
import qualified Algebra.IntegralDomain     as Integral
import qualified Algebra.Real               as Real
import qualified Algebra.Ring               as Ring
import qualified Algebra.Additive           as Additive
import qualified Algebra.ZeroTestable       as ZeroTestable
import qualified Algebra.Indexable          as Indexable

import Algebra.ZeroTestable(isZero)
import Algebra.Module((*>))
import Algebra.Algebraic((^/))

import qualified Prelude as P
import PreludeBase
import NumericPrelude hiding (signum)
import NumericPrelude.Text (showsInfixPrec, readsInfixPrec)


-- import qualified Data.Typeable as Ty

infix  6  +:, `Cons`

{- * The Complex type -}

-- | Complex numbers are an algebraic type.
data T a
  = Cons {real :: !a   -- ^ real part
         ,imag :: !a   -- ^ imaginary part
         }
  deriving (Eq)

fromReal :: Additive.C a => a -> T a
fromReal x = Cons x zero


plusPrec :: Int
plusPrec = 6

instance (Show a) => Show (T a) where
   showsPrec prec (Cons x y) = showsInfixPrec "+:" plusPrec prec x y

instance (Read a) => Read (T a) where
   readsPrec prec = readsInfixPrec "+:" plusPrec prec (+:)



{- * Functions -}

-- | Construct a complex number from real and imaginary part.
(+:) :: a -> a -> T a
(+:) = Cons

-- | Construct a complex number with negated imaginary part.
(-:) :: Additive.C a => a -> a -> T a
(-:) x y = Cons x (-y)

-- | The conjugate of a complex number.
{-# SPECIALISE conjugate :: T Double -> T Double #-}
conjugate	 :: (Additive.C a) => T a -> T a
conjugate (Cons x y) =  Cons x (-y)

-- | Scale a complex number by a real number.
{-# SPECIALISE scale :: Double -> T Double -> T Double #-}
scale		 :: (Ring.C a) => a -> T a -> T a
scale r (Cons x y) =  Cons (r * x) (r * y)

-- | Turn the point one quarter to the right.
orthoRight, orthoLeft :: (Additive.C a) => T a -> T a
orthoRight (Cons x y) = Cons   y  (-x)
orthoLeft  (Cons x y) = Cons (-y)   x

{- | Scale a complex number to magnitude 1.

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.
-}
{-# SPECIALISE signum :: T Double -> T Double #-}
signum :: (Algebraic.C a, NormedEuc.C a a, ZeroTestable.C a) => T a -> T a
signum z =
   if isZero z
     then zero
     else scale (recip (NormedEuc.norm z)) z

-- | Form a complex number from polar components of magnitude and phase.
{-# SPECIALISE fromPolar :: Double -> Double -> T Double #-}
fromPolar		 :: (Trans.C a) => a -> a -> T a
fromPolar r theta	 =  scale r (cis theta)

-- | @'cis' t@ is a complex value with magnitude @1@
-- and phase @t@ (modulo @2*'pi'@).
{-# SPECIALISE cis :: Double -> T Double #-}
cis		 :: (Trans.C a) => a -> T a
cis theta	 =  Cons (cos theta) (sin theta)

propPolar :: (Polar a, RealTrans.C a) => T a -> Bool
propPolar z  =  uncurry fromPolar (toPolar z) == z


-- | The nonnegative magnitude of a complex number.
floatMagnitude :: (P.RealFloat a, Algebraic.C a) => T a -> a
floatMagnitude (Cons x y) =  P.scaleFloat k
		     (sqrt (P.scaleFloat mk x ^ 2 +
                            P.scaleFloat mk y ^ 2))
		    where k  = max (P.exponent x) (P.exponent y)
		          mk = - k

defltMagnitude :: (Algebraic.C a) => T a -> a
defltMagnitude = sqrt . defltMagnitudeSqr

-- like NormedEuc.normSqr with lifted class constraints
defltMagnitudeSqr :: (Ring.C a) => T a -> a
defltMagnitudeSqr (Cons x y) = x^2 + y^2

-- | The phase of a complex number, in the range @(-'pi', 'pi']@.
-- If the magnitude is zero, then so is the phase.
defltPhase :: (RealTrans.C a, ZeroTestable.C a) => T a -> a
defltPhase z =
   if isZero z
     then zero   -- SLPJ July 97 from John Peterson
     else case z of (Cons x y) -> atan2 y x


{- |
Minimal implementation: toPolar or (magnitude and phase),
usually the instance definition

@
magnitude = defltMagnitude
phase     = defltPhase
@

is enough.

This class requires transcendent number types
although 'magnitude' can be computed algebraically.
-}
class RealTrans.C a => Polar a where
   {- |
   The function 'toPolar' 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 toPolar :: T Double -> (Double,Double) #--}
   toPolar   :: T a -> (a,a)
   toPolar z = (magnitude z, phase z)

   {--# SPECIALISE magnitude :: T Double -> Double #--}
   magnitude :: T a -> a
   magnitude = fst . toPolar

   {--# SPECIALISE phase :: T Double -> Double #--}
   phase     :: T a -> a
   phase     = snd . toPolar

instance Polar Float where
   magnitude = floatMagnitude
   phase     = defltPhase

instance Polar Double where
   magnitude = floatMagnitude
   phase     = defltPhase



{- * Instances of T -}

{-
complexTc = Ty.mkTyCon "Complex.T"
instance Ty.Typeable1 T where { typeOf1 _ = Ty.mkTyConApp complexTc [] }
instance Ty.Typeable a => Ty.Typeable (T a) where { typeOf = Ty.typeOfDefault }
-}

instance  (Indexable.C a) => Indexable.C (T a) where
    compare (Cons x y) (Cons x' y')  =  Indexable.compare (x,y) (x',y')

instance  (ZeroTestable.C a) => ZeroTestable.C (T a)  where
    isZero (Cons x y)  = isZero x && isZero y

instance  (Additive.C a) => Additive.C (T a)  where
    {-# SPECIALISE instance Additive.C (T Float) #-}
    {-# SPECIALISE instance Additive.C (T Double) #-}
    zero			=  Cons zero zero
    (Cons x y) + (Cons x' y')	=  Cons (x+x') (y+y')
    (Cons x y) - (Cons x' y')	=  Cons (x-x') (y-y')
    negate (Cons x y)		=  Cons (negate x) (negate y)

instance  (Ring.C a) => Ring.C (T a)  where
    {-# SPECIALISE instance Ring.C (T Float) #-}
    {-# SPECIALISE instance Ring.C (T Double) #-}
    one				=  Cons one zero
    (Cons x y) * (Cons x' y')	=  Cons (x*x'-y*y') (x*y'+y*x')
    fromInteger			=  fromReal . fromInteger

instance Vector.C T where
   zero  = zero
   (<+>) = (+)
   (*>)  = scale

-- | The '(*>)' method can't replace 'scale'
--   because it requires the Algebra.Module constraint
instance (Module.C a b) => Module.C a (T b) where
   s *> (Cons x y)  = Cons (s *> x) (s *> y)

instance (VectorSpace.C a b) => VectorSpace.C a (T b)

instance (Additive.C a, NormedSum.C a v) => NormedSum.C a (T v) where
   norm x = NormedSum.norm (real x) + NormedSum.norm (imag x)

instance (NormedEuc.Sqr a b) => NormedEuc.Sqr a (T b) where
   normSqr x = NormedEuc.normSqr (real x) + NormedEuc.normSqr (imag x)

instance (Algebraic.C a, NormedEuc.Sqr a b) => NormedEuc.C a (T b) where
   norm = NormedEuc.defltNorm

instance (Ord a, NormedMax.C a v) => NormedMax.C a (T v) where
   norm x = max (NormedMax.norm (real x)) (NormedMax.norm (imag x))


{-
  In this implementation the complex plane is structured
  as an orthogonal grid induced by the divisor z'.
  The coordinate of a cell within this grid is returned as quotient
  and the position with a cell is returned as remainder.
  The magnitude of the remainder might be larger than that of the divisor
  thus the Euclidean algorithm can fail.
-}

instance  (Integral.C a) => Integral.C (T a)  where
    divMod z z' =
       let denom = defltMagnitudeSqr z'
           zBig  = z * conjugate z'
           re    = divMod (real zBig) denom
           im    = divMod (imag zBig) denom
           q     = Cons (fst re) (fst im)
       in  (q, z-q*z')


{-
  This variant of divMod tries to come close to the origin.
  Thus the remainder has smaller magnitude than the divisor.
  This variant of divModCent can be used for Euclidean's algorithm.
-}
divModCent :: (Ord a, Integral.C a) => T a -> T a -> (T a, T a)
divModCent z z' =
   let denom = defltMagnitudeSqr z'
       zBig  = z * conjugate z'
       re    = divMod (real zBig) denom
       im    = divMod (imag zBig) denom
       q  = Cons (fst re) (fst im)
       r  = Cons (snd re) (snd im)
       q' = Cons
              (real q + if 2 * real r > denom then one else zero)
              (imag q + if 2 * imag r > denom then one else zero)
   in  (q', z-q'*z')

modCent :: (Ord a, Integral.C a) => T a -> T a -> T a
modCent z z' = snd (divModCent z z')

instance  (Ord a, Units.C a) => Units.C (T a)  where
    isUnit (Cons x y) =
       isUnit x && y==zero  ||
       isUnit y && x==zero
    stdAssociate z@(Cons x y) =
       let z' = if y<0  ||  y==0 && x<0 then negate z else z
       in  if real z'<=0 then orthoRight z' else z'
    stdUnit z@(Cons x y) =
       if z==zero
         then 1
         else
           let (x',sgn') = if y<0  ||  y==0 && x<0
                             then (negate x, -1)
                             else (x, 1)
           in  if x'<=0 then orthoLeft sgn' else sgn'


instance  (Ord a, ZeroTestable.C a, Units.C a) => PID.C (T a) where
   gcd         = euclid modCent
   extendedGCD = extendedEuclid divModCent


defltDiv :: (Field.C a) => T a -> T a -> T a
defltDiv (Cons x y) z'@(Cons x' y') =
   let d = defltMagnitudeSqr z'
   in  Cons ((x*x'+y*y') / d) ((y*x'-x*y') / d)

-- | Special implementation of @(\/)@ for floating point numbers
--   which prevent intermediate overflows.
floatDiv :: (P.RealFloat a, Field.C a) => T a -> T a -> T a
floatDiv (Cons x y) (Cons x' y') =
   let k   = - max (P.exponent x') (P.exponent y')
       x'' = P.scaleFloat k x'
       y'' = P.scaleFloat k y'
       d   = x'*x'' + y'*y''
   in  Cons ((x*x''+y*y'') / d) ((y*x''-x*y'') / d)

{-|
   In order to have an efficient implementation
   for both complex floats and exact complex numbers,
   we define the intermediate class Complex.Divisible
   which in fact implements the complex division.
   This way we avoid overlapping and undecidable instances.
   In most cases it should suffice to define
   an instance of Complex.Divisible with no method implementation
   for each instance of Fractional.
-}
class (Field.C a) => Divisible a where
    divide :: T a -> T a -> T a
    divide  =  defltDiv

instance  Divisible Float  where
    divide  =  floatDiv

instance  Divisible Double  where
    divide  =  floatDiv

instance  (PID.C a) => Divisible (Ratio.T a)


instance  (Divisible a) => Field.C (T a)  where
    (/)			=  divide
    fromRational'	=  fromReal . fromRational'

{-|
   We like to build the Complex Algebraic instance
   on top of the Algebraic instance of the scalar type.
   This poses no problem to 'sqrt'.
   However, 'Number.Complex.root' requires computing the complex argument
   which is a transcendent operation.
   In order to keep the type class dependencies clean
   for more sophisticated algebraic number types,
   we introduce a type class which actually performs the radix operation.
-}
class (Algebraic.C a) => (Power a) where
    power  ::  Rational -> T a -> T a


defltPow :: (Polar a, RealTrans.C a) =>
    Rational -> T a -> T a
defltPow r x =
    let (mag,arg) = toPolar x
    in  fromPolar (mag ^/ r)
                  (arg * fromRational' r)


instance  Power Float where
    power  =  defltPow

instance  Power Double where
    power  =  defltPow


instance  (Polar a, Real.C a, Algebraic.C a, Divisible a, Power a) =>
          Algebraic.C (T a)  where
    sqrt z@(Cons x y)  =  if z == zero
                            then zero
                            else
                              let v'    = abs y / (u'*2)
                                  u'    = sqrt ((magnitude z + abs x) / 2)
                                  (u,v) = if x < 0 then (v',u') else (u',v')
                              in  Cons u (if y < 0 then -v else v)
    (^/) = flip power


instance  (Polar a, Real.C a, RealTrans.C a, Divisible a, Power a) =>
          Trans.C (T a)  where
    {-# SPECIALISE instance Trans.C (T Float) #-}
    {-# SPECIALISE instance Trans.C (T Double) #-}
    pi                 =  fromReal pi
    exp (Cons x y)     =  scale (exp x) (cis y)
    log z              =  let (m,p) = toPolar z in Cons (log m) p

    -- use defaults for tan, tanh

    sin (Cons x y)     =  Cons (sin x * cosh y) (  cos x * sinh y)
    cos (Cons x y)     =  Cons (cos x * cosh y) (- sin x * sinh y)

    sinh (Cons x y)    =  Cons (cos y * sinh x) (sin y * cosh x)
    cosh (Cons x y)    =  Cons (cos y * cosh x) (sin y * sinh x)

    asin z             =  orthoRight (log (orthoLeft z + sqrt (1 - z^2)))
    acos z             =  orthoRight (log (z + orthoLeft (sqrt (1 - z^2))))
    atan z@(Cons x y)  =  orthoRight (log (Cons (1-y) x / sqrt (1+z^2)))

{- use the default implementation
    asinh z        =  log (z + sqrt (1+z^2))
    acosh z        =  log (z + (z+1) * sqrt ((z-1)/(z+1)))
    atanh z        =  log ((1+z) / sqrt (1-z^2))
-}


{- * legacy instances -}

legacyInstance :: a
legacyInstance =
   error "legacy Ring.C instance for simple input of numeric literals"

instance (Ring.C a, Eq a, Show a) => P.Num (T a) where
   fromInteger = fromReal . fromInteger
   negate = negate -- for unary minus
   (+)    = legacyInstance
   (*)    = legacyInstance
   abs    = legacyInstance
   signum = legacyInstance

instance (Ring.C a, Eq a, Show a, Divisible a) => P.Fractional (T a) where
   fromRational = fromRational
   (/) = legacyInstance