AERN-RnToRm-0.5: src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/ElementaryInner.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 inner-rounded
elementary functions applied to polynomials.
-}
module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.ElementaryInner
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 Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Elementary
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
{-|
-}
ienclSqrt ::
(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), Bool) ->
((ERChebPoly box b, ERChebPoly box b), Bool)
ienclSqrt maxDegree maxSize ix maxIters (e@(ln, h), isDAC) =
((lnRDown,hRDown), isDAC)
where
lnRDown = chplNeg lRUp
hRDown = chplNeg hnRUp
(_, lRUp) = enclSqrt maxDegree maxSize ix maxIters (ln,chplNeg ln)
(hnRUp, _) = enclSqrt maxDegree maxSize ix maxIters (chplNeg h,h)
--{-|
-- Approximate the pointwise exponential of a polynomial enclosure.
---}
--enclExp ::
-- (B.ERRealBase b, DomainBox box varid Int, Ord box) =>
-- 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.setGranularity 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.setGranularity coeffGr 1)
-- rNear0Encl = enclRAConst rNear0
-- _ = [rNear0Encl, pEncl] -- help the typechecker...
--
--{-|
-- Approximate the pointwise integer power of an enclosure.
---}
--enclPow ::
-- (B.ERRealBase b, Integral i, DomainBox box varid Int, Ord box) =>
-- 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) =>
-- 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.setGranularity coeffGr 1
-- pSqrEncl = enclMultiply maxDegree maxSize pEncl pEncl
-- sineErrorBound =
-- case sineErrorBoundRA of
-- ERInterval lo hi -> hi
-- ERIntervalAny -> B.plusInfinity
-- 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) =>
-- 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.setGranularity coeffGr 1
-- pSqrEncl = enclMultiply maxDegree maxSize pEncl pEncl
-- cosineErrorBound =
-- case cosineErrorBoundRA of
-- ERInterval lo hi -> hi
-- ERIntervalAny -> B.plusInfinity
-- 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) =>
-- 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) =>
-- 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