-----------------------------------------------------------------------------
-- |
-- 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
(AComplex -> AComplex -> Bool)
-> (AComplex -> AComplex -> Bool) -> Eq AComplex
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
(Int -> AComplex -> ShowS)
-> (AComplex -> String) -> ([AComplex] -> ShowS) -> Show AComplex
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
xAReal -> AReal -> AReal
forall a. Num a => a -> a -> a
*AReal
x AReal -> AReal -> AReal
forall a. Num a => a -> a -> a
+ AReal
yAReal -> AReal -> AReal
forall 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
xAReal -> AReal -> AReal
forall a. Num a => a -> a -> a
+AReal
x') AReal -> AReal -> AComplex
:+ (AReal
yAReal -> AReal -> AReal
forall a. Num a => a -> a -> a
+AReal
y')
    (AReal
x:+AReal
y) - :: AComplex -> AComplex -> AComplex
- (AReal
x':+AReal
y')   =  (AReal
xAReal -> AReal -> AReal
forall a. Num a => a -> a -> a
-AReal
x') AReal -> AReal -> AComplex
:+ (AReal
yAReal -> AReal -> AReal
forall a. Num a => a -> a -> a
-AReal
y')
    (AReal
x:+AReal
y) * :: AComplex -> AComplex -> AComplex
* (AReal
x':+AReal
y')   =  (AReal
xAReal -> AReal -> AReal
forall a. Num a => a -> a -> a
*AReal
x'AReal -> AReal -> AReal
forall a. Num a => a -> a -> a
-AReal
yAReal -> AReal -> AReal
forall a. Num a => a -> a -> a
*AReal
y') AReal -> AReal -> AComplex
:+ (AReal
xAReal -> AReal -> AReal
forall a. Num a => a -> a -> a
*AReal
y'AReal -> AReal -> AReal
forall a. Num a => a -> a -> a
+AReal
yAReal -> AReal -> AReal
forall a. Num a => a -> a -> a
*AReal
x')
    negate :: AComplex -> AComplex
negate (AReal
x:+AReal
y)       =  AReal -> AReal
forall a. Num a => a -> a
negate AReal
x AReal -> AReal -> AComplex
:+ AReal -> AReal
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
xAReal -> AReal -> AReal
forall a. Fractional a => a -> a -> a
/AReal
r AReal -> AReal -> AComplex
:+ AReal
yAReal -> AReal -> AReal
forall a. Fractional a => a -> a -> a
/AReal
r  where r :: AReal
r = AComplex -> AReal
magnitude AComplex
z
    fromInteger :: Integer -> AComplex
fromInteger Integer
n       =  Integer -> AReal
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
xAReal -> AReal -> AReal
forall a. Num a => a -> a -> a
*AReal
x'AReal -> AReal -> AReal
forall a. Num a => a -> a -> a
+AReal
yAReal -> AReal -> AReal
forall a. Num a => a -> a -> a
*AReal
y') AReal -> AReal -> AReal
forall a. Fractional a => a -> a -> a
/ AReal
d AReal -> AReal -> AComplex
:+ (AReal
yAReal -> AReal -> AReal
forall a. Num a => a -> a -> a
*AReal
x'AReal -> AReal -> AReal
forall a. Num a => a -> a -> a
-AReal
xAReal -> AReal -> AReal
forall a. Num a => a -> a -> a
*AReal
y') AReal -> AReal -> AReal
forall a. Fractional a => a -> a -> a
/ AReal
d
                           where d :: AReal
d   = AReal
x'AReal -> AReal -> AReal
forall a. Num a => a -> a -> a
*AReal
x' AReal -> AReal -> AReal
forall a. Num a => a -> a -> a
+ AReal
y'AReal -> AReal -> AReal
forall a. Num a => a -> a -> a
*AReal
y'
    fromRational :: Rational -> AComplex
fromRational Rational
a      =  Rational -> AReal
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) =
  [UPolynomial Rational] -> UPolynomial Rational
forall a. [a] -> a
head [UPolynomial Rational
q | (UPolynomial Rational
q,Integer
_) <- UPolynomial Rational -> [(UPolynomial Rational, Integer)]
forall a. Factor a => a -> [(a, Integer)]
P.factor UPolynomial Rational
p, UPolynomial Rational -> Integer
forall t. Degree t => t -> Integer
P.deg UPolynomial Rational
q Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
0, (X -> AComplex) -> Polynomial AComplex X -> AComplex
forall k v. Num k => (v -> k) -> Polynomial k v -> k
P.eval (\X
X -> AComplex
z) ((Rational -> AComplex)
-> UPolynomial Rational -> Polynomial AComplex X
forall k1 k v.
(Eq k1, Num k1) =>
(k -> k1) -> Polynomial k v -> Polynomial k1 v
P.mapCoeff Rational -> AComplex
forall a. Fractional a => Rational -> a
fromRational UPolynomial Rational
q) AComplex -> AComplex -> Bool
forall a. Eq a => a -> a -> Bool
== AComplex
0]
  where
    p :: UPolynomial Rational
p = Polynomial Rational Int
-> [Polynomial Rational Int] -> UPolynomial Rational
forall k.
(Fractional k, Ord k) =>
Polynomial k Int -> [Polynomial k Int] -> UPolynomial k
Root.findPoly (Int -> Polynomial Rational Int
forall a v. Var a v => v -> a
P.var Int
a Polynomial Rational Int
-> Polynomial Rational Int -> Polynomial Rational Int
forall a. Num a => a -> a -> a
+ Int -> Polynomial Rational Int
forall a v. Var a v => v -> a
P.var Int
b Polynomial Rational Int
-> Polynomial Rational Int -> Polynomial Rational Int
forall a. Num a => a -> a -> a
* Int -> Polynomial Rational Int
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 =
          [ (Int -> Polynomial Rational Int
forall a v. Var a v => v -> a
P.var Int
i) Polynomial Rational Int -> Integer -> Polynomial Rational Int
forall a b. (Num a, Integral b) => a -> b -> a
^ Integer
2 Polynomial Rational Int
-> Polynomial Rational Int -> Polynomial Rational Int
forall a. Num a => a -> a -> a
+ Polynomial Rational Int
1
          , UPolynomial Rational
-> (X -> Polynomial Rational Int) -> Polynomial Rational Int
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 -> Int -> Polynomial Rational Int
forall a v. Var a v => v -> a
P.var Int
a)
          , UPolynomial Rational
-> (X -> Polynomial Rational Int) -> Polynomial Rational Int
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 -> Int -> Polynomial Rational Int
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) <- [(UPolynomial (Polynomial Rational Int), [Sign])]
-> [([(Polynomial Rational Int, [Sign])],
     [Cell (Polynomial Rational Int)])]
forall v.
(Ord v, Show v, PrettyVar v) =>
[(UPolynomial (Polynomial Rational v), [Sign])]
-> [([(Polynomial Rational v, [Sign])],
     [Cell (Polynomial Rational v)])]
CAD.project' [(Polynomial Rational Int
-> Int -> UPolynomial (Polynomial Rational Int)
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) <- [(UPolynomial (Polynomial Rational Int), [Sign])]
-> [([(Polynomial Rational Int, [Sign])],
     [Cell (Polynomial Rational Int)])]
forall v.
(Ord v, Show v, PrettyVar v) =>
[(UPolynomial (Polynomial Rational v), [Sign])]
-> [([(Polynomial Rational v, [Sign])],
     [Cell (Polynomial Rational v)])]
CAD.project' [(Polynomial Rational Int
-> Int -> UPolynomial (Polynomial Rational Int)
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]
  Bool -> [()]
forall (f :: * -> *). Alternative f => Bool -> f ()
guard (Bool -> [()]) -> Bool -> [()]
forall a b. (a -> b) -> a -> b
$ [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and [Rational -> Sign
forall a. Real a => a -> Sign
Sign.signOf Rational
v Sign -> [Sign] -> Bool
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 = (Int -> Rational) -> Polynomial Rational Int -> Rational
forall k v. Num k => (v -> k) -> Polynomial k v -> k
P.eval (\Int
_ -> Rational
forall a. HasCallStack => a
undefined) Polynomial Rational Int
p]

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

  AComplex -> [AComplex]
forall (m :: * -> *) a. Monad m => a -> m a
return (AComplex -> [AComplex]) -> AComplex -> [AComplex]
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 = UPolynomial Rational
-> (X -> Polynomial Rational Int) -> Polynomial Rational Int
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 -> Int -> Polynomial Rational Int
forall a v. Var a v => v -> a
P.var Int
0 Polynomial Rational Int
-> Polynomial Rational Int -> Polynomial Rational Int
forall a. Num a => a -> a -> a
+ Int -> Polynomial Rational Int
forall a v. Var a v => v -> a
P.var Int
1 Polynomial Rational Int
-> Polynomial Rational Int -> Polynomial Rational Int
forall a. Num a => a -> a -> a
* Int -> Polynomial Rational Int
forall a v. Var a v => v -> a
P.var Int
2)
    f2 :: UPolynomial (Polynomial Rational Int)
f2 = Polynomial Rational Int
-> Int -> UPolynomial (Polynomial Rational Int)
forall k v.
(Ord k, Num k, Ord v) =>
Polynomial k v -> v -> UPolynomial (Polynomial k v)
P.toUPolynomialOf (MonomialOrder Int
-> Polynomial Rational Int
-> [Polynomial Rational Int]
-> Polynomial Rational Int
forall k v.
(Eq k, Fractional k, Ord v) =>
MonomialOrder v
-> Polynomial k v -> [Polynomial k v] -> Polynomial k v
P.reduce MonomialOrder Int
forall v. Ord v => MonomialOrder v
P.grevlex Polynomial Rational Int
f1 [Int -> Polynomial Rational Int
forall a v. Var a v => v -> a
P.var Int
2 Polynomial Rational Int
-> Polynomial Rational Int -> Polynomial Rational Int
forall a. Num a => a -> a -> a
* Int -> Polynomial Rational Int
forall a v. Var a v => v -> a
P.var Int
2 Polynomial Rational Int
-> Polynomial Rational Int -> Polynomial Rational Int
forall a. Num a => a -> a -> a
+ Polynomial Rational Int
1]) Int
2
    u :: Polynomial Rational Int
u = Monomial X
-> UPolynomial (Polynomial Rational Int) -> Polynomial Rational Int
forall k v. (Num k, Ord v) => Monomial v -> Polynomial k v -> k
P.coeff Monomial X
forall v. Monomial v
P.mone UPolynomial (Polynomial Rational Int)
f2
    v :: Polynomial Rational Int
v = Monomial X
-> UPolynomial (Polynomial Rational Int) -> Polynomial Rational Int
forall k v. (Num k, Ord v) => Monomial v -> Polynomial k v -> k
P.coeff (X -> Monomial X
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 = X -> UPolynomial Rational
forall a v. Var a v => v -> a
P.var X
X
    p :: UPolynomial Rational
p = UPolynomial Rational
x UPolynomial Rational
-> UPolynomial Rational -> UPolynomial Rational
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 = X -> UPolynomial Rational
forall a v. Var a v => v -> a
P.var X
X
    p :: UPolynomial Rational
p = UPolynomial Rational
xUPolynomial Rational -> Integer -> UPolynomial Rational
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2 UPolynomial Rational
-> UPolynomial Rational -> UPolynomial Rational
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 = X -> UPolynomial Rational
forall a v. Var a v => v -> a
P.var X
X
    p :: UPolynomial Rational
p = UPolynomial Rational
xUPolynomial Rational -> Integer -> UPolynomial Rational
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
3 UPolynomial Rational
-> UPolynomial Rational -> UPolynomial Rational
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 = X -> UPolynomial Rational
forall a v. Var a v => v -> a
P.var X
X
    p :: UPolynomial Rational
p = UPolynomial Rational
xUPolynomial Rational -> Integer -> UPolynomial Rational
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
4 UPolynomial Rational
-> UPolynomial Rational -> UPolynomial Rational
forall a. Num a => a -> a -> a
+ UPolynomial Rational
2UPolynomial Rational
-> UPolynomial Rational -> UPolynomial Rational
forall a. Num a => a -> a -> a
*UPolynomial Rational
xUPolynomial Rational -> Integer -> UPolynomial Rational
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2 UPolynomial Rational
-> UPolynomial Rational -> UPolynomial Rational
forall a. Num a => a -> a -> a
+ UPolynomial Rational
25