module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Enclosure
where
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Reduce
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Ring
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Eval
import qualified Data.Number.ER.Real.Base as B
import qualified Data.Number.ER.Real.DomainBox as DBox
import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainIntBox)
import Data.Number.ER.Real.Approx.Interval
import qualified Data.Number.ER.Real.Approx as RA
import Data.Number.ER.Misc
import qualified Data.Map as Map
enclThin ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
ERChebPoly box b ->
(ERChebPoly box b, ERChebPoly box b)
enclThin p =
(chplNeg p, p)
enclConst ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
b ->
(ERChebPoly box b, ERChebPoly box b)
enclConst c =
(chplConst (c), chplConst c)
enclBounds ix (ln, h) =
(negate $ chplUpperBound ix ln, chplUpperBound ix h)
enclEval e@(ln, h) pt
| lB > hB =
unsafePrintReturn
(
"ERChebPoly: enclEval: inverted result:"
++ "\n h = " ++ show h
++ "\n ln = " ++ show ln
++ "\n result = "
)
result
| otherwise = result
where
result = ERInterval lB hB
lB = negate $ chplEvalUp ln pt
hB = chplEvalUp h pt
enclEvalInner (ln, h) pt =
ERInterval
(negate $ chplEvalDown ln pt)
(chplEvalDown h pt)
enclRAEval e@(ln, h) pt =
result
where
result = lRA RA.\/ hRA
lRA = fst $ RA.bounds $ negate $ chplRAEval (\b -> ERInterval b b) ln pt
hRA = snd $ RA.bounds $ chplRAEval (\b -> ERInterval b b) h pt
enclRAEvalInner e@(ln, h) pt =
result
where
result =
ERInterval lB hB
lB =
case negate $ chplRAEval (\b -> ERInterval b b) ln pt of
ERInterval _ lB -> lB
hB =
case chplRAEval (\b -> ERInterval b b) h pt of
ERInterval hB _ -> hB
enclAddErr errB (pLowNeg, pHigh) =
(chplAddConstUp errB pLowNeg, chplAddConstUp errB pHigh)
enclRAConst ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
(ERInterval b) ->
(ERChebPoly box b, ERChebPoly box b)
enclRAConst (ERInterval lo hi) = (chplConst (lo), chplConst hi)
enclRAConst ERIntervalAny = (chplConst (1/0), chplConst (1/0))
enclReduceDegree maxDegree (pLowNeg, pHigh) =
(chplReduceDegreeUp maxDegree pLowNeg, chplReduceDegreeUp maxDegree pHigh)
enclReduceSize maxSize (pLowNeg, pHigh) =
(chplReduceTermCountUp maxSize pLowNeg, chplReduceTermCountUp maxSize pHigh)
enclAddConst c (pLowNeg, pHigh) =
(chplAddConstUp (c) pLowNeg, chplAddConstUp c pHigh)
enclNeg (pLowNeg, pHigh) = (pHigh, pLowNeg)
(p1LowNeg, p1High) +: (p2LowNeg, p2High) =
(p1LowNeg +^ p2LowNeg, p1High +^ p2High)
(p1LowNeg, p1High) -: (p2LowNeg, p2High) =
(p1LowNeg +^ p2High, p1High +^ p2LowNeg)
enclMultiply ::
(B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
Int ->
Int ->
(ERChebPoly box b, ERChebPoly box b) ->
(ERChebPoly box b, ERChebPoly box b) ->
(ERChebPoly box b, ERChebPoly box b)
enclMultiply maxDegr maxSize (ln1, h1) (ln2, h2) =
enclReduceSize maxSize $
enclReduceDegree maxDegr $
case (ln1UpperBound <= 0, h1UpperBound <= 0, ln2UpperBound <= 0, h2UpperBound <= 0) of
(True, _, True, _) ->
(l1l2Neg, h1h2)
(_, True, _, True) ->
(h1h2Neg, l1l2)
(True, _, _, True) ->
(h1l2Neg, l1h2)
(_, True, True, _) ->
(l1h2Neg, l1h2)
_ ->
(
(h1h2Neg `maxP` l1l2Neg) `maxP` (h1l2Neg `maxP` l1h2Neg)
,
(h1h2 `maxP` l1l2) `maxP` (h1l2 `maxP` l1h2)
)
where
ln1UpperBound = chplUpperBound ix ln1
ln2UpperBound = chplUpperBound ix ln2
h1UpperBound = chplUpperBound ix h1
h2UpperBound = chplUpperBound ix h2
ix = 10
maxP = chplMaxUp maxDegr maxSize
h1h2 = h1 *^ h2
h1h2Neg = (chplNeg h1) *^ h2
l1l2 = ln1 *^ ln2
l1l2Neg = (chplNeg ln1) *^ ln2
h1l2 = h1 *^ (chplNeg ln2)
h1l2Neg = h1 *^ ln2
l1h2 = (chplNeg ln1) *^ h2
l1h2Neg = ln1 *^ h2
enclSquare ::
(B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
Int ->
Int ->
(ERChebPoly box b, ERChebPoly box b) ->
(ERChebPoly box b, ERChebPoly box b)
enclSquare maxDegr maxSize (ln, h)
= (minZeroMaxNegSq, maxSq)
where
maxSq = maxP ln2Up h2Up
maxNegSq = maxP (chplNeg ln2Down) (chplNeg h2Down)
minZeroMaxNegSq = chplNonposUp maxDegr maxSize maxNegSq
(ln2Down, ln2Up, _) = chplMultiply ln ln
(h2Down, h2Up, _) = chplMultiply h h
maxP = chplMaxUp maxDegr maxSize
enclScaleNonneg ::
(B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
b ->
(ERChebPoly box b, ERChebPoly box b) ->
(ERChebPoly box b, ERChebPoly box b)
enclScaleNonneg ratio pEncl@(ln, h) =
(ln *^ pRatio, h *^ pRatio)
where
pRatio = chplConst ratio
enclScale ::
(B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
Int ->
Int ->
b ->
(ERChebPoly box b, ERChebPoly box b) ->
(ERChebPoly box b, ERChebPoly box b)
enclScale maxDegree maxSize ratio pEncl =
enclMultiply maxDegree maxSize pEncl (enclConst ratio)
enclRAScale ::
(B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
Int ->
Int ->
(ERInterval b) ->
(ERChebPoly box b, ERChebPoly box b) ->
(ERChebPoly box b, ERChebPoly box b)
enclRAScale maxDegree maxSize ra pEncl =
enclMultiply maxDegree maxSize pEncl (enclRAConst ra)
chplScaleRA ::
(B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
Int ->
Int ->
ERInterval b ->
ERChebPoly box b ->
(ERChebPoly box b, ERChebPoly box b)
chplScaleRA maxDegr maxSize (ERIntervalAny) p = enclRAConst ERIntervalAny
chplScaleRA maxDegr maxSize (ERInterval ratioDown ratioUp) p =
(scaledPDownNeg, scaledPUp)
where
(scaledPDownNeg, scaledPUp) =
enclMultiply maxDegr maxSize
(chplNeg p, p) (chplConst ( ratioDown), chplConst ratioUp)
chplScaleRADown m n r = chplNeg . fst . chplScaleRA m n r
chplScaleRAUp m n r = snd . chplScaleRA m n r
enclEvalTs ::
(B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
Int ->
Int ->
(ERChebPoly box b, ERChebPoly box b) ->
[(ERChebPoly box b, ERChebPoly box b)]
enclEvalTs maxDegree maxSize p1@(pLowNeg, pHigh) =
chebyIterate (enclConst 1) p1
where
chebyIterate pNm2 pNm1 =
pNm2 : (chebyIterate pNm1 pN)
where
pN =
(enclScale maxDegree maxSize 2 $
enclMultiply maxDegree maxSize p1 pNm1)
-: pNm2
enclThinTimes ::
(B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
Int ->
Int ->
ERChebPoly box b ->
(ERChebPoly box b, ERChebPoly box b) ->
(ERChebPoly box b, ERChebPoly box b)
enclThinTimes maxDegree maxSize p1 (p2LowNeg, p2High) =
(prodLowNeg, prodHigh)
where
prodHigh =
chplMaxUp maxDegree maxSize
(p1 *^ p2High)
(p1n *^ p2LowNeg)
prodLowNeg =
chplMaxUp maxDegree maxSize
(p1n *^ p2High)
(p1 *^ p2LowNeg)
p1n = chplNeg p1