-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.Data.AlgebraicNumber.Sturm
-- Copyright   :  (c) Masahiro Sakai 2012
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  portable
--
-- Reference:
--
-- * \"/Sturm's theorem/.\" Wikipedia, The Free Encyclopedia. Wikimedia Foundation, Inc.
--   2012-06-23. <http://en.wikipedia.org/wiki/Sturm%27s_theorem>
--
-- * Weisstein, Eric W. \"/Sturm Function/.\" From MathWorld--A Wolfram Web Resource.
--   <http://mathworld.wolfram.com/SturmFunction.html>
--
-----------------------------------------------------------------------------

module ToySolver.Data.AlgebraicNumber.Sturm
  ( SturmChain
  , sturmChain
  , numRoots
  , numRoots'
  , separate
  , separate'
  , halve
  , halve'
  , narrow
  , narrow'
  , approx
  , approx'
  ) where

import Data.Maybe
import qualified Data.Interval as Interval
import Data.Interval (Interval, Extended (..), (<..<=), (<=..<=))

import ToySolver.Data.Polynomial (UPolynomial, X (..))
import qualified ToySolver.Data.Polynomial as P

-- | Sturm's chain (Sturm's sequence)
type SturmChain = [UPolynomial Rational]

-- | Sturm's sequence of a polynomial
sturmChain :: UPolynomial Rational -> SturmChain
sturmChain :: UPolynomial Rational -> SturmChain
sturmChain UPolynomial Rational
p = UPolynomial Rational
p0 forall a. a -> [a] -> [a]
: UPolynomial Rational
p1 forall a. a -> [a] -> [a]
: forall {k}.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> [UPolynomial k]
go UPolynomial Rational
p0 UPolynomial Rational
p1
  where
    p0 :: UPolynomial Rational
p0 = UPolynomial Rational
p
    p1 :: UPolynomial Rational
p1 = forall k v.
(Eq k, Num k, Ord v) =>
Polynomial k v -> v -> Polynomial k v
P.deriv UPolynomial Rational
p X
P.X
    go :: UPolynomial k -> UPolynomial k -> [UPolynomial k]
go UPolynomial k
p UPolynomial k
q = if UPolynomial k
rforall a. Eq a => a -> a -> Bool
==UPolynomial k
0 then [] else UPolynomial k
r forall a. a -> [a] -> [a]
: UPolynomial k -> UPolynomial k -> [UPolynomial k]
go UPolynomial k
q UPolynomial k
r
      where
        r :: UPolynomial k
r = - (UPolynomial k
p forall k.
(Eq k, Fractional k) =>
UPolynomial k -> UPolynomial k -> UPolynomial k
`P.mod` UPolynomial k
q)

-- | The number of distinct real roots of @p@ in a given interval
numRoots
  :: UPolynomial Rational
  -> Interval Rational
  -> Int
numRoots :: UPolynomial Rational -> Interval Rational -> Int
numRoots UPolynomial Rational
p Interval Rational
ival = SturmChain -> Interval Rational -> Int
numRoots' (UPolynomial Rational -> SturmChain
sturmChain UPolynomial Rational
p) Interval Rational
ival

-- | The number of distinct real roots of @p@ in a given interval.
-- This function takes @p@'s sturm chain instead of @p@ itself.
numRoots'
  :: SturmChain
  -> Interval Rational
  -> Int
numRoots' :: SturmChain -> Interval Rational -> Int
numRoots' chain :: SturmChain
chain@(UPolynomial Rational
p:SturmChain
_) Interval Rational
ival
  | forall r. Ord r => Interval r -> Bool
Interval.null Interval Rational
ival2 = Int
0
  | Bool
otherwise =
      case (forall r. Interval r -> Extended r
Interval.lowerBound Interval Rational
ival2, forall r. Interval r -> Extended r
Interval.upperBound Interval Rational
ival2) of
        (Finite Rational
lb, Finite Rational
ub) ->
          (if Rational
lbforall a. Eq a => a -> a -> Bool
==Rational
ub then Int
0 else (Rational -> Int
n Rational
lb forall a. Num a => a -> a -> a
- Rational -> Int
n Rational
ub)) forall a. Num a => a -> a -> a
+
          (if Rational
lb forall r. Ord r => r -> Interval r -> Bool
`Interval.member` Interval Rational
ival2 Bool -> Bool -> Bool
&& Rational
lb forall k. (Eq k, Num k) => k -> UPolynomial k -> Bool
`P.isRootOf` UPolynomial Rational
p then Int
1 else Int
0) forall a. Num a => a -> a -> a
+
          (if Rational
ub forall r. Ord r => r -> Interval r -> Bool
`Interval.notMember` Interval Rational
ival2 Bool -> Bool -> Bool
&& Rational
ub forall k. (Eq k, Num k) => k -> UPolynomial k -> Bool
`P.isRootOf`  UPolynomial Rational
p then -Int
1 else Int
0)
        (Extended Rational, Extended Rational)
_ -> forall a. HasCallStack => [Char] -> a
error [Char]
"numRoots'': should not happen"
  where
    ival2 :: Interval Rational
ival2 = UPolynomial Rational -> Interval Rational -> Interval Rational
boundInterval UPolynomial Rational
p Interval Rational
ival
    n :: Rational -> Int
n Rational
x = [Rational] -> Int
countSignChanges [forall k v. Num k => (v -> k) -> Polynomial k v -> k
P.eval (\X
X -> Rational
x) UPolynomial Rational
q | UPolynomial Rational
q <- SturmChain
chain]

countSignChanges :: [Rational] -> Int
countSignChanges :: [Rational] -> Int
countSignChanges [Rational]
rs = forall a. Eq a => [a] -> Int
countChanges [Bool]
xs
  where
    xs :: [Bool]
    xs :: [Bool]
xs = forall a b. (a -> b) -> [a] -> [b]
map (Rational
0forall a. Ord a => a -> a -> Bool
<) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. (a -> Bool) -> [a] -> [a]
filter (Rational
0forall a. Eq a => a -> a -> Bool
/=) forall a b. (a -> b) -> a -> b
$ [Rational]
rs

countChanges :: Eq a => [a] -> Int
countChanges :: forall a. Eq a => [a] -> Int
countChanges [] = Int
0
countChanges (a
x:[a]
xs) = forall {t} {t}. (Eq t, Num t) => t -> [t] -> t -> t
go a
x [a]
xs Int
0
  where
    go :: t -> [t] -> t -> t
go t
x [] t
r = t
r
    go t
x1 (t
x2:[t]
xs) t
r
      | t
x1forall a. Eq a => a -> a -> Bool
==t
x2    = t -> [t] -> t -> t
go t
x1 [t]
xs t
r
      | Bool
otherwise = t -> [t] -> t -> t
go t
x2 [t]
xs (t
rforall a. Num a => a -> a -> a
+t
1)

-- | Closed interval that contains all real roots of a given polynomial.
-- 根の限界
-- <http://aozoragakuen.sakura.ne.jp/taiwa/taiwaNch02/node26.html>
bounds :: UPolynomial Rational -> (Rational, Rational)
bounds :: UPolynomial Rational -> (Rational, Rational)
bounds UPolynomial Rational
p = (-Rational
m, Rational
m)
  where
    m :: Rational
m = if UPolynomial Rational
pforall a. Eq a => a -> a -> Bool
==UPolynomial Rational
0
        then Rational
0
        else forall a. Ord a => a -> a -> a
max Rational
1 (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [forall a. Num a => a -> a
abs (Rational
cforall a. Fractional a => a -> a -> a
/Rational
s) | (Rational
c,Monomial X
_) <- forall k v. Polynomial k v -> [Term k v]
P.terms UPolynomial Rational
p] forall a. Num a => a -> a -> a
- Rational
1)
    s :: Rational
s = forall k v. Num k => MonomialOrder v -> Polynomial k v -> k
P.lc MonomialOrder X
P.nat UPolynomial Rational
p

boundInterval :: UPolynomial Rational -> Interval Rational -> Interval Rational
boundInterval :: UPolynomial Rational -> Interval Rational -> Interval Rational
boundInterval UPolynomial Rational
p Interval Rational
ival = forall r. Ord r => Interval r -> Interval r -> Interval r
Interval.intersection Interval Rational
ival (forall r. r -> Extended r
Finite Rational
lb forall r. Ord r => Extended r -> Extended r -> Interval r
<=..<= forall r. r -> Extended r
Finite Rational
ub)
  where
    (Rational
lb,Rational
ub) = UPolynomial Rational -> (Rational, Rational)
bounds UPolynomial Rational
p

-- | Disjoint intervals each of which contains exactly one real roots of the given polynoimal @p@.
-- The intervals can be further narrowed by 'narrow' or 'narrow''.
separate :: UPolynomial Rational -> [Interval Rational]
separate :: UPolynomial Rational -> [Interval Rational]
separate UPolynomial Rational
p = SturmChain -> [Interval Rational]
separate' (UPolynomial Rational -> SturmChain
sturmChain UPolynomial Rational
p)

-- | Disjoint intervals each of which contains exactly one real roots of the given polynoimal @p@.
-- The intervals can be further narrowed by 'narrow' or 'narrow''.
-- This function takes @p@'s sturm chain instead of @p@ itself.
separate' :: SturmChain -> [Interval Rational]
separate' :: SturmChain -> [Interval Rational]
separate' chain :: SturmChain
chain@(UPolynomial Rational
p:SturmChain
_) = (Rational, Rational) -> [Interval Rational]
f (UPolynomial Rational -> (Rational, Rational)
bounds UPolynomial Rational
p)
  where
    n :: Rational -> Int
n Rational
x = [Rational] -> Int
countSignChanges [forall k v. Num k => (v -> k) -> Polynomial k v -> k
P.eval (\X
X -> Rational
x) UPolynomial Rational
q | UPolynomial Rational
q <- SturmChain
chain]

    f :: (Rational, Rational) -> [Interval Rational]
f (Rational
lb,Rational
ub) =
      if Rational
lb forall k. (Eq k, Num k) => k -> UPolynomial k -> Bool
`P.isRootOf` UPolynomial Rational
p
      then forall r. Ord r => r -> Interval r
Interval.singleton Rational
lb forall a. a -> [a] -> [a]
: (Rational, Rational) -> [Interval Rational]
g (Rational
lb,Rational
ub)
      else (Rational, Rational) -> [Interval Rational]
g (Rational
lb,Rational
ub)

    g :: (Rational, Rational) -> [Interval Rational]
g (Rational
lb,Rational
ub) =
      case Rational -> Int
n Rational
lb forall a. Num a => a -> a -> a
- Rational -> Int
n Rational
ub of
        Int
0 -> []
        Int
1 -> [forall r. r -> Extended r
Finite Rational
lb forall r. Ord r => Extended r -> Extended r -> Interval r
<..<= forall r. r -> Extended r
Finite Rational
ub]
        Int
_ -> (Rational, Rational) -> [Interval Rational]
g (Rational
lb, Rational
mid) forall a. [a] -> [a] -> [a]
++ (Rational, Rational) -> [Interval Rational]
g (Rational
mid, Rational
ub)
      where
        mid :: Rational
mid = (Rational
lb forall a. Num a => a -> a -> a
+ Rational
ub) forall a. Fractional a => a -> a -> a
/ Rational
2

halve :: UPolynomial Rational -> Interval Rational -> Interval Rational
halve :: UPolynomial Rational -> Interval Rational -> Interval Rational
halve UPolynomial Rational
p Interval Rational
ival = SturmChain -> Interval Rational -> Interval Rational
halve' (UPolynomial Rational -> SturmChain
sturmChain UPolynomial Rational
p) Interval Rational
ival

halve' :: SturmChain -> Interval Rational -> Interval Rational
halve' :: SturmChain -> Interval Rational -> Interval Rational
halve' chain :: SturmChain
chain@(UPolynomial Rational
p:SturmChain
_) Interval Rational
ival
  | forall r. (Num r, Ord r) => Interval r -> r
Interval.width Interval Rational
ival forall a. Eq a => a -> a -> Bool
== Rational
0  = Interval Rational
ival
  | SturmChain -> Interval Rational -> Int
numRoots' SturmChain
chain Interval Rational
ivalL forall a. Ord a => a -> a -> Bool
> Int
0 = Interval Rational
ivalL
  | Bool
otherwise                 = Interval Rational
ivalR -- numRoots' chain ivalR > 0
  where
    Finite Rational
lb = forall r. Interval r -> Extended r
Interval.lowerBound Interval Rational
ival
    Finite Rational
ub = forall r. Interval r -> Extended r
Interval.upperBound Interval Rational
ival
    mid :: Rational
mid = (Rational
lb forall a. Num a => a -> a -> a
+ Rational
ub) forall a. Fractional a => a -> a -> a
/ Rational
2
    ivalL :: Interval Rational
ivalL = forall r.
Ord r =>
(Extended r, Boundary) -> (Extended r, Boundary) -> Interval r
Interval.interval (forall r. Interval r -> (Extended r, Boundary)
Interval.lowerBound' Interval Rational
ival) (forall r. r -> Extended r
Finite Rational
mid, Boundary
Interval.Closed)
    ivalR :: Interval Rational
ivalR = forall r.
Ord r =>
(Extended r, Boundary) -> (Extended r, Boundary) -> Interval r
Interval.interval (forall r. r -> Extended r
Finite Rational
mid, Boundary
Interval.Open) (forall r. Interval r -> (Extended r, Boundary)
Interval.upperBound' Interval Rational
ival)

narrow :: UPolynomial Rational -> Interval Rational -> Rational -> Interval Rational
narrow :: UPolynomial Rational
-> Interval Rational -> Rational -> Interval Rational
narrow UPolynomial Rational
p Interval Rational
ival Rational
size = SturmChain -> Interval Rational -> Rational -> Interval Rational
narrow' (UPolynomial Rational -> SturmChain
sturmChain UPolynomial Rational
p) Interval Rational
ival Rational
size

narrow' :: SturmChain -> Interval Rational -> Rational -> Interval Rational
narrow' :: SturmChain -> Interval Rational -> Rational -> Interval Rational
narrow' chain :: SturmChain
chain@(UPolynomial Rational
p:SturmChain
_) Interval Rational
ival Rational
size = Interval Rational -> Interval Rational
go (UPolynomial Rational -> Interval Rational -> Interval Rational
boundInterval UPolynomial Rational
p Interval Rational
ival)
  where
    go :: Interval Rational -> Interval Rational
go Interval Rational
ival
      | forall r. (Num r, Ord r) => Interval r -> r
Interval.width Interval Rational
ival forall a. Ord a => a -> a -> Bool
< Rational
size  = Interval Rational
ival
      | Bool
otherwise = Interval Rational -> Interval Rational
go (SturmChain -> Interval Rational -> Interval Rational
halve' SturmChain
chain Interval Rational
ival)

approx :: UPolynomial Rational -> Interval Rational -> Rational -> Rational
approx :: UPolynomial Rational -> Interval Rational -> Rational -> Rational
approx UPolynomial Rational
p = SturmChain -> Interval Rational -> Rational -> Rational
approx' (UPolynomial Rational -> SturmChain
sturmChain UPolynomial Rational
p)

approx' :: SturmChain -> Interval Rational -> Rational -> Rational
approx' :: SturmChain -> Interval Rational -> Rational -> Rational
approx' SturmChain
chain Interval Rational
ival Rational
epsilon = forall a. HasCallStack => Maybe a -> a
fromJust forall a b. (a -> b) -> a -> b
$ forall r. (Real r, Fractional r) => Interval r -> Maybe r
Interval.pickup forall a b. (a -> b) -> a -> b
$ SturmChain -> Interval Rational -> Rational -> Interval Rational
narrow' SturmChain
chain Interval Rational
ival Rational
epsilon