{-# 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
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

-- | Real roots of the polynomial in ascending order.
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

-- | Real roots of the polynomial in ascending order.
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]

-- p must already be factored.
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"

-- 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 (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 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'

-- | 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 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 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

-- | Same as 'properFraction'.
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

-- | Same as 'truncate'.
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'

-- | Same as 'round'.
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

-- | 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 (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

-- | 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 (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

-- | The @n@th root of @a@
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) -- select positive root
  | 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) -- 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 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

-- | 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 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)

{--------------------------------------------------------------------
  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 = forall t. Degree t => t -> Integer
P.deg 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 = forall t. Degree t => t -> Integer
P.deg AReal
x 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 = 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 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 = 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]

-- | 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 = 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 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
>)

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

-- 代数的数を係数とする多項式の根を、有理数係数多項式の根として表す
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

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

-- | Golden ratio
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)