module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Elementary
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 Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Division
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
enclSqrt ::
(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)
enclSqrt maxDegree maxSize ix maxIters p =
result
where
result
| pLowerBound >= 1 =
eSqrt p
| pLowerBound > 0 =
enclRAScale maxDegree maxSize (RAEL.sqrt ix pLowerBoundRA) $
eSqrt $
enclRAScale maxDegree maxSize (recip pLowerBoundRA) p
| otherwise =
error $ "ERChebPoly: enclSqrt: cannot confirm positivity of " ++ show p
pLowerBound = fst $ enclBounds ix p
pLowerBoundRA = ERInterval pLowerBound pLowerBound
eSqrt p@(ln,h)
| otherwise =
(chplMultiplyUp ln lRecipSqrtDown,
chplRecipUp hRecipSqrtDown)
where
lRecipSqrtDown = recipSqrtDown $ chplNeg ln
hRecipSqrtDown = recipSqrtDown $ h
chplMultiplyUp p1 p2 =
chplReduceTermCountUp maxSize $
chplReduceDegreeUp maxDegree $ p1 *^ p2
chplMultiplyDown p1 p2 =
chplReduceTermCountDown maxSize $
chplReduceDegreeDown maxDegree $ p1 *. p2
chplRecipUp p =
snd $
enclRecip maxDegree maxSize ix (maxDegree + 1) $
enclThin p
recipSqrtDown p
| chplLowerBound ix pRecipDown > 0 =
iterRecipSqrt maxIters pRecipDown
| otherwise =
chplConst $ negate $ recip $ negate $ chplUpperBound ix p
where
pRecipDown =
chplNeg $ fst $
enclRecip maxDegree maxSize ix (maxDegree + 1) $
enclThin p
iterRecipSqrt maxIters qNm1
| maxIters > 0 && qNpositive =
iterRecipSqrt (maxIters 1) qN
| otherwise = qNm1
where
qNpositive =
chplLowerBound ix qN > 0
qN =
chplMultiplyDown
(chplScaleDown (0.5) qNm1)
((chplConst 3) -. (chplMultiplyUp p $ chplMultiplyUp qNm1 qNm1))
enclExp ::
(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 ->
(ERChebPoly box b, ERChebPoly box b) ->
(ERChebPoly box b, ERChebPoly box b)
enclExp maxDegree maxSize ix pEncl =
result
where
result =
enclRAScale maxDegree maxSize expM $ enclPow maxDegree maxSize expNear0 a_int
(lowerB, upperB) = enclBounds ix pEncl
mB = (upperB + lowerB) / 2
rB = (upperB lowerB) / 2
r = ERInterval rB rB
m = ERInterval mB mB
expM = max 0 $ erExp_IR ix m
pNear0Encl =
enclRAScale maxDegree maxSize (recip a_base) (pEncl -: (enclConst mB))
rNear0 = r / a_base
a_base = fromInteger a_int
a_int = max 1 $ floor rB
expNear0 =
expTayNear0 +: (chplConst 0, chplConst (erintv_right truncCorrNear0))
expTayNear0 =
expAux pNear0Encl 1 (RA.setGranularityOuter coeffGr 1)
expAux p0Encl nextDegree thisCoeff
| nextDegree > taylorDegree =
enclRAConst thisCoeff
| otherwise =
(enclRAConst thisCoeff) +: (p0Encl *: (expAux p0Encl (nextDegree + 1) nextCoeff))
where
(*:) = enclMultiply maxDegree maxSize
nextCoeff =
thisCoeff / (fromInteger nextDegree)
taylorDegree = 1 + 2 * (ix `div` 6)
coeffGr = effIx2gran $ 10 + 3 * taylorDegree
truncCorrNear0 = expRNear0 tayRNear0
expRNear0 = erExp_R ix rNear0
tayRNear0 =
ERInterval
(negate $ getConst valueAtRNear0LowNeg)
(getConst valueAtRNear0High)
getConst p =
case chplGetConst p of Just c -> c; _ -> 0
(valueAtRNear0LowNeg, valueAtRNear0High) =
expAux rNear0Encl 1 (RA.setGranularityOuter coeffGr 1)
rNear0Encl = enclRAConst rNear0
_ = [rNear0Encl, pEncl]
enclPow ::
(Integral i, 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 ->
(ERChebPoly box b, ERChebPoly box b) ->
i ->
(ERChebPoly box b, ERChebPoly box b)
enclPow maxDegree maxSize pEncl n
| n == 0 =
enclConst 1
| n == 1 =
pEncl
| even n =
powEvenEncl
| odd n =
enclMultiply maxDegree maxSize powEvenEncl pEncl
where
powEvenEncl =
enclMultiply maxDegree maxSize powHalfEncl powHalfEncl
powHalfEncl =
enclPow maxDegree maxSize pEncl halfN
halfN = n `div` 2
enclLog ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
Int ->
Int ->
EffortIndex ->
(ERChebPoly box b, ERChebPoly box b) ->
(ERChebPoly box b, ERChebPoly box b)
enclLog maxDegree maxSize ix p =
error "ERChebPoly: chplLog: not implemented yet"
enclSine ::
(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 ->
(ERChebPoly box b, ERChebPoly box b) ->
(ERChebPoly box b, ERChebPoly box b)
enclSine maxDegree maxSize ix pEncl =
sineEncl
where
sineEncl =
enclAddErr sineErrorBound $
enclMultiply maxDegree maxSize pEncl sineTayEncl
(sineTayEncl, sineErrorTermDegree, sineErrorTermCoeffRA) =
sincosTaylorAux maxDegree maxSize True pSqrEncl taylorDegree 1 one
one = RA.setGranularityOuter coeffGr 1
pSqrEncl = enclMultiply maxDegree maxSize pEncl pEncl
sineErrorBound =
case sineErrorBoundRA of
ERInterval lo hi -> hi
where
sineErrorBoundRA =
(ranLargerEndpointRA ^ sineErrorTermDegree) * sineErrorTermCoeffHighRA
sineErrorTermCoeffHighRA =
snd $ RA.bounds $ abs sineErrorTermCoeffRA
ranLargerEndpointRA =
ERInterval ranLargerEndpoint ranLargerEndpoint
ranLargerEndpoint =
max (abs ranLowB) (abs ranHighB)
(ranLowB, ranHighB) = enclBounds ix pEncl
taylorDegree = effIx2int $ ix `div` 3
coeffGr = effIx2gran $ ix
enclCosine ::
(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 ->
(ERChebPoly box b, ERChebPoly box b) ->
(ERChebPoly box b, ERChebPoly box b)
enclCosine maxDegree maxSize ix pEncl =
cosineEncl
where
cosineEncl =
enclAddErr cosineErrorBound $
cosineTayEncl
(cosineTayEncl, cosineErrorTermDegree, cosineErrorTermCoeffRA) =
sincosTaylorAux maxDegree maxSize True pSqrEncl taylorDegree 0 one
one = RA.setGranularityOuter coeffGr 1
pSqrEncl = enclMultiply maxDegree maxSize pEncl pEncl
cosineErrorBound =
case cosineErrorBoundRA of
ERInterval lo hi -> hi
where
cosineErrorBoundRA =
(ranLargerEndpointRA ^ cosineErrorTermDegree) * cosineErrorTermCoeffHighRA
cosineErrorTermCoeffHighRA =
snd $ RA.bounds $ abs cosineErrorTermCoeffRA
ranLargerEndpointRA =
ERInterval ranLargerEndpoint ranLargerEndpoint
ranLargerEndpoint =
max (abs ranLowB) (abs ranHighB)
(ranLowB, ranHighB) = enclBounds ix pEncl
taylorDegree = effIx2int $ ix `div` 3
coeffGr = effIx2gran $ ix
sincosTaylorAux ::
(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 ->
Bool ->
(ERChebPoly box b, ERChebPoly box b) ->
Int ->
Int ->
ERInterval b ->
((ERChebPoly box b, ERChebPoly box b),
Int,
ERInterval b)
sincosTaylorAux
maxDegree maxSize resultPositive pSqrEncl tayDegree
thisDegree thisCoeffRA =
sct thisDegree thisCoeffRA
where
sct thisDegree thisCoeffRA
| nextDegree > tayDegree =
(thisCoeffEncl, nextDegree, nextCoeffRA)
| otherwise =
(resultEncl, errorTermDegree, errorTermCoeffRA)
where
thisCoeffEncl = enclRAConst thisCoeffRA
resultEncl =
thisCoeffEncl +: (enclMultiply maxDegree maxSize pSqrEncl restEncl)
(restEncl, errorTermDegree, errorTermCoeffRA) =
sct nextDegree nextCoeffRA
nextDegree = thisDegree + 2
nextCoeffRA = thisCoeffRA / nextCoeffDenominatorRA
nextCoeffDenominatorRA =
fromInteger $ toInteger $
negate $ nextDegree * (nextDegree 1)
enclAtan ::
(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 ->
(ERChebPoly box b, ERChebPoly box b) ->
(ERChebPoly box b, ERChebPoly box b)
enclAtan maxDegree maxSize ix pEncl@(pLowNeg, pHigh)
| True =
enclAtanAux maxDegree maxSize ix (Just pSquareEncl) pEncl
| otherwise =
case avoidingDivBy0 of
True ->
enclScale maxDegree maxSize 2 $
enclAtanAux maxDegree maxSize ix Nothing $
enclMultiply maxDegree maxSize pEncl $
enclRecip maxDegree maxSize ix (maxDegree + 1) $
onePlusSqrtOnePlusPSquare
where
(pLowerBound, pUpperBound) = enclBounds ix pEncl
onePlusSqrtOnePlusPSquare =
enclAddConst 1 $
enclSqrt maxDegree maxSize ix maxDegree pSquarePlus1Encl
avoidingDivBy0 =
lower1 > 0 && lower2 > 0
where
(lower1, _) = enclBounds ix pSquarePlus1Encl
(lower2, _) = enclBounds ix onePlusSqrtOnePlusPSquare
pSquareEncl =
enclSquare maxDegree maxSize pEncl
pSquarePlus1Encl =
pSquareEncl +: (enclConst 1)
enclAtanAux maxDegree maxSize ix maybePSquareEncl pEncl@(pLowNeg, pHigh)
| avoidingDivBy0 = resultEncl
| otherwise =
(piHalfUp, piHalfUp)
where
piHalfUp = chplConst $ 22/7
avoidingDivBy0 =
lower > 0
where
(lower, _) = enclBounds ix pSquarePlus1Encl
resultEncl =
enclMultiply maxDegree maxSize
pOverPSquarePlus1Encl preresEncl
preresEncl =
series termsCount 2
termsCount =
max 0 $ ix `div` 3
gran = effIx2gran ix
series termsCount coeffBase
| termsCount > 0 =
enclAddConst 1 $
enclRAScale maxDegree maxSize coeff $
enclMultiply maxDegree maxSize
pSquareOverPSquarePlus1Encl $
series (termsCount 1) (coeffBase + 2)
| otherwise =
enclAddConst 1 $
(chplConst 0, pSquareHigh)
where
coeff = coeffBase / (coeffBase + 1)
pSquareEncl@(pSquareLowNeg, pSquareHigh) =
case maybePSquareEncl of
Just pSquareEncl -> pSquareEncl
Nothing ->
enclSquare maxDegree maxSize pEncl
pSquarePlus1Encl =
pSquareEncl +: (enclConst 1)
recipPSquarePlus1Encl =
enclRecip maxDegree maxSize ix (maxDegree + 1) pSquarePlus1Encl
pSquareOverPSquarePlus1Encl =
enclMultiply maxDegree maxSize pSquareEncl recipPSquarePlus1Encl
pOverPSquarePlus1Encl =
enclMultiply maxDegree maxSize pEncl recipPSquarePlus1Encl