AERN-RnToRm-0.3.0: src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Field.hs
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE UndecidableInstances #-}
{-|
Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Field
Description : (internal) field operations 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 field arithmetic over polynomials
with rounding consistent over the whole domain.
-}
module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Field
where
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic
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, DomainIntBox)
import Data.Number.ER.Misc
import qualified Data.Map as Map
chplAffine ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
b ->
Map.Map varid b ->
ERChebPoly box b
chplAffine at0 varCoeffs =
ERChebPoly $
Map.insert chplConstTermKey at0 $
Map.mapKeys (\ i -> DBox.singleton i 1) varCoeffs
{-|
Convert a polynomial to a lower-order one that is dominated by (resp. dominates)
it closely on the domain [-1,1].
-}
chplReduceDegree ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
Int {-^ new maximal order -} ->
ERChebPoly box b ->
(ERChebPoly box b, ERChebPoly box b) {-^ lower and upper bounds with limited degree -}
chplReduceDegree maxOrder (ERChebPoly coeffs) =
(ERChebPoly newCoeffsDown, ERChebPoly newCoeffsUp)
-- errorModule "chplSetMaxOrder: not implemented yet"
where
newCoeffsUp =
Map.insertWith plusUp chplConstTermKey highOrderCompensation coeffsLowOrder
newCoeffsDown =
Map.insertWith plusDown chplConstTermKey (-highOrderCompensation) coeffsLowOrder
highOrderCompensation =
Map.fold (\ new prev -> prev + (abs new)) 0 coeffsHighOrder
(coeffsHighOrder, coeffsLowOrder) =
Map.partitionWithKey (\ k v -> chplTermOrder k > maxOrder) coeffs
chplReduceDegreeDown m = fst . chplReduceDegree m
chplReduceDegreeUp m = snd . chplReduceDegree m
instance (B.ERRealBase b, DomainBox box varid Int, Ord box) => Num (ERChebPoly box b)
where
fromInteger n =
ERChebPoly $ Map.singleton chplConstTermKey (fromInteger n)
abs (ERChebPoly coeffs) =
errorModule "abs of a polynomial not implemented, use UFB.max instead"
signum (ERChebPoly coeffs) =
errorModule "signum of a polynomial not implemented, use RA.leqReals instead"
--------- negation ----------
negate (ERChebPoly coeffs) =
ERChebPoly $ Map.map negate coeffs
--------- addition ----------
(ERChebPoly coeffs1) + (ERChebPoly coeffs2) =
ERChebPoly sumCoeffs
where
sumCoeffs =
Map.insertWith (+) chplConstTermKey maxError coeffsDown
-- point-wise sum of polynomials with coeff rounding errors
-- compensated for by enlarging the constant term
coeffsUp =
(Map.unionWith (+) coeffs1 coeffs2)
-- point-wise sum of polynomials with coeffs rounded upwards
coeffsDown =
(Map.unionWith (\c1 c2 -> - ((- c1) + (- c2))) coeffs1 coeffs2)
-- point-wise sum of polynomials with coeffs rounded upwards
maxError =
Map.fold (+) 0 $
Map.intersectionWith (-) coeffsUp coeffsDown
-- addition must round upwards on interval [-1,1]
-- non-constant terms are multiplied by quantities in [-1,1]
-- and thus can make the result drop below the exact result
-- -> to compensate add the rounding difference to the constant term
--------- multiplication ----------
(ERChebPoly coeffs1) * (ERChebPoly coeffs2) =
ERChebPoly prodCoeffs
where
prodCoeffs =
Map.insertWith (+) chplConstTermKey roundOffCompensation $
Map.map negate directProdCoeffsDown
roundOffCompensation =
Map.fold (+) 0 $
Map.unionWith (+) directProdCoeffsDown directProdCoeffsUp
(directProdCoeffsUp, directProdCoeffsDown) =
foldl addCombiCoeff (Map.empty, Map.empty) combinedCoeffs
where
addCombiCoeff
(prevCoeffsUp, prevCoeffsDown)
(coeffUp, coeffDown, (powersList, coeffCount)) =
foldl addOnce (prevCoeffsUp, prevCoeffsDown) powersList
where
addOnce (prevCoeffsUp, prevCoeffsDown) powers =
(Map.insertWith (+) powers coeffUpFrac prevCoeffsUp,
Map.insertWith (+) powers coeffDownFrac prevCoeffsDown)
coeffUpFrac = coeffUp / coeffCountB
coeffDownFrac = coeffDown / coeffCountB
coeffCountB = fromInteger coeffCount
combinedCoeffs =
[ -- (list of triples)
(
(c1 * c2) -- upwards rounded product
,
((- c1) * c2) -- downwards rounded negated product
,
combinePowers powers1 powers2
)
|
(powers1, c1) <- coeffs1List,
(powers2, c2) <- coeffs2List
]
combinePowers powers1 powers2 =
(combinedPowers, 2 ^ (length sumsDiffs))
where
combinedPowers =
map (DBox.fromAscList . (filter $ \ (k,v) -> v > 0)) $
allPairsCombinations $
sumsDiffs
sumsDiffs =
-- associative list with the sum and difference of powers for each variable
zipWith (\(k,s) (_,d) -> (k,(s,d)))
(DBox.toAscList $ DBox.unionWith (\a b -> (a + b)) powers1 powers2)
(DBox.toAscList $ DBox.unionWith (\a b -> abs (a - b)) powers1 powers2)
coeffs1List =
Map.toList coeffs1
coeffs2List =
Map.toList coeffs2
-- | multiply a polynomial by a scalar rounding downwards and upwards
chplScale ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
b ->
(ERChebPoly box b) ->
(ERChebPoly box b, ERChebPoly box b)
chplScale ratio (ERChebPoly coeffs) =
(ERChebPoly coeffsDown, ERChebPoly coeffsUp)
where
coeffsDown =
Map.insertWith plusDown chplConstTermKey (- errBound) coeffsScaled
coeffsUp =
Map.insertWith plusUp chplConstTermKey errBound coeffsScaled
(errBound, coeffsScaled) =
Map.mapAccum processTerm 0 coeffs
processTerm errBoundPrev coeff =
(errBoundPrev + errBoundHere, coeffScaledUp)
where
errBoundHere = coeffScaledUp - coeffScaledDown
coeffScaledDown = ratio `timesDown` coeff
coeffScaledUp = ratio `timesUp` coeff
chplScaleDown r = fst . chplScale r
chplScaleUp r = snd . chplScale r
-- | multiply a polynomial by a scalar interval rounding downwards and upwards
chplScaleApprox ::
(B.ERRealBase b, DomainBox box varid Int, Ord box) =>
(b, b) ->
(ERChebPoly box b) ->
(ERChebPoly box b, ERChebPoly box b)
chplScaleApprox (ratioDown, ratioUp) (ERChebPoly coeffs) =
(ERChebPoly coeffsDown, ERChebPoly coeffsUp)
where
coeffsDown =
Map.insertWith plusDown chplConstTermKey (- errBound) coeffsScaled
coeffsUp =
Map.insertWith plusUp chplConstTermKey errBound coeffsScaled
(errBound, coeffsScaled) =
Map.mapAccum processTerm 0 coeffs
processTerm errBoundPrev coeff =
(errBoundPrev + errBoundHere, coeffScaledUp)
where
errBoundHere = coeffScaledUp - coeffScaledDown
(coeffScaledDown, coeffScaledUp)
| coeff >= 0 =
(ratioDown `timesDown` coeff, ratioUp `timesUp` coeff)
| coeff < 0 =
(ratioUp `timesDown` coeff, ratioDown `timesUp` coeff)
instance (B.ERRealBase b, DomainBox box varid Int, Ord box) => Fractional (ERChebPoly box b)
where
fromRational r =
ERChebPoly $ Map.singleton chplConstTermKey (fromRational r)
--------- division ----------
_ / _ =
errorModule "for division use chplRecip from module Elementary"
instance (B.ERRealBase b, DomainBox box varid Int, Ord box) => Ord (ERChebPoly box b)
where
compare _ _ =
errorModule "cannot compare polynomials, consider using RA.leqReals instead"
--instance (B.ERRealBase b, DomainBox box varid Int, Ord box) => Real (ERChebPoly box b)
-- where
-- toRational _ =
-- errorModule "toRational: cannot convert polynomial to rational"
--
--instance (B.ERRealBase b, DomainBox box varid Int, Ord box) => RealFrac (ERChebPoly box b)
-- where
-- properFraction _ =
-- errorModule "properFraction: rounding of polynomials not implemented"