module Number.Complex
(
T(real,imag),
imaginaryUnit,
fromReal,
(+:),
(-:),
scale,
exp,
quarterLeft,
quarterRight,
fromPolar,
cis,
signum,
toPolar,
magnitude,
phase,
conjugate,
propPolar,
Power(power),
defltPow,
) where
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 Test.QuickCheck (Arbitrary, arbitrary, coarbitrary)
import Control.Monad (liftM2)
import qualified Prelude as P
import PreludeBase
import NumericPrelude hiding (signum, exp, )
import Text.Show.HT (showsInfixPrec, )
import Text.Read.HT (readsInfixPrec, )
infix 6 +:, `Cons`
data T a
= Cons {real :: !a
,imag :: !a
}
deriving (Eq)
imaginaryUnit :: Ring.C a => T a
imaginaryUnit = zero +: one
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 (+:)
instance (Arbitrary a) => Arbitrary (T a) where
arbitrary = liftM2 Cons arbitrary arbitrary
coarbitrary = undefined
(+:) :: a -> a -> T a
(+:) = Cons
(-:) :: Additive.C a => a -> a -> T a
(-:) x y = Cons x (y)
conjugate :: (Additive.C a) => T a -> T a
conjugate (Cons x y) = Cons x (y)
scale :: (Ring.C a) => a -> T a -> T a
scale r (Cons x y) = Cons (r * x) (r * y)
exp :: (Trans.C a) => T a -> T a
exp (Cons x y) = scale (Trans.exp x) (cis y)
quarterRight, quarterLeft :: (Additive.C a) => T a -> T a
quarterRight (Cons x y) = Cons y (x)
quarterLeft (Cons x y) = Cons (y) x
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
fromPolar :: (Trans.C a) => a -> a -> T a
fromPolar r theta = scale r (cis theta)
cis :: (Trans.C a) => a -> T a
cis theta = Cons (cos theta) (sin theta)
propPolar :: (RealTrans.C a) => T a -> Bool
propPolar z = uncurry fromPolar (toPolar z) == z
floatMagnitude :: (P.RealFloat a, Algebraic.C a) => T a -> a
floatMagnitude (Cons x y) =
let k = max (P.exponent x) (P.exponent y)
mk = k
in P.scaleFloat k
(sqrt (P.scaleFloat mk x ^ 2 +
P.scaleFloat mk y ^ 2))
magnitude :: (Algebraic.C a) => T a -> a
magnitude = sqrt . magnitudeSqr
magnitudeSqr :: (Ring.C a) => T a -> a
magnitudeSqr (Cons x y) = x^2 + y^2
phase :: (RealTrans.C a, ZeroTestable.C a) => T a -> a
phase z =
if isZero z
then zero
else case z of (Cons x y) -> atan2 y x
toPolar :: (RealTrans.C a) => T a -> (a,a)
toPolar z = (magnitude z, phase z)
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
zero = Cons zero zero
(Cons x y) + (Cons x' y') = Cons (x+x') (y+y')
(Cons x y) (Cons x' y') = Cons (xx') (yy')
negate (Cons x y) = Cons (negate x) (negate y)
instance (Ring.C a) => Ring.C (T a) where
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
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))
instance (Integral.C a) => Integral.C (T a) where
divMod z z' =
let denom = magnitudeSqr z'
zBig = z * conjugate z'
re = divMod (real zBig) denom
im = divMod (imag zBig) denom
q = Cons (fst re) (fst im)
in (q, zq*z')
divModCent :: (Ord a, Integral.C a) => T a -> T a -> (T a, T a)
divModCent z z' =
let denom = magnitudeSqr 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', zq'*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 quarterRight 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 quarterLeft sgn' else sgn'
instance (Ord a, ZeroTestable.C a, Units.C a) => PID.C (T a) where
gcd = euclid modCent
extendedGCD = extendedEuclid divModCent
divide :: (Field.C a) => T a -> T a -> T a
divide (Cons x y) z'@(Cons x' y') =
let d = magnitudeSqr z'
in Cons ((x*x'+y*y') / d) ((y*x'x*y') / d)
floatDivide :: (P.RealFloat a, Field.C a) => T a -> T a -> T a
floatDivide (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)
instance (Field.C a) => Field.C (T a) where
(/) = divide
fromRational' = fromReal . fromRational'
class (Algebraic.C a) => (Power a) where
power :: Rational -> T a -> T a
defltPow :: (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 (Real.C a, Algebraic.C 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 (Real.C a, RealTrans.C a, Power a) =>
Trans.C (T a) where
pi = fromReal pi
exp = exp
log z = let (m,p) = toPolar z in Cons (log m) p
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 = quarterRight (log (quarterLeft z + sqrt (1 z^2)))
acos z = quarterRight (log (z + quarterLeft (sqrt (1 z^2))))
atan z@(Cons x y) = quarterRight (log (Cons (1y) x / sqrt (1+z^2)))
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
(+) = legacyInstance
(*) = legacyInstance
abs = legacyInstance
signum = legacyInstance
instance (Field.C a, Eq a, Show a) => P.Fractional (T a) where
fromRational = fromRational
(/) = legacyInstance