{-# OPTIONS_HADDOCK show-extensions #-}
{-# LANGUAGE Rank2Types #-}
module ToySolver.Data.AlgebraicNumber.Real
(
AReal
, realRoots
, realRootsEx
, minimalPolynomial
, isolatingInterval
, isRational
, isAlgebraicInteger
, height
, rootIndex
, nthRoot
, refineIsolatingInterval
, approx
, approxInterval
, simpARealPoly
, goldenRatio
) where
import Control.Exception (assert)
import Control.Monad
import Data.List
import Data.Ratio
import qualified Data.Set as Set
import qualified Text.PrettyPrint.HughesPJClass as PP
import Text.PrettyPrint.HughesPJClass (Doc, PrettyLevel, Pretty (..), maybeParens)
import Data.Interval (Interval, Extended (..), (<=..<), (<..<=), (<..<), (<!), (>!))
import qualified Data.Interval as Interval
import ToySolver.Data.Polynomial (Polynomial, UPolynomial, X (..))
import qualified ToySolver.Data.Polynomial as P
import ToySolver.Data.AlgebraicNumber.Root
import qualified ToySolver.Data.AlgebraicNumber.Sturm as Sturm
data AReal = RealRoot (UPolynomial Rational) (Interval Rational)
deriving Int -> AReal -> ShowS
[AReal] -> ShowS
AReal -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [AReal] -> ShowS
$cshowList :: [AReal] -> ShowS
show :: AReal -> String
$cshow :: AReal -> String
showsPrec :: Int -> AReal -> ShowS
$cshowsPrec :: Int -> AReal -> ShowS
Show
realRoots :: UPolynomial Rational -> [AReal]
realRoots :: UPolynomial Rational -> [AReal]
realRoots UPolynomial Rational
p = forall a. Set a -> [a]
Set.toAscList forall a b. (a -> b) -> a -> b
$ forall a. Ord a => [a] -> Set a
Set.fromList forall a b. (a -> b) -> a -> b
$ do
(UPolynomial Rational
q,Integer
_) <- forall a. Factor a => a -> [(a, Integer)]
P.factor UPolynomial Rational
p
UPolynomial Rational -> [AReal]
realRoots' UPolynomial Rational
q
realRootsEx :: UPolynomial AReal -> [AReal]
realRootsEx :: UPolynomial AReal -> [AReal]
realRootsEx UPolynomial AReal
p
| forall (t :: * -> *). Foldable t => t Bool -> Bool
and [AReal -> Bool
isRational AReal
c | (AReal
c,Monomial X
_) <- forall k v. Polynomial k v -> [Term k v]
P.terms UPolynomial AReal
p] = UPolynomial Rational -> [AReal]
realRoots forall a b. (a -> b) -> a -> b
$ forall k1 k v.
(Eq k1, Num k1) =>
(k -> k1) -> Polynomial k v -> Polynomial k1 v
P.mapCoeff forall a. Real a => a -> Rational
toRational UPolynomial AReal
p
| Bool
otherwise = [AReal
a | AReal
a <- UPolynomial Rational -> [AReal]
realRoots (UPolynomial AReal -> UPolynomial Rational
simpARealPoly UPolynomial AReal
p), AReal
a forall k. (Eq k, Num k) => k -> UPolynomial k -> Bool
`P.isRootOf` UPolynomial AReal
p]
realRoots' :: UPolynomial Rational -> [AReal]
realRoots' :: UPolynomial Rational -> [AReal]
realRoots' UPolynomial Rational
p = do
forall (f :: * -> *). Alternative f => Bool -> f ()
guard forall a b. (a -> b) -> a -> b
$ forall t. Degree t => t -> Integer
P.deg UPolynomial Rational
p forall a. Ord a => a -> a -> Bool
> Integer
0
Interval Rational
i <- UPolynomial Rational -> [Interval Rational]
Sturm.separate UPolynomial Rational
p
forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ UPolynomial Rational -> Interval Rational -> AReal
realRoot' UPolynomial Rational
p Interval Rational
i
realRoot :: UPolynomial Rational -> Interval Rational -> AReal
realRoot :: UPolynomial Rational -> Interval Rational -> AReal
realRoot UPolynomial Rational
p Interval Rational
i =
case [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, UPolynomial Rational -> Interval Rational -> Int
Sturm.numRoots UPolynomial Rational
q Interval Rational
i forall a. Eq a => a -> a -> Bool
== Int
1] of
UPolynomial Rational
p2:[UPolynomial Rational]
_ -> UPolynomial Rational -> Interval Rational -> AReal
realRoot' UPolynomial Rational
p2 Interval Rational
i
[] -> forall a. HasCallStack => String -> a
error String
"ToySolver.Data.AlgebraicNumber.Real.realRoot: invalid interval"
realRoot' :: UPolynomial Rational -> Interval Rational -> AReal
realRoot' :: UPolynomial Rational -> Interval Rational -> AReal
realRoot' UPolynomial Rational
p Interval Rational
i = UPolynomial Rational -> Interval Rational -> AReal
RealRoot (forall k. (Fractional k, Eq k) => UPolynomial k -> UPolynomial k
normalizePoly UPolynomial Rational
p) Interval Rational
i
isZero :: AReal -> Bool
isZero :: AReal -> Bool
isZero AReal
a = Rational
0 forall r. Ord r => r -> Interval r -> Bool
`Interval.member` AReal -> Interval Rational
isolatingInterval AReal
a Bool -> Bool -> Bool
&& Rational
0 forall k. (Eq k, Num k) => k -> UPolynomial k -> Bool
`P.isRootOf` AReal -> UPolynomial Rational
minimalPolynomial AReal
a
scaleAReal :: Rational -> AReal -> AReal
scaleAReal :: Rational -> AReal -> AReal
scaleAReal Rational
r AReal
a = UPolynomial Rational -> Interval Rational -> AReal
realRoot' UPolynomial Rational
p2 Interval Rational
i2
where
p2 :: UPolynomial Rational
p2 = forall k.
(Fractional k, Eq k) =>
k -> UPolynomial k -> UPolynomial k
rootScale Rational
r (AReal -> UPolynomial Rational
minimalPolynomial AReal
a)
i2 :: Interval Rational
i2 = forall r. Ord r => r -> Interval r
Interval.singleton Rational
r forall a. Num a => a -> a -> a
* AReal -> Interval Rational
isolatingInterval AReal
a
shiftAReal :: Rational -> AReal -> AReal
shiftAReal :: Rational -> AReal -> AReal
shiftAReal Rational
r AReal
a = UPolynomial Rational -> Interval Rational -> AReal
realRoot' UPolynomial Rational
p2 Interval Rational
i2
where
p2 :: UPolynomial Rational
p2 = forall k.
(Fractional k, Eq k) =>
k -> UPolynomial k -> UPolynomial k
rootShift Rational
r (AReal -> UPolynomial Rational
minimalPolynomial AReal
a)
i2 :: Interval Rational
i2 = forall r. Ord r => r -> Interval r
Interval.singleton Rational
r forall a. Num a => a -> a -> a
+ AReal -> Interval Rational
isolatingInterval AReal
a
instance Eq AReal where
AReal
a == :: AReal -> AReal -> Bool
== AReal
b = UPolynomial Rational
p1forall a. Eq a => a -> a -> Bool
==UPolynomial Rational
p2 Bool -> Bool -> Bool
&& [UPolynomial Rational] -> Interval Rational -> Int
Sturm.numRoots' [UPolynomial Rational]
c (forall r. Ord r => Interval r -> Interval r -> Interval r
Interval.intersection Interval Rational
i1 Interval Rational
i2) forall a. Eq a => a -> a -> Bool
== Int
1
where
p1 :: UPolynomial Rational
p1 = AReal -> UPolynomial Rational
minimalPolynomial AReal
a
p2 :: UPolynomial Rational
p2 = AReal -> UPolynomial Rational
minimalPolynomial AReal
b
i1 :: Interval Rational
i1 = AReal -> Interval Rational
isolatingInterval AReal
a
i2 :: Interval Rational
i2 = AReal -> Interval Rational
isolatingInterval AReal
b
c :: [UPolynomial Rational]
c = AReal -> [UPolynomial Rational]
sturmChain AReal
a
instance Ord AReal where
compare :: AReal -> AReal -> Ordering
compare AReal
a AReal
b
| Interval Rational
i1 forall r. Ord r => Interval r -> Interval r -> Bool
>! Interval Rational
i2 = Ordering
GT
| Interval Rational
i1 forall r. Ord r => Interval r -> Interval r -> Bool
<! Interval Rational
i2 = Ordering
LT
| AReal
a forall a. Eq a => a -> a -> Bool
== AReal
b = Ordering
EQ
| Bool
otherwise = Interval Rational -> Interval Rational -> Ordering
go Interval Rational
i1 Interval Rational
i2
where
i1 :: Interval Rational
i1 = AReal -> Interval Rational
isolatingInterval AReal
a
i2 :: Interval Rational
i2 = AReal -> Interval Rational
isolatingInterval AReal
b
c1 :: [UPolynomial Rational]
c1 = AReal -> [UPolynomial Rational]
sturmChain AReal
a
c2 :: [UPolynomial Rational]
c2 = AReal -> [UPolynomial Rational]
sturmChain AReal
b
go :: Interval Rational -> Interval Rational -> Ordering
go Interval Rational
i1 Interval Rational
i2
| Interval Rational
i1 forall r. Ord r => Interval r -> Interval r -> Bool
>! Interval Rational
i2 = Ordering
GT
| Interval Rational
i1 forall r. Ord r => Interval r -> Interval r -> Bool
<! Interval Rational
i2 = Ordering
LT
| Bool
otherwise =
if forall r. (Num r, Ord r) => Interval r -> r
Interval.width Interval Rational
i1 forall a. Ord a => a -> a -> Bool
> forall r. (Num r, Ord r) => Interval r -> r
Interval.width Interval Rational
i2
then Interval Rational -> Interval Rational -> Ordering
go ([UPolynomial Rational] -> Interval Rational -> Interval Rational
Sturm.halve' [UPolynomial Rational]
c1 Interval Rational
i1) Interval Rational
i2
else Interval Rational -> Interval Rational -> Ordering
go Interval Rational
i1 ([UPolynomial Rational] -> Interval Rational -> Interval Rational
Sturm.halve' [UPolynomial Rational]
c2 Interval Rational
i2)
instance Num AReal where
AReal
a + :: AReal -> AReal -> AReal
+ AReal
b
| AReal -> Bool
isRational AReal
a = Rational -> AReal -> AReal
shiftAReal (forall a. Real a => a -> Rational
toRational AReal
a) AReal
b
| AReal -> Bool
isRational AReal
b = Rational -> AReal -> AReal
shiftAReal (forall a. Real a => a -> Rational
toRational AReal
b) AReal
a
| Bool
otherwise = UPolynomial Rational -> Interval Rational -> AReal
realRoot UPolynomial Rational
p3 Interval Rational
i3
where
p3 :: UPolynomial Rational
p3 = forall k.
(Fractional k, Ord k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
rootAdd (AReal -> UPolynomial Rational
minimalPolynomial AReal
a) (AReal -> UPolynomial Rational
minimalPolynomial AReal
b)
c1 :: [UPolynomial Rational]
c1 = AReal -> [UPolynomial Rational]
sturmChain AReal
a
c2 :: [UPolynomial Rational]
c2 = AReal -> [UPolynomial Rational]
sturmChain AReal
b
c3 :: [UPolynomial Rational]
c3 = UPolynomial Rational -> [UPolynomial Rational]
Sturm.sturmChain UPolynomial Rational
p3
i3 :: Interval Rational
i3 = Interval Rational
-> Interval Rational -> [Interval Rational] -> Interval Rational
go (AReal -> Interval Rational
isolatingInterval AReal
a) (AReal -> Interval Rational
isolatingInterval AReal
b) ([UPolynomial Rational] -> [Interval Rational]
Sturm.separate' [UPolynomial Rational]
c3)
go :: Interval Rational
-> Interval Rational -> [Interval Rational] -> Interval Rational
go Interval Rational
i1 Interval Rational
i2 [Interval Rational]
is3 =
case [Interval Rational
i5 | Interval Rational
i3 <- [Interval Rational]
is3, let i5 :: Interval Rational
i5 = forall r. Ord r => Interval r -> Interval r -> Interval r
Interval.intersection Interval Rational
i3 Interval Rational
i4, [UPolynomial Rational] -> Interval Rational -> Int
Sturm.numRoots' [UPolynomial Rational]
c3 Interval Rational
i5 forall a. Ord a => a -> a -> Bool
> Int
0] of
[] -> forall a. HasCallStack => String -> a
error String
"AReal.+: should not happen"
[Interval Rational
i5] -> Interval Rational
i5
[Interval Rational]
is5 -> Interval Rational
-> Interval Rational -> [Interval Rational] -> Interval Rational
go ([UPolynomial Rational] -> Interval Rational -> Interval Rational
Sturm.halve' [UPolynomial Rational]
c1 Interval Rational
i1) ([UPolynomial Rational] -> Interval Rational -> Interval Rational
Sturm.halve' [UPolynomial Rational]
c2 Interval Rational
i2) [[UPolynomial Rational] -> Interval Rational -> Interval Rational
Sturm.halve' [UPolynomial Rational]
c3 Interval Rational
i5 | Interval Rational
i5 <- [Interval Rational]
is5]
where
i4 :: Interval Rational
i4 = Interval Rational
i1 forall a. Num a => a -> a -> a
+ Interval Rational
i2
AReal
a * :: AReal -> AReal -> AReal
* AReal
b
| AReal -> Bool
isRational AReal
a = Rational -> AReal -> AReal
scaleAReal (forall a. Real a => a -> Rational
toRational AReal
a) AReal
b
| AReal -> Bool
isRational AReal
b = Rational -> AReal -> AReal
scaleAReal (forall a. Real a => a -> Rational
toRational AReal
b) AReal
a
| Bool
otherwise = UPolynomial Rational -> Interval Rational -> AReal
realRoot UPolynomial Rational
p3 Interval Rational
i3
where
p3 :: UPolynomial Rational
p3 = forall k.
(Fractional k, Ord k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
rootMul (AReal -> UPolynomial Rational
minimalPolynomial AReal
a) (AReal -> UPolynomial Rational
minimalPolynomial AReal
b)
c1 :: [UPolynomial Rational]
c1 = AReal -> [UPolynomial Rational]
sturmChain AReal
a
c2 :: [UPolynomial Rational]
c2 = AReal -> [UPolynomial Rational]
sturmChain AReal
b
c3 :: [UPolynomial Rational]
c3 = UPolynomial Rational -> [UPolynomial Rational]
Sturm.sturmChain UPolynomial Rational
p3
i3 :: Interval Rational
i3 = Interval Rational
-> Interval Rational -> [Interval Rational] -> Interval Rational
go (AReal -> Interval Rational
isolatingInterval AReal
a) (AReal -> Interval Rational
isolatingInterval AReal
b) ([UPolynomial Rational] -> [Interval Rational]
Sturm.separate' [UPolynomial Rational]
c3)
go :: Interval Rational
-> Interval Rational -> [Interval Rational] -> Interval Rational
go Interval Rational
i1 Interval Rational
i2 [Interval Rational]
is3 =
case [Interval Rational
i5 | Interval Rational
i3 <- [Interval Rational]
is3, let i5 :: Interval Rational
i5 = forall r. Ord r => Interval r -> Interval r -> Interval r
Interval.intersection Interval Rational
i3 Interval Rational
i4, [UPolynomial Rational] -> Interval Rational -> Int
Sturm.numRoots' [UPolynomial Rational]
c3 Interval Rational
i5 forall a. Ord a => a -> a -> Bool
> Int
0] of
[] -> forall a. HasCallStack => String -> a
error String
"AReal.*: should not happen"
[Interval Rational
i5] -> Interval Rational
i5
[Interval Rational]
is5 -> Interval Rational
-> Interval Rational -> [Interval Rational] -> Interval Rational
go ([UPolynomial Rational] -> Interval Rational -> Interval Rational
Sturm.halve' [UPolynomial Rational]
c1 Interval Rational
i1) ([UPolynomial Rational] -> Interval Rational -> Interval Rational
Sturm.halve' [UPolynomial Rational]
c2 Interval Rational
i2) [[UPolynomial Rational] -> Interval Rational -> Interval Rational
Sturm.halve' [UPolynomial Rational]
c3 Interval Rational
i5 | Interval Rational
i5 <- [Interval Rational]
is5]
where
i4 :: Interval Rational
i4 = Interval Rational
i1 forall a. Num a => a -> a -> a
* Interval Rational
i2
negate :: AReal -> AReal
negate AReal
a = Rational -> AReal -> AReal
scaleAReal (-Rational
1) AReal
a
abs :: AReal -> AReal
abs AReal
a =
case forall a. Ord a => a -> a -> Ordering
compare AReal
0 AReal
a of
Ordering
EQ -> forall a. Num a => Integer -> a
fromInteger Integer
0
Ordering
LT -> AReal
a
Ordering
GT -> forall a. Num a => a -> a
negate AReal
a
signum :: AReal -> AReal
signum AReal
a = forall a. Num a => Integer -> a
fromInteger forall a b. (a -> b) -> a -> b
$
case forall a. Ord a => a -> a -> Ordering
compare AReal
0 AReal
a of
Ordering
EQ -> Integer
0
Ordering
LT -> Integer
1
Ordering
GT -> -Integer
1
fromInteger :: Integer -> AReal
fromInteger = forall a. Fractional a => Rational -> a
fromRational forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Real a => a -> Rational
toRational
instance Fractional AReal where
fromRational :: Rational -> AReal
fromRational Rational
r = UPolynomial Rational -> Interval Rational -> AReal
realRoot' (UPolynomial Rational
x forall a. Num a => a -> a -> a
- forall k v. (Eq k, Num k) => k -> Polynomial k v
P.constant Rational
r) (forall r. Ord r => r -> Interval r
Interval.singleton Rational
r)
where
x :: UPolynomial Rational
x = forall a v. Var a v => v -> a
P.var X
X
recip :: AReal -> AReal
recip AReal
a
| AReal -> Bool
isZero AReal
a = forall a. HasCallStack => String -> a
error String
"AReal.recip: zero division"
| Bool
otherwise = UPolynomial Rational -> Interval Rational -> AReal
realRoot' UPolynomial Rational
p2 Interval Rational
i2
where
p2 :: UPolynomial Rational
p2 = forall k. (Fractional k, Eq k) => UPolynomial k -> UPolynomial k
rootRecip (AReal -> UPolynomial Rational
minimalPolynomial AReal
a)
c1 :: [UPolynomial Rational]
c1 = AReal -> [UPolynomial Rational]
sturmChain AReal
a
c2 :: [UPolynomial Rational]
c2 = UPolynomial Rational -> [UPolynomial Rational]
Sturm.sturmChain UPolynomial Rational
p2
i2 :: Interval Rational
i2 = Interval Rational -> [Interval Rational] -> Interval Rational
go (AReal -> Interval Rational
isolatingInterval AReal
a) ([UPolynomial Rational] -> [Interval Rational]
Sturm.separate' [UPolynomial Rational]
c2)
go :: Interval Rational -> [Interval Rational] -> Interval Rational
go Interval Rational
i1 [Interval Rational]
is2 =
case [Interval Rational
i2 | Interval Rational
i2 <- [Interval Rational]
is2, forall r. Ord r => r -> Interval r -> Bool
Interval.member Rational
1 (Interval Rational
i1 forall a. Num a => a -> a -> a
* Interval Rational
i2)] of
[] -> forall a. HasCallStack => String -> a
error String
"AReal.recip: should not happen"
[Interval Rational
i2] -> Interval Rational
i2
[Interval Rational]
is2' -> Interval Rational -> [Interval Rational] -> Interval Rational
go ([UPolynomial Rational] -> Interval Rational -> Interval Rational
Sturm.halve' [UPolynomial Rational]
c1 Interval Rational
i1) [[UPolynomial Rational] -> Interval Rational -> Interval Rational
Sturm.halve' [UPolynomial Rational]
c2 Interval Rational
i2 | Interval Rational
i2 <- [Interval Rational]
is2']
instance Real AReal where
toRational :: AReal -> Rational
toRational AReal
x
| AReal -> Bool
isRational AReal
x =
let p :: UPolynomial Rational
p = AReal -> UPolynomial Rational
minimalPolynomial AReal
x
a :: Rational
a = 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 Rational
p
b :: Rational
b = forall k v. (Num k, Ord v) => Monomial v -> Polynomial k v -> k
P.coeff forall v. Monomial v
P.mone UPolynomial Rational
p
in - Rational
b forall a. Fractional a => a -> a -> a
/ Rational
a
| Bool
otherwise = forall a. HasCallStack => String -> a
error String
"AReal.toRational: proper algebraic number"
instance RealFrac AReal where
properFraction :: forall b. Integral b => AReal -> (b, AReal)
properFraction = forall b. Integral b => AReal -> (b, AReal)
properFraction'
truncate :: forall b. Integral b => AReal -> b
truncate = forall b. Integral b => AReal -> b
truncate'
round :: forall b. Integral b => AReal -> b
round = forall b. Integral b => AReal -> b
round'
ceiling :: forall b. Integral b => AReal -> b
ceiling = forall b. Integral b => AReal -> b
ceiling'
floor :: forall b. Integral b => AReal -> b
floor = forall b. Integral b => AReal -> b
floor'
approx
:: AReal
-> Rational
-> Rational
approx :: AReal -> Rational -> Rational
approx AReal
a Rational
epsilon =
if AReal -> Bool
isRational AReal
a
then forall a. Real a => a -> Rational
toRational AReal
a
else [UPolynomial Rational] -> Interval Rational -> Rational -> Rational
Sturm.approx' (AReal -> [UPolynomial Rational]
sturmChain AReal
a) (AReal -> Interval Rational
isolatingInterval AReal
a) Rational
epsilon
approxInterval
:: AReal
-> Rational
-> Interval Rational
approxInterval :: AReal -> Rational -> Interval Rational
approxInterval AReal
a Rational
epsilon =
if AReal -> Bool
isRational AReal
a
then forall r. Ord r => r -> Interval r
Interval.singleton (forall a. Real a => a -> Rational
toRational AReal
a)
else [UPolynomial Rational]
-> Interval Rational -> Rational -> Interval Rational
Sturm.narrow' (AReal -> [UPolynomial Rational]
sturmChain AReal
a) (AReal -> Interval Rational
isolatingInterval AReal
a) Rational
epsilon
properFraction' :: Integral b => AReal -> (b, AReal)
properFraction' :: forall b. Integral b => AReal -> (b, AReal)
properFraction' AReal
x =
case forall a. Ord a => a -> a -> Ordering
compare AReal
x AReal
0 of
Ordering
EQ -> (b
0, AReal
0)
Ordering
GT -> (forall a. Num a => Integer -> a
fromInteger Integer
floor_x, AReal
x forall a. Num a => a -> a -> a
- forall a. Num a => Integer -> a
fromInteger Integer
floor_x)
Ordering
LT -> (forall a. Num a => Integer -> a
fromInteger Integer
ceiling_x, AReal
x forall a. Num a => a -> a -> a
- forall a. Num a => Integer -> a
fromInteger Integer
ceiling_x)
where
floor_x :: Integer
floor_x = forall b. Integral b => AReal -> b
floor' AReal
x
ceiling_x :: Integer
ceiling_x = forall b. Integral b => AReal -> b
ceiling' AReal
x
truncate' :: Integral b => AReal -> b
truncate' :: forall b. Integral b => AReal -> b
truncate' = forall a b. (a, b) -> a
fst forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall b. Integral b => AReal -> (b, AReal)
properFraction'
round' :: Integral b => AReal -> b
round' :: forall b. Integral b => AReal -> b
round' AReal
x =
case forall a. Num a => a -> a
signum (forall a. Num a => a -> a
abs AReal
r forall a. Num a => a -> a -> a
- AReal
0.5) of
-1 -> b
n
AReal
0 -> if forall a. Integral a => a -> Bool
even b
n then b
n else b
m
AReal
1 -> b
m
AReal
_ -> forall a. HasCallStack => String -> a
error String
"round default defn: Bad value"
where
(b
n,AReal
r) = forall b. Integral b => AReal -> (b, AReal)
properFraction' AReal
x
m :: b
m = if AReal
r forall a. Ord a => a -> a -> Bool
< AReal
0 then b
n forall a. Num a => a -> a -> a
- b
1 else b
n forall a. Num a => a -> a -> a
+ b
1
ceiling' :: Integral b => AReal -> b
ceiling' :: forall b. Integral b => AReal -> b
ceiling' AReal
a =
if [UPolynomial Rational] -> Interval Rational -> Int
Sturm.numRoots' [UPolynomial Rational]
chain (forall r. Ord r => Interval r -> Interval r -> Interval r
Interval.intersection Interval Rational
i2 Interval Rational
i3) forall a. Ord a => a -> a -> Bool
>= Int
1
then forall a. Num a => Integer -> a
fromInteger Integer
ceiling_lb
else forall a. Num a => Integer -> a
fromInteger Integer
ceiling_ub
where
chain :: [UPolynomial Rational]
chain = AReal -> [UPolynomial Rational]
sturmChain AReal
a
i2 :: Interval Rational
i2 = [UPolynomial Rational]
-> Interval Rational -> Rational -> Interval Rational
Sturm.narrow' [UPolynomial Rational]
chain (AReal -> Interval Rational
isolatingInterval AReal
a) (Rational
1forall a. Fractional a => a -> a -> a
/Rational
2)
(Finite Rational
lb, Boundary
inLB) = forall r. Interval r -> (Extended r, Boundary)
Interval.lowerBound' Interval Rational
i2
(Finite Rational
ub, Boundary
inUB) = forall r. Interval r -> (Extended r, Boundary)
Interval.upperBound' Interval Rational
i2
ceiling_lb :: Integer
ceiling_lb = forall a b. (RealFrac a, Integral b) => a -> b
ceiling Rational
lb
ceiling_ub :: Integer
ceiling_ub = forall a b. (RealFrac a, Integral b) => a -> b
ceiling Rational
ub
i3 :: Interval Rational
i3 = forall r. Extended r
NegInf forall r. Ord r => Extended r -> Extended r -> Interval r
<..<= forall a. Num a => Integer -> a
fromInteger Integer
ceiling_lb
floor' :: Integral b => AReal -> b
floor' :: forall b. Integral b => AReal -> b
floor' AReal
a =
if [UPolynomial Rational] -> Interval Rational -> Int
Sturm.numRoots' [UPolynomial Rational]
chain (forall r. Ord r => Interval r -> Interval r -> Interval r
Interval.intersection Interval Rational
i2 Interval Rational
i3) forall a. Ord a => a -> a -> Bool
>= Int
1
then forall a. Num a => Integer -> a
fromInteger Integer
floor_ub
else forall a. Num a => Integer -> a
fromInteger Integer
floor_lb
where
chain :: [UPolynomial Rational]
chain = AReal -> [UPolynomial Rational]
sturmChain AReal
a
i2 :: Interval Rational
i2 = [UPolynomial Rational]
-> Interval Rational -> Rational -> Interval Rational
Sturm.narrow' [UPolynomial Rational]
chain (AReal -> Interval Rational
isolatingInterval AReal
a) (Rational
1forall a. Fractional a => a -> a -> a
/Rational
2)
(Finite Rational
lb, Boundary
inLB) = forall r. Interval r -> (Extended r, Boundary)
Interval.lowerBound' Interval Rational
i2
(Finite Rational
ub, Boundary
inUB) = forall r. Interval r -> (Extended r, Boundary)
Interval.upperBound' Interval Rational
i2
floor_lb :: Integer
floor_lb = forall a b. (RealFrac a, Integral b) => a -> b
floor Rational
lb
floor_ub :: Integer
floor_ub = forall a b. (RealFrac a, Integral b) => a -> b
floor Rational
ub
i3 :: Interval Rational
i3 = forall a. Num a => Integer -> a
fromInteger Integer
floor_ub forall r. Ord r => Extended r -> Extended r -> Interval r
<=..< forall r. Extended r
PosInf
nthRoot :: Integer -> AReal -> AReal
nthRoot :: Integer -> AReal -> AReal
nthRoot Integer
n AReal
a
| Integer
n forall a. Ord a => a -> a -> Bool
<= Integer
0 = forall a. HasCallStack => String -> a
error String
"ToySolver.Data.AlgebraicNumver.Root.nthRoot"
| forall a. Integral a => a -> Bool
even Integer
n =
if AReal
a forall a. Ord a => a -> a -> Bool
< AReal
0
then forall a. HasCallStack => String -> a
error String
"ToySolver.Data.AlgebraicNumver.Root.nthRoot: no real roots"
else forall a. HasCallStack => Bool -> a -> a
assert (forall (t :: * -> *) a. Foldable t => t a -> Int
length [AReal]
bs forall a. Eq a => a -> a -> Bool
== Int
2) (forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [AReal]
bs)
| Bool
otherwise =
forall a. HasCallStack => Bool -> a -> a
assert (forall (t :: * -> *) a. Foldable t => t a -> Int
length [AReal]
bs forall a. Eq a => a -> a -> Bool
== Int
1) (forall a. [a] -> a
head [AReal]
bs)
where
bs :: [AReal]
bs = Integer -> AReal -> [AReal]
nthRoots Integer
n AReal
a
nthRoots :: Integer -> AReal -> [AReal]
nthRoots :: Integer -> AReal -> [AReal]
nthRoots Integer
n AReal
_ | Integer
n forall a. Ord a => a -> a -> Bool
<= Integer
0 = []
nthRoots Integer
n AReal
a | forall a. Integral a => a -> Bool
even Integer
n Bool -> Bool -> Bool
&& AReal
a forall a. Ord a => a -> a -> Bool
< AReal
0 = []
nthRoots Integer
n AReal
a = forall a. (a -> Bool) -> [a] -> [a]
filter AReal -> Bool
check (UPolynomial Rational -> [AReal]
realRoots UPolynomial Rational
p2)
where
p1 :: UPolynomial Rational
p1 = AReal -> UPolynomial Rational
minimalPolynomial AReal
a
p2 :: UPolynomial Rational
p2 = forall k.
(Fractional k, Ord k) =>
Integer -> UPolynomial k -> UPolynomial k
rootNthRoot Integer
n UPolynomial Rational
p1
c1 :: [UPolynomial Rational]
c1 = AReal -> [UPolynomial Rational]
sturmChain AReal
a
ok0 :: Interval Rational
ok0 = AReal -> Interval Rational
isolatingInterval AReal
a
ng0 :: [Interval Rational]
ng0 = forall a b. (a -> b) -> [a] -> [b]
map AReal -> Interval Rational
isolatingInterval forall a b. (a -> b) -> a -> b
$ forall a. Eq a => a -> [a] -> [a]
delete AReal
a forall a b. (a -> b) -> a -> b
$ UPolynomial Rational -> [AReal]
realRoots UPolynomial Rational
p1
check :: AReal -> Bool
check :: AReal -> Bool
check AReal
b = Interval Rational
-> [Interval Rational] -> Interval Rational -> Bool
loop Interval Rational
ok0 [Interval Rational]
ng0 (AReal -> Interval Rational
isolatingInterval AReal
b)
where
c2 :: [UPolynomial Rational]
c2 = AReal -> [UPolynomial Rational]
sturmChain AReal
b
loop :: Interval Rational
-> [Interval Rational] -> Interval Rational -> Bool
loop Interval Rational
ok [Interval Rational]
ng Interval Rational
i
| [UPolynomial Rational] -> Interval Rational -> Int
Sturm.numRoots' [UPolynomial Rational]
c1 Interval Rational
ok' forall a. Eq a => a -> a -> Bool
== Int
0 = Bool
False
| forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Interval Rational]
ng' = Bool
True
| Bool
otherwise =
Interval Rational
-> [Interval Rational] -> Interval Rational -> Bool
loop ([UPolynomial Rational] -> Interval Rational -> Interval Rational
Sturm.halve' [UPolynomial Rational]
c1 Interval Rational
ok')
(forall a b. (a -> b) -> [a] -> [b]
map (\Interval Rational
i3 -> [UPolynomial Rational] -> Interval Rational -> Interval Rational
Sturm.halve' [UPolynomial Rational]
c1 Interval Rational
i3) [Interval Rational]
ng')
([UPolynomial Rational] -> Interval Rational -> Interval Rational
Sturm.halve' [UPolynomial Rational]
c2 Interval Rational
i)
where
i2 :: Interval Rational
i2 = Interval Rational
i forall a b. (Num a, Integral b) => a -> b -> a
^ Integer
n
ok' :: Interval Rational
ok' = forall r. Ord r => Interval r -> Interval r -> Interval r
Interval.intersection Interval Rational
i2 Interval Rational
ok
ng' :: [Interval Rational]
ng' = forall a. (a -> Bool) -> [a] -> [a]
filter (\Interval Rational
i3 -> [UPolynomial Rational] -> Interval Rational -> Int
Sturm.numRoots' [UPolynomial Rational]
c1 Interval Rational
i3 forall a. Eq a => a -> a -> Bool
/= Int
0) forall a b. (a -> b) -> a -> b
$
forall a b. (a -> b) -> [a] -> [b]
map (forall r. Ord r => Interval r -> Interval r -> Interval r
Interval.intersection Interval Rational
i2) [Interval Rational]
ng
refineIsolatingInterval :: AReal -> AReal
refineIsolatingInterval :: AReal -> AReal
refineIsolatingInterval a :: AReal
a@(RealRoot UPolynomial Rational
p Interval Rational
i) =
if forall r. (Num r, Ord r) => Interval r -> r
Interval.width Interval Rational
i forall a. Ord a => a -> a -> Bool
<= Rational
0
then AReal
a
else UPolynomial Rational -> Interval Rational -> AReal
RealRoot UPolynomial Rational
p (UPolynomial Rational -> Interval Rational -> Interval Rational
Sturm.halve UPolynomial Rational
p Interval Rational
i)
minimalPolynomial :: AReal -> UPolynomial Rational
minimalPolynomial :: AReal -> UPolynomial Rational
minimalPolynomial (RealRoot UPolynomial Rational
p Interval Rational
_) = UPolynomial Rational
p
sturmChain :: AReal -> Sturm.SturmChain
sturmChain :: AReal -> [UPolynomial Rational]
sturmChain AReal
a = UPolynomial Rational -> [UPolynomial Rational]
Sturm.sturmChain (AReal -> UPolynomial Rational
minimalPolynomial AReal
a)
isolatingInterval :: AReal -> Interval Rational
isolatingInterval :: AReal -> Interval Rational
isolatingInterval (RealRoot UPolynomial Rational
_ Interval Rational
i) = Interval Rational
i
instance P.Degree AReal where
deg :: AReal -> Integer
deg AReal
a = forall t. Degree t => t -> Integer
P.deg forall a b. (a -> b) -> a -> b
$ AReal -> UPolynomial Rational
minimalPolynomial AReal
a
isRational :: AReal -> Bool
isRational :: AReal -> Bool
isRational AReal
x = forall t. Degree t => t -> Integer
P.deg AReal
x forall a. Eq a => a -> a -> Bool
== Integer
1
isAlgebraicInteger :: AReal -> Bool
isAlgebraicInteger :: AReal -> Bool
isAlgebraicInteger AReal
x = forall a. Num a => a -> a
abs (forall k v. Num k => MonomialOrder v -> Polynomial k v -> k
P.lc MonomialOrder X
P.nat (forall k v.
(ContPP k, Ord v) =>
Polynomial k v -> Polynomial (PPCoeff k) v
P.pp (AReal -> UPolynomial Rational
minimalPolynomial AReal
x))) forall a. Eq a => a -> a -> Bool
== Integer
1
height :: AReal -> Integer
height :: AReal -> Integer
height AReal
x = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [forall a. Num a => a -> a
abs Integer
c | (Integer
c,Monomial X
_) <- forall k v. Polynomial k v -> [Term k v]
P.terms forall a b. (a -> b) -> a -> b
$ forall k v.
(ContPP k, Ord v) =>
Polynomial k v -> Polynomial (PPCoeff k) v
P.pp forall a b. (a -> b) -> a -> b
$ AReal -> UPolynomial Rational
minimalPolynomial AReal
x]
rootIndex :: AReal -> Int
rootIndex :: AReal -> Int
rootIndex AReal
a = Int
idx
where
as :: [AReal]
as = UPolynomial Rational -> [AReal]
realRoots' (AReal -> UPolynomial Rational
minimalPolynomial AReal
a)
Just Int
idx = forall a. Eq a => a -> [a] -> Maybe Int
elemIndex AReal
a [AReal]
as
instance Pretty AReal where
pPrintPrec :: PrettyLevel -> Rational -> AReal -> Doc
pPrintPrec PrettyLevel
lv Rational
prec AReal
r =
Bool -> Doc -> Doc
maybeParens (Rational
prec forall a. Ord a => a -> a -> Bool
> Rational
appPrec) forall a b. (a -> b) -> a -> b
$
[Doc] -> Doc
PP.hsep [String -> Doc
PP.text String
"RealRoot", forall a. Pretty a => PrettyLevel -> Rational -> a -> Doc
pPrintPrec PrettyLevel
lv (Rational
appPrecforall a. Num a => a -> a -> a
+Rational
1) UPolynomial Rational
p, Int -> Doc
PP.int (AReal -> Int
rootIndex AReal
r)]
where
p :: UPolynomial Rational
p = AReal -> UPolynomial Rational
minimalPolynomial AReal
r
appPrec :: Rational
appPrec = Rational
10
instance P.PrettyCoeff AReal where
pPrintCoeff :: PrettyLevel -> Rational -> AReal -> Doc
pPrintCoeff = forall a. Pretty a => PrettyLevel -> Rational -> a -> Doc
pPrintPrec
isNegativeCoeff :: AReal -> Bool
isNegativeCoeff = (AReal
0forall a. Ord a => a -> a -> Bool
>)
simpARealPoly :: UPolynomial AReal -> UPolynomial Rational
simpARealPoly :: UPolynomial AReal -> UPolynomial Rational
simpARealPoly UPolynomial AReal
p = forall k l.
(Fractional k, Ord k) =>
(l -> UPolynomial k) -> UPolynomial l -> UPolynomial k
rootSimpPoly AReal -> UPolynomial Rational
minimalPolynomial UPolynomial AReal
p
goldenRatio :: AReal
goldenRatio :: AReal
goldenRatio = (AReal
1 forall a. Num a => a -> a -> a
+ AReal
root5) forall a. Fractional a => a -> a -> a
/ AReal
2
where
[AReal
_, AReal
root5] = forall a. Ord a => [a] -> [a]
sort forall a b. (a -> b) -> a -> b
$ UPolynomial Rational -> [AReal]
realRoots' ((forall a v. Var a v => v -> a
P.var X
X)forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2 forall a. Num a => a -> a -> a
- UPolynomial Rational
5)