{-# 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 plusDown 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) | otherwise = error $ "chplScaleApprox: " ++ show 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 leqReals or compareApprox 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"