module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Division
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.Eval
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.Enclosure
import qualified Data.Number.ER.Real.Approx as RA
import qualified Data.Number.ER.Real.Approx.Elementary as RAEL
import qualified Data.Number.ER.Real.Base as B
import Data.Number.ER.Real.Approx.Interval
import Data.Number.ER.Real.Arithmetic.Elementary
import qualified Data.Number.ER.BasicTypes.DomainBox as DBox
import Data.Number.ER.BasicTypes.DomainBox (VariableID(..), DomainBox, DomainIntBox, DomainBoxMappable)
import Data.Number.ER.BasicTypes
import Data.Number.ER.Misc
import qualified Data.Map as Map
enclRecip ::
(B.ERRealBase b,
DomainBox box varid Int, Ord box, Show varid,
DomainIntBox boxra varid (ERInterval b),
DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) =>
Int ->
Int ->
EffortIndex ->
Int ->
(ERChebPoly box b, ERChebPoly box b) ->
(ERChebPoly box b, ERChebPoly box b)
enclRecip maxDegree maxSize ix tauDegr pEncl@(pLowNeg, pHigh)
| pIsConst =
enclRAConst (recip pConst)
| upperB < 0 =
enclNeg $ enclRecip maxDegree maxSize ix tauDegr (enclNeg pEncl)
| lowerB > 0 =
case allInterimsBounded of
True -> resEncl
False -> (chplConst 0, chplConst (B.plusInfinity))
| otherwise =
error $
"ERChebPoly: enclRecip: "
++ "cannot deal with estimated range " ++ show ranp
++ "of polynomial enclosure: \n" ++ show pEncl
where
ranp = ERInterval lowerB upperB
(lowerB, upperB) = enclBounds ix pEncl
(pIsConst, pConst) =
case (chplGetConst pLowNeg, chplGetConst pHigh) of
(Just pConstLowNeg, Just pConstHigh) ->
(True, ERInterval ( pConstLowNeg) pConstHigh)
_ ->
(False, error "ERChebPoly: chplRecip: internal error")
tauDegree = max 2 tauDegr
coeffGr = effIx2gran $ ix
k = intLogUp 2 $ ceiling (recip lowerB)
upperBtr = upperB * 2^k
pAbove1Encl =
enclScale maxDegree maxSize (2^k) pEncl
trT1Encl =
enclAddConst (1) (enclRAScale maxDegree maxSize nu (enclAddConst (1) pAbove1Encl))
nu = recip nuInv
nuInv = (RA.setMinGranularityOuter coeffGr (ERInterval upperBtr upperBtr) 1) / 2
nuPlus1 = nu + 1
nuInvPlus1 = nuInv + 1
nuInvDiv2 = nuInv / 2
trTis =
enclEvalTs maxDegree maxSize trT1Encl
resEncl = (resLowNeg, resHigh)
resLowNeg =
chplScaleUp (2^k) $
chplScaleUp errScaleDownB $
scaledTrTisSumLowNeg
resHigh
| errScaleUpB > 0 =
chplScaleUp (2^k) $
chplScaleUp errScaleUpB $
scaledTrTisSumHigh
| otherwise =
chplScaleUp (2^k) $
chplAddConstUp errAddUpB scaledTrTisSumHigh
ERInterval errScaleDownB _ = nuOverNuPlusTauAns
nuOverNuPlusTauAns = (nu / (nu + tauAbs))
ERInterval _ errScaleUpB = nuOverNuMinusTauAns
nuOverNuMinusTauAns = (nu / (nu tauAbs))
ERInterval _ errAddUpB = tauAbsTimesNuInv
tauAbsTimesNuInv = tauAbs * nuInv
allInterimsBounded =
and $ map RA.isBounded [nuOverNuPlusTauAns, nuOverNuMinusTauAns, nuOverNuMinusTauAns]
tauAbs = abs tau
tau = recip tauInv
(scaledTrTisSumLowNeg, scaledTrTisSumHigh) =
foldl1 (+:) $ zipWith scaleTerm c0n trTis
scaleTerm c trTEncl =
enclRAScale maxDegree maxSize (c * tau) trTEncl
c0n = c0 : c1n
tauInv = c0 * nuInvPlus1 + c1 * nuInvDiv2
c0 = c1 * nuPlus1 c2/2
(c1 : c2 : _) = c1n
c1n = reverse $ take n $ csRev
n = tauDegree
csRev =
cn : cnM1 : (csAux cn cnM1)
where
cn = 1
cnM1 = 2 * nuPlus1
csAux cn cnM1 =
cnM2 : (csAux cnM1 cnM2)
where
cnM2 = cn 2 * nuPlus1 * cnM1