module Number.ComplexSquareRoot where

-- import qualified Algebra.Algebraic as Algebraic
import qualified Algebra.RealField as RealField
import qualified Algebra.RealRing as RealRing
-- import qualified Algebra.Field as Field
import qualified Algebra.Ring as Ring
import qualified Algebra.Additive as Additive
import qualified Algebra.ZeroTestable as ZeroTestable

import qualified Number.Complex as Complex

import Test.QuickCheck (Arbitrary, arbitrary, )

import Control.Monad (liftM2, )

import qualified NumericPrelude.Numeric as NP
import NumericPrelude.Numeric hiding (recip, )
import NumericPrelude.Base
import Prelude ()

{- |
Represent the square root of a complex number
without actually having to compute a square root.
If the Bool is False,
then the square root is represented with positive real part
or zero real part and positive imaginary part.
If the Bool is True the square root is negated.
-}
data T a = Cons Bool (Complex.T a)
   deriving (Show)

{- |
You must use @fmap@ only for number type conversion.
-}
instance Functor T where
   fmap f (Cons n x) = Cons n (fmap f x)

instance (ZeroTestable.C a) => ZeroTestable.C (T a) where
   isZero (Cons _b s) = isZero s

instance (ZeroTestable.C a, Eq a) => Eq (T a) where
   (Cons xb xs) == (Cons yb ys) =
      isZero xs && isZero ys  ||
      xb==yb && xs==ys

instance (Arbitrary a) => Arbitrary (T a) where
   arbitrary = liftM2 Cons arbitrary arbitrary


fromNumber :: (RealRing.C a) => Complex.T a -> T a
fromNumber x =
   Cons
      (case compare zero (Complex.real x) of
         LT -> False
         GT -> True
         EQ -> Complex.imag x < zero)
      (x^2)

-- htam:Wavelet.DyadicResultant.parityFlip
toNumber :: (RealRing.C a, Complex.Power a) => T a -> Complex.T a
toNumber (Cons n x) =
   case sqrt x of y -> if n then NP.negate y else y


one :: (Ring.C a) => T a
one = Cons False NP.one

inUpperHalfplane :: (Additive.C a, Ord a) => Complex.T a -> Bool
inUpperHalfplane x =
   case compare (Complex.imag x) zero of
      GT -> True
      LT -> False
      EQ -> Complex.real x < zero

mul, mulAlt, mulAlt2 :: (RealRing.C a) => T a -> T a -> T a
mul (Cons xb xs) (Cons yb ys) =
   let zs = xs*ys
   in  Cons
          ((xb /= yb) /=
             case (inUpperHalfplane xs,
                   inUpperHalfplane ys,
                   inUpperHalfplane zs) of
                (True,True,False) -> True
                (False,False,True) -> True
                _ -> False)
          zs

mulAlt (Cons xb xs) (Cons yb ys) =
   let zs = xs*ys
   in  Cons
          ((xb /= yb) /=
             let xi = Complex.imag xs
                 yi = Complex.imag ys
                 zi = Complex.imag zs
             in  (xi>=zero) /= (yi>=zero) &&
                 (xi>=zero) /= (zi>=zero))
          zs

mulAlt2 (Cons xb xs) (Cons yb ys) =
   let zs = xs*ys
   in  Cons
          ((xb /= yb) /=
             let xi = Complex.imag xs
                 yi = Complex.imag ys
                 zi = Complex.imag zs
             in  xi*yi<zero && xi*zi<zero)
          zs

div :: (RealField.C a) => T a -> T a -> T a
div x y = mul x (recip y)

recip :: (RealField.C a) => T a -> T a
recip (Cons b s) =
   Cons
      (b /= (Complex.imag s == zero && Complex.real s < zero))
      (NP.recip s)