AERN-RnToRm-0.5: src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs
{-# LANGUAGE FlexibleContexts #-}
{-|
Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Elementary
Description : (internal) elementary functions applied to polynomials
Copyright : (c) 2007-2008 Michal Konecny
License : BSD3
Maintainer : mik@konecny.aow.cz
Stability : experimental
Portability : portable
Internal module for "Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom".
Implementation of elementary functions applied to polynomials.
-}
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
{-|
Approximate the pointwise square root of a polynomial enclosure.
-}
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 {-^ maximum polynomial degree -} ->
Int {-^ maximum term count -} ->
EffortIndex {-^ for calls to other ER functions -} ->
Int {-^ how many times to iterate -} ->
(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)
-- | chplUpperBound ix ln >= 0 =
-- error $ "ERChebPoly: enclSqrt: internal error at eSqrt: cannot show l is positive"
-- | chplLowerBound ix h <= 0 =
-- error $ "ERChebPoly: enclSqrt: internal error at eSqrt: cannot show h is positive"
| otherwise =
-- assuming p >= 1, which implies 1 <= sqrt p <= p and 0 < 1/p <= 1/(sqrt p)
(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))
{-|
Approximate the pointwise exponential of a polynomial enclosure.
-}
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 {-^ maximum polynomial degree -} ->
Int {-^ maximum term count -} ->
EffortIndex {-^ used to derive minimum approx Taylor degree -} ->
(ERChebPoly box b, ERChebPoly box b) ->
(ERChebPoly box b, ERChebPoly box b)
enclExp maxDegree maxSize ix pEncl =
-- unsafePrintReturn
-- (
-- "chplExp:" ++
-- "\n pEncl = " ++ show pEncl ++
-- "\n upperB = " ++ show upperB ++
-- "\n lowerB = " ++ show lowerB ++
-- "\n m = " ++ show m ++
-- "\n expM = " ++ show expM ++
-- "\n r = " ++ show r ++
-- "\n a_int = " ++ show a_int ++
-- "\n a_base = " ++ show a_base ++
-- "\n pNear0Encl = " ++ show (pNear0Encl) ++
-- "\n expNear0 = " ++ show (expNear0) ++
---- "\n chplPow maxDegree (expNear0Up pNear0Up) a_int = " ++ show (chplPow maxDegree (expNear0Up pNear0Up) a_int)
-- "\n result = "
-- )
-- $
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
-- scale the problem down for polynomials whose value is always near zero:
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 -- could this be too high?
expNear0 =
expTayNear0 +: (chplConst 0, chplConst (erintv_right truncCorrNear0))
-- the difference between exact exp and finite Taylor expanstion is an increasing function
-- therefore it is enough to compensate the error at the right-most point
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
-- correction of truncation error (highest at the right-most point):
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] -- help the typechecker...
{-|
Approximate the pointwise integer power of an enclosure.
-}
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 {-^ maximum polynomial degree -} ->
Int {-^ maximum term count -} ->
(ERChebPoly box b, ERChebPoly box b) ->
i ->
(ERChebPoly box b, ERChebPoly box b)
{-^ lower (negated) and upper bound -}
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
{-|
Approximate the pointwise natural logarithm of an enclosure.
-}
enclLog ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
Int {-^ maximum polynomial degree -} ->
Int {-^ maximum term count -} ->
EffortIndex {-^ ?? -} ->
(ERChebPoly box b, ERChebPoly box b) ->
(ERChebPoly box b, ERChebPoly box b)
enclLog maxDegree maxSize ix p =
error "ERChebPoly: chplLog: not implemented yet"
{-|
Approximate the pointwise sine of an enclosure.
Assuming the polynomial range is [-pi/2, pi/2].
-}
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 {-^ maximum polynomial degree -} ->
Int {-^ maximum term count -} ->
EffortIndex {-^ how hard to try (determines Taylor degree and granularity) -} ->
(ERChebPoly box b, ERChebPoly box b) ->
(ERChebPoly box b, ERChebPoly box b)
enclSine maxDegree maxSize ix pEncl =
-- unsafePrint
-- (
-- "ERChebPoly: enclSine: "
-- ++ "\n pEncl = " ++ show pEncl
-- ++ "\n ranLargerEndpoint = " ++ show ranLargerEndpoint
-- ++ "\n sineEncl = " ++ show sineEncl
-- ) $
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
{-|
Approximate the pointwise cosine of an enclosure.
Assuming the polynomial range is [-pi/2, pi/2].
-}
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 {-^ maximum polynomial degree -} ->
Int {-^ maximum term count -} ->
EffortIndex {-^ how hard to try (determines Taylor degree and granularity) -} ->
(ERChebPoly box b, ERChebPoly box b) ->
(ERChebPoly box b, ERChebPoly box b)
enclCosine maxDegree maxSize ix pEncl =
-- unsafePrint
-- (
-- "ERChebPoly: chplCosine: "
-- ++ "\n pEncl = " ++ show pEncl
-- ++ "\n ranLargerEndpoint = " ++ show ranLargerEndpoint
-- ++ "\n cosineEncl = " ++ show cosineEncl
-- ) $
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 {-^ maximum polynomial degree -} ->
Int {-^ maximum term count -} ->
Bool {-^ is sine ? -} ->
(ERChebPoly box b, ERChebPoly box b) ->
Int {-^ how far to go in the Taylor series -} ->
Int {-^ degree of the term now being constructed -} ->
ERInterval b {-^ the coefficient of the term now being constructed -} ->
((ERChebPoly box b, ERChebPoly box b),
Int,
ERInterval b)
{-^
Bounds for the series result and information about the first discarded term,
from which some bound on the uniform error can be deduced.
-}
sincosTaylorAux
maxDegree maxSize resultPositive pSqrEncl tayDegree
thisDegree thisCoeffRA =
sct thisDegree thisCoeffRA
where
sct thisDegree thisCoeffRA
| nextDegree > tayDegree =
-- unsafePrint
-- (
-- "ERChebPoly: sincosTaylorAux: "
-- ++ "\n thisCoeffRA = " ++ show thisCoeffRA
-- ++ "\n nextDegree = " ++ show nextDegree
-- )
(thisCoeffEncl, nextDegree, nextCoeffRA)
| otherwise =
-- unsafePrint
-- (
-- "ERChebPoly: chplSine: taylorAux: "
-- ++ "\n thisCoeffRA = " ++ show thisCoeffRA
-- ++ "\n nextDegree = " ++ show nextDegree
-- ++ "\n errorTermCoeffRA = " ++ show errorTermCoeffRA
-- ++ "\n errorTermDegree = " ++ show errorTermDegree
-- )
(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)
{-|
Approximate the pointwise arcus tangens of an enclosure.
-}
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 {-^ maximum polynomial degree -} ->
Int {-^ maximum term count -} ->
EffortIndex {-^ how far to go in the Euler's series -} ->
(ERChebPoly box b, ERChebPoly box b) ->
(ERChebPoly box b, ERChebPoly box b)
{- arctan using Euler's series:
(http://en.wikipedia.org/wiki/Inverse_trigonometric_function#Infinite_series)
(x / (1 + x^2)) * (1 + t*2*1/(2*1 + 1)*(1 + t*2*2/(2*2 + 1)*(1 + ... (1 + t*2*n/(2*n+1)*(1 + ...)))))
where
t = x^2/(1 + x^2)
where the tail (1 + t*2*n/(2*n+1)*(1 + ...)) is inside the interval:
[1, 1 + x^2]
-}
enclAtan maxDegree maxSize ix pEncl@(pLowNeg, pHigh)
| True = -- pLowerBound >= (-3) && pUpperBound <= 3 =
enclAtanAux maxDegree maxSize ix (Just pSquareEncl) pEncl
| otherwise = -- too far from 0, needs atan(x) = 2*atan(x/(1+sqrt(1+x^2)))
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) -- [-22/14,22/14] is always safe...
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