-----------------------------------------------------------------------------
-- |
-- 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 UPolynomial Rational -> SturmChain -> SturmChain
forall a. a -> [a] -> [a]
: UPolynomial Rational
p1 UPolynomial Rational -> SturmChain -> SturmChain
forall a. a -> [a] -> [a]
: UPolynomial Rational -> UPolynomial Rational -> SturmChain
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 = UPolynomial Rational -> X -> UPolynomial Rational
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
rUPolynomial k -> UPolynomial k -> Bool
forall a. Eq a => a -> a -> Bool
==UPolynomial k
0 then [] else UPolynomial k
r UPolynomial k -> [UPolynomial k] -> [UPolynomial k]
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 UPolynomial k -> UPolynomial k -> UPolynomial k
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
  | Interval Rational -> Bool
forall r. Ord r => Interval r -> Bool
Interval.null Interval Rational
ival2 = Int
0
  | Bool
otherwise =
      case (Interval Rational -> Extended Rational
forall r. Interval r -> Extended r
Interval.lowerBound Interval Rational
ival2, Interval Rational -> Extended Rational
forall r. Interval r -> Extended r
Interval.upperBound Interval Rational
ival2) of
        (Finite Rational
lb, Finite Rational
ub) ->
          (if Rational
lbRational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
==Rational
ub then Int
0 else (Rational -> Int
n Rational
lb Int -> Int -> Int
forall a. Num a => a -> a -> a
- Rational -> Int
n Rational
ub)) Int -> Int -> Int
forall a. Num a => a -> a -> a
+
          (if Rational
lb Rational -> Interval Rational -> Bool
forall r. Ord r => r -> Interval r -> Bool
`Interval.member` Interval Rational
ival2 Bool -> Bool -> Bool
&& Rational
lb Rational -> UPolynomial Rational -> Bool
forall k. (Eq k, Num k) => k -> UPolynomial k -> Bool
`P.isRootOf` UPolynomial Rational
p then Int
1 else Int
0) Int -> Int -> Int
forall a. Num a => a -> a -> a
+
          (if Rational
ub Rational -> Interval Rational -> Bool
forall r. Ord r => r -> Interval r -> Bool
`Interval.notMember` Interval Rational
ival2 Bool -> Bool -> Bool
&& Rational
ub Rational -> UPolynomial Rational -> Bool
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)
_ -> [Char] -> Int
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 [(X -> Rational) -> UPolynomial Rational -> Rational
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 = [Bool] -> Int
forall a. Eq a => [a] -> Int
countChanges [Bool]
xs
  where
    xs :: [Bool]
    xs :: [Bool]
xs = (Rational -> Bool) -> [Rational] -> [Bool]
forall a b. (a -> b) -> [a] -> [b]
map (Rational
0Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
<) ([Rational] -> [Bool])
-> ([Rational] -> [Rational]) -> [Rational] -> [Bool]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Rational -> Bool) -> [Rational] -> [Rational]
forall a. (a -> Bool) -> [a] -> [a]
filter (Rational
0Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
/=) ([Rational] -> [Bool]) -> [Rational] -> [Bool]
forall a b. (a -> b) -> a -> b
$ [Rational]
rs

countChanges :: Eq a => [a] -> Int
countChanges :: [a] -> Int
countChanges [] = Int
0
countChanges (a
x:[a]
xs) = a -> [a] -> Int -> Int
forall t a. (Eq t, Num a) => t -> [t] -> a -> a
go a
x [a]
xs Int
0
  where
    go :: t -> [t] -> a -> a
go t
x [] a
r = a
r
    go t
x1 (t
x2:[t]
xs) a
r
      | t
x1t -> t -> Bool
forall a. Eq a => a -> a -> Bool
==t
x2    = t -> [t] -> a -> a
go t
x1 [t]
xs a
r
      | Bool
otherwise = t -> [t] -> a -> a
go t
x2 [t]
xs (a
ra -> a -> a
forall a. Num a => a -> a -> a
+a
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
pUPolynomial Rational -> UPolynomial Rational -> Bool
forall a. Eq a => a -> a -> Bool
==UPolynomial Rational
0
        then Rational
0
        else Rational -> Rational -> Rational
forall a. Ord a => a -> a -> a
max Rational
1 ([Rational] -> Rational
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Rational -> Rational
forall a. Num a => a -> a
abs (Rational
cRational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/Rational
s) | (Rational
c,Monomial X
_) <- UPolynomial Rational -> [(Rational, Monomial X)]
forall k v. Polynomial k v -> [Term k v]
P.terms UPolynomial Rational
p] Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Rational
1)
    s :: Rational
s = MonomialOrder X -> UPolynomial Rational -> Rational
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 = Interval Rational -> Interval Rational -> Interval Rational
forall r. Ord r => Interval r -> Interval r -> Interval r
Interval.intersection Interval Rational
ival (Rational -> Extended Rational
forall r. r -> Extended r
Finite Rational
lb Extended Rational -> Extended Rational -> Interval Rational
forall r. Ord r => Extended r -> Extended r -> Interval r
<=..<= Rational -> Extended Rational
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 [(X -> Rational) -> UPolynomial Rational -> Rational
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 Rational -> UPolynomial Rational -> Bool
forall k. (Eq k, Num k) => k -> UPolynomial k -> Bool
`P.isRootOf` UPolynomial Rational
p
      then Rational -> Interval Rational
forall r. Ord r => r -> Interval r
Interval.singleton Rational
lb Interval Rational -> [Interval Rational] -> [Interval Rational]
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 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Rational -> Int
n Rational
ub of
        Int
0 -> []
        Int
1 -> [Rational -> Extended Rational
forall r. r -> Extended r
Finite Rational
lb Extended Rational -> Extended Rational -> Interval Rational
forall r. Ord r => Extended r -> Extended r -> Interval r
<..<= Rational -> Extended Rational
forall r. r -> Extended r
Finite Rational
ub]
        Int
_ -> (Rational, Rational) -> [Interval Rational]
g (Rational
lb, Rational
mid) [Interval Rational] -> [Interval Rational] -> [Interval Rational]
forall a. [a] -> [a] -> [a]
++ (Rational, Rational) -> [Interval Rational]
g (Rational
mid, Rational
ub)
      where
        mid :: Rational
mid = (Rational
lb Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
ub) Rational -> Rational -> Rational
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
  | Interval Rational -> Rational
forall r. (Num r, Ord r) => Interval r -> r
Interval.width Interval Rational
ival Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational
0  = Interval Rational
ival
  | SturmChain -> Interval Rational -> Int
numRoots' SturmChain
chain Interval Rational
ivalL Int -> Int -> Bool
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 = Interval Rational -> Extended Rational
forall r. Interval r -> Extended r
Interval.lowerBound Interval Rational
ival
    Finite Rational
ub = Interval Rational -> Extended Rational
forall r. Interval r -> Extended r
Interval.upperBound Interval Rational
ival
    mid :: Rational
mid = (Rational
lb Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
ub) Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
2
    ivalL :: Interval Rational
ivalL = (Extended Rational, Boundary)
-> (Extended Rational, Boundary) -> Interval Rational
forall r.
Ord r =>
(Extended r, Boundary) -> (Extended r, Boundary) -> Interval r
Interval.interval (Interval Rational -> (Extended Rational, Boundary)
forall r. Interval r -> (Extended r, Boundary)
Interval.lowerBound' Interval Rational
ival) (Rational -> Extended Rational
forall r. r -> Extended r
Finite Rational
mid, Boundary
Interval.Closed)
    ivalR :: Interval Rational
ivalR = (Extended Rational, Boundary)
-> (Extended Rational, Boundary) -> Interval Rational
forall r.
Ord r =>
(Extended r, Boundary) -> (Extended r, Boundary) -> Interval r
Interval.interval (Rational -> Extended Rational
forall r. r -> Extended r
Finite Rational
mid, Boundary
Interval.Open) (Interval Rational -> (Extended Rational, Boundary)
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
      | Interval Rational -> Rational
forall r. (Num r, Ord r) => Interval r -> r
Interval.width Interval Rational
ival Rational -> Rational -> Bool
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 = Maybe Rational -> Rational
forall a. HasCallStack => Maybe a -> a
fromJust (Maybe Rational -> Rational) -> Maybe Rational -> Rational
forall a b. (a -> b) -> a -> b
$ Interval Rational -> Maybe Rational
forall r. (Real r, Fractional r) => Interval r -> Maybe r
Interval.pickup (Interval Rational -> Maybe Rational)
-> Interval Rational -> Maybe Rational
forall a b. (a -> b) -> a -> b
$ SturmChain -> Interval Rational -> Rational -> Interval Rational
narrow' SturmChain
chain Interval Rational
ival Rational
epsilon