{-# OPTIONS_HADDOCK show-extensions #-}
{-# LANGUAGE Rank2Types #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.Data.AlgebraicNumber.Real
-- Copyright   :  (c) Masahiro Sakai 2012-2013
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  non-portable
--
-- Algebraic reals
--
-- Reference:
--
-- * Why the concept of a field extension is a natural one
--   <http://www.dpmms.cam.ac.uk/~wtg10/galois.html>
--
-----------------------------------------------------------------------------
module ToySolver.Data.AlgebraicNumber.Real
  (
  -- * Algebraic real type
    AReal

  -- * Construction
  , realRoots
  , realRootsEx

  -- * Properties
  , minimalPolynomial
  , isolatingInterval
  , isRational
  , isAlgebraicInteger
  , height
  , rootIndex

  -- * Operations
  , nthRoot
  , refineIsolatingInterval

  -- * Approximation
  , approx
  , approxInterval

  -- * Misc
  , 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

{--------------------------------------------------------------------
  Construction
--------------------------------------------------------------------}

-- | Algebraic real numbers.
data AReal = RealRoot (UPolynomial Rational) (Interval Rational)
  deriving Int -> AReal -> ShowS
[AReal] -> ShowS
AReal -> String
(Int -> AReal -> ShowS)
-> (AReal -> String) -> ([AReal] -> ShowS) -> Show AReal
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> AReal -> ShowS
showsPrec :: Int -> AReal -> ShowS
$cshow :: AReal -> String
show :: AReal -> String
$cshowList :: [AReal] -> ShowS
showList :: [AReal] -> ShowS
Show

-- | Real roots of the polynomial in ascending order.
realRoots :: UPolynomial Rational -> [AReal]
realRoots :: UPolynomial Rational -> [AReal]
realRoots UPolynomial Rational
p = Set AReal -> [AReal]
forall a. Set a -> [a]
Set.toAscList (Set AReal -> [AReal]) -> Set AReal -> [AReal]
forall a b. (a -> b) -> a -> b
$ [AReal] -> Set AReal
forall a. Ord a => [a] -> Set a
Set.fromList ([AReal] -> Set AReal) -> [AReal] -> Set AReal
forall a b. (a -> b) -> a -> b
$ do
  (UPolynomial Rational
q,Integer
_) <- UPolynomial Rational -> [(UPolynomial Rational, Integer)]
forall a. Factor a => a -> [(a, Integer)]
P.factor UPolynomial Rational
p
  UPolynomial Rational -> [AReal]
realRoots' UPolynomial Rational
q

-- | Real roots of the polynomial in ascending order.
realRootsEx :: UPolynomial AReal -> [AReal]
realRootsEx :: UPolynomial AReal -> [AReal]
realRootsEx UPolynomial AReal
p
  | [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and [AReal -> Bool
isRational AReal
c | (AReal
c,Monomial X
_) <- UPolynomial AReal -> [(AReal, Monomial X)]
forall k v. Polynomial k v -> [Term k v]
P.terms UPolynomial AReal
p] = UPolynomial Rational -> [AReal]
realRoots (UPolynomial Rational -> [AReal])
-> UPolynomial Rational -> [AReal]
forall a b. (a -> b) -> a -> b
$ (AReal -> Rational) -> UPolynomial AReal -> UPolynomial Rational
forall k1 k v.
(Eq k1, Num k1) =>
(k -> k1) -> Polynomial k v -> Polynomial k1 v
P.mapCoeff AReal -> Rational
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 AReal -> UPolynomial AReal -> Bool
forall k. (Eq k, Num k) => k -> UPolynomial k -> Bool
`P.isRootOf` UPolynomial AReal
p]

-- p must already be factored.
realRoots' :: UPolynomial Rational -> [AReal]
realRoots' :: UPolynomial Rational -> [AReal]
realRoots' UPolynomial Rational
p = do
  Bool -> [()]
forall (f :: * -> *). Alternative f => Bool -> f ()
guard (Bool -> [()]) -> Bool -> [()]
forall a b. (a -> b) -> a -> b
$ UPolynomial Rational -> Integer
forall t. Degree t => t -> Integer
P.deg UPolynomial Rational
p Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
0
  Interval Rational
i <- UPolynomial Rational -> [Interval Rational]
Sturm.separate UPolynomial Rational
p
  AReal -> [AReal]
forall a. a -> [a]
forall (m :: * -> *) a. Monad m => a -> m a
return (AReal -> [AReal]) -> AReal -> [AReal]
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
_) <- 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, UPolynomial Rational -> Interval Rational -> Int
Sturm.numRoots UPolynomial Rational
q Interval Rational
i Int -> Int -> Bool
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
    []   -> String -> AReal
forall a. HasCallStack => String -> a
error String
"ToySolver.Data.AlgebraicNumber.Real.realRoot: invalid interval"

-- p must already be factored.
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 (UPolynomial Rational -> UPolynomial Rational
forall k. (Fractional k, Eq k) => UPolynomial k -> UPolynomial k
normalizePoly UPolynomial Rational
p) Interval Rational
i

{--------------------------------------------------------------------
  Operations
--------------------------------------------------------------------}

isZero :: AReal -> Bool
isZero :: AReal -> Bool
isZero AReal
a = Rational
0 Rational -> Interval Rational -> Bool
forall r. Ord r => r -> Interval r -> Bool
`Interval.member` AReal -> Interval Rational
isolatingInterval AReal
a Bool -> Bool -> Bool
&& Rational
0 Rational -> UPolynomial Rational -> Bool
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 = Rational -> UPolynomial Rational -> UPolynomial Rational
forall k.
(Fractional k, Eq k) =>
k -> UPolynomial k -> UPolynomial k
rootScale Rational
r (AReal -> UPolynomial Rational
minimalPolynomial AReal
a)
    i2 :: Interval Rational
i2 = Rational -> Interval Rational
forall r. Ord r => r -> Interval r
Interval.singleton Rational
r Interval Rational -> Interval Rational -> Interval Rational
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 = Rational -> UPolynomial Rational -> UPolynomial Rational
forall k.
(Fractional k, Eq k) =>
k -> UPolynomial k -> UPolynomial k
rootShift Rational
r (AReal -> UPolynomial Rational
minimalPolynomial AReal
a)
    i2 :: Interval Rational
i2 = Rational -> Interval Rational
forall r. Ord r => r -> Interval r
Interval.singleton Rational
r Interval Rational -> Interval Rational -> Interval Rational
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
p1UPolynomial Rational -> UPolynomial Rational -> Bool
forall a. Eq a => a -> a -> Bool
==UPolynomial Rational
p2 Bool -> Bool -> Bool
&& [UPolynomial Rational] -> Interval Rational -> Int
Sturm.numRoots' [UPolynomial Rational]
c (Interval Rational -> Interval Rational -> Interval Rational
forall r. Ord r => Interval r -> Interval r -> Interval r
Interval.intersection Interval Rational
i1 Interval Rational
i2) Int -> Int -> Bool
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 Interval Rational -> Interval Rational -> Bool
forall r. Ord r => Interval r -> Interval r -> Bool
>! Interval Rational
i2  = Ordering
GT
    | Interval Rational
i1 Interval Rational -> Interval Rational -> Bool
forall r. Ord r => Interval r -> Interval r -> Bool
<! Interval Rational
i2  = Ordering
LT
    | AReal
a AReal -> AReal -> Bool
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 Interval Rational -> Interval Rational -> Bool
forall r. Ord r => Interval r -> Interval r -> Bool
>! Interval Rational
i2 = Ordering
GT
        | Interval Rational
i1 Interval Rational -> Interval Rational -> Bool
forall r. Ord r => Interval r -> Interval r -> Bool
<! Interval Rational
i2 = Ordering
LT
        | Bool
otherwise =
            if Interval Rational -> Rational
forall r. (Num r, Ord r) => Interval r -> r
Interval.width Interval Rational
i1 Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
> Interval Rational -> Rational
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 (AReal -> Rational
forall a. Real a => a -> Rational
toRational AReal
a) AReal
b
    | AReal -> Bool
isRational AReal
b = Rational -> AReal -> AReal
shiftAReal (AReal -> Rational
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 = UPolynomial Rational
-> UPolynomial Rational -> UPolynomial Rational
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 = Interval Rational -> Interval Rational -> Interval Rational
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 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0] of
          []   -> String -> Interval Rational
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 Interval Rational -> Interval Rational -> Interval Rational
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 (AReal -> Rational
forall a. Real a => a -> Rational
toRational AReal
a) AReal
b
    | AReal -> Bool
isRational AReal
b = Rational -> AReal -> AReal
scaleAReal (AReal -> Rational
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 = UPolynomial Rational
-> UPolynomial Rational -> UPolynomial Rational
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 = Interval Rational -> Interval Rational -> Interval Rational
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 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0] of
          []   -> String -> Interval Rational
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 Interval Rational -> Interval Rational -> Interval Rational
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 AReal -> AReal -> Ordering
forall a. Ord a => a -> a -> Ordering
compare AReal
0 AReal
a of
      Ordering
EQ -> Integer -> AReal
forall a. Num a => Integer -> a
fromInteger Integer
0
      Ordering
LT -> AReal
a
      Ordering
GT -> AReal -> AReal
forall a. Num a => a -> a
negate AReal
a

  signum :: AReal -> AReal
signum AReal
a = Integer -> AReal
forall a. Num a => Integer -> a
fromInteger (Integer -> AReal) -> Integer -> AReal
forall a b. (a -> b) -> a -> b
$
    case AReal -> AReal -> Ordering
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 = Rational -> AReal
forall a. Fractional a => Rational -> a
fromRational (Rational -> AReal) -> (Integer -> Rational) -> Integer -> AReal
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Rational
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 UPolynomial Rational
-> UPolynomial Rational -> UPolynomial Rational
forall a. Num a => a -> a -> a
- Rational -> UPolynomial Rational
forall k v. (Eq k, Num k) => k -> Polynomial k v
P.constant Rational
r) (Rational -> Interval Rational
forall r. Ord r => r -> Interval r
Interval.singleton Rational
r)
    where
      x :: UPolynomial Rational
x = X -> UPolynomial Rational
forall a v. Var a v => v -> a
P.var X
X

  recip :: AReal -> AReal
recip AReal
a
    | AReal -> Bool
isZero AReal
a  = String -> AReal
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 = UPolynomial Rational -> UPolynomial Rational
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, Rational -> Interval Rational -> Bool
forall r. Ord r => r -> Interval r -> Bool
Interval.member Rational
1 (Interval Rational
i1 Interval Rational -> Interval Rational -> Interval Rational
forall a. Num a => a -> a -> a
* Interval Rational
i2)] of
            [] -> String -> Interval Rational
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 = Monomial X -> UPolynomial Rational -> Rational
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 Rational
p
            b :: Rational
b = Monomial X -> UPolynomial Rational -> Rational
forall k v. (Num k, Ord v) => Monomial v -> Polynomial k v -> k
P.coeff Monomial X
forall v. Monomial v
P.mone UPolynomial Rational
p
        in - Rational
b Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
a
    | Bool
otherwise  = String -> Rational
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 = AReal -> (b, AReal)
forall b. Integral b => AReal -> (b, AReal)
properFraction'
  truncate :: forall b. Integral b => AReal -> b
truncate       = AReal -> b
forall b. Integral b => AReal -> b
truncate'
  round :: forall b. Integral b => AReal -> b
round          = AReal -> b
forall b. Integral b => AReal -> b
round'
  ceiling :: forall b. Integral b => AReal -> b
ceiling        = AReal -> b
forall b. Integral b => AReal -> b
ceiling'
  floor :: forall b. Integral b => AReal -> b
floor          = AReal -> b
forall b. Integral b => AReal -> b
floor'

-- | Returns approximate rational value such that @abs (a - approx a epsilon) <= epsilon@.
approx
  :: AReal    -- ^ a
  -> Rational -- ^ epsilon
  -> Rational
approx :: AReal -> Rational -> Rational
approx AReal
a Rational
epsilon =
  if AReal -> Bool
isRational AReal
a
    then AReal -> Rational
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

-- | Returns approximate interval such that @width (approxInterval a epsilon) <= epsilon@.
approxInterval
  :: AReal    -- ^ a
  -> Rational -- ^ epsilon
  -> Interval Rational
approxInterval :: AReal -> Rational -> Interval Rational
approxInterval AReal
a Rational
epsilon =
  if AReal -> Bool
isRational AReal
a
    then Rational -> Interval Rational
forall r. Ord r => r -> Interval r
Interval.singleton (AReal -> Rational
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

-- | Same as 'properFraction'.
properFraction' :: Integral b => AReal -> (b, AReal)
properFraction' :: forall b. Integral b => AReal -> (b, AReal)
properFraction' AReal
x =
  case AReal -> AReal -> Ordering
forall a. Ord a => a -> a -> Ordering
compare AReal
x AReal
0 of
    Ordering
EQ -> (b
0, AReal
0)
    Ordering
GT -> (Integer -> b
forall a. Num a => Integer -> a
fromInteger Integer
floor_x, AReal
x AReal -> AReal -> AReal
forall a. Num a => a -> a -> a
- Integer -> AReal
forall a. Num a => Integer -> a
fromInteger Integer
floor_x)
    Ordering
LT -> (Integer -> b
forall a. Num a => Integer -> a
fromInteger Integer
ceiling_x, AReal
x AReal -> AReal -> AReal
forall a. Num a => a -> a -> a
- Integer -> AReal
forall a. Num a => Integer -> a
fromInteger Integer
ceiling_x)
  where
    floor_x :: Integer
floor_x   = AReal -> Integer
forall b. Integral b => AReal -> b
floor' AReal
x
    ceiling_x :: Integer
ceiling_x = AReal -> Integer
forall b. Integral b => AReal -> b
ceiling' AReal
x

-- | Same as 'truncate'.
truncate' :: Integral b => AReal -> b
truncate' :: forall b. Integral b => AReal -> b
truncate' = (b, AReal) -> b
forall a b. (a, b) -> a
fst ((b, AReal) -> b) -> (AReal -> (b, AReal)) -> AReal -> b
forall b c a. (b -> c) -> (a -> b) -> a -> c
. AReal -> (b, AReal)
forall b. Integral b => AReal -> (b, AReal)
properFraction'

-- | Same as 'round'.
round' :: Integral b => AReal -> b
round' :: forall b. Integral b => AReal -> b
round' AReal
x =
  case AReal -> AReal
forall a. Num a => a -> a
signum (AReal -> AReal
forall a. Num a => a -> a
abs AReal
r AReal -> AReal -> AReal
forall a. Num a => a -> a -> a
- AReal
0.5) of
    -1 -> b
n
    AReal
0  -> if b -> Bool
forall a. Integral a => a -> Bool
even b
n then b
n else b
m
    AReal
1  -> b
m
    AReal
_  -> String -> b
forall a. HasCallStack => String -> a
error String
"round default defn: Bad value"
  where
    (b
n,AReal
r) = AReal -> (b, AReal)
forall b. Integral b => AReal -> (b, AReal)
properFraction' AReal
x
    m :: b
m = if AReal
r AReal -> AReal -> Bool
forall a. Ord a => a -> a -> Bool
< AReal
0 then b
n b -> b -> b
forall a. Num a => a -> a -> a
- b
1 else b
n b -> b -> b
forall a. Num a => a -> a -> a
+ b
1

-- | Same as 'ceiling'.
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 (Interval Rational -> Interval Rational -> Interval Rational
forall r. Ord r => Interval r -> Interval r -> Interval r
Interval.intersection Interval Rational
i2 Interval Rational
i3) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
1
    then Integer -> b
forall a. Num a => Integer -> a
fromInteger Integer
ceiling_lb
    else Integer -> b
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
1Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/Rational
2)
    (Finite Rational
lb, Boundary
inLB) = Interval Rational -> (Extended Rational, Boundary)
forall r. Interval r -> (Extended r, Boundary)
Interval.lowerBound' Interval Rational
i2
    (Finite Rational
ub, Boundary
inUB) = Interval Rational -> (Extended Rational, Boundary)
forall r. Interval r -> (Extended r, Boundary)
Interval.upperBound' Interval Rational
i2
    ceiling_lb :: Integer
ceiling_lb = Rational -> Integer
forall b. Integral b => Rational -> b
forall a b. (RealFrac a, Integral b) => a -> b
ceiling Rational
lb
    ceiling_ub :: Integer
ceiling_ub = Rational -> Integer
forall b. Integral b => Rational -> b
forall a b. (RealFrac a, Integral b) => a -> b
ceiling Rational
ub
    i3 :: Interval Rational
i3 = Extended Rational
forall r. Extended r
NegInf Extended Rational -> Extended Rational -> Interval Rational
forall r. Ord r => Extended r -> Extended r -> Interval r
<..<= Integer -> Extended Rational
forall a. Num a => Integer -> a
fromInteger Integer
ceiling_lb

-- | Same as 'floor'.
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 (Interval Rational -> Interval Rational -> Interval Rational
forall r. Ord r => Interval r -> Interval r -> Interval r
Interval.intersection Interval Rational
i2 Interval Rational
i3) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
1
    then Integer -> b
forall a. Num a => Integer -> a
fromInteger Integer
floor_ub
    else Integer -> b
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
1Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/Rational
2)
    (Finite Rational
lb, Boundary
inLB) = Interval Rational -> (Extended Rational, Boundary)
forall r. Interval r -> (Extended r, Boundary)
Interval.lowerBound' Interval Rational
i2
    (Finite Rational
ub, Boundary
inUB) = Interval Rational -> (Extended Rational, Boundary)
forall r. Interval r -> (Extended r, Boundary)
Interval.upperBound' Interval Rational
i2
    floor_lb :: Integer
floor_lb = Rational -> Integer
forall b. Integral b => Rational -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor Rational
lb
    floor_ub :: Integer
floor_ub = Rational -> Integer
forall b. Integral b => Rational -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor Rational
ub
    i3 :: Interval Rational
i3 = Integer -> Extended Rational
forall a. Num a => Integer -> a
fromInteger Integer
floor_ub Extended Rational -> Extended Rational -> Interval Rational
forall r. Ord r => Extended r -> Extended r -> Interval r
<=..< Extended Rational
forall r. Extended r
PosInf

-- | The @n@th root of @a@
nthRoot :: Integer -> AReal -> AReal
nthRoot :: Integer -> AReal -> AReal
nthRoot Integer
n AReal
a
  | Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
0 = String -> AReal
forall a. HasCallStack => String -> a
error String
"ToySolver.Data.AlgebraicNumver.Root.nthRoot"
  | Integer -> Bool
forall a. Integral a => a -> Bool
even Integer
n =
      if AReal
a AReal -> AReal -> Bool
forall a. Ord a => a -> a -> Bool
< AReal
0
      then String -> AReal
forall a. HasCallStack => String -> a
error String
"ToySolver.Data.AlgebraicNumver.Root.nthRoot: no real roots"
      else Bool -> AReal -> AReal
forall a. HasCallStack => Bool -> a -> a
assert ([AReal] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [AReal]
bs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
2) ([AReal] -> AReal
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [AReal]
bs) -- select positive root
  | Bool
otherwise =
      Bool -> AReal -> AReal
forall a. HasCallStack => Bool -> a -> a
assert ([AReal] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [AReal]
bs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1) ([AReal] -> AReal
forall a. HasCallStack => [a] -> a
head [AReal]
bs) -- select unique root
  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 Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
0 = []
nthRoots Integer
n AReal
a | Integer -> Bool
forall a. Integral a => a -> Bool
even Integer
n Bool -> Bool -> Bool
&& AReal
a AReal -> AReal -> Bool
forall a. Ord a => a -> a -> Bool
< AReal
0 = []
nthRoots Integer
n AReal
a = (AReal -> Bool) -> [AReal] -> [AReal]
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 = Integer -> UPolynomial Rational -> UPolynomial Rational
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 = (AReal -> Interval Rational) -> [AReal] -> [Interval Rational]
forall a b. (a -> b) -> [a] -> [b]
map AReal -> Interval Rational
isolatingInterval ([AReal] -> [Interval Rational]) -> [AReal] -> [Interval Rational]
forall a b. (a -> b) -> a -> b
$ AReal -> [AReal] -> [AReal]
forall a. Eq a => a -> [a] -> [a]
delete AReal
a ([AReal] -> [AReal]) -> [AReal] -> [AReal]
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' Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = Bool
False
          | [Interval Rational] -> Bool
forall a. [a] -> Bool
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')
                   ((Interval Rational -> Interval Rational)
-> [Interval Rational] -> [Interval Rational]
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 Interval Rational -> Integer -> Interval Rational
forall a b. (Num a, Integral b) => a -> b -> a
^ Integer
n
            ok' :: Interval Rational
ok' = Interval Rational -> Interval Rational -> Interval Rational
forall r. Ord r => Interval r -> Interval r -> Interval r
Interval.intersection Interval Rational
i2 Interval Rational
ok
            ng' :: [Interval Rational]
ng' = (Interval Rational -> Bool)
-> [Interval Rational] -> [Interval Rational]
forall a. (a -> Bool) -> [a] -> [a]
filter (\Interval Rational
i3 -> [UPolynomial Rational] -> Interval Rational -> Int
Sturm.numRoots' [UPolynomial Rational]
c1 Interval Rational
i3 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
0) ([Interval Rational] -> [Interval Rational])
-> [Interval Rational] -> [Interval Rational]
forall a b. (a -> b) -> a -> b
$
                    (Interval Rational -> Interval Rational)
-> [Interval Rational] -> [Interval Rational]
forall a b. (a -> b) -> [a] -> [b]
map (Interval Rational -> Interval Rational -> Interval Rational
forall r. Ord r => Interval r -> Interval r -> Interval r
Interval.intersection Interval Rational
i2) [Interval Rational]
ng

-- | Same algebraic real, but represented using finer grained 'isolatingInterval'.
refineIsolatingInterval :: AReal -> AReal
refineIsolatingInterval :: AReal -> AReal
refineIsolatingInterval a :: AReal
a@(RealRoot UPolynomial Rational
p Interval Rational
i) =
  if Interval Rational -> Rational
forall r. (Num r, Ord r) => Interval r -> r
Interval.width Interval Rational
i Rational -> Rational -> Bool
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)

{--------------------------------------------------------------------
  Properties
--------------------------------------------------------------------}

-- | The polynomial of which the algebraic number is root.
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)

-- | Isolating interval that separate the number from other roots of 'minimalPolynomial' of it.
isolatingInterval :: AReal -> Interval Rational
isolatingInterval :: AReal -> Interval Rational
isolatingInterval (RealRoot UPolynomial Rational
_ Interval Rational
i) = Interval Rational
i

-- | Degree of the algebraic number.
--
-- If the algebraic number's 'minimalPolynomial' has degree @n@,
-- then the algebraic number is said to be degree @n@.
instance P.Degree AReal where
  deg :: AReal -> Integer
deg AReal
a = UPolynomial Rational -> Integer
forall t. Degree t => t -> Integer
P.deg (UPolynomial Rational -> Integer)
-> UPolynomial Rational -> Integer
forall a b. (a -> b) -> a -> b
$ AReal -> UPolynomial Rational
minimalPolynomial AReal
a

-- | Whether the algebraic number is a rational.
isRational :: AReal -> Bool
isRational :: AReal -> Bool
isRational AReal
x = AReal -> Integer
forall t. Degree t => t -> Integer
P.deg AReal
x Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
1

-- | Whether the algebraic number is a root of a polynomial with integer
-- coefficients with leading coefficient @1@ (a monic polynomial).
isAlgebraicInteger :: AReal -> Bool
isAlgebraicInteger :: AReal -> Bool
isAlgebraicInteger AReal
x = Integer -> Integer
forall a. Num a => a -> a
abs (MonomialOrder X -> Polynomial Integer X -> Integer
forall k v. Num k => MonomialOrder v -> Polynomial k v -> k
P.lc MonomialOrder X
P.nat (UPolynomial Rational -> Polynomial (PPCoeff Rational) X
forall v.
Ord v =>
Polynomial Rational v -> Polynomial (PPCoeff Rational) v
forall k v.
(ContPP k, Ord v) =>
Polynomial k v -> Polynomial (PPCoeff k) v
P.pp (AReal -> UPolynomial Rational
minimalPolynomial AReal
x))) Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
1

-- | Height of the algebraic number.
--
-- The height of an algebraic number is the greatest absolute value of the
-- coefficients of the irreducible and primitive polynomial with integral
-- rational coefficients.
height :: AReal -> Integer
height :: AReal -> Integer
height AReal
x = [Integer] -> Integer
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Integer -> Integer
forall a. Num a => a -> a
abs Integer
c | (Integer
c,Monomial X
_) <- Polynomial Integer X -> [(Integer, Monomial X)]
forall k v. Polynomial k v -> [Term k v]
P.terms (Polynomial Integer X -> [(Integer, Monomial X)])
-> Polynomial Integer X -> [(Integer, Monomial X)]
forall a b. (a -> b) -> a -> b
$ UPolynomial Rational -> Polynomial (PPCoeff Rational) X
forall v.
Ord v =>
Polynomial Rational v -> Polynomial (PPCoeff Rational) v
forall k v.
(ContPP k, Ord v) =>
Polynomial k v -> Polynomial (PPCoeff k) v
P.pp (UPolynomial Rational -> Polynomial (PPCoeff Rational) X)
-> UPolynomial Rational -> Polynomial (PPCoeff Rational) X
forall a b. (a -> b) -> a -> b
$ AReal -> UPolynomial Rational
minimalPolynomial AReal
x]

-- | root index, satisfying
--
-- @
-- 'realRoots' ('minimalPolynomial' a) !! rootIndex a == a
-- @
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 = AReal -> [AReal] -> Maybe Int
forall a. Eq a => a -> [a] -> Maybe Int
elemIndex AReal
a [AReal]
as

{--------------------------------------------------------------------
  Pretty printing
--------------------------------------------------------------------}

instance Pretty AReal where
  pPrintPrec :: PrettyLevel -> Rational -> AReal -> Doc
pPrintPrec PrettyLevel
lv Rational
prec AReal
r =
    Bool -> Doc -> Doc
maybeParens (Rational
prec Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
> Rational
appPrec) (Doc -> Doc) -> Doc -> Doc
forall a b. (a -> b) -> a -> b
$
      [Doc] -> Doc
PP.hsep [String -> Doc
PP.text String
"RealRoot", PrettyLevel -> Rational -> UPolynomial Rational -> Doc
forall a. Pretty a => PrettyLevel -> Rational -> a -> Doc
pPrintPrec PrettyLevel
lv (Rational
appPrecRational -> Rational -> Rational
forall 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 = PrettyLevel -> Rational -> AReal -> Doc
forall a. Pretty a => PrettyLevel -> Rational -> a -> Doc
pPrintPrec
  isNegativeCoeff :: AReal -> Bool
isNegativeCoeff = (AReal
0AReal -> AReal -> Bool
forall a. Ord a => a -> a -> Bool
>)

{--------------------------------------------------------------------
  Manipulation of polynomials
--------------------------------------------------------------------}

-- 代数的数を係数とする多項式の根を、有理数係数多項式の根として表す
simpARealPoly :: UPolynomial AReal -> UPolynomial Rational
simpARealPoly :: UPolynomial AReal -> UPolynomial Rational
simpARealPoly UPolynomial AReal
p = (AReal -> UPolynomial Rational)
-> UPolynomial AReal -> UPolynomial Rational
forall k l.
(Fractional k, Ord k) =>
(l -> UPolynomial k) -> UPolynomial l -> UPolynomial k
rootSimpPoly AReal -> UPolynomial Rational
minimalPolynomial UPolynomial AReal
p

{--------------------------------------------------------------------
  Misc
--------------------------------------------------------------------}

-- | Golden ratio
goldenRatio :: AReal
goldenRatio :: AReal
goldenRatio = (AReal
1 AReal -> AReal -> AReal
forall a. Num a => a -> a -> a
+ AReal
root5) AReal -> AReal -> AReal
forall a. Fractional a => a -> a -> a
/ AReal
2
  where
    [AReal
_, AReal
root5] = [AReal] -> [AReal]
forall a. Ord a => [a] -> [a]
sort ([AReal] -> [AReal]) -> [AReal] -> [AReal]
forall a b. (a -> b) -> a -> b
$ UPolynomial Rational -> [AReal]
realRoots' ((X -> UPolynomial Rational
forall a v. Var a v => v -> a
P.var X
X)UPolynomial 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
5)