{-# LANGUAGE FlexibleContexts #-}
{-|
Module      :  Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Elementary
Description :  (internal) elementary functions applied to polynomials
Copyright   :  (c) 2007-2008 Michal Konecny

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 =
(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.
-}
chplSineCosine ::
(B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
Bool {-^ True iff sine, False iff cosine -} ->
Int {-^ maximum polynomial degree -} ->
EffortIndex {-^ minimum approx Taylor degree -} ->
ERChebPoly box b ->
(ERChebPoly box b, ERChebPoly box b)
chplSineCosine isSine maxDegree ix p
-- p - 2k*pi range within [-pi/2, pi/2]:
| ranfNear0 `RA.refines` plusMinusPiHalf =
--        unsafePrint
--        (
--            "ERChebPoly: chplSineCosine: [-pi/2, pi/2]: "
--            ++ "\n p = " ++ show p
--            ++ "\n ranf = " ++ show ranf
--            ++ "\n k = " ++ show k
--            ++ "\n ranfNear0 = " ++ show ranfNear0
--        ) \$
case isSine of
True -> sineShifted (- k2pi)
False -> cosineShifted (- k2pi)
-- p - 2k*pi range within [0, pi]:
| (ranfNear0 - piHalf) `RA.refines` plusMinusPiHalf =
--        unsafePrint
--        (
--            "ERChebPoly: chplSineCosine: [0, pi]: "
--            ++ "\n p = " ++ show p
--            ++ "\n ranf = " ++ show ranf
--            ++ "\n k = " ++ show k
--            ++ "\n ranfNear0 = " ++ show ranfNear0
--        ) \$
case isSine of
-- use sin(x) = cos(x - pi/2) and cos(x) = - sin(x - pi/2):
True -> cosineShifted (- k2pi - piHalf)
False -> sineShiftedNegated (- k2pi - piHalf)
-- p - 2k*pi range within [-pi, 0]:
| (ranfNear0 + piHalf) `RA.refines` plusMinusPiHalf =
case isSine of
-- use sin(x) = - cos(x + pi/2) and cos(x) = sin(x + pi/2):
True -> cosineShiftedNegated (-k2pi + piHalf)
False -> sineShifted (-k2pi + piHalf)
-- p - 2k*pi range within [pi/2, 3pi/2]:
| (ranfNear0 - pi) `RA.refines` plusMinusPiHalf =
-- use sin(x) = - sin(x - pi) and cos(x) = - cos(x - pi)
case isSine of
True -> sineShiftedNegated (- k2pi - pi)
False -> cosineShiftedNegated (- k2pi - pi)
| otherwise = (chplConst (-1), chplConst 1)
--    (expDownwards, expUpwards + valueAtRDnNeg + (chplConst expRUp))
where
ranfNear0 = ranf - k2pi
k2pi = k * 2 * pi
plusMinusPiHalf = (-piHalfLO) RA.\/ piHalfLO
pi = RAEL.pi ix
piHalf = pi / 2
(piHalfLO, piHalfHI) = RA.bounds piHalf
ranf =
ERInterval
(negate \$ chplUpperBoundAffine 10 (-p))
(chplUpperBoundAffine 10 p)
k =
fromInteger \$ floor \$
case (pi + ranf) / (2 * pi) of ERInterval lo hi -> lo

sineShiftedNegated shift =
boundsNegate \$ sineShifted shift

cosineShiftedNegated shift =
boundsNegate \$ cosineShifted shift

boundsNegate (pLO, pHI) = (- pHI, - pLO)

sineShifted shift =
boundsAddErr shiftWidthB \$ sineTaylor (p + shiftPoly) (ranf + shift)
where
shiftPoly = chplConst shiftLOB
ERInterval shiftLOB shiftHIB = shift
shiftWidthB = shiftHIB - shiftLOB

cosineShifted shift =
boundsAddErr shiftWidthB \$ cosineTaylor (p + shiftPoly) (ranf + shift)
where
shiftPoly = chplConst shiftLOB
ERInterval shiftLOB shiftHIB = shift
shiftWidthB = shiftHIB - shiftLOB

(pLO `plusDown` (- errPoly), pHI + errPoly)
where
errPoly = chplConst errB

sineTaylor x xran =
(sineDown, sineUp)
where
sineUp =
chplReduceDegreeUp maxDegree \$
x * sineUpTaylor + (chplConst sineUpErrorBound)
(sineUpTaylor, sineUpErrorTermDegree, sineUpErrorTermCoeff) =
taylorAux x 1 (B.setGranularity coeffGr 1)
sineUpErrorBound =
case sineUpErrorBoundRA of ERInterval lo hi -> hi
where
sineUpErrorBoundRA =
(xranLargerEndpoint ^ (1 + sineUpErrorTermDegree)) * sineUpErrorTermCoeffRA
sineUpErrorTermCoeffRA =
abs \$
ERInterval sineUpErrorTermCoeff sineUpErrorTermCoeff
sineDown =
negate \$ chplReduceDegreeUp maxDegree \$
x * sineDownTaylorNeg + (chplConst \$ sineDownErrorBound)
(sineDownTaylorNeg, sineDownErrorTermDegree, sineDownErrorTermCoeff) =
taylorAux x 1 (B.setGranularity coeffGr (-1))
sineDownErrorBound =
case sineDownErrorBoundRA of ERInterval lo hi -> hi
where
sineDownErrorBoundRA =
(xranLargerEndpoint ^ (1 + sineDownErrorTermDegree)) * sineDownErrorTermCoeffRA
sineDownErrorTermCoeffRA =
abs \$
ERInterval sineDownErrorTermCoeff sineDownErrorTermCoeff
xranLargerEndpoint =
max (abs xranLO) (abs xranHI)
(xranLO, xranHI) = RA.bounds xran

cosineTaylor x xran =
--        unsafePrint
--        (
--            "ERChebPoly.Elementary: chplSineCosine: cosineTaylor: "
--            ++ "\n xran = " ++ show xran
--            ++ "\n cosineUpErrorBound = " ++ show cosineUpErrorBound
--            ++ "\n cosineUpErrorTermDegree = " ++ show cosineUpErrorTermDegree
--            ++ "\n cosineUpErrorTermCoeff = " ++ show cosineUpErrorTermCoeff
--            ++ "\n xranLargerEndpoint = " ++ show xranLargerEndpoint
--        )
(cosineDown, cosineUp)
where
cosineUp =
chplReduceDegreeUp maxDegree \$
cosineUpTaylor + (chplConst cosineUpErrorBound)
(cosineUpTaylor, cosineUpErrorTermDegree, cosineUpErrorTermCoeff) =
taylorAux x 0 (B.setGranularity coeffGr 1)
cosineUpErrorBound
| odd (cosineUpErrorTermDegree `div` 2) = 0
| otherwise =
case cosineUpErrorBoundRA of ERInterval lo hi -> hi
where
cosineUpErrorBoundRA =
(xranLargerEndpoint ^ (cosineUpErrorTermDegree)) * cosineUpErrorTermCoeffRA
cosineUpErrorTermCoeffRA =
abs \$
ERInterval cosineUpErrorTermCoeff cosineUpErrorTermCoeff
cosineDown =
negate \$ chplReduceDegreeUp maxDegree \$
cosineDownTaylorNeg + (chplConst \$ cosineDownErrorBound)
(cosineDownTaylorNeg, cosineDownErrorTermDegree, cosineDownErrorTermCoeff) =
taylorAux x 0 (B.setGranularity coeffGr (-1))
cosineDownErrorBound
| even (cosineDownErrorTermDegree `div` 2) = 0
| otherwise =
case cosineDownErrorBoundRA of ERInterval lo hi -> hi
where
cosineDownErrorBoundRA =
(xranLargerEndpoint ^ (cosineDownErrorTermDegree)) * cosineDownErrorTermCoeffRA
cosineDownErrorTermCoeffRA =
abs \$
ERInterval cosineDownErrorTermCoeff cosineDownErrorTermCoeff
xranLargerEndpoint =
max (abs xranLO) (abs xranHI)
(xranLO, xranHI) = RA.bounds xran

taylorAux p0 thisDegree thisCoeff
| nextDegree > taylorDegree =
--                unsafePrint
--                (
--                    "ERChebPoly: chplSine: taylorAux: "
--                    ++ "\n thisCoeff = " ++ show thisCoeff
--                    ++ "\n nextDegree = " ++ show nextDegree
--                )
(chplConst thisCoeff, nextDegree, nextCoeff)
| otherwise =
--                unsafePrint
--                (
--                    "ERChebPoly: chplSine: taylorAux: "
--                    ++ "\n thisCoeff = " ++ show thisCoeff
--                    ++ "\n nextDegree = " ++ show nextDegree
--                    ++ "\n errorTermCoeff = " ++ show errorTermCoeff
--                    ++ "\n errorTermDegree = " ++ show errorTermDegree
--                )
(chplReduceDegreeUp maxDegree \$
(chplConst thisCoeff) + p0 * p0 * rest,
errorTermDegree, errorTermCoeff)
where
(rest, errorTermDegree, errorTermCoeff) =
taylorAux p0 nextDegree nextCoeff
nextDegree = thisDegree + 2
nextCoeff =
thisCoeff / (fromInteger \$ negate \$ nextDegree * (nextDegree - 1))
taylorDegree = ix `div` 3
coeffGr = effIx2gran \$ ix

{-|
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
| 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

tauDegree = effIx2int (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