-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.Data.AlgebraicNumber.Complex
-- Copyright   :  (c) Masahiro Sakai 2013
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  portable
--
-- Algebraic Numbers
--
-----------------------------------------------------------------------------
module ToySolver.Data.AlgebraicNumber.Complex
    (
    -- * Rectangular form
      AComplex((:+))

    -- * Properties
    , realPart          -- :: AComplex -> AReal
    , imagPart          -- :: AComplex -> AReal
    , magnitude         -- :: AComplex -> AReal
    , minimalPolynomial -- :: AComplex -> UPolynomial Rational

    -- * Operations
    , conjugate         -- :: AComplex -> AComplex
    ) where

import Control.Monad
import qualified Data.Sign as Sign
import qualified Data.Map as Map
import Data.Maybe

import qualified ToySolver.Arith.CAD as CAD
import qualified ToySolver.Data.AlgebraicNumber.Real as AReal
import ToySolver.Data.AlgebraicNumber.Real (AReal)
import qualified ToySolver.Data.AlgebraicNumber.Root as Root
import qualified ToySolver.Data.Polynomial as P
import ToySolver.Data.Polynomial (Polynomial, UPolynomial, X (..))

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.
data AComplex = !AReal :+ !AReal
    deriving (AComplex -> AComplex -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: AComplex -> AComplex -> Bool
$c/= :: AComplex -> AComplex -> Bool
== :: AComplex -> AComplex -> Bool
$c== :: AComplex -> AComplex -> Bool
Eq, Int -> AComplex -> ShowS
[AComplex] -> ShowS
AComplex -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [AComplex] -> ShowS
$cshowList :: [AComplex] -> ShowS
show :: AComplex -> String
$cshow :: AComplex -> String
showsPrec :: Int -> AComplex -> ShowS
$cshowsPrec :: Int -> AComplex -> ShowS
Show)

-- -----------------------------------------------------------------------------
-- Functions over Complex

-- | Extracts the real part of a complex number.
realPart :: AComplex -> AReal
realPart :: AComplex -> AReal
realPart (AReal
x :+ AReal
_) = AReal
x

-- | Extracts the imaginary part of a complex number.
imagPart :: AComplex -> AReal
imagPart :: AComplex -> AReal
imagPart (AReal
_ :+ AReal
y) = AReal
y

-- | The conjugate of a complex number.
conjugate :: AComplex -> AComplex
conjugate :: AComplex -> AComplex
conjugate (AReal
x :+ AReal
y) =  AReal
x AReal -> AReal -> AComplex
:+ (-AReal
y)

-- | The nonnegative magnitude of a complex number.
magnitude :: AComplex -> AReal
magnitude :: AComplex -> AReal
magnitude (AReal
x :+ AReal
y) = Integer -> AReal -> AReal
AReal.nthRoot Integer
2 (AReal
xforall a. Num a => a -> a -> a
*AReal
x forall a. Num a => a -> a -> a
+ AReal
yforall a. Num a => a -> a -> a
*AReal
y)

-- -----------------------------------------------------------------------------
-- Instances of Complex

instance  Num AComplex where
    (AReal
x:+AReal
y) + :: AComplex -> AComplex -> AComplex
+ (AReal
x':+AReal
y')   =  (AReal
xforall a. Num a => a -> a -> a
+AReal
x') AReal -> AReal -> AComplex
:+ (AReal
yforall a. Num a => a -> a -> a
+AReal
y')
    (AReal
x:+AReal
y) - :: AComplex -> AComplex -> AComplex
- (AReal
x':+AReal
y')   =  (AReal
xforall a. Num a => a -> a -> a
-AReal
x') AReal -> AReal -> AComplex
:+ (AReal
yforall a. Num a => a -> a -> a
-AReal
y')
    (AReal
x:+AReal
y) * :: AComplex -> AComplex -> AComplex
* (AReal
x':+AReal
y')   =  (AReal
xforall a. Num a => a -> a -> a
*AReal
x'forall a. Num a => a -> a -> a
-AReal
yforall a. Num a => a -> a -> a
*AReal
y') AReal -> AReal -> AComplex
:+ (AReal
xforall a. Num a => a -> a -> a
*AReal
y'forall a. Num a => a -> a -> a
+AReal
yforall a. Num a => a -> a -> a
*AReal
x')
    negate :: AComplex -> AComplex
negate (AReal
x:+AReal
y)       =  forall a. Num a => a -> a
negate AReal
x AReal -> AReal -> AComplex
:+ forall a. Num a => a -> a
negate AReal
y
    abs :: AComplex -> AComplex
abs AComplex
z               =  AComplex -> AReal
magnitude AComplex
z AReal -> AReal -> AComplex
:+ AReal
0
    signum :: AComplex -> AComplex
signum (AReal
0:+AReal
0)       =  AComplex
0
    signum z :: AComplex
z@(AReal
x:+AReal
y)     =  AReal
xforall a. Fractional a => a -> a -> a
/AReal
r AReal -> AReal -> AComplex
:+ AReal
yforall a. Fractional a => a -> a -> a
/AReal
r  where r :: AReal
r = AComplex -> AReal
magnitude AComplex
z
    fromInteger :: Integer -> AComplex
fromInteger Integer
n       =  forall a. Num a => Integer -> a
fromInteger Integer
n AReal -> AReal -> AComplex
:+ AReal
0

instance  Fractional AComplex  where
    (AReal
x:+AReal
y) / :: AComplex -> AComplex -> AComplex
/ (AReal
x':+AReal
y')   =  (AReal
xforall a. Num a => a -> a -> a
*AReal
x'forall a. Num a => a -> a -> a
+AReal
yforall a. Num a => a -> a -> a
*AReal
y') forall a. Fractional a => a -> a -> a
/ AReal
d AReal -> AReal -> AComplex
:+ (AReal
yforall a. Num a => a -> a -> a
*AReal
x'forall a. Num a => a -> a -> a
-AReal
xforall a. Num a => a -> a -> a
*AReal
y') forall a. Fractional a => a -> a -> a
/ AReal
d
                           where d :: AReal
d   = AReal
x'forall a. Num a => a -> a -> a
*AReal
x' forall a. Num a => a -> a -> a
+ AReal
y'forall a. Num a => a -> a -> a
*AReal
y'
    fromRational :: Rational -> AComplex
fromRational Rational
a      =  forall a. Fractional a => Rational -> a
fromRational Rational
a AReal -> AReal -> AComplex
:+ AReal
0

-- -----------------------------------------------------------------------------

-- | The polynomial of which the algebraic number is root.
minimalPolynomial :: AComplex -> UPolynomial Rational
minimalPolynomial :: AComplex -> UPolynomial Rational
minimalPolynomial z :: AComplex
z@(AReal
x :+ AReal
y) =
  forall a. [a] -> a
head [UPolynomial Rational
q | (UPolynomial Rational
q,Integer
_) <- forall a. Factor a => a -> [(a, Integer)]
P.factor UPolynomial Rational
p, forall t. Degree t => t -> Integer
P.deg UPolynomial Rational
q forall a. Ord a => a -> a -> Bool
> Integer
0, forall k v. Num k => (v -> k) -> Polynomial k v -> k
P.eval (\X
X -> AComplex
z) (forall k1 k v.
(Eq k1, Num k1) =>
(k -> k1) -> Polynomial k v -> Polynomial k1 v
P.mapCoeff forall a. Fractional a => Rational -> a
fromRational UPolynomial Rational
q) forall a. Eq a => a -> a -> Bool
== AComplex
0]
  where
    p :: UPolynomial Rational
p = forall k.
(Fractional k, Ord k) =>
Polynomial k Int -> [Polynomial k Int] -> UPolynomial k
Root.findPoly (forall a v. Var a v => v -> a
P.var Int
a forall a. Num a => a -> a -> a
+ forall a v. Var a v => v -> a
P.var Int
b forall a. Num a => a -> a -> a
* forall a v. Var a v => v -> a
P.var Int
i) [Polynomial Rational Int]
ps
      where
        a, b, i :: Int
        i :: Int
i = Int
0
        a :: Int
a = Int
1
        b :: Int
b = Int
2
        ps :: [Polynomial Rational Int]
ps =
          [ (forall a v. Var a v => v -> a
P.var Int
i) forall a b. (Num a, Integral b) => a -> b -> a
^ Integer
2 forall a. Num a => a -> a -> a
+ Polynomial Rational Int
1
          , forall k v2 v1.
(Eq k, Num k, Ord v2) =>
Polynomial k v1 -> (v1 -> Polynomial k v2) -> Polynomial k v2
P.subst (AReal -> UPolynomial Rational
AReal.minimalPolynomial AReal
x) (\X
X -> forall a v. Var a v => v -> a
P.var Int
a)
          , forall k v2 v1.
(Eq k, Num k, Ord v2) =>
Polynomial k v1 -> (v1 -> Polynomial k v2) -> Polynomial k v2
P.subst (AReal -> UPolynomial Rational
AReal.minimalPolynomial AReal
y) (\X
X -> forall a v. Var a v => v -> a
P.var Int
b)
          ]

-- -----------------------------------------------------------------------------

-- | Roots of the polynomial
roots :: UPolynomial Rational -> [AComplex]
roots :: UPolynomial Rational -> [AComplex]
roots UPolynomial Rational
f = do
  let cs1 :: [(Polynomial Rational Int, [Sign])]
cs1 = [ (Polynomial Rational Int
u, [Sign
Sign.Zero]), (Polynomial Rational Int
v, [Sign
Sign.Zero]) ]
  ([(Polynomial Rational Int, [Sign])]
cs2, [Cell (Polynomial Rational Int)]
cells2) <- forall v.
(Ord v, Show v, PrettyVar v) =>
[(UPolynomial (Polynomial Rational v), [Sign])]
-> [([(Polynomial Rational v, [Sign])],
     [Cell (Polynomial Rational v)])]
CAD.project' [(forall k v.
(Ord k, Num k, Ord v) =>
Polynomial k v -> v -> UPolynomial (Polynomial k v)
P.toUPolynomialOf Polynomial Rational Int
p Int
0, [Sign]
ss) | (Polynomial Rational Int
p,[Sign]
ss) <- [(Polynomial Rational Int, [Sign])]
cs1]
  ([(Polynomial Rational Int, [Sign])]
cs3, [Cell (Polynomial Rational Int)]
cells3) <- forall v.
(Ord v, Show v, PrettyVar v) =>
[(UPolynomial (Polynomial Rational v), [Sign])]
-> [([(Polynomial Rational v, [Sign])],
     [Cell (Polynomial Rational v)])]
CAD.project' [(forall k v.
(Ord k, Num k, Ord v) =>
Polynomial k v -> v -> UPolynomial (Polynomial k v)
P.toUPolynomialOf Polynomial Rational Int
p Int
1, [Sign]
ss) | (Polynomial Rational Int
p,[Sign]
ss) <- [(Polynomial Rational Int, [Sign])]
cs2]
  forall (f :: * -> *). Alternative f => Bool -> f ()
guard forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *). Foldable t => t Bool -> Bool
and [forall a. Real a => a -> Sign
Sign.signOf Rational
v forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [Sign]
ss | (Polynomial Rational Int
p,[Sign]
ss) <- [(Polynomial Rational Int, [Sign])]
cs3, let v :: Rational
v = forall k v. Num k => (v -> k) -> Polynomial k v -> k
P.eval (\Int
_ -> forall a. HasCallStack => a
undefined) Polynomial Rational Int
p]

  let m3 :: Map k a
m3 = forall k a. Map k a
Map.empty
  AReal
yval <- forall a. [Maybe a] -> [a]
catMaybes [forall v.
Ord v =>
Model v -> Cell (Polynomial Rational v) -> Maybe AReal
CAD.findSample forall k a. Map k a
m3 Cell (Polynomial Rational Int)
cell | Cell (Polynomial Rational Int)
cell <- [Cell (Polynomial Rational Int)]
cells3]
  let m2 :: Map Int AReal
m2 = forall k a. Ord k => k -> a -> Map k a -> Map k a
Map.insert Int
1 AReal
yval forall k a. Map k a
m3
  AReal
xval <- forall a. [Maybe a] -> [a]
catMaybes [forall v.
Ord v =>
Model v -> Cell (Polynomial Rational v) -> Maybe AReal
CAD.findSample Map Int AReal
m2 Cell (Polynomial Rational Int)
cell | Cell (Polynomial Rational Int)
cell <- [Cell (Polynomial Rational Int)]
cells2]

  forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ AReal
xval AReal -> AReal -> AComplex
:+ AReal
yval

  where

    -- f(x + yi) = u(x,y) + v(x,y)i
    f1 :: P.Polynomial Rational Int
    f1 :: Polynomial Rational Int
f1 = forall k v2 v1.
(Eq k, Num k, Ord v2) =>
Polynomial k v1 -> (v1 -> Polynomial k v2) -> Polynomial k v2
P.subst UPolynomial Rational
f (\X
X -> forall a v. Var a v => v -> a
P.var Int
0 forall a. Num a => a -> a -> a
+ forall a v. Var a v => v -> a
P.var Int
1 forall a. Num a => a -> a -> a
* forall a v. Var a v => v -> a
P.var Int
2)
    f2 :: UPolynomial (Polynomial Rational Int)
f2 = forall k v.
(Ord k, Num k, Ord v) =>
Polynomial k v -> v -> UPolynomial (Polynomial k v)
P.toUPolynomialOf (forall k v.
(Eq k, Fractional k, Ord v) =>
MonomialOrder v
-> Polynomial k v -> [Polynomial k v] -> Polynomial k v
P.reduce forall v. Ord v => MonomialOrder v
P.grevlex Polynomial Rational Int
f1 [forall a v. Var a v => v -> a
P.var Int
2 forall a. Num a => a -> a -> a
* forall a v. Var a v => v -> a
P.var Int
2 forall a. Num a => a -> a -> a
+ Polynomial Rational Int
1]) Int
2
    u :: Polynomial Rational Int
u = forall k v. (Num k, Ord v) => Monomial v -> Polynomial k v -> k
P.coeff forall v. Monomial v
P.mone UPolynomial (Polynomial Rational Int)
f2
    v :: Polynomial Rational Int
v = forall k v. (Num k, Ord v) => Monomial v -> Polynomial k v -> k
P.coeff (forall a v. Var a v => v -> a
P.var X
X) UPolynomial (Polynomial Rational Int)
f2

-- -----------------------------------------------------------------------------

test1 :: UPolynomial Rational
test1 = AComplex -> UPolynomial Rational
minimalPolynomial (AReal
sqrt2 AReal -> AReal -> AComplex
:+ AReal
sqrt3)
  where
    sqrt2 :: AReal
sqrt2 = Integer -> AReal -> AReal
AReal.nthRoot Integer
2 AReal
2
    sqrt3 :: AReal
sqrt3 = Integer -> AReal -> AReal
AReal.nthRoot Integer
2 AReal
3

test2 :: AReal
test2 = AComplex -> AReal
magnitude (AReal
sqrt2 AReal -> AReal -> AComplex
:+ AReal
sqrt3)
  where
    sqrt2 :: AReal
sqrt2 = Integer -> AReal -> AReal
AReal.nthRoot Integer
2 AReal
2
    sqrt3 :: AReal
sqrt3 = Integer -> AReal -> AReal
AReal.nthRoot Integer
2 AReal
3

test3 :: [AComplex]
test3 = UPolynomial Rational -> [AComplex]
roots UPolynomial Rational
p
  where
    x :: UPolynomial Rational
x = forall a v. Var a v => v -> a
P.var X
X
    p :: UPolynomial Rational
p = UPolynomial Rational
x forall a. Num a => a -> a -> a
- UPolynomial Rational
2

test4 :: [AComplex]
test4 = UPolynomial Rational -> [AComplex]
roots UPolynomial Rational
p
  where
    x :: UPolynomial Rational
x = forall a v. Var a v => v -> a
P.var X
X
    p :: UPolynomial Rational
p = UPolynomial Rational
xforall a b. (Num a, Integral b) => a -> b -> a
^Integer
2 forall a. Num a => a -> a -> a
- UPolynomial Rational
2

test5 :: [AComplex]
test5 = UPolynomial Rational -> [AComplex]
roots UPolynomial Rational
p
  where
    x :: UPolynomial Rational
x = forall a v. Var a v => v -> a
P.var X
X
    p :: UPolynomial Rational
p = UPolynomial Rational
xforall a b. (Num a, Integral b) => a -> b -> a
^Integer
3 forall a. Num a => a -> a -> a
- UPolynomial Rational
1

test6 :: [AComplex]
test6 = UPolynomial Rational -> [AComplex]
roots UPolynomial Rational
p
  where
    x :: UPolynomial Rational
x = forall a v. Var a v => v -> a
P.var X
X
    p :: UPolynomial Rational
p = UPolynomial Rational
xforall a b. (Num a, Integral b) => a -> b -> a
^Integer
4 forall a. Num a => a -> a -> a
+ UPolynomial Rational
2forall a. Num a => a -> a -> a
*UPolynomial Rational
xforall a b. (Num a, Integral b) => a -> b -> a
^Integer
2 forall a. Num a => a -> a -> a
+ UPolynomial Rational
25