module ToySolver.Data.AlgebraicNumber.Complex
(
AComplex((:+))
, realPart
, imagPart
, magnitude
, minimalPolynomial
, conjugate
) 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 :+
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)
realPart :: AComplex -> AReal
realPart :: AComplex -> AReal
realPart (AReal
x :+ AReal
_) = AReal
x
imagPart :: AComplex -> AReal
imagPart :: AComplex -> AReal
imagPart (AReal
_ :+ AReal
y) = AReal
y
conjugate :: AComplex -> AComplex
conjugate :: AComplex -> AComplex
conjugate (AReal
x :+ AReal
y) = AReal
x AReal -> AReal -> AComplex
:+ (-AReal
y)
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)
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
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 :: 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
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