```{-# LANGUAGE FlexibleContexts #-}
{-|
Module      :  Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Compose
Description :  (internal) composition of 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 pointwise consistently rounded polynomial composition.
-}
module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Compose
where

import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic
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 qualified Data.Number.ER.Real.Approx as RA
import qualified Data.Number.ER.Real.Base as B
import qualified Data.Number.ER.Real.DomainBox as DBox
import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox)
import Data.Number.ER.Misc

import qualified Data.Map as Map

{-|
Compose a polynomial and an enclosure, producing a correcly rounded enclosure,
assuming the second polynomial maps [-1,1] into [-1,1].
-}
enclCompose ::
(B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
Int {-^ max degree for result -} ->
Int {-^ max approx size for result -} ->
ERChebPoly box b {-^ @f@ -} ->
varid {-^ variable @v@ to substitute in @f@ -} ->
(ERChebPoly box b, ERChebPoly box b)
{-^ enclosure of a function @f_v@ to substitute for @v@
that maps @[-1,1]@ into @[-1,1]@  -} ->
(ERChebPoly box b, ERChebPoly box b)
{-^ lower bound and upper bound -}

enclCompose maxDegree maxSize p@(ERChebPoly coeffs) substVar substEncl =
result
{------------------------------
The algorithm: separate from the polynomial
all terms for each degree of the substituted variable,
giving rise to a number of polynomials.
These polynomials are then used as coefficients multiplying
the enclosure evaluations of the Chebyshev polynomials
over the substituted enclosure.
-------------------------------}
where
result =
Map.fold (+:) (enclConst 0) \$ Map.mapWithKey evalDegree degreePolynomialMap
degreePolynomialMap =
Map.foldWithKey extractTerm Map.empty coeffs
extractTerm term c prevPolynomMap =
Map.insertWith Map.union substVarDegree (Map.singleton termNoSubstVar c) prevPolynomMap
where
substVarDegree = DBox.findWithDefault 0 substVar term
termNoSubstVar = DBox.delete substVar term
evalDegree degree degreeCoeffs =
enclMultiply maxDegree maxSize (substPolyDegrees !! degree) (chplNeg degreePoly, degreePoly)
where
degreePoly = ERChebPoly degreeCoeffs
substPolyDegrees =
enclEvalTs maxSize maxDegree substEncl

{------------------------------
The following algorithm is quite wasteful when the polynomial
contains other variables besides the one being substituted.
-------------------------------}
--chplComposeWithEncl maxDegree maxSize p@(ERChebPoly coeffs) substVar substEncl =
--    result
--    where
--    result =
--        foldl (+:) (enclConst 0) \$ map evalTerm \$ Map.toList coeffs
--    evalTerm (term, c) =
--        enclScale c \$
--            foldl (enclMultiply maxDegree maxSize) (enclConst 1) \$
--                map evalVar \$ DBox.toList term
--    evalVar (var, degree) =
--        case var == substVar of
--            True ->
--                substPolyDegrees !! degree
--            False ->
--                (chplNeg varPoly, varPoly)
--        where
--        varPoly =
--            ERChebPoly \$ Map.singleton (DBox.singleton var degree) 1
--    substPolyDegrees =
--        enclEvalTs maxSize maxDegree substEncl

{-|
Compose two polynomials, rounding upwards
provided the second polynomial maps [-1,1] into [-1,1].
-}
enclComposeMany ::
(B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>
Int {-^ max degree for result -} ->
Int {-^ max approx size for result -} ->
ERChebPoly box b ->
Map.Map varid (ERChebPoly box b, ERChebPoly box b)
{-^ variables to substitute and the enclosures to substitute for each of them respectively  -} ->
(ERChebPoly box b, ERChebPoly box b)
{-^ lower bound (negated) and upper bound -}
enclComposeMany maxDegree maxSize p@(ERChebPoly coeffs) substitutions =
result
where
result =
foldl (+:) (enclConst 0) \$ map evalTerm \$ Map.toList coeffs
evalTerm (term, c) =
enclScale maxDegree maxSize c \$
foldl (enclMultiply maxDegree maxSize) (enclConst 1) \$
map evalVar \$ DBox.toList term
evalVar (varID, degree) =
case Map.lookup varID substDegrees of
Nothing ->
(chplNeg varPoly, varPoly)
Just pvDegrees ->
pvDegrees !! degree
where
varPoly =
ERChebPoly \$ Map.singleton (DBox.singleton varID degree) 1
substDegrees =
Map.map mkPVDegrees substitutions
mkPVDegrees pvEncl =
enclEvalTs maxSize maxDegree pvEncl

```