AERN-RnToRm-0.4.2: 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.Eval
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Field
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds
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.Real.DomainBox as DBox
import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainIntBox)
import Data.Number.ER.BasicTypes
import Data.Number.ER.Misc
import qualified Data.Map as Map
{-|
Approximate the pointwise square root of a polynomial
by another polynomial from below and from above.
-}
chplSqrt ::
(B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
Int {-^ maximum polynomial degree -} ->
EffortIndex {-^ ?? -} ->
ERChebPoly box b ->
(ERChebPoly box b, ERChebPoly box b)
chplSqrt maxDegree ix p =
error "ERChebPoly: chplSqrt: not implemented yet"
{-|
Approximate the pointwise exponential of a polynomial
by another polynomial from below and from above.
-}
chplExp ::
(B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
Int {-^ maximum polynomial degree -} ->
EffortIndex {-^ minimum approx Taylor degree -} ->
ERChebPoly box b ->
(ERChebPoly box b, ERChebPoly box b)
chplExp maxDegree ix p =
-- unsafePrint
-- (
-- "chplExp:" ++
-- "\n expM = " ++ show expM ++
-- "\n upperB = " ++ show upperB ++
-- "\n lowerB = " ++ show lowerB ++
-- "\n a_int = " ++ show a_int ++
-- "\n expNear0Dn pNear0Dn = " ++ show (expNear0Dn pNear0Dn) ++
-- "\n chplPow maxDegree (expNear0Up pNear0Up) 2000 = " ++ show (chplPow maxDegree (expNear0Up pNear0Up) 2000)
-- )
-- $
(expDownwards, expUpwards + valueAtRDnNeg + (chplConst expRUp))
where
expUpwards =
(chplConst expMUp) * (chplPow maxDegree (expNear0Up pNear0Up) a_int)
expDownwards =
(chplConst expMDn) * (chplPow maxDegree (expNear0Dn pNear0Dn) a_int)
upperB = chplUpperBoundAffine ix p
lowerB = - (chplUpperBoundAffine ix (- p))
m = (upperB + lowerB) / 2
r = (upperB - lowerB) / 2
expMUp = erintv_right expM
expMDn = erintv_left expM
expM =
erExp_R ix (ERInterval m m)
pNear0Up = (p - (chplConst m)) * (chplConst $ recip a_base)
pNear0Dn = - (((-p) + (chplConst m)) * (chplConst $ recip a_base))
a_base = fromInteger a_int
a_int = max 1 $ floor r -- could this be too high?
expNear0Up p0 =
expAux p0 1 (B.setGranularity coeffGr 1)
expNear0Dn p0 =
negate $ expAux p0 1 (B.setGranularity coeffGr (-1))
expAux p0 nextDegree thisCoeff
| nextDegree > taylorDegree =
chplConst thisCoeff
| otherwise =
snd $ chplReduceDegree maxDegree $
(chplConst thisCoeff) + p0 * (expAux p0 (nextDegree + 1) nextCoeff)
where
nextCoeff =
thisCoeff / (fromInteger nextDegree)
taylorDegree = 1 + 2 * (ix `div` 6)
coeffGr = effIx2gran $ 10 + 3 * taylorDegree
expRUp = erintv_right expR
expR = erExp_R ix (ERInterval r r)
valueAtRDnNeg =
expAux (chplConst r) 1 (B.setGranularity coeffGr (-1))
{-|
Approximate the pointwise integer power of a polynomial by another polynomial from above.
-}
chplPow ::
(B.ERRealBase b, Integral i, DomainBox box varid Int, Ord box) =>
Int {-^ maximum polynomial degree -} ->
ERChebPoly box b ->
i ->
ERChebPoly box b
chplPow maxDegree p n
| n == 0 =
chplConst 1
| n == 1 =
p
| even n =
snd $ chplReduceDegree maxDegree $ powHalfN * powHalfN
| odd n =
snd $ chplReduceDegree maxDegree $
p *
(snd $ chplReduceDegree maxDegree $
powHalfN * powHalfN)
where
powHalfN =
chplPow maxDegree p halfN
halfN = n `div` 2
{-|
Approximate the pointwise natural logarithm of a polynomial
by another polynomial from below and from above.
-}
chplLog ::
(B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
Int {-^ maximum polynomial degree -} ->
EffortIndex {-^ ?? -} ->
ERChebPoly box b ->
(ERChebPoly box b, ERChebPoly box b)
chplLog maxDegree ix p =
error "ERChebPoly: chplLog: not implemented yet"
{-|
Approximate the pointwise sine of a polynomial
by another polynomial from below and from above.
Assuming the polynomial range is [-pi/2, pi/2].
-}
chplSine ::
(B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
Int {-^ maximum polynomial degree -} ->
EffortIndex {-^ how hard to try (determines Taylor degree and granularity) -} ->
ERChebPoly box b ->
(ERChebPoly box b, ERChebPoly box b)
chplSine maxDegree ix p =
-- unsafePrint
-- (
-- "ERChebPoly: sineTaylor: "
-- ++ "\n p = " ++ show p
-- ++ "\n ranLargerEndpoint = " ++ show ranLargerEndpoint
-- ++ "\n sineUp = " ++ show sineUp
-- ++ "\n sineDown = " ++ show sineDown
-- ) $
(sineDown, sineUp)
where
(sineDown, sineUp) =
boundsAddErr sineErrorBound $
chplThinTimesEncl maxDegree p (sineDownTaylor, sineUpTaylor)
((sineDownTaylor, sineUpTaylor),
sineErrorTermDegree,
(sineErrorTermCoeffDown, sineErrorTermCoeffUp)) =
sincosTaylorAux True (chplSquare maxDegree p) taylorDegree 1 (one, one)
one = B.setGranularity coeffGr 1
sineErrorBound =
case sineErrorBoundRA of ERInterval lo hi -> hi
where
sineErrorBoundRA =
(ranLargerEndpointRA ^ (sineErrorTermDegree)) * sineErrorTermCoeffRA
sineErrorTermCoeffRA =
ERInterval sineErrorTermCoeff sineErrorTermCoeff
sineErrorTermCoeff =
max (abs sineErrorTermCoeffDown) (abs sineErrorTermCoeffUp)
ranLargerEndpointRA =
ERInterval ranLargerEndpoint ranLargerEndpoint
ranLargerEndpoint =
max (abs ranLO) (abs ranHI)
ranLO = negate $ chplUpperBoundAffine ix (-p)
ranHI = chplUpperBoundAffine ix p
taylorDegree = effIx2int $ ix `div` 3
coeffGr = effIx2gran $ ix
boundsAddErr errB (pLO, pHI) =
(pLO `plusDown` (- errPoly), pHI + errPoly)
where
errPoly = chplConst errB
{-|
Approximate the pointwise sine of a polynomial
by another polynomial from below and from above.
Assuming the polynomial range is [-pi/2, pi/2].
-}
chplCosine ::
(B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
Int {-^ maximum polynomial degree -} ->
EffortIndex {-^ how hard to try (determines Taylor degree and granularity) -} ->
ERChebPoly box b ->
(ERChebPoly box b, ERChebPoly box b)
chplCosine maxDegree ix p =
-- unsafePrint
-- (
-- "ERChebPoly: chplCosine: "
-- ++ "\n p = " ++ show p
-- ++ "\n ranLargerEndpoint = " ++ show ranLargerEndpoint
-- ++ "\n cosineUp = " ++ show cosineUp
-- ++ "\n cosineDown = " ++ show cosineDown
-- ) $
(cosineDown, cosineUp)
where
(cosineDown, cosineUp) =
boundsAddErr cosineErrorBound $
(cosineDownTaylor, cosineUpTaylor)
((cosineDownTaylor, cosineUpTaylor),
cosineErrorTermDegree,
(cosineErrorTermCoeffDown, cosineErrorTermCoeffUp)) =
sincosTaylorAux True (chplSquare maxDegree p) taylorDegree 0 (one, one)
one = B.setGranularity coeffGr 1
cosineErrorBound =
case cosineErrorBoundRA of ERInterval lo hi -> hi
where
cosineErrorBoundRA =
(ranLargerEndpointRA ^ (cosineErrorTermDegree)) * cosineErrorTermCoeffRA
cosineErrorTermCoeffRA =
ERInterval cosineErrorTermCoeff cosineErrorTermCoeff
cosineErrorTermCoeff =
max (abs cosineErrorTermCoeffDown) (abs cosineErrorTermCoeffUp)
ranLargerEndpointRA =
ERInterval ranLargerEndpoint ranLargerEndpoint
ranLargerEndpoint =
max (abs ranLO) (abs ranHI)
ranLO = negate $ chplUpperBoundAffine ix (-p)
ranHI = chplUpperBoundAffine ix p
taylorDegree = effIx2int $ ix `div` 3
coeffGr = effIx2gran $ ix
sincosTaylorAux ::
(B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
Bool ->
(ERChebPoly box b, ERChebPoly box b) ->
Int {-^ how far to go in the Taylor series -} ->
Int {-^ degree of the term now being constructed -} ->
(b,b) ->
((ERChebPoly box b, ERChebPoly box b),
Int,
(b,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 resultPositive pSquares@(pSquareDown, pSquareUp)
maxDegree thisDegree (thisCoeffDown, thisCoeffUp)
| nextDegree > maxDegree =
-- unsafePrint
-- (
-- "ERChebPoly: chplSine: taylorAux: "
-- ++ "\n thisCoeff = " ++ show thisCoeff
-- ++ "\n nextDegree = " ++ show nextDegree
-- )
((thisCoeffDownP, thisCoeffUpP), nextDegree, (nextCoeffDown, nextCoeffUp))
| otherwise =
-- unsafePrint
-- (
-- "ERChebPoly: chplSine: taylorAux: "
-- ++ "\n thisCoeff = " ++ show thisCoeff
-- ++ "\n nextDegree = " ++ show nextDegree
-- ++ "\n errorTermCoeff = " ++ show errorTermCoeff
-- ++ "\n errorTermDegree = " ++ show errorTermDegree
-- )
((resultDown, resultUp), errorTermDegree, errorTermCoeffs)
where
thisCoeffDownP = chplConst thisCoeffDown
thisCoeffUpP = chplConst thisCoeffUp
resultDown
| resultPositive =
-- ie rest's ideal value is negative and thisCoeff is positive
chplReduceDegreeDown maxDegree $
thisCoeffDownP `plusDown` (pSquareUp `timesDown` restDown)
| otherwise =
-- ie rest's ideal value is positive and thisCoeff is negative
chplReduceDegreeDown maxDegree $
thisCoeffDownP `plusDown` (pSquareDown `timesDown` restDown)
resultUp
| resultPositive =
-- ie rest's ideal value is negative and thisCoeff is positive
chplReduceDegreeUp maxDegree $
thisCoeffUpP `plusUp` (pSquareDown `timesUp` restUp)
| otherwise =
-- ie rest's ideal value is positive and thisCoeff is negative
chplReduceDegreeUp maxDegree $
thisCoeffUpP `plusUp` (pSquareUp `timesUp` restUp)
((restDown, restUp), errorTermDegree, errorTermCoeffs) =
sincosTaylorAux (not resultPositive) pSquares maxDegree nextDegree (nextCoeffDown, nextCoeffUp)
nextDegree = thisDegree + 2
nextCoeffUp
| resultPositive =
thisCoeffDown / nextCoeffDenominator -- positive / negative
| otherwise =
thisCoeffUp / nextCoeffDenominator -- negative / negative
nextCoeffDown
| resultPositive =
thisCoeffUp `divDown` nextCoeffDenominator -- positive / negative
| otherwise =
thisCoeffDown `divDown` nextCoeffDenominator -- negative / negative
nextCoeffDenominator =
fromInteger $ toInteger $ negate $ nextDegree * (nextDegree - 1)
divDown a b = negate $ a / (negate b)
{-|
Approximate the pointwise natural logarithm of a polynomial
by another polynomial from below and from above.
-}
chplAtan ::
(B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
Int {-^ maximum polynomial degree -} ->
EffortIndex {-^ ?? -} ->
ERChebPoly box b ->
(ERChebPoly box b, ERChebPoly box b)
{- atan using Euler's 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 + (x^2*2n/(2n + 1)), 1 + x^2]
-}
chplAtan maxDegree ix p
| avoidingDivBy0 =
-- unsafePrint
-- (
-- "ERChebPoly.Elementary: chplAtan: "
-- ++ "\n maxDegree = " ++ show maxDegree
-- ++ "\n p = " ++ show p
-- ++ "\n pSquareDn = " ++ show pSquareDn
-- ++ "\n pSquareUp = " ++ show pSquareUp
-- ++ "\n pOverPSquarePlus1Dn = " ++ show pOverPSquarePlus1Dn
-- ++ "\n pOverPSquarePlus1Up = " ++ show pOverPSquarePlus1Up
-- ++ "\n preresDn = " ++ show preresDn
-- ++ "\n preresUp = " ++ show preresUp
-- ++ "\n resDn = " ++ show resDn
-- ++ "\n resUp = " ++ show resUp
-- )
(resDn, resUp)
| otherwise =
(chplConst (-2), chplConst 2) -- this is always safe...
where
avoidingDivBy0 =
(chplUpperBoundAffine ix (- pSquarePlus1Dn) < 0)
&&
(chplUpperBoundAffine ix (- pSquarePlus1Up) < 0)
resDn =
negate $
chplMaxUp maxDegree
(chplReduceDegreeUp maxDegree $
pOverPSquarePlus1DnNeg `timesUp` preresDn) -- beware: pOverPSquarePlus1Dn can be negative
(chplReduceDegreeUp maxDegree $
pOverPSquarePlus1DnNeg `timesUp` preresUp)
where
pOverPSquarePlus1DnNeg = negate pOverPSquarePlus1Dn
resUp =
chplMaxUp maxDegree
(chplReduceDegreeUp maxDegree $
pOverPSquarePlus1Up `timesUp` preresDn) -- beware: pOverPSquarePlus1Up can be negative
(chplReduceDegreeUp maxDegree $
pOverPSquarePlus1Up `timesUp` preresUp)
(preresDn, preresUp) = seriesDnUp termsCount 2
termsCount = max 0 $ ix `div` 3
gran = effIx2gran ix
seriesDnUp termsCount coeffBase
| termsCount > 0 =
(
chplReduceDegreeDown maxDegree $
1 `plusDown`
(pSquareOverPSquarePlus1Dn -- >=0
`timesDown` (chplConst coeffDn) -- >=0
`timesDown` restDn -- >=0
)
,
chplReduceDegreeUp maxDegree $
1 `plusUp`
(pSquareOverPSquarePlus1Up -- >=0
`timesUp` (chplConst coeffUp) -- >=0
`timesUp` restUp -- >=0
)
)
| otherwise =
(
1 `plusDown` (pSquareDn `timesDown` (chplConst coeffDn)) -- both >=0
,
1 `plusUp` pSquareUp
)
where
(restDn, restUp) = seriesDnUp (termsCount - 1) (coeffBase + 2)
coeffUp = coeffBaseB / (coeffBaseB `plusDown` 1)
coeffDn = negate $ coeffBaseB / (negate $ coeffBaseB `plusUp` 1)
coeffBaseB = B.setMinGranularity gran $ fromInteger coeffBase
(pSquareDn, pSquareUp) = chplSquare maxDegree p
pSquarePlus1Dn = pSquareDn `plusDown` 1
pSquarePlus1Up = pSquareUp `plusUp` 1
recipPSquarePlus1Dn = chplRecipDn maxDegree ix pSquarePlus1Up
recipPSquarePlus1Up = chplRecipUp maxDegree ix pSquarePlus1Dn
-- -- safely compute the square of an enclosure:
-- pSquareDn = chplMinDn m pUpTDnpUp (chplMinDn m pDnTDnpUp pDnTDnpDn)
-- pSquareUp = chplMaxUp m pUpTUppUp (chplMaxUp m pDnTUppUp pDnTUppDn)
-- pUpTDnpUp = pUp `timesDown` pUp
-- pDnTDnpUp = pDn `timesDown` pUp
-- pDnTDnpDn = pDn `timesDown` pDn
-- pUpTUppUp = pUp `timesUp` pUp
-- pDnTUppUp = pDn `timesUp` pUp
-- pDnTUppDn = pDn `timesUp` pDn
-- mMinus1 = m - 1
pSquareOverPSquarePlus1Up =
pSquareUp `timesUp` recipPSquarePlus1Up -- both >=0
pSquareOverPSquarePlus1Dn =
pSquareDn `timesDown` recipPSquarePlus1Dn -- both >=0 (one enclosure may dip below 0, not a problem)
-- negate $
-- chplMaxUp maxDegree
-- (pSquareDnNeg `timesUp` recipPSquarePlus1Up) -- beware: pSquareDn may dip below 0
-- (pSquareDnNeg `timesUp` recipPSquarePlus1Dn)
-- where
-- pSquareDnNeg = negate pSquareDn
pOverPSquarePlus1Up =
chplMaxUp maxDegree
(p `timesUp` recipPSquarePlus1Up)
(p `timesUp` recipPSquarePlus1Dn) -- beware: p can be negative
pOverPSquarePlus1Dn =
negate $
chplMaxUp maxDegree
(pn `timesUp` recipPSquarePlus1Up) -- beware: pn can be positive
(pn `timesUp` recipPSquarePlus1Dn)
where
pn = negate p
chplRecipDn m i = fst . chplRecip m i
chplRecipUp m i = snd . chplRecip m i
{-|
Approximate the pointwise cosine of a polynomial
by another polynomial from below and from above
using the tau method
as described in [Mason & Handscomb 2003, p 62].
-}
chplRecip ::
(B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
Int {-^ maximum polynomial degree -} ->
EffortIndex {-^ minimum approx degree -} ->
ERChebPoly box b ->
(ERChebPoly box b, ERChebPoly box b)
chplRecip maxDegree ix p@(ERChebPoly pCoeffs)
| pIsConst =
(chplConst $ - (recip (- pConst)), chplConst $ recip pConst)
| upperB < 0 = -- range negative
(\(lo, hi) -> (-hi, -lo)) $ chplRecip maxDegree ix (negate p)
| lowerB > 0 = -- range positive
-- unsafePrint
-- (
-- "ERChebPoly: chplRecip: "
-- ++ "\n k = " ++ show k
-- ++ "\n lowerB = " ++ show lowerB
-- ++ "\n tau = " ++ (show $ recip tauInv)
-- )
(resDn, resUp)
| otherwise = -- cannot establish 0 freedom
error $
"ERChebPoly: chplRecip: "
++ "cannot deal with estimated range " ++ show ranp
++ "of polynomial: \n" ++ show p
where
ranp = ERInterval lowerB upperB
lowerB = - (chplUpperBoundAffine ix (- p))
upperB = chplUpperBoundAffine ix p
(pIsConst, pConst) =
case chplGetConst p of
Nothing -> (False, error "ChebyshevBase.Polynom.Elementary.chplRecip")
Just coeff -> (True, coeff)
tauDegree = effIx2int (max 2 $ ix `div` 3)
coeffGr = effIx2gran $ ix
-- translate p to have range above 1:
k = intLogUp 2 $ ceiling (1/lowerB) -- 2^k * lowerB >= 1
upperBtr = upperB * 2^k -- upper bound of translated poly
(pAbove1Dn, pAbove1Up) = -- p multiplied by 2^k; range in [1,upperBtr]
chplScale (2^k) p
-- translate T_1 to domain [0, upperBtr] and apply it to (pAbove1 - 1):
-- T'_1 = nu * (p - 1) - 1
trT1Dn =
(chplScaleDown nuLOB (pAbove1Dn - 1)) - 1
trT1Up =
(chplScaleUp nuHIB (pAbove1Up - 1)) - 1
nu = recip nuInv -- auxiliary constant
ERInterval nuLOB nuHIB = nu
nuInv = (RA.setMinGranularity coeffGr (ERInterval upperBtr upperBtr)) / 2
nuPlus1 = nu + 1
nuInvPlus1 = nuInv + 1
nuInvDiv2 = nuInv / 2
-- define such translated T_i's for all i >= 0:
trTis =
map (mapPair (chplReduceDegreeDown maxDegree, chplReduceDegreeUp maxDegree)) $
chebyEvalTsRoundDownUp trT1Dn
-- construct the result from interval coefficients:
resDn =
chplScaleDown (2^k) $
(-tauAbsUpPoly) `plusDown`
(chplScaleUp tauAbsDnB $
sumDown $
(- errPoly) : (zipWith scaleDn cis trTis))
resUp =
chplScaleUp (2^k) $
(tauAbsUpPoly) `plusUp`
(chplScaleUp tauAbsUpB $
sumUp $
(errPoly) : (zipWith scaleUp cis trTis))
scaleDn c (trTDn, trTUp)
| r >= 0 = chplScaleDown r trTDn
| otherwise = chplScaleDown r trTUp
where
r = c * tauSign
scaleUp c (trTDn, trTUp)
| r >= 0 = chplScaleUp r trTUp
| otherwise = chplScaleUp r trTDn
where
r = c * tauSign
tauAbsUpPoly = chplConst $ tauAbsUpB
tauSign =
case RA.compareReals tauInv 0 of
Just GT -> 1
Just LT -> -1
ERInterval tauAbsDnB tauAbsUpB = abs $ recip tauInv
cis =
map (\(ERInterval lo hi) -> hi) c0n
errPoly = chplConst err
err =
foldl1 plusUp $
map (\(ERInterval lo hi) -> hi - lo) c0n
-- work out the coefficients in interval arithmetic using the tau method:
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