AERN-RnToRm 0.4.2 → 0.4.9
raw patch · 29 files changed
+3893/−1504 lines, 29 filesdep +QuickCheckdep +directorydep +filepathdep ~AERN-Realdep ~base
Dependencies added: QuickCheck, directory, filepath, time
Dependency ranges changed: AERN-Real, base
Files
- AERN-RnToRm.cabal +28/−22
- ChangeLog +13/−1
- src/Data/Number/ER/RnToRm/Approx.hs +113/−16
- src/Data/Number/ER/RnToRm/Approx/DomTransl.hs +77/−12
- src/Data/Number/ER/RnToRm/TestingDefs.hs +25/−6
- src/Data/Number/ER/RnToRm/UnitDom/Approx.hs +82/−3
- src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs +341/−218
- src/Data/Number/ER/RnToRm/UnitDom/Base.hs +292/−247
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs +63/−31
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Basic.hs +92/−34
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Bounds.hs +217/−131
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Compose.hs +138/−0
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Division.hs +168/−0
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs +249/−418
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Enclosure.hs +311/−0
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Eval.hs +23/−76
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Field.hs +0/−228
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Integration.hs +59/−61
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Reduce.hs +85/−0
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Ring.hs +218/−0
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Bounds.hs +46/−0
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Compose.hs +114/−0
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Division.hs +78/−0
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Elementary.hs +120/−0
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Enclosure.hs +106/−0
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Generate.hs +592/−0
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Reduce.hs +37/−0
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Ring.hs +47/−0
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Run.hs +159/−0
AERN-RnToRm.cabal view
@@ -1,5 +1,5 @@ Name: AERN-RnToRm-Version: 0.4.2+Version: 0.4.9 Cabal-Version: >= 1.2 Build-Type: Simple License: BSD3@@ -10,7 +10,7 @@ Stability: experimental Category: Data, Math Synopsis: polynomial function enclosures (PFEs) approximating exact real functions-Tested-with: GHC ==6.8.3+Tested-with: GHC ==6.10.1 Description: AERN-RnToRm provides datatypes and abstractions for approximating functions@@ -32,43 +32,49 @@ with Taylor Models is included in the paper <http://www-users.aston.ac.uk/~konecnym/papers/cfv08.html>. .- Simple examples of usage can be found in module @Demo.hs@ in folder @tests@.+ Simple examples of usage can be found in folder @tests@. Extra-source-files: tests/Demo.hs tests/ISin3.hs Data-files: ChangeLog -Flag containers-in-base- Default: False- Library hs-source-dirs: src- if flag(containers-in-base)- Build-Depends:- base < 3, binary >= 0.4, html >= 1.0, AERN-Real >= 0.9.7- else- Build-Depends:- base >= 3, containers, binary >= 0.4, html >= 1.0, AERN-Real >= 0.9.7+ Build-Depends:+ AERN-Real >= 0.9.9, base >= 3, base < 4, containers, binary >= 0.4, html >= 1.0, QuickCheck >= 1.2, QuickCheck < 2, time, filepath, directory Exposed-modules: Data.Number.ER.RnToRm,- Data.Number.ER.RnToRm.BisectionTree.Path,- Data.Number.ER.RnToRm.BisectionTree.Integration,- Data.Number.ER.RnToRm.BisectionTree,- Data.Number.ER.RnToRm.DefaultRepr,- Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic,- Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Elementary,- Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Field, Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Integration,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Compose, Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Eval, Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Ring,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Enclosure,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Elementary,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Reduce,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Division,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Compose,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Bounds,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Ring,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Enclosure,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Elementary,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Generate,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Reduce,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Division,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Run, Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom, Data.Number.ER.RnToRm.UnitDom.Base, Data.Number.ER.RnToRm.UnitDom.Approx.Interval, Data.Number.ER.RnToRm.UnitDom.Approx,+ Data.Number.ER.RnToRm.TestingDefs,+ Data.Number.ER.RnToRm.DefaultRepr,+ Data.Number.ER.RnToRm.BisectionTree.Integration,+ Data.Number.ER.RnToRm.BisectionTree.Path,+ Data.Number.ER.RnToRm.BisectionTree,+ Data.Number.ER.RnToRm.Approx.DomEdges, Data.Number.ER.RnToRm.Approx.DomTransl, Data.Number.ER.RnToRm.Approx.PieceWise,- Data.Number.ER.RnToRm.Approx.DomEdges, Data.Number.ER.RnToRm.Approx.Tuple,- Data.Number.ER.RnToRm.Approx,- Data.Number.ER.RnToRm.TestingDefs + Data.Number.ER.RnToRm.Approx
ChangeLog view
@@ -1,5 +1,17 @@+0.4.9:+ * Added a quickcheck testing harness for the polynomial arithmetic core.+ * Rewritten polynomial arithmetic core.+ * Fixed many rounding errors affecting almost all operations.+ * New operation: substitution into an enclosure of a *monotone* function.+ * In enclosure arithmetic, now can set a limit on the size of each enclosure representation.+ This is important for many-variate polynomials that tend to have very many terms.+ +0.4.3:+ * fixed two serious errors in exponentiation of PFEs+ * added composition of a monotone function approx with another function approx+ and implemented it for PFEs on individual domain boxes 0.4.2: 1 December 2008- * substantially improved division by a constant PFE+ * substantially improved division by a constant PFE (polynomial function enclosure) * added proper handling of overflown coefficients 0.4.1: 30 September 2008 * updated to work with AERN-Real 0.9.7
src/Data/Number/ER/RnToRm/Approx.hs view
@@ -18,7 +18,9 @@ ( ERFnApprox(..), ERFnDomApprox(..),- bisectUnbisectDepth+ bisectUnbisectDepth,+ keyPointsConsistencyCheck,+ keyPointsPointwiseConsistencyCheck ) where @@ -29,6 +31,8 @@ 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 {-|@@ -48,7 +52,7 @@ parts of the function's domain. -} class - (RA.ERApprox fa, RA.ERIntApprox domra, RA.ERIntApprox ranra, + (RA.ERApprox fa, RA.ERIntApprox fa, RA.ERIntApprox domra, RA.ERIntApprox ranra, DomainBox box varid domra) => ERFnApprox box varid domra ranra fa | fa -> box varid domra ranra@@ -73,6 +77,9 @@ This reduces the degree immediately if necessary and also affects all operations performed with this value later.++ May also set the maximum size of the approximations to a default+ based on the degree and the dimension of this enclosure. -} setMaxDegree :: Int -> fa -> fa {-| @@ -81,6 +88,27 @@ -} getMaxDegree :: fa -> Int {-| + Get the internal size of the approximation + (usually number of polynomial terms). + -}+ getSize :: fa -> Int+ {-| + Set an upper bound on the size of this function approximation.+ + This reduces the size immediately if necessary and also+ affects all operations performed with this value later.+ -}+ setMaxSize :: Int -> fa -> fa+ {-| + Get the current uppend bound on the size associated + with this function approximation. + -}+ getMaxSize :: fa -> Int+ {-| + List all variables that are actually used in the approximation.+ -}+ getVariables :: fa -> [varid]+ {-| Give a close upper bound of the precision of the range at the best approximated point in the domain. -}@@ -115,7 +143,7 @@ -} scale :: ranra -> fa -> fa {-|- Intersect one enclosure by another but only on a box within its domain.+ Intersect one approximation by another but only on a box within its domain. -} partialIntersect :: EffortIndex -> @@ -144,22 +172,34 @@ Fix some variables in the function to the given exact values. -} partialEval :: box -> fa -> fa- {-| - A simple and limited composition of functions.- - It is primarily intended to be used for precomposition with affine functions.+ {-|+ A simple and limited composition of functions applicable+ only when the range-defining function is non-decreasing. -} - composeThin ::- fa {-^ enclosure of @f@ -} ->- Map.Map varid fa- {-^ specifies the variables to substitute and for each such variable @v@, - gives an /exact/ enclosure of a function @f_v@ to substitute for @v@ -} ->- fa - {-^ enclosure of @f[v |-> f_v]@ + composeNonDecreasing ::+ fa {-^ enclosure of @f@, @f@ is non-decreasing in variable @var@ -} ->+ varid {-^ variable @var@ to get substituted in @f@ -} ->+ fa {-^ enclosure of @f_var@, to be substituted for @var@ -} -> + fa+ {-^ enclosure of @f[var |-> f_var]@ - BEWARE: Enclosure is probably incorrect where values of @f_v@ are outside the domain of @v@ in @f@.+ BEWARE: Enclosure is probably incorrect where values of + @f_v@ are outside the domain of @v@ in @f@. -}-+ {-|+ A simple and limited composition of functions applicable+ only when the range-defining function is non-increasing. + -} + composeNonIncreasing ::+ fa {-^ enclosure of @f@, @f@ is non-increasing in variable @var@ -} ->+ varid {-^ variable @var@ to get substituted in @f@ -} ->+ fa {-^ enclosure of @f_var@, to be substituted for @var@ -} -> + fa+ {-^ enclosure of @f[var |-> f_var]@ + + BEWARE: Enclosure is probably incorrect where values of + @f_v@ are outside the domain of @v@ in @f@.+ -} {-| This class extends 'ERFnApprox' by:@@ -301,3 +341,60 @@ fRDone = aux restVars depthsToGoM1 fR (fL, fR) = bisect var Nothing f depthsToGoM1 = depthsToGo - 1++{-|+ Check that a pointwise operation previously performed on function approximations is + consistent with the same operation performed on+ selected points in the domain of these functions.+ The selected points are the centres of all faces of all dimensions,+ which includes the corners.+ + The result of this function is the list of points in which + the consistency check failed. The result of the operation+ is also included both for the real number version and the+ function approximation version.+-} +keyPointsPointwiseConsistencyCheck ::+ (ERFnDomApprox box varid domra ranra fa) =>+ ([ranra] -> ranra) {-^ function @G@ acting on real numbers -} ->+ [fa] {-^ approximations of input functions -} ->+ fa {-^ alleged approximation of @G@ applied pointwise to the above function approximations -} ->+ [(box, ranra, ranra)]+keyPointsPointwiseConsistencyCheck g fIns fRes =+ keyPointsConsistencyCheck gPointwise fRes+ where+ gPointwise ptB =+ g $ map ((\[x] -> x) . eval ptB) fIns+ +{-|+ Check that a function approximation is consistent with+ a real function that is meant to compute the same function.+ + The result of this function is the list of points in which + the consistency check failed. The result of the operation+ is also included both for the real number version and the+ function approximation version.+-} +keyPointsConsistencyCheck ::+ (ERFnDomApprox box varid domra ranra fa) =>+ (box -> ranra) {-^ function @G@ acting on tuples of real numbers -} ->+ fa {-^ alleged approximation of @G@ over a domain box -} ->+ [(box, ranra, ranra)]+keyPointsConsistencyCheck g fRes =+ filter (isInConsistent) $ map testPoint points+ where+ points = map DBox.fromList $ allCombinations $ map getVarPoints varDoms+ varDoms = DBox.toList $ dom fRes+ getVarPoints (var, dom) =+ (var, [domL, domM, domR])+ where+ (domL, domR) = RA.bounds dom+ (domM, _) = RA.bounds $ (domL + domR)/2+ testPoint ptB =+ (ptB, gResPt, fResPt)+ where+ gResPt = g ptB+ [fResPt] = eval ptB fRes+ isInConsistent (_, gResPt, fResPt) =+ RA.isDisjoint gResPt fResPt+
src/Data/Number/ER/RnToRm/Approx/DomTransl.hs view
@@ -25,7 +25,6 @@ ERFnDomTranslApprox(..), DomTransl(..) ) where- import qualified Data.Number.ER.RnToRm.Approx as FA import qualified Data.Number.ER.RnToRm.UnitDom.Approx as UFA import qualified Data.Number.ER.Real.Approx as RA@@ -35,6 +34,8 @@ import Data.Number.ER.BasicTypes import Data.Number.ER.Misc +import Data.Number.ER.RnToRm.UnitDom.Approx.Interval+ import qualified Text.Html as H import Data.Typeable@@ -218,8 +219,8 @@ translateUfaToDom ufa dtrB -- gr = 20 + (RA.getGranularity ufa) -translateUfaToDom ufa dtrB =- FA.composeThin ufa $ +translateUfaToDom ufa dtrB = -- this is unsafe, use only for printing!+ UFA.composeWithThin ufa $ Map.fromAscList $ map mkToUnitUFA $ DBox.toAscList dtrB@@ -289,7 +290,8 @@ Fractional (ERFnDomTranslApprox dtrbox varid ufa domra) where fromRational r = ERFnDomTranslApprox (fromRational r) DBox.noinfo- recip (ERFnDomTranslApprox ufa dtrB) =+ recip f@(ERFnDomTranslApprox ufa dtrB) =+-- unsafePrintReturn ("DomTransl: recip of " ++ show f ++ "\n = ") $ ERFnDomTranslApprox (recip ufa) dtrB instance @@ -361,7 +363,8 @@ where abs ix (ERFnDomTranslApprox ufa dtrB) = ERFnDomTranslApprox (RAEL.abs ix ufa) dtrB- exp ix (ERFnDomTranslApprox ufa dtrB) =+ exp ix f@(ERFnDomTranslApprox ufa dtrB) =+-- unsafePrintReturn ("DomTransl: exp of " ++ show f ++ "\n = ") $ ERFnDomTranslApprox (RAEL.exp ix ufa) dtrB log ix (ERFnDomTranslApprox ufa dtrB) = ERFnDomTranslApprox (RAEL.log ix ufa) dtrB@@ -376,6 +379,7 @@ (UFA.ERUnitFnApprox box varid domra ranra ufa, DomainBoxMappable dtrbox box varid (DomTransl domra) domra, DomainIntBox box varid domra, + Show varid, Show box, DomainBoxMappable box dtrbox varid domra (DomTransl domra), Eq dtrbox, Ord dtrbox) => FA.ERFnApprox box varid domra ranra (ERFnDomTranslApprox dtrbox varid ufa domra)@@ -386,10 +390,14 @@ FA.domra2ranra (erfnUnitApprox fa) d ranra2domra fa r = FA.ranra2domra (erfnUnitApprox fa) r- setMaxDegree maxDegree (ERFnDomTranslApprox ufa dtrB) =- ERFnDomTranslApprox (FA.setMaxDegree maxDegree ufa) dtrB getMaxDegree (ERFnDomTranslApprox ufa _) = FA.getMaxDegree ufa+ setMaxDegree maxDegree (ERFnDomTranslApprox ufa dtrB) =+ ERFnDomTranslApprox (FA.setMaxDegree maxDegree ufa) dtrB+ getMaxSize (ERFnDomTranslApprox ufa _) =+ FA.getMaxSize ufa+ setMaxSize maxSize (ERFnDomTranslApprox ufa dtrB) =+ ERFnDomTranslApprox (FA.setMaxSize maxSize ufa) dtrB getRangeApprox (ERFnDomTranslApprox ufa dtrB) = FA.getRangeApprox ufa volume (ERFnDomTranslApprox ufa dtrB) =@@ -415,7 +423,46 @@ where dtrBNoVars = DBox.difference dtrB substitutions- + composeNonDecreasing+ fOuter@(ERFnDomTranslApprox ufaOuter dtrBOuter)+ varid+ fInner@(ERFnDomTranslApprox ufaInner dtrBInner)+ =+-- unsafePrintReturn+-- (+-- "ER.RnToRm.DomTransl: composeNonDecreasing: "+-- ++ "\n fOuter = " ++ show fOuter+-- ++ "\n varid = " ++ show varid+-- ++ "\n fInner = " ++ show fInner+-- ++ "\n inconsistencies = " ++ show (FA.keyPointsConsistencyCheck resultReals result)+-- ++ "\n result = "+-- )+-- $+ result+ where+ resultReals ptB = -- this is only used for consistency checking...+ (\[x] -> x) $ FA.eval ptBOuter fOuter+ where+ ptBOuter =+ DBox.insert varid fInnerVal ptB+ fInnerVal =+ FA.ranra2domra fInner $+ (\[x] -> x) $ FA.eval ptB fInner+ result = ERFnDomTranslApprox ufaComp dtrBComp + dtrBComp = + DBox.union (DBox.delete varid dtrBOuter) dtrBInner+ ufaComp = + FA.composeNonDecreasing ufaOuter varid ufaInnerUnitDom+ ufaInnerUnitDom =+ UFA.const [dtrVarConst]+ ++ (FA.scale dtrVarSlope ufaInner)+ dtrVarSlope =+ FA.domra2ranra ufaInner $ dtrToUnitSlope dtrVar+ dtrVarConst =+ FA.domra2ranra ufaInner $ dtrToUnitConst dtrVar+ dtrVar =+ DBox.lookup "ER.RnToRm.DomTransl: composeNonDecreasing: " varid dtrBOuter --instance -- (UFA.ERUnitFnApprox box varid domra ranra ufa, @@ -439,6 +486,7 @@ instance (UFA.ERUnitFnApprox box varid domra ranra ufa, DomainIntBox box varid domra,+ Show varid, Show box, DomainBoxMappable dtrbox box varid (DomTransl domra) domra, DomainBoxMappable box dtrbox varid domra (DomTransl domra), Eq dtrbox, Ord dtrbox) =>@@ -496,8 +544,8 @@ errMsg = "ERFnDomTranslApprox: FA.bisect: var " ++ showVar var ++ " not in the domain of " ++ show f- ufaLeft = FA.composeThin ufa $ Map.singleton var toLeft - ufaRight = FA.composeThin ufa $ Map.singleton var toRight+ ufaLeft = UFA.composeWithThin ufa $ Map.singleton var toLeft + ufaRight = UFA.composeWithThin ufa $ Map.singleton var toRight dtrLeft = DBox.insert var (makeDomTransl domLeft) dtrB dtrRight = DBox.insert var (makeDomTransl domRight) dtrB domLeft = domL RA.\/ pt@@ -527,17 +575,34 @@ ptGr = RA.setMinGranularity gran $ FA.domra2ranra ufa pt integrate ix fD@(ERFnDomTranslApprox ufaD dtrBD) x integdomBox- origin (ERFnDomTranslApprox ufaInit dtrBInit) =+ origin fI@(ERFnDomTranslApprox ufaInit dtrBInit) =+-- unsafePrintReturn+-- (+-- "ER.RnToRm.DomTransl: integrate: "+-- ++ "\n fD = " ++ show fD+-- ++ "\n variable = " ++ show x+-- ++ "\n origin = " ++ show origin+-- ++ "\n fI = " ++ show fI+-- ++ "\n ufaD = " ++ show ufaD+-- ++ "\n ufaDadj = " ++ show ufaDadj+-- ++ "\n originAdj = " ++ show originAdj+-- ++ "\n ufaI = " ++ show ufaI+-- ++ "\n ufaI(originAdj) = " ++ show (FA.eval (DBox.singleton x originAdj) ufaI)+-- ++ "\n result = "+-- )+-- $ ERFnDomTranslApprox ufaI dtrBD where ufaI = UFA.integrate ix ufaDadj x - (dtrToUnit trX origin) + originAdj ufaInit ufaDadj = FA.scale (FA.domra2ranra ufaD $ dtrFromUnitSlope trX) $ ufaD+ originAdj = + dtrToUnit trX origin trX = DBox.findWithDefault err x dtrBD err =
src/Data/Number/ER/RnToRm/TestingDefs.hs view
@@ -31,9 +31,14 @@ fapd04X0 = (FA.proj (DBox.fromAscList [(0,0 RA.\/ 4)]) 0) :: (FAPD B) fapd13X0 = (FA.proj (DBox.fromAscList [(0,1 RA.\/ 3)]) 0) :: (FAPD B) fapd12X1 = (FA.proj (DBox.fromAscList [(1,1 RA.\/ 2)]) 1) :: (FAPD B)-fapdUX0 = (FA.proj (DBox.fromAscList [(0,(-1) RA.\/ 1)]) 0) :: (FAPD B)-fapdUX1 = (FA.proj (DBox.fromAscList [(1,(-1) RA.\/ 1)]) 1) :: (FAPD B)+fapdUX0 = FA.setMaxDegree 2 (FA.proj (DBox.fromAscList [(0,(-1) RA.\/ 1)]) 0) :: (FAPD B)+fapdUX1 = FA.setMaxDegree 2 (FA.proj (DBox.fromAscList [(1,(-1) RA.\/ 1)]) 1) :: (FAPD B)+fapdUX2 = FA.setMaxDegree 2 (FA.proj (DBox.fromAscList [(2,(-1) RA.\/ 1)]) 2) :: (FAPD B) +fapdT1 = (1 + fapdUX2) * (1 + fapdUX2)+fapdT2 = fapdUX0 * fapdUX1 +fapdT3 = FA.composeNonDecreasing fapdT1 2 fapdT2+ fapeConst1 = (FA.const DBox.noinfo [1]) :: (FAPE B) fapeConstU = (FA.const DBox.noinfo [(-1) RA.\/ 1]) :: (FAPE B) fapeConst01 = (FA.const DBox.noinfo [0 RA.\/ 1]) :: (FAPE B)@@ -70,16 +75,30 @@ testIntegrP = FA.integrateMeasureImprovement 1 (FA.setMaxDegree 0 fapwUConst13InitPt) 0 (DBox.unary $ 0 RA.\/ 0.5) 0 fapwUConst13InitPt ++jas1 =+ FA.integrate+ 0+ f+ 0+ DBox.noinfo+ 1+ 0++f =+ RAEL.exp 100 x+ x = -- FA.bisectUnbisectDepth 1 $- FA.setMaxDegree 4 - fapwUUX0+ FA.setMaxDegree 10+-- fapwUUX10+ fapd13X0 y = -- FA.bisectUnbisectDepth 1 $ FA.setMaxDegree 4 - fapwUUX1- +-- fapwUUX1+ fapd12X1 xLR = snd $ FA.bisect 0 Nothing $ fst $ FA.bisect 0 Nothing $ x
src/Data/Number/ER/RnToRm/UnitDom/Approx.hs view
@@ -17,15 +17,20 @@ -} module Data.Number.ER.RnToRm.UnitDom.Approx (- ERUnitFnApprox(..)+ ERUnitFnApprox(..),+ keyPointsConsistencyCheck,+ keyPointsPointwiseConsistencyCheck ) where -import Data.Number.ER.RnToRm.Approx+import qualified Data.Number.ER.Real.Approx as RA+import qualified Data.Number.ER.RnToRm.Approx as FA 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 {-|@@ -37,7 +42,7 @@ where the domain has to be known. -} -class (ERFnApprox box varid domra ranra fa) => +class (FA.ERFnApprox box varid domra ranra fa) => ERUnitFnApprox box varid domra ranra fa | fa -> box varid domra ranra where@@ -56,6 +61,21 @@ [ranra] {-^ values at 0 -} -> Map.Map varid ([ranra]) {-^ ascents of each base vector -} -> fa+ {-|+ A simple and limited composition of functions.+ + It is primarily intended to be used for precomposition with affine functions.+ -} + composeWithThin ::+ fa {-^ enclosure of @f@ -} ->+ Map.Map varid fa+ {-^ specifies the variables to substitute and for each such variable @v@, + gives an /exact/ enclosure of a function @f_v@ to substitute for @v@ -} ->+ fa + {-^ enclosure of @f[v |-> f_v]@ + + BEWARE: Enclosure is probably incorrect where values of @f_v@ are outside the domain of @v@ in @f@.+ -} {-| Find close upper and lower bounds of the volume of the entire enclosure. A negative volume means that the enclosure is certainly inconsistent.@@ -90,3 +110,62 @@ domra {-^ origin in terms of @x@; this has to be exact! -} -> fa {-^ values at origin -} -> fa++ +{-|+ Check that a pointwise operation previously performed on function approximations is + consistent with the same operation performed on+ selected points in the domain of these functions.+ The selected points are the centres of all faces of all dimensions,+ which includes the corners.+ + The result of this function is the list of points in which + the consistency check failed. The result of the operation+ is also included both for the real number version and the+ function approximation version.+-} +keyPointsPointwiseConsistencyCheck ::+ (ERUnitFnApprox box varid domra ranra fa) =>+ ([ranra] -> ranra) {-^ function @G@ acting on real numbers -} ->+ [fa] {-^ approximations of input functions -} ->+ fa {-^ alleged approximation of @G@ applied pointwise to the above function approximations -} ->+ [(box, ranra, ranra)]+keyPointsPointwiseConsistencyCheck g fIns fRes =+ keyPointsConsistencyCheck gPointwise fRes+ where+ gPointwise ptB =+ g $ map ((\[x] -> x) . FA.eval ptB) fIns+ +{-|+ Check that a function approximation is consistent with+ a real function that is meant to compute the same function.+ + The result of this function is the list of points in which + the consistency check failed. The result of the operation+ is also included both for the real number version and the+ function approximation version.+-} +keyPointsConsistencyCheck ::+ (ERUnitFnApprox box varid domra ranra fa) =>+ (box -> ranra) {-^ function @G@ acting on tuples of real numbers -} ->+ fa {-^ alleged approximation of @G@ over a domain box -} ->+ [(box, ranra, ranra)]+keyPointsConsistencyCheck g fRes =+ filter (isInConsistent) $ map testPoint points+ where+ points = map DBox.fromList $ allCombinations $ map getVarPoints varDoms+ varDoms = map (\v -> (v,unitInterval)) $ FA.getVariables fRes+ unitInterval = (-1) RA.\/ 1+ getVarPoints (var, dom) =+ (var, [domL, domM, domR])+ where+ (domL, domR) = RA.bounds dom+ (domM, _) = RA.bounds $ (domL + domR)/2+ testPoint ptB =+ (ptB, gResPt, fResPt)+ where+ gResPt = g ptB+ [fResPt] = FA.eval ptB fRes+ isInConsistent (_, gResPt, fResPt) =+ RA.isDisjoint gResPt fResPt+
src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs view
@@ -32,6 +32,7 @@ import qualified Data.Number.ER.RnToRm.Approx as FA import qualified Data.Number.ER.RnToRm.UnitDom.Approx as UFA import qualified Data.Number.ER.RnToRm.UnitDom.Base as UFB+import Data.Number.ER.RnToRm.UnitDom.Base ((+^),(-^),(*^),multiplyEncl) import qualified Data.Number.ER.Real.Approx as RA import qualified Data.Number.ER.Real.Approx.Elementary as RAEL @@ -67,8 +68,8 @@ | ERFnInterval {- erfnUpper :: fb, erfnLowerNeg :: fb,+ erfnUpper :: fb, erfnContext :: ERFnContext, erfnGlobal :: ra }@@ -89,24 +90,26 @@ ERFnContext { erfnMaxDegree :: Int,+ erfnMaxSize :: Int, erfnCoeffGranularity :: Granularity } deriving (Show, Typeable, Data) instance Binary ERFnContext where- put (ERFnContext a b) = put a >> put b- get = get >>= \a -> get >>= \b -> return (ERFnContext a b)+ put (ERFnContext a b c) = put a >> put b >> put c+ get = get >>= \a -> get >>= \b -> get >>= \c -> return (ERFnContext a b c) erfnContextDefault = ERFnContext { erfnMaxDegree = 2,+ erfnMaxSize = 20, erfnCoeffGranularity = 10 } -erfnContextUnify (ERFnContext dg1 gr1) (ERFnContext dg2 gr2) =- ERFnContext (max dg1 dg2) (max gr1 gr2)+erfnContextUnify (ERFnContext dg1 sz1 gr1) (ERFnContext dg2 sz2 gr2) =+ ERFnContext (max dg1 dg2) (max sz1 sz2) (max gr1 gr2) instance @@ -114,37 +117,41 @@ Show (ERFnInterval fb ra) where show (ERFnIntervalAny _) = "ERFnIntervalAny"- show (ERFnInterval h ln ctxt gl) =+ show (ERFnInterval ln h ctxt gl) = "\nERFnInterval"- ++ "\n upper = " ++ show h- ++ "\n lower = " ++ show (-ln)+ ++ "\n upper = " ++ ufbShow h+ ++ "\n lower = " ++ ufbShow (UFB.neg ln) -- ++ " global = " ++ show gl ++ "\n" -- ++ " context = " ++ show ctxt ++ "\n"+ where+ ufbShow = UFB.showDiGrCmp 10 False False instance (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>- H.HTML (ERFnInterval fb ra) + H.HTML (ERFnInterval fb ra) where toHtml (ERFnIntervalAny ctxt) = H.toHtml "ERFnIntervalAny"- toHtml (ERFnInterval h ln ctxt gl) =+ toHtml (ERFnInterval ln h ctxt gl) = -- H.toHtml $ -- abovesTable -- [ -- H.toHtml "ERFnInterval", H.toHtml $ H.simpleTable [H.border 2] [] [- [H.toHtml "upper = ", H.toHtml $ show h],- [H.toHtml "lower = ", H.toHtml $ show (- ln)]+ [H.toHtml "upper = ", H.toHtml $ ufbShow h],+ [H.toHtml "lower = ", H.toHtml $ ufbShow (UFB.neg ln)] ] -- ]+ where+ ufbShow = UFB.showDiGrCmp 10 False False instance (UFB.ERUnitFnBase boxb boxra varid b ra fb) => Eq (ERFnInterval fb ra) where- (ERFnInterval h1 ln1 ctxt1 gl1) - == (ERFnInterval h2 ln2 ctxt2 gl2) =+ (ERFnInterval ln1 h1 ctxt1 gl1) + == (ERFnInterval ln2 h2 ctxt2 gl2) = error "ERFnInterval: equality not implemented" _ == _ = error "ERFnInterval: equality not implemented"@@ -154,144 +161,134 @@ Ord (ERFnInterval fb ra) where compare - (ERFnInterval h1 ln1 ctxt1 gl1) - (ERFnInterval h2 ln2 ctxt2 gl2) =+ (ERFnInterval ln1 h1 ctxt1 gl1) + (ERFnInterval ln2 h2 ctxt2 gl2) = error "ERFnInterval: comparison not implemented; consider leqReals or compareApprox from class ERApprox instead" compare _ _ = error "ERFnInterval: comparison not implemented; consider leqReals or compareApprox from class ERApprox instead" instance - (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>+ (UFB.ERUnitFnBase boxb boxra varid b ra fb, Show varid, Show boxra) => Num (ERFnInterval fb ra) where fromInteger n = UFA.const [fromInteger n] negate f@(ERFnIntervalAny _) = f- negate (ERFnInterval h ln ctxt gl) =- (ERFnInterval ln h ctxt (negate gl))- (ERFnInterval h1 ln1 ctxt1 gl1) + (ERFnInterval h2 ln2 ctxt2 gl2) =+ negate (ERFnInterval ln h ctxt gl) =+ (ERFnInterval h ln ctxt (negate gl))+ (ERFnInterval ln1 h1 ctxt1 gl1) + (ERFnInterval ln2 h2 ctxt2 gl2) = normalise $- ERFnInterval (h1 + h2) (ln1 + ln2) ctxt (gl1 + gl2)+ ERFnInterval (reduceSzUp ln) (reduceSzUp h) ctxt (gl1 + gl2) where+ ln = ln1 +^ ln2+ h = h1 +^ h2+ reduceSzUp = UFB.reduceSizeUp maxSize+ maxSize = erfnMaxSize ctxt ctxt = erfnContextUnify ctxt1 ctxt2 f1 + f2 = ERFnIntervalAny ctxt where ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2)- (ERFnInterval h1 ln1 ctxt1 gl1) * (ERFnInterval h2 ln2 ctxt2 gl2) =+ (ERFnInterval ln1 h1 ctxt1 gl1) * (ERFnInterval ln2 h2 ctxt2 gl2) = normalise $- ERFnInterval h ln ctxt (gl1 * gl2)+ ERFnInterval ln h ctxt (gl1 * gl2) where- (h, ln) =- case (RA.leqReals 0 gl1, RA.leqReals gl1 0, RA.leqReals 0 gl2, RA.leqReals gl2 0) of- (Just True, _, Just True, _) -> -- both non-negative- (h1h2, l1l2Neg)- (_, Just True, _, Just True) -> -- both non-positive- (l1l2, h1h2Neg)- (Just True, _, _, Just True) -> -- first non-negative, second non-positive- (l1h2, h1l2Neg)- (_, Just True, Just True, _) -> -- first non-positive, second non-negative- (h1l2, l1h2Neg)- _ -> -- one of both may be crossing zero- ((h1h2 `maxP` l1l2) `maxP` (h1l2 `maxP` l1h2),- (h1h2Neg `maxP` l1l2Neg) `maxP` (h1l2Neg `maxP` l1h2Neg))- where- h1h2 = UFB.reduceDegreeUp maxDegr $ h1 * h2- h1h2Neg = UFB.reduceDegreeUp maxDegr $ (negate h1) * h2- l1l2 = UFB.reduceDegreeUp maxDegr $ ln1 * ln2- l1l2Neg = UFB.reduceDegreeUp maxDegr $ (negate ln1) * ln2- h1l2 = UFB.reduceDegreeUp maxDegr $ h1 * (negate ln2)- h1l2Neg = UFB.reduceDegreeUp maxDegr $ h1 * ln2- l1h2 = UFB.reduceDegreeUp maxDegr $ (negate ln1) * h2- l1h2Neg = UFB.reduceDegreeUp maxDegr $ ln1 * h2- maxP p1 p2 = fst $ UFB.max maxDegr p1 p2- - ctxt = erfnContextUnify ctxt1 ctxt2+ (ln, h) = multiplyEncl maxDegr maxSize (ln1, h1) (ln2, h2) maxDegr = erfnMaxDegree ctxt+ maxSize = erfnMaxSize ctxt+ ctxt = erfnContextUnify ctxt1 ctxt2 f1 * f2 = ERFnIntervalAny ctxt where ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2) instance - (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>- Fractional (ERFnInterval fb ra) + (UFB.ERUnitFnBase boxb boxra varid b ra fb, Show varid, Show boxra) =>+ Fractional (ERFnInterval fb ra) where fromRational r = UFA.const [fromRational r] recip f@(ERFnIntervalAny _) = f- recip (ERFnInterval h ln ctxt gl) + recip (ERFnInterval ln h ctxt gl) | certainNoZero = normalise $- ERFnInterval lRecipUp hnRecipUp ctxt (recip gl)+ ERFnInterval lnR hR ctxt (recip gl) | otherwise = ERFnIntervalAny ctxt where+ (hR, lnR) = UFB.recipEncl maxDegr maxSize ix (h,ln) certainNoZero = certainAboveZero || certainBelowZero certainAboveZero = UFB.upperBound ix ln < 0 certainBelowZero = UFB.upperBound ix h < 0 - hnRecipUp =- UFB.recipUp maxDegr ix (negate h) - lRecipUp =- UFB.recipUp maxDegr ix (negate ln)+-- hnRecipUp =+-- UFB.recipUp maxDegr maxSize ix (negate h) +-- lRecipUp =+-- UFB.recipUp maxDegr maxSize ix (negate ln) maxDegr = erfnMaxDegree ctxt+ maxSize = erfnMaxSize ctxt ix = int2effIx $ 3 * maxDegr instance- (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>+ (UFB.ERUnitFnBase boxb boxra varid b ra fb, Show varid, Show boxra) => RA.ERApprox (ERFnInterval fb ra) where initialiseBaseArithmetic _ =- UFB.initialiseBaseArithmetic (0 :: fb)+ UFB.initialiseBaseArithmetic (UFB.const 0 :: fb) getGranularity (ERFnIntervalAny ctxt) = erfnCoeffGranularity ctxt- getGranularity (ERFnInterval h ln ctxt gl) =+ getGranularity (ERFnInterval ln h ctxt gl) = max (erfnCoeffGranularity ctxt) $ - max (UFB.getGranularity h) (UFB.getGranularity ln)+ max (UFB.getGranularity ln) (UFB.getGranularity h) setGranularity gran (ERFnIntervalAny ctxt) = ERFnIntervalAny $ ctxt { erfnCoeffGranularity = gran }- setGranularity gran (ERFnInterval h ln ctxt gl) =+ setGranularity gran (ERFnInterval ln h ctxt gl) = ERFnInterval - (UFB.setGranularity gran h) (UFB.setGranularity gran ln) + (UFB.setGranularity gran ln) (UFB.setGranularity gran h) (ctxt { erfnCoeffGranularity = gran }) gl setMinGranularity gran (ERFnIntervalAny ctxt) = ERFnIntervalAny (ctxt { erfnCoeffGranularity = max gran (erfnCoeffGranularity ctxt) })- setMinGranularity gran (ERFnInterval h ln ctxt gl) =+ setMinGranularity gran (ERFnInterval ln h ctxt gl) = ERFnInterval - (UFB.setMinGranularity gran h) (UFB.setMinGranularity gran ln) + (UFB.setMinGranularity gran ln) (UFB.setMinGranularity gran h) (ctxt { erfnCoeffGranularity = max gran (erfnCoeffGranularity ctxt) }) gl -- getPrecision (ERFnIntervalAny _) = 0 -- getPrecision f = intLog 2 (1 + (fst $ RA.integerBounds (FA.volume f))) -- wrong! - f1@(ERFnInterval h1 ln1 ctxt1 gl1) /\ f2@(ERFnInterval h2 ln2 ctxt2 gl2) =+ f1@(ERFnInterval ln1 h1 ctxt1 gl1) /\ f2@(ERFnInterval ln2 h2 ctxt2 gl2) = ---- #ifdef RUNTIME_CHECKS ---- check ("ERFnInterval: /\\:\n f1:\n" ++ show f1 ++ " f2:\n" ++ show f2 ++ "\n result:\n") $ ---- #endif normalise $- ERFnInterval (snd $ UFB.min maxDegr h1 h2) (snd $ UFB.min maxDegr ln1 ln2) ctxt (gl1 RA./\ gl2)+ ERFnInterval + (UFB.minUp maxDegr maxSize ln1 ln2) + (UFB.minUp maxDegr maxSize h1 h2) + ctxt (gl1 RA./\ gl2) where ctxt = erfnContextUnify ctxt1 ctxt2 maxDegr = erfnMaxDegree ctxt- (ERFnIntervalAny ctxt1) /\ (ERFnInterval h2 ln2 ctxt2 gl2) =- ERFnInterval h2 ln2 ctxt gl2+ maxSize = erfnMaxSize ctxt+ (ERFnIntervalAny ctxt1) /\ (ERFnInterval ln2 h2 ctxt2 gl2) =+ ERFnInterval ln2 h2 ctxt gl2 where ctxt = erfnContextUnify ctxt1 ctxt2- (ERFnInterval h1 ln1 ctxt1 gl1) /\ (ERFnIntervalAny ctxt2) =- ERFnInterval h1 ln1 ctxt gl1+ (ERFnInterval ln1 h1 ctxt1 gl1) /\ (ERFnIntervalAny ctxt2) =+ ERFnInterval ln1 h1 ctxt gl1 where ctxt = erfnContextUnify ctxt1 ctxt2 f1 /\ f2 = ERFnIntervalAny ctxt where ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2)- leqReals = erfnintLeq+ leqReals f1 f2 = +-- unsafePrint ("ERInterval: leqReals: sizes: " ++ show (FA.getSize f1) ++ ", " ++ show (FA.getSize f2)) $ + erfnintLeq f1 f2 refines _ (ERFnIntervalAny _) = True refines (ERFnIntervalAny _) _ = False- refines (ERFnInterval h1 ln1 _ _) (ERFnInterval h2 ln2 _ _) = - (UFB.upperBound 10 (h2 - h1) >= 0)+ refines (ERFnInterval ln1 h1 _ _) (ERFnInterval ln2 h2 _ _) = + (UFB.upperBound 10 (ln2 -^ ln1) >= 0) &&- (UFB.upperBound 10 (ln2 - ln1) >= 0)+ (UFB.upperBound 10 (h2 -^ h1) >= 0) compareApprox (ERFnIntervalAny _) (ERFnIntervalAny _) = EQ compareApprox (ERFnIntervalAny _) _ = LT compareApprox _ (ERFnIntervalAny _) = GT- compareApprox (ERFnInterval h1 ln1 _ _) (ERFnInterval h2 ln2 _ _) =+ compareApprox (ERFnInterval ln1 h1 _ _) (ERFnInterval ln2 h2 _ _) = compareComposeMany [ UFB.compareApprox h1 h2,@@ -306,16 +303,16 @@ isClearlyBelow (ERFnIntervalAny _) _ = False isClearlyBelow _ (ERFnIntervalAny _) = False isClearlyBelow f g- | UFB.upperBound 10 (erfnUpper f + erfnLowerNeg g) <= 0 = True+ | UFB.upperBound 10 (erfnUpper f +^ erfnLowerNeg g) <= 0 = True | otherwise = False isClearlyStrictlyBelow (ERFnIntervalAny _) _ = False isClearlyStrictlyBelow _ (ERFnIntervalAny _) = False isClearlyStrictlyBelow f g- | UFB.upperBound 10 (erfnUpper f + erfnLowerNeg g) < 0 = True+ | UFB.upperBound 10 (erfnUpper f +^ erfnLowerNeg g) < 0 = True | otherwise = False instance- (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>+ (UFB.ERUnitFnBase boxb boxra varid b ra fb, Show varid, Show boxra) => RA.ERIntApprox (ERFnInterval fb ra) where -- doubleBounds = :: ira -> (Double, Double) @@ -323,37 +320,38 @@ -- integerBounds :: ira -> (ExtendedInteger, ExtendedInteger) bisectDomain maybePt (ERFnIntervalAny c) = error "ERFnInterval: RA.bisectDomain: cannot bisect ERFnIntervalAny"- bisectDomain maybePt (ERFnInterval u ln c g) =- (ERFnInterval midUp ln c g,- ERFnInterval u (negate midDown) c g)+ bisectDomain maybePt (ERFnInterval ln h c g) =+ (ERFnInterval ln midUp c g,+ ERFnInterval midDownNeg h c g) where- (midDown, midUp) =+ (midDownNeg, midUp) = case maybePt of Nothing ->- (negate $ (ln - u) / 2, (u - ln) / 2)- Just (ERFnInterval uPt lnPt _ _) ->- (negate lnPt, uPt)+ (UFB.scaleUp (1/2) $ ln -^ h, UFB.scaleUp (1/2) $ h -^ ln)+ Just (ERFnInterval lnPt hPt _ _) ->+ (lnPt, hPt) bounds (ERFnIntervalAny c) = error "ERFnInterval: RA.bounds: cannot get bounds for ERFnIntervalAny"- bounds (ERFnInterval u ln c g) =- (ERFnInterval (negate ln) ln c g,- ERFnInterval u (negate u) c g) - f1@(ERFnInterval u1 ln1 c1 g1) \/ f2@(ERFnInterval u2 ln2 c2 g2) =+ bounds (ERFnInterval ln h c g) =+ (ERFnInterval ln (UFB.neg ln) c g,+ ERFnInterval (UFB.neg h) h c g) + f1@(ERFnInterval ln1 h1 c1 g1) \/ f2@(ERFnInterval ln2 h2 c2 g2) = ---- #ifdef RUNTIME_CHECKS ---- check ("ERFnInterval: abs:\n f1:\n" ++ show f1 ++ " f2:\n" ++ show f2 ++ "\n result:\n") $ ---- #endif normalise $- ERFnInterval u ln c (g1 RA.\/ g2)+ ERFnInterval ln h c (g1 RA.\/ g2) where- u = UFB.maxUp maxDegree u1 u2- ln = UFB.maxUp maxDegree ln1 ln2+ ln = UFB.maxUp maxDegree maxSize ln1 ln2+ h = UFB.maxUp maxDegree maxSize h1 h2 c = erfnContextUnify c1 c2 maxDegree = erfnMaxDegree c- (ERFnIntervalAny ctxt1) \/ (ERFnInterval h2 ln2 ctxt2 gl2) =+ maxSize = erfnMaxSize c+ (ERFnIntervalAny ctxt1) \/ (ERFnInterval ln2 h2 ctxt2 gl2) = ERFnIntervalAny ctxt where ctxt = erfnContextUnify ctxt1 ctxt2- (ERFnInterval h1 ln1 ctxt1 gl1) \/ (ERFnIntervalAny ctxt2) =+ (ERFnInterval ln1 h1 ctxt1 gl1) \/ (ERFnIntervalAny ctxt2) = ERFnIntervalAny ctxt where ctxt = erfnContextUnify ctxt1 ctxt2@@ -362,70 +360,88 @@ ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2) instance- (UFB.ERUnitFnBase boxb boxra varid b ra fb, RAEL.ERApproxElementary ra, RealFrac b) =>+ (UFB.ERUnitFnBase boxb boxra varid b ra fb, + RAEL.ERApproxElementary ra, RealFrac b, + Show varid, Show boxra) => RAEL.ERApproxElementary (ERFnInterval fb ra) where -- default abs does not work because we do not have Prelude.abs abs _ f@(ERFnIntervalAny _) = f- abs _ f@(ERFnInterval u ln c g) =+ abs _ f@(ERFnInterval ln h c g) = ---- #ifdef RUNTIME_CHECKS ---- check ("ERFnInterval: abs:\n f:\n" ++ show f ++ "\n result:\n") $ ---- #endif normalise $- ERFnInterval maxulnUp maxunl0Dn c (abs g)+ ERFnInterval minhln0Up maxhlnUp c (abs g) where+ maxhlnUp = UFB.maxUp maxDegree maxSize h ln + minhln0Up =+ UFB.minUp maxDegree maxSize (UFB.const 0) $+ UFB.minUp maxDegree maxSize h ln maxDegree = erfnMaxDegree c- maxulnUp = snd $ UFB.max maxDegree u ln - maxunl0Dn =- fst $ UFB.max maxDegree 0 $- fst $ UFB.max maxDegree (- u) (- ln)+ maxSize = erfnMaxSize c exp ix f@(ERFnIntervalAny _) = f- exp ix f@(ERFnInterval u ln c g) = + exp ix f@(ERFnInterval ln h c g) = normalise $- ERFnInterval uExp lExpNeg c (RAEL.exp ix g)+ ERFnInterval lExpNeg hExp c (RAEL.exp ix g) where maxDegree = erfnMaxDegree c- uExp = snd $ UFB.exp maxDegree ix u- lExpNeg = - negate $ fst $ UFB.exp maxDegree ix (negate ln) + maxSize = erfnMaxSize c+ (lExpNeg, hExp) =+ case (UFB.upperBound ix (h +^ ln) <= 1) of+ True -> + UFB.expEncl maxDegree maxSize ix (ln, h)+ False ->+ (lExpNeg, hExp)+ where+ (lExpNeg, _) = UFB.expEncl maxDegree maxSize ix (ln, UFB.neg ln)+ (_, hExp) = UFB.expEncl maxDegree maxSize ix (UFB.neg h,h) sin ix f@(ERFnIntervalAny c) = - ERFnInterval 1 1 c ((-1) RA.\/ 1)- sin ix f@(ERFnInterval u ln c g) =+ ERFnInterval one one c ((-1) RA.\/ 1)+ where+ one = UFB.const 1+ sin ix f@(ERFnInterval ln h c g) = -- unsafePrint -- ( -- "ERFnInterval: RAEL.sin: "--- ++ "\n u = " ++ show u+-- ++ "\n h = " ++ show h -- ++ "\n ln = " ++ show ln--- ++ "\n uSin = " ++ show uSin+-- ++ "\n hSin = " ++ show hSin -- ++ "\n lSinNeg = " ++ show lSinNeg -- ) $ ---- #ifdef RUNTIME_CHECKS ---- check ("ERFnInterval: sin:\n f:\n" ++ show f ++ "\n result:\n") $ ---- #endif normalise $- ERFnInterval uSin (- lSin) c (RAEL.sin ix g)+ ERFnInterval lSinNeg hSin c (RAEL.sin ix g) where- (lSin, uSin) = sincos True maxDegree ix u (-ln) + (lSinNeg, hSin) = sincos True maxDegree maxSize ix (ln, h) maxDegree = erfnMaxDegree c+ maxSize = erfnMaxSize c cos ix f@(ERFnIntervalAny c) =- ERFnInterval 1 1 c ((-1) RA.\/ 1)- cos ix f@(ERFnInterval u ln c g) =+ ERFnInterval one one c ((-1) RA.\/ 1)+ where+ one = UFB.const 1+ cos ix f@(ERFnInterval ln h c g) = -- unsafePrint -- ( -- "ERFnInterval: RAEL.cos: "--- ++ "\n u = " ++ show u+-- ++ "\n h = " ++ show h -- ++ "\n ln = " ++ show ln -- ++ "\n uCos = " ++ show uCos -- ++ "\n lCosNeg = " ++ show lCosNeg -- ) $ normalise $- ERFnInterval uCos (- lCos) c (RAEL.cos ix g)+ ERFnInterval lCosNeg hCos c (RAEL.cos ix g) where- (lCos, uCos) = sincos False maxDegree ix u (-ln) + (lCosNeg, hCos) = sincos False maxDegree maxSize ix (ln,h) maxDegree = erfnMaxDegree c+ maxSize = erfnMaxSize c atan ix f@(ERFnIntervalAny c) =- ERFnInterval 1 1 c ((-1) RA.\/ 1)- atan ix f@(ERFnInterval u ln c g) =+ ERFnInterval one one c ((-1) RA.\/ 1)+ where+ one = UFB.const 1+ atan ix f@(ERFnInterval ln h c g) = -- unsafePrint -- ( -- "ERFnInterval: RAEL.atan: "@@ -435,23 +451,30 @@ -- ++ "\n lAtanNeg = " ++ show lAtanNeg -- ) $ normalise $- ERFnInterval uAtan lAtanNeg c (RAEL.atan ix g)+ ERFnInterval lAtanNeg hAtan c (RAEL.atan ix g) where maxDegree = erfnMaxDegree c+ maxSize = erfnMaxSize c -- ix = int2effIx maxDegree- uAtan = snd $ UFB.atan maxDegree ix u- lAtanNeg = - negate $ fst $ UFB.atan maxDegree ix (negate ln) + (lAtanNeg, hAtan) = + case (UFB.upperBound ix (h +^ ln) <= 1) of+ True ->+ UFB.atanEncl maxDegree maxSize ix (ln, h)+ False ->+ (lAtanNeg, hAtan)+ where+ (lAtanNeg, _) = UFB.atanEncl maxDegree maxSize ix (ln, UFB.neg ln)+ (_, hAtan) = UFB.atanEncl maxDegree maxSize ix (UFB.neg h,h) sincos :: (UFB.ERUnitFnBase boxb boxra varid b ra fb, RAEL.ERApproxElementary ra, RealFrac b) => Bool {-^ True iff sine, False iff cosine -} -> Int {-^ maximum representation degree -} -> + Int {-^ maximum approx size -} -> EffortIndex {-^ how hard to try to eliminate truncation errors -} -> - fb ->- fb ->+ (fb, fb) -> (fb, fb)-sincos isSine maxDegree ix u l+sincos isSine maxDegree maxSize ix (ln,h) -- p - 2k*pi range within [-pi/2, pi/2]: | ranfNear0 `RA.refines` plusMinusPiHalf = -- unsafePrint@@ -524,6 +547,7 @@ (UFB.const (-1), UFB.const 1) -- (expDownwards, expUpwards + valueAtRDnNeg + (UFB.const expRUp)) where+-- l = UFB.neg ln ranfNear0 = ranf - k2pi k2pi = k * 2 * pi plusMinusPiHalf = (-piHalfLO) RA.\/ piHalfLO@@ -532,11 +556,10 @@ (piHalfLO, piHalfHI) = RA.bounds piHalf ranf = ERInterval - (UFB.lowerBound 10 l) - (UFB.upperBound 10 u)- k = - fromInteger $ floor $ - case (pi + ranf) / (2 * pi) of ERInterval lo hi -> lo+ (negate $ UFB.upperBound 10 ln) + (UFB.upperBound 10 h)+ k = fromInteger $ toInteger kEI+ (kEI,_) = RA.integerBounds $ 0.5 + (ranf / (2*pi)) sineShiftedNegated shift = boundsNegate $ sineShifted shift@@ -544,120 +567,180 @@ cosineShiftedNegated shift = boundsNegate $ cosineShifted shift - boundsNegate (pLO, pHI) = (- pHI, - pLO)+ boundsNegate (pLONeg, pHI) = (pHI, pLONeg) - sineShifted shift =- boundsAddErr shiftWidthB (lSinDown, uSinUp)+ sineShifted shift = -- moving to domain where sinus is non-decreasing+ case (UFB.upperBound ix (h +^ ln) <= 0.25) of+ True -> + UFB.sinEncl maxDegree maxSize ix (lnShifted, hShifted)+ False ->+ (lSinNeg, hSin)+ where+ (lSinNeg, _) = UFB.sinEncl maxDegree maxSize ix (ln, UFB.neg ln)+ (_, hSin) = UFB.sinEncl maxDegree maxSize ix (UFB.neg h,h) where- lSinDown = fst $ UFB.sin maxDegree ix (l `plusUp` shiftPoly)- uSinUp = snd $ UFB.sin maxDegree ix (u `plusDown` shiftPoly) - shiftPoly = UFB.const shiftLOB+ lnShifted = ln +^ (UFB.const (- shiftLOB))+ hShifted = h +^ (UFB.const shiftHIB) ERInterval shiftLOB shiftHIB = shift- shiftWidthB = shiftHIB - shiftLOB++ - cosineShifted shift =- boundsAddErr shiftWidthB $ - (UFB.minDown maxDegree lCosDown uCosDown,- UFB.maxUp maxDegree lCosUp uCosUp - + (snd $ UFB.scale 0.5 (u-l))) -- important near 0+ cosineShifted shift = -- moving to domain where cosinus is non-decreasing+ case (UFB.upperBound ix (h +^ ln) <= 0.25) of+ True -> + UFB.cosEncl maxDegree maxSize ix (lnShifted, hShifted)+ False ->+ (UFB.minUp maxDegree maxSize lCosDownNeg hCosDownNeg,+ UFB.maxUp maxDegree maxSize lCosUp hCosUp + +^ (UFB.scaleUp 0.5 (h +^ ln))) + -- this term is important when enclosure hits 0;+ -- without it, the result could miss cosine's maximum at 0 where- (lCosDown, lCosUp) = UFB.cos maxDegree ix (l `plusUp` shiftPoly)- (uCosDown, uCosUp) = UFB.cos maxDegree ix (u `plusDown` shiftPoly) - shiftPoly = UFB.const shiftLOB+ (lCosDownNeg, lCosUp) = UFB.cosEncl maxDegree maxSize ix (ln, UFB.neg ln)+ (hCosDownNeg, hCosUp) = UFB.cosEncl maxDegree maxSize ix (UFB.neg h,h)+ lnShifted = ln +^ (UFB.const (- shiftLOB))+ hShifted = h +^ (UFB.const shiftHIB) ERInterval shiftLOB shiftHIB = shift- shiftWidthB = shiftHIB - shiftLOB - boundsAddErr errB (pLO, pHI) =- (pLO `plusDown` (- errPoly), pHI + errPoly)+ boundsAddErr errB (pLONeg, pHI) =+ (pLONeg +^ errPoly, pHI +^ errPoly) where errPoly = UFB.const errB normalise f@(ERFnIntervalAny c) = f-normalise f@(ERFnInterval u ln c g)- | UFB.isValid u && UFB.isValid ln = f+normalise f@(ERFnInterval ln h c g)+ | UFB.isValid h && UFB.isValid ln = f | otherwise = ERFnIntervalAny c check callerLocation f@(ERFnIntervalAny c) = f-check callerLocation f@(ERFnInterval u ln c g) =+check callerLocation f@(ERFnInterval ln h c g) = ERFnInterval - (UFB.check (callerLocation ++ "upper: ") u) + (UFB.check (callerLocation ++ "upper: ") h) (UFB.check (callerLocation ++ "neg lower: ") ln) c g instance - (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>+ (UFB.ERUnitFnBase boxb boxra varid b ra fb, Show varid, Show boxra) => FA.ERFnApprox boxra varid ra ra (ERFnInterval fb ra) where check = check domra2ranra _ = id ranra2domra _ = id+ getMaxDegree (ERFnIntervalAny c) =+ erfnMaxDegree c+ getMaxDegree (ERFnInterval _ _ c _) =+ erfnMaxDegree c setMaxDegree maxDegr (ERFnIntervalAny c) = ERFnIntervalAny (c { erfnMaxDegree = maxDegr } )- setMaxDegree maxDegr (ERFnInterval u ln c g) =+ setMaxDegree maxDegr (ERFnInterval ln h c g) = ERFnInterval - (UFB.reduceDegreeUp maxDegr u) (UFB.reduceDegreeUp maxDegr ln)+ (UFB.reduceDegreeUp maxDegr h) (c { erfnMaxDegree = maxDegr } ) g- getMaxDegree (ERFnIntervalAny c) =- erfnMaxDegree c- getMaxDegree (ERFnInterval _ _ c _) =- erfnMaxDegree c+ getSize (ERFnIntervalAny c) = 0+ getSize (ERFnInterval ln h c g) =+ max (UFB.getSize ln) (UFB.getSize h)+ getMaxSize (ERFnIntervalAny c) =+ erfnMaxSize c+ getMaxSize (ERFnInterval _ _ c _) =+ erfnMaxSize c+ setMaxSize maxSize (ERFnIntervalAny c) =+ ERFnIntervalAny (c { erfnMaxDegree = maxSize } )+ setMaxSize maxSize (ERFnInterval ln h c g) =+ ERFnInterval + (UFB.reduceSizeUp maxSize ln)+ (UFB.reduceSizeUp maxSize h)+ (c { erfnMaxSize = maxSize } )+ g+ getVariables (ERFnIntervalAny _) = []+ getVariables (ERFnInterval ln h _ _) = UFB.getVariables h getRangeApprox (ERFnIntervalAny _) = RA.bottomApprox - getRangeApprox (ERFnInterval u ln c g) =- UFB.raFromEndpoints u+ getRangeApprox (ERFnInterval ln h c g) =+ UFB.raFromEndpoints h ( (- (UFB.upperBound 10 ln)) ,- (UFB.upperBound 10 u)+ (UFB.upperBound 10 h) ) scale ratio f@(ERFnIntervalAny c) = f- scale ratio f@(ERFnInterval u ln c g) = + scale ratio f@(ERFnInterval ln h c g) = ---- #ifdef RUNTIME_CHECKS ---- FA.check ("ERFnInterval: scale:\n before:\n" ++ show f ++ "\n after:\n") $ ---- #endif normalise $ case RA.compareReals ratio 0 of Just GT -> - ERFnInterval (UFB.scaleApproxUp ratio u) (UFB.scaleApproxUp ratio ln) c g+ ERFnInterval (scaleUp ratio ln) (scaleUp ratio h) c g Just LT -> - ERFnInterval (UFB.scaleApproxUp (- ratio) ln) (UFB.scaleApproxUp (- ratio) u) c g+ ERFnInterval (scaleUp (- ratio) h) (scaleUp (- ratio) ln) c g _ -> (UFA.const [ratio]) * f+ where+ scaleUp = UFB.scaleApproxUp maxDegree maxSize+ maxDegree = erfnMaxDegree c+ maxSize = erfnMaxSize c eval ptBox (ERFnIntervalAny c) = [RA.bottomApprox]- eval ptBox (ERFnInterval u ln c g) =+ eval ptBox (ERFnInterval ln h c g) = [lo RA.\/ up] where- up = UFB.evalApprox ptBox u+ up = UFB.evalApprox ptBox h lo = negate $ UFB.evalApprox ptBox ln partialEval substitutions f@(ERFnIntervalAny c) = f- partialEval substitutions f@(ERFnInterval u ln c g) =+ partialEval substitutions f@(ERFnInterval ln h c g) = normalise $- (ERFnInterval uP lnP c g)+ (ERFnInterval lnP hP c g) where- uP = UFB.partialEvalApproxUp substitutions u+ hP = UFB.partialEvalApproxUp substitutions h lnP = UFB.partialEvalApproxUp substitutions ln-- composeThin- f@(ERFnIntervalAny ctxt)- substitutions =- f- composeThin- f@(ERFnInterval h1 ln1 ctxt1 gl1)- substitutions =- (ERFnInterval h ln ctxt1 gl1)+ composeNonDecreasing+ fOuter@(ERFnInterval lnOuter hOuter cOuter gOuter)+ varid+ fInner@(ERFnInterval lnInner hInner cInner gInner) =+-- unsafePrintReturn+-- (+-- "ER.RnToRm.UnitDom.Interval: composeNonDecreasing: "+-- ++ "\n fOuter = " ++ show fOuter+-- ++ "\n varid = " ++ show varid+-- ++ "\n fInner = " ++ show fInner+-- ++ "\n inconsistencies = " ++ show (UFA.keyPointsConsistencyCheck resultReals result)+-- ++ "\n result = "+-- )+-- $+ result where- h = UFB.composeUp maxDegree h1 ufbSubstitutions - ln = UFB.composeUp maxDegree ln1 ufbSubstitutions- ufbSubstitutions = Map.map erfnUpper substitutions- maxDegree = erfnMaxDegree ctxt1 --- ctxt = erfnContextUnify ctxt1 ctxt2+ resultReals ptB = -- this is only used for consistency checking...+ (\[x] -> x) $ FA.eval ptBOuter fOuter+ where+ ptBOuter =+ DBox.insert varid fInnerVal ptB+ fInnerVal =+ FA.ranra2domra fInner $+ (\[x] -> x) $ FA.eval ptB fInner+ + result = ERFnInterval ln h c gOuter+ h =+ erfnUpper $ + UFA.composeWithThin fOuter $+ Map.singleton varid+ (ERFnInterval (UFB.neg hInner) hInner cInner gInner)+ ln =+ erfnLowerNeg $+ UFA.composeWithThin fOuter $+ Map.singleton varid $+ (ERFnInterval lnInner (UFB.neg lnInner) cInner gInner)+ c = erfnContextUnify cOuter cInner+ + composeNonDecreasing fOuter varid fInner = + ERFnIntervalAny c+ where+ c = erfnContextUnify (erfnContext fOuter) (erfnContext fInner) instance - (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>+ (UFB.ERUnitFnBase boxb boxra varid b ra fb, Show varid, Show boxra) => UFA.ERUnitFnApprox boxra varid ra ra (ERFnInterval fb ra) where bottomApprox =@@ -670,8 +753,8 @@ normalise $ ERFnInterval {- erfnUpper = fbH, erfnLowerNeg = fbLNeg,+ erfnUpper = fbH, erfnContext = context, erfnGlobal = val }@@ -694,8 +777,8 @@ normalise $ ERFnInterval {- erfnUpper = fbH, erfnLowerNeg = fbLNeg,+ erfnUpper = fbH, erfnContext = context, erfnGlobal = UFB.raFromEndpoints fbH@@ -723,6 +806,42 @@ { erfnCoeffGranularity = coeffGranularity }+ composeWithThin+ f@(ERFnIntervalAny ctxt)+ substitutions =+ f+ composeWithThin+ f@(ERFnInterval ln1 h1 ctxt1 gl1)+ substitutions =+-- unsafePrintReturn+-- (+-- "ER.RnToRm.UnitDom.Interval: composeWithThin: "+-- ++ "\n f = " ++ show f+-- ++ "\n substitutions = " ++ show substitutions+-- ++ "\n inconsistencies = " ++ show (UFA.keyPointsConsistencyCheck resultReals result)+-- ++ "\n result = "+-- )+-- $+ result+ where+ resultReals ptB = -- this is only used for consistency checking...+ (\[x] -> x) $+ FA.eval ptBOuter f+ where+ ptBOuter =+ foldl insertVal ptB $ Map.toList substitutions+ insertVal ptB (varid, fInner) =+ DBox.insert varid (evalPtB fInner) ptB+ evalPtB fInner =+ FA.ranra2domra fInner $ (\[x] -> x) $ FA.eval ptB fInner+ + result = ERFnInterval ln h ctxt1 gl1 + ln = UFB.composeManyUp maxDegree maxSize ln1 ufbSubstitutions+ h = UFB.composeManyUp maxDegree maxSize h1 ufbSubstitutions + ufbSubstitutions = Map.map erfnUpper substitutions+ maxDegree = erfnMaxDegree ctxt1 + maxSize = erfnMaxSize ctxt1 +-- ctxt = erfnContextUnify ctxt1 ctxt2 intersectMeasureImprovement ix vars f1@(ERFnIntervalAny ctxt1) f2@(ERFnIntervalAny ctxt2) =@@ -731,19 +850,19 @@ ctxt = erfnContextUnify ctxt1 ctxt2 intersectMeasureImprovement ix vars f1@(ERFnIntervalAny ctxt1) - f2@(ERFnInterval h2 ln2 ctxt2 gl2) =- (ERFnInterval h2 ln2 ctxt gl2, 1 / 0)+ f2@(ERFnInterval ln2 h2 ctxt2 gl2) =+ (ERFnInterval ln2 h2 ctxt gl2, 1 / 0) where ctxt = erfnContextUnify ctxt1 ctxt2 intersectMeasureImprovement ix vars- f1@(ERFnInterval h1 ln1 ctxt1 gl1) + f1@(ERFnInterval ln1 h1 ctxt1 gl1) f2@(ERFnIntervalAny ctxt2) = - (ERFnInterval h1 ln1 ctxt gl1, 1)+ (ERFnInterval ln1 h1 ctxt gl1, 1) where ctxt = erfnContextUnify ctxt1 ctxt2 intersectMeasureImprovement ix vars- f1@(ERFnInterval h1 ln1 ctxt1 gl1) - f2@(ERFnInterval h2 ln2 ctxt2 gl2) =+ f1@(ERFnInterval ln1 h1 ctxt1 gl1) + f2@(ERFnInterval ln2 h2 ctxt2 gl2) = case RA.compareReals improvementRA 1 of Just LT -> (f1, 1) -- intersection made it worse, keep original _ -> (intersection, improvementRA)@@ -765,16 +884,18 @@ f1Volume = UFA.volume vars f1 ctxt = erfnContextUnify ctxt1 ctxt2 volume vars (ERFnIntervalAny c) = 1/0- volume vars (ERFnInterval u ln c g) =--- unsafePrint ("ERFnInterval: volume: result = " ++ show result) $ result--- where--- result =- UFB.raFromEndpoints u $ UFB.volumeAboveZero vars (u + ln)+ volume vars (ERFnInterval ln h c g) =+ UFB.raFromEndpoints h (volL, volH)+ where + volH = UFB.volumeAboveZeroUp vars (ln +^ h)+ volL = negate $ UFB.volumeAboveZeroUp vars (l +^ hn)+ l = UFB.neg ln+ hn = UFB.neg h integrate _ f@(ERFnIntervalAny c) _ _ _ = f integrate - ix fD@(ERFnInterval u ln c g) x - origin fI@(ERFnInterval uInit lnInit cInit gInit) =--- unsafePrint+ ix fD@(ERFnInterval ln h c g) x + origin fI@(ERFnInterval lnInit hInit cInit gInit) =+-- unsafePrintReturn -- ( -- "ERFnInterval: integrate: " -- ++ "\n u = " ++ show u@@ -792,35 +913,37 @@ -- ++ "\n lnIuOriginU = " ++ show lnIuOriginU -- ++ "\n uIov = " ++ show uIov -- ++ "\n lnIov = " ++ show lnIov+-- ++ "\n result = " -- )+-- $ ---- #ifdef RUNTIME_CHECKS ---- check ("ERFnInterval: integrate:\n fD:\n" ++ show fD ++ "\n fI:\n" ++ show fI ++ "\n result:\n") $ ---- #endif normalise $- (ERFnInterval uIov lnIov c gIov)+ (ERFnInterval lnIov hIov c gIov) where -- perform raw integration of both bounds:- (uIuL, uIuU) = + (hIuL, hIuH) = -- mapPair (UFB.reduceDegreeDown maxDegree, UFB.reduceDegreeUp maxDegree) $ - UFB.integrate x u- (lnIuL, lnIuU) = + UFB.integrate x h+ (lnIuL, lnIuH) = -- mapPair (UFB.reduceDegreeDown maxDegree, UFB.reduceDegreeUp maxDegree) $ UFB.integrate x ln maxDegree = erfnMaxDegree c+ maxSize = erfnMaxSize c -- constrain the raw integrals to the origin:- uIuOriginL = UFB.composeDown maxDegree uIuL substXOrigin- uIuOriginU = UFB.composeUp maxDegree uIuU substXOrigin- lnIuOriginL = UFB.composeDown maxDegree lnIuL substXOrigin- lnIuOriginU = UFB.composeUp maxDegree lnIuU substXOrigin- substXOrigin = Map.singleton x originUFB- originUFB = UFB.const $ fst $ UFB.raEndpoints u origin- -- adjust the raw integrated functions enclose the initial condition function: - uIov = - UFB.reduceDegreeUp maxDegree $- uIuU + uInit - uIuOriginL + (uIuOriginU - uIuOriginL)+ (hIuOriginLNeg, hIuOriginH) =+ UFB.composeEncl maxDegree maxSize hIuL x originEncl+ (lnIuOriginLNeg, lnIuOriginH) = + UFB.composeEncl maxDegree maxSize lnIuL x originEncl+ originEncl = UFB.constEncl $ UFB.raEndpoints h origin+ -- adjust the raw integrated functions to enclose the initial condition function: + hIov = + UFB.reduceSizeUp maxSize $+ hIuH +^ hInit +^ hIuOriginLNeg +^ (hIuOriginH +^ hIuOriginLNeg) lnIov = - UFB.reduceDegreeUp maxDegree $- lnIuU + lnInit - lnIuOriginL + (lnIuOriginU - lnIuOriginL)+ UFB.reduceSizeUp maxSize $+ lnIuH +^ lnInit +^ lnIuOriginLNeg +^ (lnIuOriginH +^ lnIuOriginLNeg) gIov = gInit + g * ((1 - origin) RA.\/ (-1 - origin))
src/Data/Number/ER/RnToRm/UnitDom/Base.hs view
@@ -19,7 +19,7 @@ -} module Data.Number.ER.RnToRm.UnitDom.Base where -import Prelude hiding (min, max, recip)+import Prelude hiding (min, max, recip, const) import qualified Data.Number.ER.Real.DomainBox as DBox import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainIntBox)@@ -27,344 +27,389 @@ import qualified Data.Number.ER.Real.Base as B import qualified Data.Number.ER.Real.Approx as RA +import Data.Number.ER.Misc+ import qualified Data.Map as Map import Data.Typeable class - (B.ERRealBase b, RA.ERIntApprox ra, Fractional ufb, Ord ufb,+ (B.ERRealBase b, RA.ERIntApprox ra, Ord ufb, DomainBox boxb varid b, DomainIntBox boxra varid ra) => ERUnitFnBase boxb boxra varid b ra ufb | ufb -> boxb boxra varid b ra where++ {--------------} + {----- Miscellaneous associated operations -----}+ {--------------} ++ {-| This should be evaluated before using any of the following operations. -} initialiseBaseArithmetic :: ufb -> IO () initialiseBaseArithmetic _ = B.initialiseBaseArithmetic (0 :: b)+ {-|- Check internal consistency, typically absence of NaN.+ Convert from the associated interval type to the base type.+ (The types are determined by the given example function.) -}- isValid :: ufb -> Bool+ raEndpoints :: + ufb {-^ this parameter is not used except for type checking -} -> + ra -> + (b,b) {-|- A linear ordering, which can be syntactic and rather arbitrary. + Convert from the base type to the associated interval type. + (The types are determined by the given example function.) -}+ raFromEndpoints :: + ufb {-^ this parameter is not used except for type checking -} -> + (b,b) ->+ ra++ {-|+ A linear ordering on basic functions, which can be syntactic and rather arbitrary. + -} compareApprox :: ufb -> ufb -> Ordering++ showDiGrCmp :: + Int {- ^ number of decimal digits to show -} ->+ Bool {-^ whether to show granularity -} ->+ Bool {-^ whether to show internal structure -} ->+ ufb -> String+ + {--------------} + {----- Structural analysis and update of functions -----}+ {--------------} ++ {-|+ Check internal consistency of the basic function, typically absence of NaN.+ -}+ isValid :: ufb -> Bool {-| - Check internal consistency of the function and report problem if any.+ Check internal consistency of the basic function and report problem if any. -} check :: String {-^ indentification of caller location for easier debugging -} -> ufb -> ufb+ + {-| + Get the granularity of the coefficients inside this basic function.+ -} getGranularity :: ufb -> Granularity setMinGranularity :: Granularity -> ufb -> ufb setGranularity :: Granularity -> ufb -> ufb- {-| Construct a constant function. -}- const :: b -> ufb- {-| Construct an affine function. -}- affine :: - b {-^ value at 0 -} ->- Map.Map varid b {-^ ascent of each base vector -} -> - ufb- {-| - Multiply a function by a scalar, - rounding downwards and upwards. - -} - scale :: b -> ufb -> (ufb, ufb) - {-| - Multiply a function by an approximation of a scalar, - rounding downwards and upwards. - -} - scaleApprox :: ra -> ufb -> (ufb, ufb) - {-| - Multiply a function by an approximation of a scalar, - rounding downwards. - -} - scaleApproxDown :: ra -> ufb -> ufb- scaleApproxDown ratio = fst . scaleApprox ratio - {-| - Multiply a function by an approximation of a scalar, - rounding upwards. - -} - scaleApproxUp :: ra -> ufb -> ufb- scaleApproxUp ratio = snd . scaleApprox ratio + {-| - Get the degree of this particular function.+ Get the degree of this basic function. If the function is a polynomial, this function should return its degree. -} getDegree :: ufb -> Int {-| - Decrease the degree of function approximation, - rounding pointwise downwards and upwards.+ Decrease the degree of a basic function, rounding pointwise upwards. -}- reduceDegree :: Int -> ufb -> (ufb, ufb)- {-| - Decrease the degree of function approximation, rounding pointwise downwards.+ reduceDegreeUp :: Int -> ufb -> ufb+ + {-|+ Get the term size of this basic function.+ + If the function is a polynomial, this function should+ return the number of terms in the polynomial. -}- reduceDegreeDown :: Int -> ufb -> ufb- reduceDegreeDown maxDegr = fst . reduceDegree maxDegr+ getSize :: ufb -> Int {-| - Decrease the degree of function approximation, rounding pointwise upwards.+ Decrease the size of this basic function, rounding pointwise upwards. -}- reduceDegreeUp :: Int -> ufb -> ufb- reduceDegreeUp maxDegr = snd . reduceDegree maxDegr- {-| - Approximate the integral of p (with 0 at 0) from below and from above.+ reduceSizeUp :: Int -> ufb -> ufb+ + {-|+ Get a list of all variables featured in this basic function. -}- integrate :: - varid {-^ variable to integrate by -} -> - ufb {-^ p(x) -} -> - (ufb, ufb)- {-| Approximate the integral of p (with 0 at 0) from below. -}- integrateDown :: - varid {-^ variable to integrate by -} -> - ufb {-^ p(x) -} -> - ufb- integrateDown x = fst . integrate x- {-| Approximate the integral of p (with 0 at 0) from above. -}- integrateUp :: - varid {-^ variable to integrate by -} -> - ufb {-^ p(x) -} -> + getVariables :: ufb -> [varid]+ + {--------------} + {----- Construction of basic functions -----}+ {--------------} + + {-| Construct a constant basic function. -}+ const :: b -> ufb+ + {-| Construct a constant basic enclosure (negated lower bound, upper bound). -}+ constEncl :: (b,b) -> (ufb, ufb)+ + {-| Construct an affine basic function. -}+ affine :: + b {-^ value at 0 -} ->+ Map.Map varid b {-^ ascent of each base vector -} -> ufb- integrateUp x = snd . integrate x- {-| - Measure the volume between a function - and the zero hyperplane on the domain @[-1,1]^n@.- -}- volumeAboveZero :: - [varid] {-^ axes to include in the measuring domain -} -> - ufb -> (b,b)++ {--------------}+ {----- Pointwise order operations ----------} + {--------------}+ {-|- Find an upper bound of the function over @[-1,1]^n@.+ Find an upper bound of a basic function over @[-1,1]^n@. -} upperBound :: EffortIndex -> ufb -> b+ {-|- Find a lower bound of the function over @[-1,1]^n@.+ Approximate the function @max(f1,f2)@ from above. -}- lowerBound :: EffortIndex -> ufb -> b- lowerBound ix f = negate $ upperBound ix (negate f)- {-| - Approximate the function max(0,p(x)) from below and from above.+ maxUp :: + Int {-^ max degree for result -} -> + Int {-^ max approx size for result -} ->+ ufb {-^ @f1@ -} -> + ufb {-^ @f2@ -} -> + ufb+ {-|+ Approximate the function @min(f1,f2)@ from above. -}- nonneg ::+ minUp :: Int {-^ max degree for result -} -> - ufb {-^ p(x) -} -> - (ufb, ufb)- {-| - Approximate the function 1/p(x) from below and from above.+ Int {-^ max approx size for result -} ->+ ufb {-^ @f1@ -} -> + ufb {-^ @f2@ -} -> + ufb+ + {--------------} + {----- Field operations ----------}+ {--------------} + + {-| Pointwise exact negation of a basic function -}+ neg :: ufb -> ufb++ {-|+ Multiply a basic function by a scalar, rounding upwards. -}- recip :: - Int {-^ max degree for result -} ->- EffortIndex -> - ufb {-^ p(x) -} -> - (ufb, ufb)+ scaleUp :: b -> ufb -> ufb+ {-| - Approximate the function 1/p(x) from below.- -}- recipDown :: Int -> EffortIndex -> ufb -> ufb- recipDown maxDegr ix a = fst $ recip maxDegr ix a+ Multiply a basic function by an approximation of a scalar, + rounding upwards. + -} + scaleApproxUp :: + Int {-^ maximum polynomial degree -} -> + Int {-^ maximum term count -} -> + ra -> ufb -> ufb+ + {-| Pointwise upwards rounded addition -}+ (+^) :: ufb -> ufb -> ufb+ {-| Pointwise upwards rounded subtraction -}+ (-^) :: ufb -> ufb -> ufb+ {-| Pointwise upwards rounded multiplication -}+ (*^) :: ufb -> ufb -> ufb+ + {-| Enclosure multiplication ++ IMPORTANT: enclosure = (negated lower bound, upper bound) + -}+ multiplyEncl :: + Int {-^ maximum polynomial degree -} -> + Int {-^ maximum term count -} -> + (ufb,ufb) -> (ufb,ufb) -> (ufb, ufb)+ {-| - Approximate the function 1/p(x) from above.+ Approximate the function @1/f@ from above, assuming+ @f@ does not hit zero in the unit domain. -}- recipUp :: Int -> EffortIndex -> ufb -> ufb- recipUp maxDegr ix a = snd $ recip maxDegr ix a+ recipUp :: Int -> Int -> EffortIndex -> ufb -> ufb+ {-|- Approximate the function max(p_1(x),p_2(x)) from below and from above.+ Approximate the reciprocal of an enclosure, assuming+ @f@ does not hit zero in the unit domain.+ + IMPORTANT: enclosure = (negated lower bound, upper bound) -}- max :: - Int {-^ max degree for result -} -> - ufb {-^ p_1(x) -} -> - ufb {-^ p_2(x) -} -> - (ufb, ufb)+ recipEncl :: + Int {-^ max degree for result -} ->+ Int {-^ max approx size for result -} ->+ EffortIndex -> + (ufb,ufb) {-^ enclosure of @f@ -} -> + (ufb,ufb)++ {--------------}+ {----- Evaluation and composition of functions -----}+ {--------------}+ {-|- Approximate the function max(p_1(x),p_2(x)) from below.+ Evaluate a basic function at a point rounding upwards + using a basic number for both the point and the result. -}- maxDown :: - Int {-^ max degree for result -} -> - ufb {-^ p_1(x) -} -> - ufb {-^ p_2(x) -} -> - ufb- maxDown maxDegr a b = fst $ max maxDegr a b+ evalUp :: boxb -> ufb -> b+ {-|- Approximate the function max(p_1(x),p_2(x)) from above.+ Safely evaluate a basic function at a point using a real number approximation+ for both the point and the result. -}- maxUp :: - Int {-^ max degree for result -} -> - ufb {-^ p_1(x) -} -> - ufb {-^ p_2(x) -} -> - ufb- maxUp maxDegr a b = snd $ max maxDegr a b+ evalApprox :: boxra -> ufb -> ra+ {-|- Approximate the function min(p_1(x),p_2(x)) from below and from above.+ Partially evaluate a basic function at a lower-dimensional point + given using a real number approximation.+ Approximate the resulting function from above. -}- min :: - Int {-^ max degree for result -} -> - ufb {-^ p_1(x) -} -> - ufb {-^ p_2(x) -} -> - (ufb, ufb)- min maxDegr p1 p2 = -- default implementation using symmetry with ufbMax- (negate hi, negate lo)- where- (lo, hi) = max maxDegr (negate p1) (negate p2)- {-|- Approximate the function min(p_1(x),p_2(x)) from below.+ partialEvalApproxUp :: boxra -> ufb -> ufb++ {-| + Compose two basic functions, rounding downwards and upwards, + assuming @f_v@ ranges within the domain @[-1,1]@. -}- minDown :: + composeUp :: Int {-^ max degree for result -} -> - ufb {-^ p_1(x) -} -> - ufb {-^ p_2(x) -} -> - ufb- minDown maxDegr a b = fst $ min maxDegr a b- {-|- Approximate the function min(p_1(x),p_2(x)) from above.+ Int {-^ max approx size for result -} ->+ ufb {-^ function @f@ -} -> + varid {-^ variable @v@ to substitute in @f@ -} -> + ufb + {-^ function @f_v@ to substitute for @v@ + that maps @[-1,1]@ into @[-1,1]@ -} ->+ ufb {-^ pointwise upper bound of @f[v |-> f_v]@ -}++ {-| + Compose two basic functions, rounding downwards and upwards, + assuming @f_v@ ranges within the domain @[-1,1]@. -}- minUp :: + composeEncl :: Int {-^ max degree for result -} -> - ufb {-^ p_1(x) -} -> - ufb {-^ p_2(x) -} -> - ufb- minUp maxDegr a b = snd $ min maxDegr a b+ Int {-^ max approx size for result -} ->+ ufb {-^ function @f@ -} -> + varid {-^ variable @v@ to substitute in @f@ -} -> + (ufb, ufb) + {-^ enclosure of a function @f_v@ to substitute for @v@ + that maps @[-1,1]@ into @[-1,1]@ -} ->+ (ufb, ufb) {-^ enclosure of @f[v |-> f_v]@ -}++ {-| + Substitute several variables in a basic function with other basic functions, + rounding downwards and upwards, assuming each @f_v@ ranges + within the domain @[-1,1]@. + -} + composeManyUp ::+ Int {-^ max degree for result -} -> + Int {-^ max approx size for result -} ->+ ufb {-^ function @f@ -} -> + Map.Map varid ufb + {-^ variables to substitute and for each variable @v@, + function @f_v@ to substitute for @v@ + that maps @[-1,1]@ into @[-1,1]@ -} ->+ ufb {-^ pointwise upper bound of @f[v |-> f_v]@ -}++ {-| + Substitute several variables in a basic function with other basic functions, + rounding downwards and upwards, assuming each @f_v@ ranges + within the domain @[-1,1]@. + -} + composeManyEncls ::+ Int {-^ max degree for result -} -> + Int {-^ max approx size for result -} ->+ ufb {-^ function @f@ -} -> + Map.Map varid (ufb, ufb) + {-^ variables to substitute and for each variable @v@, + enclosure of a function @f_v@ to substitute for @v@ + that maps @[-1,1]@ into @[-1,1]@ -} ->+ (ufb, ufb) {-^ enclosure of @f[v |-> f_v]@ -}++ {--------------}+ {----- Selected elementary operations ----------} + {--------------}+ {-|- Approximate @sqrt(p(x))@ from below and from above.+ Approximate @sqrt(f)@ for enclosures. -}- sqrt :: + sqrtEncl :: Int {-^ max degree for result -} -> + Int {-^ max approx size for result -} -> EffortIndex {-^ how hard to try when approximating exp as a polynomial -} -> - ufb {-^ p(x) -} -> + (ufb, ufb) {-^ @f@ -} -> (ufb, ufb) {-|- Approximate @exp(p(x))@ from below and from above.+ Approximate @exp(f)@ for enclosures. -}- exp :: + expEncl :: Int {-^ max degree for result -} -> + Int {-^ max approx size for result -} -> EffortIndex {-^ how hard to try when approximating exp as a polynomial -} -> - ufb {-^ p(x) -} -> + (ufb, ufb) {-^ @f@ -} -> (ufb, ufb) {-| - Approximate @log(p(x))@ from below and from above.+ Approximate @log(f)@ for enclosures. -}- log :: + logEncl :: Int {-^ max degree for result -} -> + Int {-^ max approx size for result -} -> EffortIndex {-^ how hard to try when approximating log as a polynomial -} -> - ufb {-^ p(x) -} -> + (ufb, ufb) {-^ @f@ -} -> (ufb, ufb) {-| - Approximate @sin(p(x))@ from below and from above,- assuming the range of p is within [-pi/2,pi/2].+ Approximate @sin(f)@ for enclosures,+ assuming the range of @f@ is within @[-pi/2,pi/2]@. -}- sin :: + sinEncl :: Int {-^ max degree for result -} -> + Int {-^ max approx size for result -} -> EffortIndex {-^ how hard to try when approximating sin as a polynomial -} -> - ufb {-^ p(x) -} -> - (ufb, ufb)+ (ufb, ufb) {-^ @f@ -} -> + (ufb, ufb) {-|- Approximate @cos(p(x))@ from below and from above,- assuming the range of p is within [-pi/2,pi/2].+ Approximate @cos(f)@ for enclosures,+ assuming the range of @f@ is within @[-pi/2,pi/2]@. -}- cos :: + cosEncl :: Int {-^ max degree for result -} -> + Int {-^ max approx size for result -} -> EffortIndex {-^ how hard to try when approximating cos as a polynomial -} -> - ufb {-^ p(x) -} -> + (ufb, ufb) {-^ @f@ -} -> (ufb, ufb) {-|- Approximate @atan(p(x))@ from below and from above.+ Approximate @atan(f)@ for enclosures. -}- atan :: + atanEncl :: Int {-^ max degree for result -} -> + Int {-^ max approx size for result -} -> EffortIndex {-^ how hard to try when approximating cos as a polynomial -} -> - ufb {-^ p(x) -} -> + (ufb, ufb) {-^ @f@ -} -> (ufb, ufb)- {-| - Evaluate at a point, rounding upwards and downwards.- -}- eval :: boxb -> ufb -> (b, b)- {-| - Evaluate at a point, rounding downwards.- -}- evalDown :: boxb -> ufb -> b- evalDown pt = fst . eval pt- {-| - Evaluate at a point, rounding downwards.- -}- evalUp :: boxb -> ufb -> b- evalUp pt = snd . eval pt- {-|- Safely evaluate at a point using a real number approximation- for both the point and the result.- -}- evalApprox :: boxra -> ufb -> ra- {-|- Partially evaluate at a lower-dimensional point - given using a real number approximation.- Approximate the resulting function from below and from above.- -}- partialEvalApprox :: boxra -> ufb -> (ufb, ufb)- {-|- Partially evaluate at a lower-dimensional point - given using a real number approximation.- Approximate the resulting function from below.- -}- partialEvalApproxDown :: boxra -> ufb -> ufb- partialEvalApproxDown substitutions = fst . partialEvalApprox substitutions+ + {--------------}+ {----- Approximate symbolic integration ----------} + {--------------}+ {-|- Partially evaluate at a lower-dimensional point - given using a real number approximation.- Approximate the resulting function from above.+ Approximate the primitive function of @f@ from below and from above. -}- partialEvalApproxUp :: boxra -> ufb -> ufb- partialEvalApproxUp substitutions = snd . partialEvalApprox substitutions- {-| - Compose two functions, rounding upwards and downwards- provided each @f_v@ ranges within the domain @[-1,1]@. - -} - compose ::- Int {-^ max degree for result -} -> - ufb {-^ function @f@ -} -> - Map.Map varid ufb - {-^ variables to substitute and for each variable @v@, - function @f_v@ to substitute for @v@ - that maps @[-1,1]@ into @[-1,1]@ -} ->- (ufb, ufb) {-^ upper and lower bounds of @f[v |-> f_v]@ -}- {-| - Compose two functions, rounding downwards- provided each @f_v@ ranges within the domain @[-1,1]@. - -} - composeDown ::- Int {-^ max degree for result -} -> - ufb {-^ function @f1@ -} -> - Map.Map varid ufb - {-^ variables to substitute and for each variable @v@, - function @f_v@ to substitute for @v@ - that maps @[-1,1]@ into @[-1,1]@ -} ->- ufb {-^ a lower bound of @f1.f2@ -}- composeDown maxDegr f = fst . compose maxDegr f + integrate ::+ varid {-^ variable to integrate by -} -> + ufb {-^ @f@ -} -> + (ufb, ufb)+ {-| - Compose two functions, rounding upwards- provided each @f_v@ ranges within the domain @[-1,1]@. - -} - composeUp ::- Int {-^ max degree for result -} -> - ufb {-^ function @f1@ -} -> - Map.Map varid ufb - {-^ variables to substitute and for each variable @v@, - function @f_v@ to substitute for @v@ - that maps @[-1,1]@ into @[-1,1]@ -} ->- ufb {-^ an upper bound of @f1.f2@ -}- composeUp maxDegr f = snd . compose maxDegr f - {-|- Convert from the interval type to the base type.- (The types are determined by the given example function.)- -}- raEndpoints :: - ufb {-^ this parameter is not used except for type checking -} -> - ra -> - (b,b)- {-|- Convert from the base type to the interval type. - (The types are determined by the given example function.)+ Measure the volume between a function + and the zero hyperplane on the domain @[-1,1]^n@. -}- raFromEndpoints :: - ufb {-^ this parameter is not used except for type checking -} -> - (b,b) ->- ra+ volumeAboveZeroUp :: + [varid] + {-^ dimensions to include in the measuring domain; + have to include all those present in @f@ -} -> + ufb {-^ @f@ -} -> + b+ volumeAboveZeroUp vars p =+-- unsafePrint ("chplVolumeAboveZero: returning:" ++ show result) $+-- unsafePrint ("chplVolumeAboveZero: vars = " ++ show vars) $+ result+ where+ result = integUpAtEvenCorners - integDownAtOddCorners+ integUpAtEvenCorners = sumUp $ map (\pt -> evalUp pt integUp) evenCorners+ integDownAtOddCorners = sumUp $ map (\pt -> evalUp pt integDownNeg) oddCorners+ evenCorners = map (DBox.fromList) evenCornersL+ oddCorners = map (DBox.fromList) oddCornersL+ (evenCornersL, oddCornersL) =+ allPairsCombinationsEvenOdd $ zip vars $ repeat (1,-1)+ integUp = integrateByAllVars snd p vars+ integDownNeg = neg $ integrateByAllVars fst p vars+ integrateByAllVars pick p [] = p+ integrateByAllVars pick p (x : xs) =+ integrateByAllVars pick ip xs+ where+ ip = pick $ integrate x p+
src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs view
@@ -1,5 +1,4 @@ {-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE UndecidableInstances #-} {-|@@ -27,10 +26,14 @@ where import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic-import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Field import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Eval+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Reduce+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.Compose import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Integration+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Division import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Elementary import qualified Data.Number.ER.RnToRm.UnitDom.Base as UFB@@ -38,6 +41,8 @@ import Data.Number.ER.Real.Approx.Interval import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox) +import qualified Data.Map as Map+ {- code for testing purpose, to be deleted later -} import Data.Number.ER.Real.DefaultRepr import Data.Number.ER.Real.DomainBox.IntMap@@ -47,50 +52,77 @@ x2 = chplVar 2 :: P x3 = chplVar 3 :: P x4 = chplVar 4 :: P-p1 = x1 * x1 * x1 + x1 * (x2 + 2) * (x3 - 3)+p1 = x1 *^ x1 *^ x1 +^ x1 *^ (x2 +^ (chplConst 2)) *^ (x3 -^ (chplConst 3)) {- end of code for testing purposes -} - instance (B.ERRealBase rb, RealFrac rb, DomainBox box varid Int, Ord box,- DomainBoxMappable boxb boxbb varid rb [(rb,rb)],+ DomainBoxMappable boxb boxras varid rb [ERInterval rb], DomainBoxMappable boxra boxras varid (ERInterval rb) [ERInterval rb], DomainIntBox boxra varid (ERInterval rb)) => (UFB.ERUnitFnBase boxb boxra varid rb (ERInterval rb) (ERChebPoly box rb)) where+ {----- Miscellaneous associated operations -----}+ raEndpoints _ (ERInterval l h) = (l,h)+ raEndpoints _ ERIntervalAny = (- B.plusInfinity, B.plusInfinity)+ raFromEndpoints _ (l,h) = normaliseERInterval (ERInterval l h)+ compareApprox = chplCompareApprox+ showDiGrCmp = chplShow + + {----- Structural analysis and update of functions -----} isValid = chplHasNoNaNOrInfty check = chplCheck- compareApprox = chplCompareApprox getGranularity = chplGetGranularity setMinGranularity = chplSetMinGranularity setGranularity = chplSetGranularity+ getDegree = chplGetDegree+ reduceDegreeUp = chplReduceDegreeUp+ getSize = chplCountTerms+ reduceSizeUp = chplReduceTermCountUp+ getVariables = chplGetVars+ + {----- Construction of basic functions -----} const = chplConst+ constEncl (low, high) = (chplConst (-low), chplConst high) affine = chplAffine- scale = chplScale- scaleApprox (ERInterval ratioDown ratioUp) = chplScaleApprox (ratioDown, ratioUp) --- Arity = chplGetArity- getDegree = chplGetDegree- reduceDegree = chplReduceDegree- volumeAboveZero = chplVolumeAboveZero+ + {----- Pointwise order operations ----------} + upperBound = chplUpperBound+ maxUp = chplMaxUp+ minUp = chplMinUp+ + {----- Field operations ----------}+ neg = chplNeg+ scaleUp = chplScaleUp+ scaleApproxUp = chplScaleRAUp+ (+^) = (+^)+ (-^) = (-^)+ (*^) = (*^)+ multiplyEncl = enclMultiply+ recipUp md mt ix f = snd $ enclRecip md mt ix (md + 1) (chplNeg f, f)+ recipEncl md mt ix = enclRecip md mt ix (md + 1)+ + {----- Evaluation and composition of functions -----}+ evalUp pt f = chplEvalUp f pt+ evalApprox x ufb = chplRAEval (\ b -> ERInterval b b) ufb x+ + partialEvalApproxUp substitutions ufb =+ snd $ + chplPartialRAEval (UFB.raEndpoints ufb) ufb substitutions+ composeUp m n f v fv = snd $ enclCompose m n f v (enclThin fv) + composeEncl = enclCompose+ composeManyUp m n f subst = snd $ enclComposeMany m n f (Map.map enclThin subst)+ composeManyEncls = enclComposeMany++ {----- Selected elementary operations ----------}+ sqrtEncl = enclSqrt + expEncl = enclExp+ logEncl = enclLog+ sinEncl = enclSine+ cosEncl = enclCosine+ atanEncl = enclAtan+ integrate = chplIntegrate- upperBound = chplUpperBoundAffine--- upperBound = chplUpperBoundQuadr- nonneg = chplNonneg- recip = chplRecip- max = chplMax- sqrt = chplSqrt- exp = chplExp- log = chplLog- sin = chplSine- cos = chplCosine- atan = chplAtan- eval = chplEval- evalApprox ufb x = chplEvalApprox (\ b -> ERInterval b b) ufb x- partialEvalApprox substitutions ufb = - chplPartialEvalApprox (UFB.raEndpoints ufb) substitutions ufb- raEndpoints _ (ERInterval l h) = (l,h)- raEndpoints _ ERIntervalAny = (- B.plusInfinity, B.plusInfinity)- raFromEndpoints _ (l,h) = normaliseERInterval (ERInterval l h)- compose = chplCompose+
src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Basic.hs view
@@ -46,7 +46,7 @@ { chplCoeffs :: (Map.Map (TermKey box) b) }- deriving (Eq, Typeable, Data)+ deriving (Eq, Ord, Typeable, Data) type TermKey box = box @@ -138,7 +138,7 @@ (ERChebPoly $ Map.singleton chplConstTermKey val) {-|- make a basic "x" polynomial for a given variable number + Make a basic "x" polynomial for a given variable number. -} chplVar :: (B.ERRealBase b, DomainBox box varid Int, Ord box) => @@ -147,41 +147,66 @@ chplVar varName = ERChebPoly $ Map.singleton (DBox.singleton varName 1) 1 ---{-|--- Make a univariate polynomial given by a series of coefficients--- in the Chebyshev basis. ----}---chplMakeUnivariate ::+{-|+ Construct an affine polynomial.+-}+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+++--chplRemoveZeroTermsDown, chplRemoveZeroTermsUp :: -- (B.ERRealBase b, DomainBox box varid Int, Ord box) =>--- varid ->--- [(Int, b)] {-^ list of pairs: degree of Chebyshev polynomial + coefficient -} ->--- ERChebPoly box b---chplMakeUnivariate varName powCoeffPairs =--- ERChebPoly $ Map.fromList $ map encodePow powCoeffPairs+-- ERChebPoly box b -> ERChebPoly box b+--chplRemoveZeroTermsDown = chplNeg . fst . chplRemoveZeroTerms+--chplRemoveZeroTermsUp = snd . chplRemoveZeroTerms++--chplRemoveZeroTerms ::+-- (B.ERRealBase b, DomainBox box varid Int, Ord box) =>+-- ERChebPoly box b -> (ERChebPoly box b, ERChebPoly box b)+--chplRemoveZeroTerms (ERChebPoly coeffs) =+-- (chplNeg $ ERChebPoly $ coeffsNo0T0Down,+-- ERChebPoly $ coeffsNo0T0Up) -- where--- encodePow (pow, coeff) =--- (DBox.singleton varName pow, coeff)+-- coeffsNo0T0Down =+-- Map.insertWith plusDown chplConstTermKey (- err) coeffsNo0T0+-- coeffsNo0T0Up =+-- Map.insertWith plusUp chplConstTermKey err coeffsNo0T0+-- (coeffsNo0T0, err) = +-- foldl addTermNo0T0 (Map.empty, 0) $ Map.toList coeffs+-- addTermNo0T0 (prevCoeffs, prevErr) (term, coeff) +-- | coeff == 0 =+-- (prevCoeffs, prevErr)+-- | otherwise =+-- (newCoeffs, newErr)+-- where+-- newTerm =+-- DBox.filter (> 0) term+-- newCoeffs = +-- Map.insert newTerm newCoeffUp prevCoeffs+-- newCoeffUp = prevCoeff + coeff+-- newCoeffDown = prevCoeff `plusDown` coeff+-- prevCoeff =+-- Map.findWithDefault 0 newTerm prevCoeffs+-- newErr = prevErr + newCoeffUp - newCoeffDown -chplNormaliseDown, chplNormaliseUp ::+chplRemoveZeroTermsUp :: (B.ERRealBase b, DomainBox box varid Int, Ord box) => ERChebPoly box b -> ERChebPoly box b-chplNormaliseUp = snd . chplNormalise-chplNormaliseDown = fst . chplNormalise--chplNormalise ::- (B.ERRealBase b, DomainBox box varid Int, Ord box) =>- ERChebPoly box b -> (ERChebPoly box b, ERChebPoly box b)-chplNormalise (ERChebPoly coeffs) =- (ERChebPoly $ coeffsNo0T0Down,- ERChebPoly $ coeffsNo0T0Up)+chplRemoveZeroTermsUp (ERChebPoly coeffs) =+ ERChebPoly coeffsNo0T0Up where- coeffsNo0T0Down =- Map.insertWith plusDown chplConstTermKey err coeffsNo0T0 coeffsNo0T0Up = Map.insertWith plusUp chplConstTermKey err coeffsNo0T0 (coeffsNo0T0, err) = foldl addTermNo0T0 (Map.empty, 0) $ Map.toList coeffs- addTermNo0T0 (prevCoeffs, prevErr) (term, coeff) + addTermNo0T0 (prevCoeffs, prevErr) (term, coeff) | coeff == 0 = (prevCoeffs, prevErr) | otherwise =@@ -195,12 +220,42 @@ newCoeffDown = prevCoeff `plusDown` coeff prevCoeff = Map.findWithDefault 0 newTerm prevCoeffs- newErr = newCoeffUp - newCoeffDown+ newErr = prevErr + newCoeffUp - newCoeffDown +--chplRemoveLowCoeffsDown, chplRemoveLowCoeffsUp ::+-- (B.ERRealBase b, DomainBox box varid Int, Ord box) =>+-- b -> ERChebPoly box b -> ERChebPoly box b+--chplRemoveLowCoeffsDown maxCoeff = chplNeg . fst . chplRemoveLowCoeffs maxCoeff+--chplRemoveLowCoeffsUp maxCoeff = snd . chplRemoveLowCoeffs maxCoeff++--chplRemoveLowCoeffs ::+-- (B.ERRealBase b, DomainBox box varid Int, Ord box) =>+-- b -> ERChebPoly box b -> (ERChebPoly box b, ERChebPoly box b)+--chplRemoveLowCoeffs maxCoeff (ERChebPoly coeffs) =+-- (chplNeg $ ERChebPoly $ coeffsNoLowDown,+-- ERChebPoly $ coeffsNoLowUp)+-- where+-- coeffsNoLowDown =+-- Map.insertWith plusDown chplConstTermKey (- err) coeffsNoLow+-- coeffsNoLowUp =+-- Map.insertWith plusUp chplConstTermKey err coeffsNoLow+-- err = sum $ map abs $ Map.elems coeffsLow+-- (coeffsLow, coeffsNoLow) = +-- Map.partition (\ c -> abs c < maxCoeff) coeffs++chplCountTerms ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) =>+ ERChebPoly box b -> Int+chplCountTerms (ERChebPoly coeffs) =+ Map.size coeffs+++{------------------ Formatting ------------------------}+ instance (B.ERRealBase b, DomainBox box varid Int, Ord box) => Show (ERChebPoly box b) where--- show = chplShow True- show = chplShow False+-- show = chplShow 8 False True+ show = chplShow 8 False False {-| Convert a polynomial to a string representation,@@ -208,15 +263,17 @@ -} chplShow :: (B.ERRealBase b, DomainBox box varid Int, Ord box) =>+ Int {- ^ number of decimal digits to show -} ->+ Bool {-^ whether to show granularity -} -> Bool {-^ show the polynomial also in its native Chebyshev basis -} -> ERChebPoly box b -> String-chplShow showChebyshevBasis (ERChebPoly coeffs) +chplShow digitsToShow showGranularity showChebyshevBasis (ERChebPoly coeffs) | showChebyshevBasis = "\n" ++ inChebBasis ++ " = \n" ++ inXBasis | otherwise = inXBasis where inChebBasis = - showCoeffs showTermT $ coeffs+ showCoeffs showTermT $ Map.filter (\c -> c /= 0) $ coeffs inXBasis = showCoeffs showTermX $ chebToXBasis coeffs showCoeffs showTerm coeffs =@@ -231,7 +288,7 @@ showC coeff ++ "*" ++ (concat $ map showX $ DBox.toList term) showT (var, deg) = "T" ++ show deg ++ "(" ++ showVar var ++ ")" showX (var, deg) = showVar var ++ "^" ++ show deg- showC = B.showDiGrCmp 8 False False+ showC = B.showDiGrCmp digitsToShow showGranularity False {-| conversion of polynomials from Chebyshev basis to the X^n basis@@ -243,7 +300,8 @@ (Map.Map (TermKey box) b) {-^ polynomial in Chebyshev basis -} -> (Map.Map (TermKey box) b) {-^ approxition of the equivalent polynomial in X^n basis -} chebToXBasis coeffs =- Map.foldWithKey addTerm Map.empty coeffs+ Map.filter (\c -> c /= 0) $+ Map.foldWithKey addTerm Map.empty coeffs where addTerm term coeff prevXCoeffs = Map.unionWith (+) prevXCoeffs $
src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Bounds.hs view
@@ -17,8 +17,9 @@ 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.Reduce import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Eval-import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Field import qualified Data.Number.ER.Real.Approx as RA import qualified Data.Number.ER.Real.Base as B@@ -37,60 +38,71 @@ Find an upper bound on a polynomial over the unit domain [-1,1]^n. -}-chplUpperBoundAffine ::+chplUpperBound :: (B.ERRealBase b, DomainBox box varid Int, Ord box) => EffortIndex {-^ how hard to try -} -> ERChebPoly box b -> b-chplUpperBoundAffine ix (ERChebPoly coeffs) =- affiBound coeffs- where- affiBound coeffs =- Map.fold (+) constTerm absCoeffs- where- absCoeffs = Map.map abs $ Map.delete chplConstTermKey coeffs- constTerm = Map.findWithDefault 0 chplConstTermKey coeffs-+chplUpperBound ix p = snd $ chplBounds ix p {-|- Find a close upper bound on an affine polynomial over the + Find a lower bound on a polynomial over the unit domain [-1,1]^n. -}-chplUpperBoundAffineCorners ::- (B.ERRealBase b, DomainBox box varid Int, Ord box,- DomainBoxMappable boxb boxbb varid b [(b,b)], Num varid, Enum varid) => +chplLowerBound ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) => EffortIndex {-^ how hard to try -} -> ERChebPoly box b -> b-chplUpperBoundAffineCorners ix p@(ERChebPoly coeffs) =- affiBound (coeffs, vars)+chplLowerBound ix p = fst $ chplBounds ix p++{-|+ Find both lower and upper bounds on a polynomial over the + unit domain [-1,1]^n. +-}+chplBounds ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) => + EffortIndex {-^ how hard to try -} ->+ ERChebPoly box b ->+ (b,b)+chplBounds = chplBoundsAffine++{-|+ Find bounds on a polynomial over the unit domain [-1,1]^n.+ + Fast but inaccurate method, in essence+ taking the maximum of the upper affine reduction.+-}+chplBoundsAffine ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) => + EffortIndex {-^ how hard to try -} ->+ ERChebPoly box b ->+ (b,b)+chplBoundsAffine ix p@(ERChebPoly coeffs) =+-- unsafePrintReturn+-- (+-- "chplBoundsAffine:"+-- ++ "\n p = " ++ show p+-- ++ "\n noConstCoeffAbsSum = " ++ show noConstCoeffAbsSum+-- ++ "\n result = "+-- ) + result where- vars = chplGetVars p- affiBound (coeffs, vars)- | null vars =- Map.findWithDefault 0 chplConstTermKey coeffs- | otherwise =- foldl1 max cornerValues- where- cornerValues =- map (\pt -> chplEvalUp pt p) corners- where--- corners :: [boxb]- corners = - map (DBox.fromList . (zip [1..n])) $ prod n- where- n = fromInteger $ toInteger $ length vars- -- n-fold product list of [-1,1]- prod n - | n == 1 = [[-1],[1]]- | otherwise =- (map ((-1):) prodNm1) ++ (map (1:) $ prodNm1)- where- prodNm1 = prod (n-1)+ result =+ (constTerm `plusDown` (- noConstCoeffAbsSum),+ constTerm `plusUp` noConstCoeffAbsSum)+ noConstCoeffAbsSum = Map.fold plusUp 0 absCoeffs+ absCoeffs = Map.map abs $ Map.delete chplConstTermKey coeffs+ constTerm = Map.findWithDefault 0 chplConstTermKey coeffs {-| Find a close upper bound on a quadratic polynomial over the unit domain [-1,1]^n. ++ Much slower and somewhat more accurate method, in essence+ taking the maximum of the upper quadratic reduction.+ + !!! Not yet properly tested !!! -} chplUpperBoundQuadr :: (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box,@@ -101,9 +113,10 @@ ERChebPoly box b -> b chplUpperBoundQuadr ix p@(ERChebPoly coeffs) =- quadBound (coeffs, vars)+ quadBound (coeffsQ, vars) where- vars = chplGetVars p+ pQ@(ERChebPoly coeffsQ) = chplReduceDegreeUp 2 p+ vars = chplGetVars pQ quadBound (coeffs, vars) | null vars = Map.findWithDefault 0 chplConstTermKey coeffs@@ -122,7 +135,7 @@ (and $ map maybeInUnit $ DBox.elems peak) , erintv_right $- chplEvalApprox makeInterval peak p + chplRAEval makeInterval p peak ) Nothing -> (False, undefined) where@@ -167,7 +180,7 @@ newVars = var `delete` vars substVar isOne = chplCoeffs $- sum $ + foldl (+^) (chplConst 0) $ map (makeMonomial isOne) $ Map.toList coeffs makeMonomial isOne (term, coeff) =@@ -187,37 +200,61 @@ _ -> [(term, coeff)] -chplMaxDn m a b = fst $ chplMax m a b-chplMaxUp m a b = snd $ chplMax m a b-chplMinDn m a b = fst $ chplMin m a b-chplMinUp m a b = snd $ chplMin m a b--chplMin m a b =- (-u,-l)- where- (l,u) = chplMax m (-a) (-b)- {-| Approximate from below and from above the pointwise maximum of two polynomials -} chplMax ::- (B.ERRealBase b, DomainBox box varid Int, Ord box) => + (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) => Int {-^ maximum polynomial degree -} -> + Int {-^ maximum term count -} -> ERChebPoly box b -> ERChebPoly box b -> (ERChebPoly box b, ERChebPoly box b)-chplMax maxDegree p1 p2 =- (p1 `plusDown` differenceDown, p1 `plusUp` differenceUp)+chplMax maxDegree maxSize p1 p2 =+ (p1 +. differenceDown, p1 +^ differenceUp) where- (differenceDown, differenceUp) = chplNonneg maxDegree $ p2 - p1+ (differenceDown, _) = chplNonneg maxDegree maxSize p2MinusP1Down+ (_, differenceUp) = chplNonneg maxDegree maxSize $ p2MinusP1Up+ (p2MinusP1Down, p2MinusP1Up, _) = chplAdd p2 (chplNeg p1) +chplMaxDn m s a b = fst $ chplMax m s a b+chplMaxUp m s a b = snd $ chplMax m s a b+chplMinDn m s a b = fst $ chplMin m s a b+chplMinUp m s a b = snd $ chplMin m s a b+ {-|+ Approximate from below and from above the pointwise minimum of two polynomials+-}+chplMin ::+ (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) => + Int {-^ maximum polynomial degree -} -> + Int {-^ maximum term count -} -> + ERChebPoly box b ->+ ERChebPoly box b ->+ (ERChebPoly box b, ERChebPoly box b)+chplMin m s a b =+ (chplNeg u,chplNeg l)+ where+ (l,u) = chplMax m s (chplNeg a) (chplNeg b)++chplNonnegDown m s p = fst $ chplNonneg m s p+chplNonnegUp m s p = snd $ chplNonneg m s p +chplNonposDown m s p = fst $ chplNonpos m s p+chplNonposUp m s p = snd $ chplNonpos m s p ++chplNonpos m s p =+ (chplNeg h, chplNeg l)+ where+ (l,h) = chplNonneg m s (chplNeg p)++{-| Approximate the function max(0,p(x)) by a polynomial from below and from above. -} chplNonneg ::- (B.ERRealBase b, DomainBox box varid Int, Ord box) => + (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) => Int {-^ maximum polynomial degree -} -> + Int {-^ maximum term count -} -> ERChebPoly box b -> (ERChebPoly box b, ERChebPoly box b) chplNonneg = chplNonnegCubic@@ -226,112 +263,161 @@ A version of 'chplNonneg' using a cubic approximation. -} chplNonnegCubic ::- (B.ERRealBase b, DomainBox box varid Int, Ord box) => + (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) => Int {-^ maximum polynomial degree -} -> + Int {-^ maximum term count -} -> ERChebPoly box b -> (ERChebPoly box b, ERChebPoly box b)-chplNonnegCubic maxDegree p+chplNonnegCubic maxDegree maxSize p | upperB <= 0 = (chplConst 0, chplConst 0) | lowerB >= 0 = (p, p)+ | not allInterimsBounded = (chplConst (1/0), chplConst (1/0)) | otherwise = -- ie lowerB < 0 < upperB: polynomial may be crossing 0...+-- unsafePrintReturn+-- (+-- "chplNonnegCubic:"+-- ++ "\n p = " ++ show p+-- ++ "\n maxDegree = " ++ show maxDegree+-- ++ "\n maxSize = " ++ show maxSize+-- ++ "\n upperB = " ++ show upperB+-- ++ "\n lowerB = " ++ show lowerB+-- ++ "\n a0 = " ++ show a0+-- ++ "\n a1 = " ++ show a1+-- ++ "\n a2 = " ++ show a2+-- ++ "\n a3 = " ++ show a3+-- ++ "\n b = " ++ show b+-- ++ "\n rb = " ++ show rb+-- ++ "\n correctionB = " ++ show correctionB+-- ++ "\n valueAt0B = " ++ show valueAt0B+-- ++ "\n result = "+-- ) -- work out the cubic polynomial (a3*x^3 + a2*x^2 + a1*x + a0) / b -- that hits 0 at lowerB with derivative 0 -- and hits upperB at upperB with derivative 1 - (cubicAppliedOnPDown - valueAt0, cubicAppliedOnPUp + (chplConst correction))+ (chplAddConstDown (- valueAt0B) cubicAppliedOnPDown, + chplAddConstUp correctionB cubicAppliedOnPUp) where - upperB = chplUpperBoundAffine 10 p - lowerB = - (chplUpperBoundAffine 10 (- p))- cubicAppliedOnPUp = evalCubic multiplyByPUp- cubicAppliedOnPDown = evalCubic multiplyByPDown- evalCubic multiplyByP =- p0 * (chplConst $ recip b)+ (lowerB, upperB) = chplBounds 10 p+ (cubicAppliedOnPDown, cubicAppliedOnPUp, width) =+ p0 `scaleByPositiveConsts` (rbLo, rbHi) where- p0 = multiplyByP p1 + (chplConst a0) -- ie p*(p*(p * a3 + a2) + a1) + a0- p1 = multiplyByP p2 + (chplConst a1) -- ie p*(p * a3 + a2) + a1- p2 = multiplyByP p3 + (chplConst a2) -- ie p * a3 + a2- p3 = chplConst a3- multiplyByPUp =- chplReduceDegreeUp maxDegree . (p *)- multiplyByPDown =- chplReduceDegreeDown maxDegree . (p *)+ p0 = (multiplyByP p1) `addConsts` (a0Lo, a0Hi) -- ie p*(p*(p * a3 + a2) + a1) + a0 enclosure+ p1 = (multiplyByP p2) `addConsts` (a1Lo, a1Hi) -- ie p*(p * a3 + a2) + a1 enclosure+ p2 = (multiplyByP p3) `addConsts` (a2Lo, a2Hi) -- ie p * a3 + a2 enclosure+ p3 = (chplConst a3Lo, chplConst a3Hi, a3Hi - a3Lo) -- ie a3 enclosure+ multiplyByP (lo,hi,wd) =+ (ploRed, phiRed, pwd)+ where+ ploRed = reduceDgSzDown plo+ phiRed = reduceDgSzUp phi + pwd = chplUpperBound 10 $ phiRed -^ ploRed + (plo, phi, _) = chplTimesLoHi p (lo,hi,wd)+ reduceDgSzUp =+ chplReduceTermCountUp maxSize . chplReduceDegreeUp maxDegree+ reduceDgSzDown =+ chplReduceTermCountDown maxSize . chplReduceDegreeDown maxDegree+ addConsts (lo, hi, wd) (cLo, cHi) =+ (alo, ahi, wd + wdlo + wdhi)+ where+ (alo, _, wdlo) = chplAddConst cLo lo + (_, ahi, wdhi) = chplAddConst cHi hi + scaleByPositiveConsts (lo, hi, wd) (cLo, cHi) =+ (slo, shi, wd + wdlo + wdhi)+ where+ (slo, _, wdlo) = chplScale cLo lo + (_, shi, wdhi) = chplScale cHi hi + + -- convert interval coefficients to pairs of bounds:+ ERInterval rbLo rbHi = rb+ ERInterval a3Lo a3Hi = a3+ ERInterval a2Lo a2Hi = a2+ ERInterval a1Lo a1Hi = a1+ ERInterval a0Lo a0Hi = a0+ allInterimsBounded = + and $ map RA.isBounded [w, s, rb, a0, a1, a2, a3, correction] {- The cubic polynomial's coefficients are calculated by solving a system of 4 linear eqs. The generic solution is as follows:- b = (r - l)^3+ b = (r - l)^3 always positive a3 = -(r + l) a2 = 2*(r^2 + r*l + l^2) a1 = -l*(4*r^2 + r*l + l^2) a0 = 2*r^2*l^2 -}- r = upperB- l = lowerB- b = - ((r - l) * ((r - l) * (l - r))) - -- this one has to round downwards because it is a denominator- a3 = (- r) + (- l) -- remember to round upwards!- a2 = 2*(r2rll2Up)- a1 = (- l) * (r2rll2Up + 3*rSqUp) -- since l < 0, the other argument is rounded upwards- a0 = 2 * rSqUp * lSqUp- r2rll2Up = rSqUp + r*l + lSqUp - rSqUp = r*r- lSqUp = l*l- rSqDown = -((-r)*r)- lSqDown = -((-l)*l)+ rb = recip b+ b = w3 -- = w^3 -- see below+ w = r - l+ r = ERInterval upperB upperB+ l = ERInterval lowerB lowerB+ --+ a3 = - s+ s = r + l+ --+ a2 = 2 * (r2PrlPl2)+ r2PrlPl2 = s2 - rl+ rl = r * l+ --+ a1 = (- l) * (r2PrlPl2 + 3*r2)+ a0 = 2*r2*l2+ -- interval arithmetic tricks to speed it up and reduce dependency errors:+ w3 = ERInterval (wLo * wLo * wLo) (wHi * wHi * wHi) -- x^3 is monotone + ERInterval wLo wHi = w+ s2 = ERInterval (max 0 s2Lo) s2Hi+ s2Lo = min sLo2 sHi2 + s2Hi = max sLo2 sHi2+ sLo2 = sLo * sLo+ sHi2 = sHi * sHi + ERInterval sLo sHi = s + r2 = ERInterval (upperB `timesDown` upperB) (upperB `timesUp` upperB) + l2 = ERInterval (lowerB `timesDown` lowerB) (lowerB `timesUp` lowerB) {- The cubic polynomial may sometimes fail to dominate x or sometimes it dips below 0. Work out the amount by which it has to be lifted up to fix these problems. -}- correction- | 2*rSqDown < l*(r + l) =- erintv_right $- (peak0 * (peak0 * (peak0 * (-a3I) - a2I) - a1I) - a0I) / bI- | 2*lSqDown < r*(r + l) =- erintv_right $- ((peakP * (peakP * (peakP * (-a3I) - a2I) - a1I) - a0I) / bI) + peakP- | otherwise = 0+ ERInterval _ correctionB = correction+ correction =+ case (RA.compareReals (2 * r2) (l*s), RA.compareReals (2 * l2) (r*s)) of+ (Just LT, _) ->+ (peak0 * (peak0 * (peak0 * (-a3) - a2) - a1) - a0) / b+ (_, Just LT) ->+ ((peakP * (peakP * (peakP * (-a3) - a2) - a1) - a0) / b) + peakP+ _ -> 0 where- -- these have to be computed interval-based:- [a0I, a1I, a2I, a3I, bI, lI, rI] = - map (\x -> ERInterval x x) [a0,a1,a2,a3,b,l,r]- peak0 = (lI + 4*rI*rI/(lI+rI)) / 3 - peakP = (rI + 4*lI*lI/(lI+rI)) / 3+ peak0 = (l + 4*r2/s) / 3 + peakP = (r + 4*l2/s) / 3 {- The same cubic polynomial can be used as a lower bound when we subtract its value at 0 rounded upwards. -}- valueAt0 = chplConst $ a0 / b+ valueAt0B = + case a0 / b of+ ERInterval lo hi -> hi+ ERIntervalAny -> 1/0 {-|- Multiply a thin enclosure by a non-thin enclosure+ Multiply a polynomial by an enclosure (with non-negated lower bound). -}-chplThinTimesEncl ::+chplTimesLoHi :: (B.ERRealBase b, DomainBox box varid Int, Ord box) => - Int {-^ maximum polynomial degree -} -> ERChebPoly box b ->- (ERChebPoly box b, ERChebPoly box b) ->- (ERChebPoly box b, ERChebPoly box b)-chplThinTimesEncl maxDegree p1 (p2LO, p2HI) =- (prodLO, prodHI)+ (ERChebPoly box b, ERChebPoly box b, b) ->+ (ERChebPoly box b, ERChebPoly box b, b)+chplTimesLoHi p1 (p2Low, p2High, p2Width) =+ (prodMid -. (chplConst width), + prodMid +^ (chplConst width), + 2 * width) where- prodHI =- chplMaxUp maxDegree - (p1 `timesUp` p2HI)- (p1 `timesUp` p2LO) -- beware: p1 can be negative- prodLO =- negate $- chplMaxUp maxDegree - (p1n `timesUp` p2HI)- (p1n `timesUp` p2LO)- p1n = negate p1--{-|- Safely multiply a polynomial by itself.--}-chplSquare ::- (B.ERRealBase b, DomainBox box varid Int, Ord box) => - Int {-^ maximum polynomial degree -} -> - ERChebPoly box b ->- (ERChebPoly box b, ERChebPoly box b)-chplSquare maxDegree p =- (p `timesDown` p, p `timesUp` p)+ prodMid = prodLowUp+ (prodLowDown, prodLowUp, prodLowWidth) = + chplMultiply p1 p2Low+ (prodHighDown, prodHighUp, prodHighWidth) = + chplMultiply p1 p2High+ width = + p1Norm `timesUp` p2Width `plusUp` prodLowWidth `plusUp` prodHighWidth+ p1Norm = + max (abs $ p1LowerBound) (abs $ p1UpperBound)+ (p1LowerBound, p1UpperBound) = + chplBounds ix p1+ ix = 10
+ src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Compose.hs view
@@ -0,0 +1,138 @@+{-# LANGUAGE FlexibleContexts #-}+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Compose+ Description : (internal) composition of 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 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+
+ src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Division.hs view
@@ -0,0 +1,168 @@+{-# LANGUAGE FlexibleContexts #-}+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Division+ Description : (internal) division of 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 elementary functions applied to polynomials.+-}+module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Division +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 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 cosine of a polynomial + by another polynomial from below and from above+ using the tau method + as described in [Mason & Handscomb 2003, p 62]. +-}+enclRecip ::+ (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) => + Int {-^ maximum polynomial degree -} -> + Int {-^ maximum term count -} -> + EffortIndex {-^ minimum approx degree -} -> + Int {-^ degree of tau expansion -} -> + (ERChebPoly box b, ERChebPoly box b) ->+ (ERChebPoly box b, ERChebPoly box b)+enclRecip maxDegree maxSize ix tauDegr pEncl@(pLowNeg, pHigh)+ | pIsConst =+ enclRAConst (recip pConst)+ | upperB < 0 = -- range negative+ enclNeg $ enclRecip maxDegree maxSize ix tauDegr (enclNeg pEncl)+ | lowerB > 0 = -- range positive+-- unsafePrintReturn+-- (+-- "ERChebPoly: enclRecip: "+-- ++ "\n pEncl = " ++ show pEncl+-- ++ "\n lowerB = " ++ show lowerB+-- ++ "\n upperB = " ++ show upperB+-- ++ "\n k = " ++ show k+-- ++ "\n pAbove1Encl = " ++ show pAbove1Encl+-- ++ "\n trT1Encl = " ++ show trT1Encl+-- ++ "\n nu = " ++ show nu+-- ++ "\n c0n = " ++ show c0n+--+-- ++ "\n tauDegree = " ++ (show $ tauDegree)+-- ++ "\n tauInv = " ++ (show $ tauInv)+-- ++ "\n tau = " ++ (show $ recip tauInv)+-- ++ "\n errScaleUp = " ++ (show $ errScaleUp)+-- ++ "\n errScaleDown = " ++ (show $ errScaleDown)+-- ++ "\n resEncl = "+-- ) $+ case allInterimsBounded of+ True -> resEncl+ False -> (chplConst 0, chplConst (1/0))+ | otherwise = -- cannot establish 0 freedom+ error $+ "ERChebPoly: enclRecip: "+ ++ "cannot deal with estimated range " ++ show ranp+ ++ "of polynomial enclosure: \n" ++ show pEncl+ where+ ranp = ERInterval lowerB upperB+ (lowerB, upperB) = enclBounds ix pEncl+ + (pIsConst, pConst) = + case (chplGetConst pLowNeg, chplGetConst pHigh) of+ (Just pConstLowNeg, Just pConstHigh) ->+ (True, ERInterval (- pConstLowNeg) pConstHigh)+ _ ->+ (False, error "ERChebPoly: chplRecip: internal error")+ + tauDegree = max 2 tauDegr+ coeffGr = effIx2gran $ ix+ + -- translate p to have range above 1:+ k = intLogUp 2 $ ceiling (recip lowerB) -- 2^k * lowerB >= 1+ upperBtr = upperB * 2^k -- upper bound of translated poly+ pAbove1Encl = -- p multiplied by 2^k; range in [1,upperBtr] + enclScale maxDegree maxSize (2^k) pEncl+ + -- translate T_1 to domain [0, upperBtr-1] and apply it to x = (pAbove1 - 1):+ -- T'_1(x) = nu * x - 1 where nu = 2/(upperBtr - 1)+ trT1Encl = + enclAddConst (-1) (enclRAScale maxDegree maxSize nu (enclAddConst (-1) pAbove1Encl))+ nu = recip nuInv -- auxiliary constant+ nuInv = (RA.setMinGranularity coeffGr (ERInterval upperBtr upperBtr) - 1) / 2+ + nuPlus1 = nu + 1+ nuInvPlus1 = nuInv + 1+ nuInvDiv2 = nuInv / 2+ + -- define such translated T_i's for all i >= 0:+ trTis =+ enclEvalTs maxDegree maxSize trT1Encl+ + -- construct the result from interval coefficients:+ resEncl = (resLowNeg, resHigh)+ resLowNeg =+ chplScaleUp (2^k) $+ chplScaleUp errScaleDownB $+ scaledTrTisSumLowNeg+ resHigh+ | errScaleUpB > 0 =+ chplScaleUp (2^k) $+ chplScaleUp errScaleUpB $+ scaledTrTisSumHigh+ | otherwise =+ chplScaleUp (2^k) $+ chplAddConstUp errAddUpB scaledTrTisSumHigh++ ERInterval errScaleDownB _ = nuOverNuPlusTauAns + nuOverNuPlusTauAns = (nu / (nu + tauAbs))+ ERInterval _ errScaleUpB = nuOverNuMinusTauAns + nuOverNuMinusTauAns = (nu / (nu - tauAbs)) + ERInterval _ errAddUpB = tauAbsTimesNuInv + tauAbsTimesNuInv = tauAbs * nuInv+ + allInterimsBounded =+ and $ map RA.isBounded [nuOverNuPlusTauAns, nuOverNuMinusTauAns, nuOverNuMinusTauAns]+ + tauAbs = abs tau+ tau = recip tauInv+ + (scaledTrTisSumLowNeg, scaledTrTisSumHigh) =+ foldl1 (+:) $ zipWith scaleTerm c0n trTis+ scaleTerm c trTEncl =+ enclRAScale maxDegree maxSize (c * tau) trTEncl + + -- 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
src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs view
@@ -17,9 +17,12 @@ 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.Field+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 qualified Data.Number.ER.Real.Approx as RA import qualified Data.Number.ER.Real.Approx.Elementary as RAEL@@ -34,540 +37,368 @@ import qualified Data.Map as Map {-|- Approximate the pointwise square root of a polynomial - by another polynomial from below and from above. + Approximate the pointwise exponential of a square root enclosure. -}-chplSqrt ::+enclSqrt :: (B.ERRealBase b, RealFrac 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, ERChebPoly box b)-chplSqrt maxDegree ix p =+enclSqrt maxDegree maxSize ix p = error "ERChebPoly: chplSqrt: not implemented yet" {-|- Approximate the pointwise exponential of a polynomial - by another polynomial from below and from above. + Approximate the pointwise exponential of a polynomial enclosure. -}-chplExp ::+enclExp :: (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) => Int {-^ maximum polynomial degree -} -> - EffortIndex {-^ minimum approx Taylor degree -} -> - ERChebPoly box b ->+ Int {-^ maximum term count -} -> + EffortIndex {-^ used to derive minimum approx Taylor degree -} -> + (ERChebPoly box b, ERChebPoly box b) -> (ERChebPoly box b, ERChebPoly box b)-chplExp maxDegree ix p =--- unsafePrint+enclExp maxDegree maxSize ix pEncl =+-- unsafePrintReturn -- ( --- "chplExp:" ++ --- "\n expM = " ++ show expM +++-- "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 expNear0Dn pNear0Dn = " ++ show (expNear0Dn pNear0Dn) ++--- "\n chplPow maxDegree (expNear0Up pNear0Up) 2000 = " ++ show (chplPow maxDegree (expNear0Up pNear0Up) 2000)--- ) +-- "\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 = "+-- ) -- $ - (expDownwards, expUpwards + valueAtRDnNeg + (chplConst expRUp))+ result 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))+ 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 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+ 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 =- chplConst thisCoeff+ enclRAConst thisCoeff | otherwise =- snd $ chplReduceDegree maxDegree $- (chplConst thisCoeff) + p0 * (expAux p0 (nextDegree + 1) nextCoeff)+ (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- expRUp = erintv_right expR- expR = erExp_R ix (ERInterval r r)- valueAtRDnNeg = - expAux (chplConst r) 1 (B.setGranularity coeffGr (-1))+ -- 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 a polynomial by another polynomial from above. + Approximate the pointwise integer power of an enclosure. -}-chplPow ::- (B.ERRealBase b, Integral i, DomainBox box varid Int, Ord box) => +enclPow ::+ (B.ERRealBase b, RealFrac b, Integral i, DomainBox box varid Int, Ord box) => Int {-^ maximum polynomial degree -} -> - ERChebPoly box b ->+ Int {-^ maximum term count -} -> + (ERChebPoly box b, ERChebPoly box b) -> i ->- ERChebPoly box b-chplPow maxDegree p n+ (ERChebPoly box b, ERChebPoly box b)+ {-^ lower (negated) and upper bound -}+enclPow maxDegree maxSize pEncl n | n == 0 =- chplConst 1+ enclConst 1 | n == 1 =- p + pEncl | even n =- snd $ chplReduceDegree maxDegree $ powHalfN * powHalfN+ powEvenEncl | odd n =- snd $ chplReduceDegree maxDegree $ - p * - (snd $ chplReduceDegree maxDegree $- powHalfN * powHalfN)+ enclMultiply maxDegree maxSize powEvenEncl pEncl where- powHalfN =- chplPow maxDegree p halfN+ powEvenEncl =+ enclMultiply maxDegree maxSize powHalfEncl powHalfEncl + powHalfEncl = + enclPow maxDegree maxSize pEncl halfN halfN = n `div` 2 {-|- Approximate the pointwise natural logarithm of a polynomial - by another polynomial from below and from above. + Approximate the pointwise natural logarithm of an enclosure. -}-chplLog ::+enclLog :: (B.ERRealBase b, RealFrac 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, ERChebPoly box b)-chplLog maxDegree ix p =+enclLog maxDegree maxSize ix p = error "ERChebPoly: chplLog: not implemented yet" {-|- Approximate the pointwise sine of a polynomial - by another polynomial from below and from above.+ Approximate the pointwise sine of an enclosure. Assuming the polynomial range is [-pi/2, pi/2]. -}-chplSine ::+enclSine :: (B.ERRealBase b, RealFrac 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, ERChebPoly box b)-chplSine maxDegree ix p =+enclSine maxDegree maxSize ix pEncl = -- unsafePrint -- (--- "ERChebPoly: sineTaylor: "--- ++ "\n p = " ++ show p+-- "ERChebPoly: enclSine: "+-- ++ "\n pEncl = " ++ show pEncl -- ++ "\n ranLargerEndpoint = " ++ show ranLargerEndpoint--- ++ "\n sineUp = " ++ show sineUp--- ++ "\n sineDown = " ++ show sineDown+-- ++ "\n sineEncl = " ++ show sineEncl -- ) $- (sineDown, sineUp)+ sineEncl where- (sineDown, sineUp) =- boundsAddErr sineErrorBound $- chplThinTimesEncl maxDegree p (sineDownTaylor, sineUpTaylor) - ((sineDownTaylor, sineUpTaylor), - sineErrorTermDegree, - (sineErrorTermCoeffDown, sineErrorTermCoeffUp)) =- sincosTaylorAux True (chplSquare maxDegree p) taylorDegree 1 (one, one)- one = B.setGranularity coeffGr 1+ 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+ case sineErrorBoundRA of + ERInterval lo hi -> hi+ ERIntervalAny -> 1/0 where sineErrorBoundRA = - (ranLargerEndpointRA ^ (sineErrorTermDegree)) * sineErrorTermCoeffRA- sineErrorTermCoeffRA =- ERInterval sineErrorTermCoeff sineErrorTermCoeff- sineErrorTermCoeff =- max (abs sineErrorTermCoeffDown) (abs sineErrorTermCoeffUp)+ (ranLargerEndpointRA ^ sineErrorTermDegree) * sineErrorTermCoeffHighRA+ sineErrorTermCoeffHighRA =+ snd $ RA.bounds $ abs sineErrorTermCoeffRA ranLargerEndpointRA = ERInterval ranLargerEndpoint ranLargerEndpoint ranLargerEndpoint =- max (abs ranLO) (abs ranHI)- ranLO = negate $ chplUpperBoundAffine ix (-p)- ranHI = chplUpperBoundAffine ix p+ max (abs ranLowB) (abs ranHighB)+ (ranLowB, ranHighB) = enclBounds ix pEncl taylorDegree = effIx2int $ ix `div` 3 coeffGr = effIx2gran $ ix -boundsAddErr errB (pLO, pHI) =- (pLO `plusDown` (- errPoly), pHI + errPoly)- where- errPoly = chplConst errB- {-|- Approximate the pointwise sine of a polynomial - by another polynomial from below and from above.+ Approximate the pointwise cosine of an enclosure. Assuming the polynomial range is [-pi/2, pi/2]. -}-chplCosine ::+enclCosine :: (B.ERRealBase b, RealFrac 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, ERChebPoly box b)-chplCosine maxDegree ix p =+enclCosine maxDegree maxSize ix pEncl = -- unsafePrint -- ( -- "ERChebPoly: chplCosine: "--- ++ "\n p = " ++ show p+-- ++ "\n pEncl = " ++ show pEncl -- ++ "\n ranLargerEndpoint = " ++ show ranLargerEndpoint--- ++ "\n cosineUp = " ++ show cosineUp--- ++ "\n cosineDown = " ++ show cosineDown+-- ++ "\n cosineEncl = " ++ show cosineEncl -- ) $- (cosineDown, cosineUp)+ cosineEncl where- (cosineDown, cosineUp) =- boundsAddErr cosineErrorBound $- (cosineDownTaylor, cosineUpTaylor) - ((cosineDownTaylor, cosineUpTaylor), - cosineErrorTermDegree, - (cosineErrorTermCoeffDown, cosineErrorTermCoeffUp)) =- sincosTaylorAux True (chplSquare maxDegree p) taylorDegree 0 (one, one)- one = B.setGranularity coeffGr 1+ 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+ case cosineErrorBoundRA of + ERInterval lo hi -> hi+ ERIntervalAny -> 1/0 where- cosineErrorBoundRA =- (ranLargerEndpointRA ^ (cosineErrorTermDegree)) * cosineErrorTermCoeffRA- cosineErrorTermCoeffRA =- ERInterval cosineErrorTermCoeff cosineErrorTermCoeff- cosineErrorTermCoeff =- max (abs cosineErrorTermCoeffDown) (abs cosineErrorTermCoeffUp)+ cosineErrorBoundRA = + (ranLargerEndpointRA ^ cosineErrorTermDegree) * cosineErrorTermCoeffHighRA+ cosineErrorTermCoeffHighRA =+ snd $ RA.bounds $ abs cosineErrorTermCoeffRA ranLargerEndpointRA = ERInterval ranLargerEndpoint ranLargerEndpoint ranLargerEndpoint =- max (abs ranLO) (abs ranHI)- ranLO = negate $ chplUpperBoundAffine ix (-p)- ranHI = chplUpperBoundAffine ix p+ max (abs ranLowB) (abs ranHighB)+ (ranLowB, ranHighB) = enclBounds ix pEncl taylorDegree = effIx2int $ ix `div` 3 coeffGr = effIx2gran $ ix sincosTaylorAux :: (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>- Bool -> + 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 -} ->- (b,b) -> + ERInterval b {-^ the coefficient of the term now being constructed -} -> ((ERChebPoly box b, ERChebPoly box b), Int,- (b,b))+ 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 resultPositive pSquares@(pSquareDown, pSquareUp) - maxDegree thisDegree (thisCoeffDown, thisCoeffUp)- | nextDegree > maxDegree =--- unsafePrint--- (--- "ERChebPoly: chplSine: taylorAux: "--- ++ "\n thisCoeff = " ++ show thisCoeff--- ++ "\n nextDegree = " ++ show nextDegree--- )- ((thisCoeffDownP, thisCoeffUpP), nextDegree, (nextCoeffDown, nextCoeffUp))- | otherwise =--- unsafePrint--- (--- "ERChebPoly: chplSine: taylorAux: "--- ++ "\n thisCoeff = " ++ show thisCoeff--- ++ "\n nextDegree = " ++ show nextDegree--- ++ "\n errorTermCoeff = " ++ show errorTermCoeff--- ++ "\n errorTermDegree = " ++ show errorTermDegree--- )- ((resultDown, resultUp), errorTermDegree, errorTermCoeffs) +sincosTaylorAux + maxDegree maxSize resultPositive pSqrEncl tayDegree + thisDegree thisCoeffRA =+ sct thisDegree thisCoeffRA where- thisCoeffDownP = chplConst thisCoeffDown- thisCoeffUpP = chplConst thisCoeffUp- resultDown- | resultPositive = - -- ie rest's ideal value is negative and thisCoeff is positive- chplReduceDegreeDown maxDegree $- thisCoeffDownP `plusDown` (pSquareUp `timesDown` restDown)- | otherwise =- -- ie rest's ideal value is positive and thisCoeff is negative- chplReduceDegreeDown maxDegree $- thisCoeffDownP `plusDown` (pSquareDown `timesDown` restDown)- resultUp- | resultPositive = - -- ie rest's ideal value is negative and thisCoeff is positive- chplReduceDegreeUp maxDegree $- thisCoeffUpP `plusUp` (pSquareDown `timesUp` restUp)- | otherwise =- -- ie rest's ideal value is positive and thisCoeff is negative- chplReduceDegreeUp maxDegree $- thisCoeffUpP `plusUp` (pSquareUp `timesUp` restUp)- ((restDown, restUp), errorTermDegree, errorTermCoeffs) =- sincosTaylorAux (not resultPositive) pSquares maxDegree nextDegree (nextCoeffDown, nextCoeffUp)- nextDegree = thisDegree + 2- nextCoeffUp- | resultPositive = - thisCoeffDown / nextCoeffDenominator -- positive / negative- | otherwise = - thisCoeffUp / nextCoeffDenominator -- negative / negative- nextCoeffDown - | resultPositive = - thisCoeffUp `divDown` nextCoeffDenominator -- positive / negative- | otherwise = - thisCoeffDown `divDown` nextCoeffDenominator -- negative / negative- nextCoeffDenominator =- fromInteger $ toInteger $ negate $ nextDegree * (nextDegree - 1)- divDown a b = negate $ a / (negate b) + 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 natural logarithm of a polynomial - by another polynomial from below and from above. + Approximate the pointwise arcus tangens of an enclosure. -}-chplAtan ::+enclAtan :: (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) => Int {-^ maximum polynomial degree -} -> - EffortIndex {-^ ?? -} -> - ERChebPoly box b ->+ 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)-{- atan using Euler's 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 + ...)))))+{- 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 + (x^2*2n/(2n + 1)), 1 + x^2]+ [1, 1 + x^2] -}-chplAtan maxDegree ix p - | avoidingDivBy0 = --- unsafePrint--- (--- "ERChebPoly.Elementary: chplAtan: "--- ++ "\n maxDegree = " ++ show maxDegree--- ++ "\n p = " ++ show p--- ++ "\n pSquareDn = " ++ show pSquareDn--- ++ "\n pSquareUp = " ++ show pSquareUp--- ++ "\n pOverPSquarePlus1Dn = " ++ show pOverPSquarePlus1Dn--- ++ "\n pOverPSquarePlus1Up = " ++ show pOverPSquarePlus1Up--- ++ "\n preresDn = " ++ show preresDn--- ++ "\n preresUp = " ++ show preresUp--- ++ "\n resDn = " ++ show resDn--- ++ "\n resUp = " ++ show resUp--- )- (resDn, resUp)- | otherwise =- (chplConst (-2), chplConst 2) -- this is always safe... +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- avoidingDivBy0 = - (chplUpperBoundAffine ix (- pSquarePlus1Dn) < 0)- &&- (chplUpperBoundAffine ix (- pSquarePlus1Up) < 0)- resDn = - negate $- chplMaxUp maxDegree - (chplReduceDegreeUp maxDegree $ - pOverPSquarePlus1DnNeg `timesUp` preresDn) -- beware: pOverPSquarePlus1Dn can be negative- (chplReduceDegreeUp maxDegree $- pOverPSquarePlus1DnNeg `timesUp` preresUp)+ (pLowerBound, pUpperBound) = enclBounds ix pEncl+ onePlusSqrtOnePlusPSquare =+ enclAddConst 1 $+ enclSqrt maxDegree maxSize ix pSquarePlus1Encl+ avoidingDivBy0 =+ lower1 > 0 && lower2 > 0 where- pOverPSquarePlus1DnNeg = negate pOverPSquarePlus1Dn- resUp = - chplMaxUp maxDegree - (chplReduceDegreeUp maxDegree $- pOverPSquarePlus1Up `timesUp` preresDn) -- beware: pOverPSquarePlus1Up can be negative- (chplReduceDegreeUp maxDegree $- pOverPSquarePlus1Up `timesUp` preresUp)- (preresDn, preresUp) = seriesDnUp termsCount 2- termsCount = max 0 $ ix `div` 3+ (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- seriesDnUp termsCount coeffBase + series termsCount coeffBase | termsCount > 0 =- (- chplReduceDegreeDown maxDegree $- 1 `plusDown` - (pSquareOverPSquarePlus1Dn -- >=0 - `timesDown` (chplConst coeffDn) -- >=0 - `timesDown` restDn -- >=0- )- ,- chplReduceDegreeUp maxDegree $- 1 `plusUp`- (pSquareOverPSquarePlus1Up -- >=0 - `timesUp` (chplConst coeffUp) -- >=0 - `timesUp` restUp -- >=0- )- )+ enclAddConst 1 $+ enclRAScale maxDegree maxSize coeff $+ enclMultiply maxDegree maxSize + pSquareOverPSquarePlus1Encl $+ series (termsCount - 1) (coeffBase + 2) | otherwise =- (- 1 `plusDown` (pSquareDn `timesDown` (chplConst coeffDn)) -- both >=0- ,- 1 `plusUp` pSquareUp- )- where- (restDn, restUp) = seriesDnUp (termsCount - 1) (coeffBase + 2)- coeffUp = coeffBaseB / (coeffBaseB `plusDown` 1)- coeffDn = negate $ coeffBaseB / (negate $ coeffBaseB `plusUp` 1)- coeffBaseB = B.setMinGranularity gran $ fromInteger coeffBase- (pSquareDn, pSquareUp) = chplSquare maxDegree p- pSquarePlus1Dn = pSquareDn `plusDown` 1- pSquarePlus1Up = pSquareUp `plusUp` 1- recipPSquarePlus1Dn = chplRecipDn maxDegree ix pSquarePlus1Up- recipPSquarePlus1Up = chplRecipUp maxDegree ix pSquarePlus1Dn--- -- safely compute the square of an enclosure:--- pSquareDn = chplMinDn m pUpTDnpUp (chplMinDn m pDnTDnpUp pDnTDnpDn)--- pSquareUp = chplMaxUp m pUpTUppUp (chplMaxUp m pDnTUppUp pDnTUppDn) --- pUpTDnpUp = pUp `timesDown` pUp--- pDnTDnpUp = pDn `timesDown` pUp--- pDnTDnpDn = pDn `timesDown` pDn--- pUpTUppUp = pUp `timesUp` pUp--- pDnTUppUp = pDn `timesUp` pUp--- pDnTUppDn = pDn `timesUp` pDn--- mMinus1 = m - 1- pSquareOverPSquarePlus1Up = - pSquareUp `timesUp` recipPSquarePlus1Up -- both >=0- pSquareOverPSquarePlus1Dn = - pSquareDn `timesDown` recipPSquarePlus1Dn -- both >=0 (one enclosure may dip below 0, not a problem)--- negate $--- chplMaxUp maxDegree--- (pSquareDnNeg `timesUp` recipPSquarePlus1Up) -- beware: pSquareDn may dip below 0--- (pSquareDnNeg `timesUp` recipPSquarePlus1Dn)--- where--- pSquareDnNeg = negate pSquareDn- pOverPSquarePlus1Up =- chplMaxUp maxDegree - (p `timesUp` recipPSquarePlus1Up)- (p `timesUp` recipPSquarePlus1Dn) -- beware: p can be negative- pOverPSquarePlus1Dn =- negate $- chplMaxUp maxDegree- (pn `timesUp` recipPSquarePlus1Up) -- beware: pn can be positive- (pn `timesUp` recipPSquarePlus1Dn)+ enclAddConst 1 $+ (chplConst 0, pSquareHigh) where- pn = negate p--chplRecipDn m i = fst . chplRecip m i-chplRecipUp m i = snd . chplRecip m i--{-|- 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@(ERChebPoly pCoeffs)- | pIsConst = - (chplConst $ - (recip (- pConst)), chplConst $ recip pConst)- | 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- - (pIsConst, pConst) = - case chplGetConst p of- Nothing -> (False, error "ChebyshevBase.Polynom.Elementary.chplRecip")- Just coeff -> (True, coeff)- - tauDegree = effIx2int (max 2 $ 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 + coeff = coeffBase / (coeffBase + 1) - -- 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+ 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
+ src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Enclosure.hs view
@@ -0,0 +1,311 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE UndecidableInstances #-}+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Enclosure+ 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 selected operations working on pairs+ of polynomials understood as function enclosures.+ These are needed to define composition and some elementary operations.+-}+module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Enclosure++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.Ring+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Eval++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.Real.Approx.Interval+import qualified Data.Number.ER.Real.Approx as RA+import Data.Number.ER.Misc++import qualified Data.Map as Map++enclThin ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) =>+ ERChebPoly box b ->+ (ERChebPoly box b, ERChebPoly box b)+enclThin p =+ (chplNeg p, p)++enclConst ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) =>+ b ->+ (ERChebPoly box b, ERChebPoly box b)+enclConst c =+ (chplConst (-c), chplConst c)++enclBounds ix (ln, h) =+ (negate $ chplUpperBound ix ln, chplUpperBound ix h)++enclEval e@(ln, h) pt + | lB > hB =+ unsafePrintReturn+ (+ "ERChebPoly: enclEval: inverted result:"+ ++ "\n h = " ++ show h + ++ "\n ln = " ++ show ln + ++ "\n result = "+ )+ result+ | otherwise = result+ where+ result = ERInterval lB hB+ lB = negate $ chplEvalUp ln pt+ hB = chplEvalUp h pt++enclEvalInner (ln, h) pt =+-- normaliseERInterval $+ ERInterval + (negate $ chplEvalDown ln pt)+ (chplEvalDown h pt)++enclRAEval e@(ln, h) pt =+ result + where+ result = lRA RA.\/ hRA+ lRA = fst $ RA.bounds $ negate $ chplRAEval (\b -> ERInterval b b) ln pt+ hRA = snd $ RA.bounds $ chplRAEval (\b -> ERInterval b b) h pt++enclRAEvalInner e@(ln, h) pt =+-- unsafePrintReturn+-- (+-- "ERChebPoly: enclRAEvalInner: "+-- ++ "\n lB = " ++ show lB+-- ++ "\n hB = " ++ show hB+-- ++ "\n result = "+-- )+ result + where+ result =+-- normaliseERInterval $ + ERInterval lB hB+ lB = + case negate $ chplRAEval (\b -> ERInterval b b) ln pt of+ ERInterval _ lB -> lB+ hB = + case chplRAEval (\b -> ERInterval b b) h pt of+ ERInterval hB _ -> hB++enclAddErr errB (pLowNeg, pHigh) =+ (chplAddConstUp errB pLowNeg, chplAddConstUp errB pHigh)+++enclRAConst ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) =>+ (ERInterval b) ->+ (ERChebPoly box b, ERChebPoly box b)+enclRAConst (ERInterval lo hi) = (chplConst (-lo), chplConst hi)+enclRAConst ERIntervalAny = (chplConst (-1/0), chplConst (1/0))++enclReduceDegree maxDegree (pLowNeg, pHigh) =+ (chplReduceDegreeUp maxDegree pLowNeg, chplReduceDegreeUp maxDegree pHigh) + +enclReduceSize maxSize (pLowNeg, pHigh) =+ (chplReduceTermCountUp maxSize pLowNeg, chplReduceTermCountUp maxSize pHigh) + +enclAddConst c (pLowNeg, pHigh) =+ (chplAddConstUp (-c) pLowNeg, chplAddConstUp c pHigh)++enclNeg (pLowNeg, pHigh) = (pHigh, pLowNeg)++(p1LowNeg, p1High) +: (p2LowNeg, p2High) = + (p1LowNeg +^ p2LowNeg, p1High +^ p2High)+ +(p1LowNeg, p1High) -: (p2LowNeg, p2High) =+ (p1LowNeg +^ p2High, p1High +^ p2LowNeg)++enclMultiply ::+ (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>+ Int {-^ maximum polynomial degree -} -> + Int {-^ maximum term count -} -> + (ERChebPoly box b, ERChebPoly box b) -> + (ERChebPoly box b, ERChebPoly box b) ->+ (ERChebPoly box b, ERChebPoly box b)+enclMultiply maxDegr maxSize (ln1, h1) (ln2, h2) =+ enclReduceSize maxSize $+ enclReduceDegree maxDegr $+ case (ln1UpperBound <= 0, h1UpperBound <= 0, ln2UpperBound <= 0, h2UpperBound <= 0) of+ (True, _, True, _) -> -- both non-negative+ (l1l2Neg, h1h2)+ (_, True, _, True) -> -- both non-positive+ (h1h2Neg, l1l2)+ (True, _, _, True) -> -- first non-negative, second non-positive+ (h1l2Neg, l1h2)+ (_, True, True, _) -> -- first non-positive, second non-negative+ (l1h2Neg, l1h2)+ _ -> -- one of both may be crossing zero+ (+ (h1h2Neg `maxP` l1l2Neg) `maxP` (h1l2Neg `maxP` l1h2Neg)+ ,+ (h1h2 `maxP` l1l2) `maxP` (h1l2 `maxP` l1h2)+ )+ where+ ln1UpperBound = chplUpperBound ix ln1+ ln2UpperBound = chplUpperBound ix ln2+ h1UpperBound = chplUpperBound ix h1+ h2UpperBound = chplUpperBound ix h2+ ix = 10+ maxP = chplMaxUp maxDegr maxSize+ + h1h2 = h1 *^ h2+ h1h2Neg = (chplNeg h1) *^ h2+ l1l2 = ln1 *^ ln2+ l1l2Neg = (chplNeg ln1) *^ ln2+ h1l2 = h1 *^ (chplNeg ln2)+ h1l2Neg = h1 *^ ln2+ l1h2 = (chplNeg ln1) *^ h2+ l1h2Neg = ln1 *^ h2+++enclSquare ::+ (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>+ Int {-^ maximum polynomial degree -} -> + Int {-^ maximum term count -} -> + (ERChebPoly box b, ERChebPoly box b) ->+ (ERChebPoly box b, ERChebPoly box b)+enclSquare maxDegr maxSize (ln, h)+ {-+ formula:+ (ln, h)^2 =+ ( minUp( 0, maxUp( - ln *. ln, - h *. h)), maxUp(ln *^ ln, h *^ h) )+ -}+-- | minZeroHelps = + = (minZeroMaxNegSq, maxSq)+-- | otherwise =+-- (maxNegSq, maxSq)+ where+ maxSq = maxP ln2Up h2Up+ maxNegSq = maxP (chplNeg ln2Down) (chplNeg h2Down)+ minZeroMaxNegSq = chplNonposUp maxDegr maxSize maxNegSq +-- minZeroHelps =+-- (maxNegSqUpperB > 0) && (minZeroMaxNegSqUpperB / maxNegSqUpperB < 1/2)+-- maxNegSqUpperB = chplUpperBound 10 maxNegSq+-- minZeroMaxNegSqUpperB = chplUpperBound 10 minZeroMaxNegSq+ + (ln2Down, ln2Up, _) = chplMultiply ln ln+ (h2Down, h2Up, _) = chplMultiply h h+ +-- reduceDegrSize = reduceSize maxSize . reduceDegree maxDegr+ maxP = chplMaxUp maxDegr maxSize+ + ++ +{-| + Multiply an enclosure by a scalar + assuming the enclosure is non-negative on the whole unit domain.+-} +enclScaleNonneg ::+ (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>+ b {-^ ratio to scale by -} -> + (ERChebPoly box b, ERChebPoly box b) -> + (ERChebPoly box b, ERChebPoly box b)+enclScaleNonneg ratio pEncl@(ln, h) =+ (ln *^ pRatio, h *^ pRatio)+ where+ pRatio = chplConst ratio++{-| + Multiply an enclosure by a scalar.+-} +enclScale ::+ (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>+ Int {-^ maximum polynomial degree -} -> + Int {-^ maximum term count -} ->+ b {-^ ratio to scale by -} -> + (ERChebPoly box b, ERChebPoly box b) -> + (ERChebPoly box b, ERChebPoly box b)+enclScale maxDegree maxSize ratio pEncl =+ enclMultiply maxDegree maxSize pEncl (enclConst ratio) ++enclRAScale ::+ (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) => + Int {-^ maximum polynomial degree -} -> + Int {-^ maximum term count -} ->+ (ERInterval b) -> + (ERChebPoly box b, ERChebPoly box b) ->+ (ERChebPoly box b, ERChebPoly box b)+enclRAScale maxDegree maxSize ra pEncl =+ enclMultiply maxDegree maxSize pEncl (enclRAConst ra) ++{-|+ Multiply a polynomial by a scalar interval, returning an enclosure.+-} +chplScaleRA :: + (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>+ Int {-^ maximum polynomial degree -} -> + Int {-^ maximum term count -} -> + ERInterval b {-^ lower and upper bounds on the ratio to scale by -} -> + ERChebPoly box b ->+ (ERChebPoly box b, ERChebPoly box b)+chplScaleRA maxDegr maxSize (ERIntervalAny) p = enclRAConst ERIntervalAny+chplScaleRA maxDegr maxSize (ERInterval ratioDown ratioUp) p =+ (scaledPDownNeg, scaledPUp)+ where+ (scaledPDownNeg, scaledPUp) =+ enclMultiply maxDegr maxSize + (chplNeg p, p) (chplConst (- ratioDown), chplConst ratioUp)++chplScaleRADown m n r = chplNeg . fst . chplScaleRA m n r+chplScaleRAUp m n r = snd . chplScaleRA m n r++{-|+ Evaluate the Chebyshev polynomials of the first kind+ applied to a given polynomial, yielding a list of polynomial enclosures. +-}+enclEvalTs ::+ (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, ERChebPoly box b) {-^ bounds of a polynomial enclosure to evaluate -} ->+ [(ERChebPoly box b, ERChebPoly box b)]+enclEvalTs maxDegree maxSize p1@(pLowNeg, pHigh) =+ chebyIterate (enclConst 1) p1+ where+ chebyIterate pNm2 pNm1 =+ pNm2 : (chebyIterate pNm1 pN)+ where+ pN = + (enclScale maxDegree maxSize 2 $ + enclMultiply maxDegree maxSize p1 pNm1) + -: pNm2++{-|+ Multiply a polynomial by an enclosure using min/max+-}+enclThinTimes ::+ (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) => + Int {-^ maximum polynomial degree -} -> + Int {-^ maximum term count -} -> + ERChebPoly box b ->+ (ERChebPoly box b, ERChebPoly box b) ->+ (ERChebPoly box b, ERChebPoly box b)+enclThinTimes maxDegree maxSize p1 (p2LowNeg, p2High) =+ (prodLowNeg, prodHigh)+ where+ prodHigh =+ chplMaxUp maxDegree maxSize+ (p1 *^ p2High)+ (p1n *^ p2LowNeg) -- beware: p1 can cross zero+ prodLowNeg =+ chplMaxUp maxDegree maxSize+ (p1n *^ p2High)+ (p1 *^ p2LowNeg)+ p1n = chplNeg p1++
src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Eval.hs view
@@ -1,7 +1,7 @@ {-# LANGUAGE FlexibleContexts #-} {-| Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Eval- Description : (internal) evaluation of polynomials+ Description : (internal) evaluation of polynomials at a point Copyright : (c) 2007-2008 Michal Konecny License : BSD3 @@ -17,12 +17,12 @@ where import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic-import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Field 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.Real.Approx.Interval import Data.Number.ER.Misc import qualified Data.Map as Map@@ -32,46 +32,34 @@ -} chplEval :: (B.ERRealBase b, DomainBox box varid Int, Ord box, - DomainBoxMappable boxb boxbb varid b [(b,b)]) => - boxb -> + DomainBoxMappable boxb boxbb varid b [ERInterval b]) => ERChebPoly box b ->+ boxb -> (b, b)-chplEval vals (ERChebPoly coeffs) =- (foldl plusDown 0 termValsLo, foldl plusUp 0 termValsHi)+chplEval (ERChebPoly coeffs) vals =+ case resultRA of+ ERInterval low high -> (low, high)+ ERIntervalAny -> (-1/0,1/0)+ ERIntervalEmpty -> (1/0, -1/0) where- (termValsLo, termValsHi) =- unzip $ map evalTerm $ Map.toList coeffs+ resultRA =+ sum $ map evalTerm $ Map.toList coeffs evalTerm (term, c) =- (foldl timesDown c valsLo, foldl timesUp c valsHi)- where- (valsLo, valsHi) = - unzip $ map evalVar $ DBox.toList term+ foldl (*) (ERInterval c c) $ map evalVar $ DBox.toList term evalVar (varID, degree) =- (DBox.lookup "ERChebPoly.Eval: chplEval" varID valsDegrees) !! degree+ (DBox.lookup "ERChebPoly.Eval: chplEval: " varID valsDegrees) !! degree valsDegrees =- DBox.map chebyEvalTsRoundDownUp vals+ DBox.map (chebyEvalTsExact . \a->(ERInterval a a)) $ vals chplEvalDown, chplEvalUp :: (B.ERRealBase b, DomainBox box varid Int, Ord box, - DomainBoxMappable boxb boxbb varid b [(b,b)]) => - boxb -> + DomainBoxMappable boxb boxbb varid b [ERInterval b]) => ERChebPoly box b ->+ boxb -> b chplEvalUp pt = snd . chplEval pt chplEvalDown pt = fst . chplEval pt -chebyEvalTsRoundDownUp ::- (Num v) =>- v -> [(v,v)]-chebyEvalTsRoundDownUp val =- chebyIterate (1,1) (val, val)- where- chebyIterate tNm2@(tNm2Down, tNm2Up) tNm1@(tNm1Down, tNm1Up) =- tNm2 : (chebyIterate tNm1 (tNDown, tNUp))- where- tNUp = 2 * val * tNm1Up - tNm2Down - tNDown = ((2 * val) `timesDown` tNm1Down) - tNm2Up - chebyEvalTsExact :: (Num v) => v -> [v] @@ -86,16 +74,16 @@ {-| Evaluate a polynomial at a real number approximation -}-chplEvalApprox ::+chplRAEval :: (B.ERRealBase b, RA.ERApprox ra, DomainBox box varid Int, Ord box, DomainBoxMappable boxra boxras varid ra [ra], DomainIntBox boxra varid ra) => (b -> ra) -> - boxra -> ERChebPoly box b ->+ boxra -> ra-chplEvalApprox b2ra vals (ERChebPoly coeffs) =+chplRAEval b2ra (ERChebPoly coeffs) vals = sum $ map evalTerm $ Map.toList coeffs where evalTerm (term, c) =@@ -109,16 +97,16 @@ Substitute several variables in a polynomial with real number approximations, rounding downwards and upwards. -}-chplPartialEvalApprox ::+chplPartialRAEval :: (B.ERRealBase b, RA.ERApprox ra, DomainBox box varid Int, Ord box,- DomainBoxMappable boxra boxras varid ra [ra], + DomainBoxMappable boxra boxras varid ra [ra], DomainIntBox boxra varid ra) => (ra -> (b,b)) ->- boxra -> ERChebPoly box b ->+ boxra -> (ERChebPoly box b, ERChebPoly box b)-chplPartialEvalApprox ra2endpts substitutions (ERChebPoly coeffs) =+chplPartialRAEval ra2endpts (ERChebPoly coeffs) substitutions = (ERChebPoly $ Map.insertWith plusDown chplConstTermKey (- corr) coeffsSubstDown, ERChebPoly $ Map.insertWith plusUp chplConstTermKey corr coeffsSubstUp) where@@ -147,44 +135,3 @@ (DBox.lookup "ERChebPoly.Eval: chplPartialEvalApprox: " varID valsDegrees) !! degree valsDegrees = DBox.map chebyEvalTsExact substitutions- --{-|- Compose two polynomials, rounding upwards- provided the second polynomial maps [-1,1] into [-1,1].--}-chplCompose ::- (B.ERRealBase b, DomainBox box varid Int, Ord box) =>- Int ->- ERChebPoly box b ->- Map.Map varid (ERChebPoly box b) - {-^ variable to substitute, polynomial to substitute -} ->- (ERChebPoly box b, ERChebPoly box b)-chplCompose maxDegree p@(ERChebPoly coeffs) substitutions =- (foldl plusDown 0 termValsLo, foldl plusUp 0 termValsHi)- where- (termValsLo, termValsHi) =- unzip $ map evalTerm $ Map.toList coeffs- evalTerm (term, c) =- (foldl timesDown cPoly valsLo, foldl timesUp cPoly valsHi)- where- cPoly = chplConst c- (valsLo, valsHi) = - unzip $ map evalVar $ DBox.toList term- evalVar (varID, degree) =- case Map.lookup varID substDegrees of- Nothing ->- (varPoly, varPoly)- Just pvDegrees ->- pvDegrees !! degree- where- varPoly = - ERChebPoly $ Map.singleton (DBox.singleton varID degree) 1- substDegrees =- Map.map mkPVDegrees substitutions- mkPVDegrees pv =- map - (mapPair - (chplReduceDegreeDown maxDegree, - chplReduceDegreeUp maxDegree)) $ - chebyEvalTsRoundDownUp pv
− src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Field.hs
@@ -1,228 +0,0 @@-{-# 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" -
src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Integration.hs view
@@ -19,12 +19,13 @@ 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.Ring import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds 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.Real.Approx.Interval import Data.Number.ER.Misc import qualified Data.Map as Map@@ -49,15 +50,17 @@ varid {-^ variable to integrate by -} -> ERChebPoly box b -> (ERChebPoly box b, ERChebPoly box b)-chplIntegrate x (ERChebPoly coeffs) =--- unsafePrint+chplIntegrate x p@(ERChebPoly coeffs) =+-- unsafePrintReturn -- ( -- "ERChebPoly: integrate:"--- ++ "\n pNp1Down = " ++ chplShow True pNp1Down --- ++ "\n pNm1Up = " ++ chplShow True pNm1Up +-- ++ "\n p = " ++ show p+-- ++ "\n result = " -- )- (chplNormaliseDown $ pNp1Down - pNm1Up, - chplNormaliseUp $ pNp1Up - pNm1Down)+ (pNp1Down -. pNm1Up, + pNp1Up -^ pNm1Down)+-- (chplRemoveZeroTermsDown $ pNp1Down - pNm1Up, +-- chplRemoveZeroTermsUp $ pNp1Up - pNm1Down) where pNp1Up = ERChebPoly $ @@ -84,12 +87,13 @@ | n == 0 = ((termKeyNp1, coeff):prevTerms, prevErr) | n == 1 =- ((termKeyNm1, coeff0Up):(termKeyNp1, coeffNp1Up):prevTerms, prevErr + coeff0Err + coeffNp1Err)+ ((termKeyN0, coeff0Up):(termKeyNp1, coeffNp1Up):prevTerms, prevErr + coeff0Err + coeffNp1Err) | otherwise = ((termKeyNp1, coeffNp1Up):prevTerms, prevErr + coeffNp1Err) where termKeyNp1 = DBox.insert x (n + 1) termKey termKeyNm1 = DBox.insert x (n - 1) termKey + termKeyN0 = DBox.delete x termKey n = DBox.findWithDefault 0 x termKey coeffNp1Err = coeffNp1Up - coeffNp1Down coeffNp1Up = coeff / (2*nB + 2)@@ -100,7 +104,7 @@ coeff0Err = coeff0Up - coeff0Down cfNm1 (prevTerms, prevErr) (termKey, coeff) | n == 0 || n == 1 =- ((chplConstTermKey, 0):prevTerms, prevErr)+ (prevTerms, prevErr) | otherwise = ((termKeyNm1, coeffNm1Up):prevTerms, prevErr + coeffNm1Err) where@@ -111,59 +115,53 @@ nB = fromInteger $ toInteger n coeffNm1Err = coeffNm1Up - coeffNm1Down -{-|- measure the volume between a polynomial and the zero axis on [-1,1]^n--}-chplVolumeAboveZero ::- (B.ERRealBase b, DomainBox box varid Int, Ord box, - DomainBoxMappable boxb boxbb varid b [(b,b)]) =>- [varid] ->- ERChebPoly box b ->- (b,b)-chplVolumeAboveZero vars p@(ERChebPoly coeffs) =--- unsafePrint ("chplVolumeAboveZero: returning:" ++ show result) $--- unsafePrint ("chplVolumeAboveZero: vars = " ++ show vars) $- result- where- result = - (- (integUpAtOddCorners - integDownAtEvenCorners), integUpAtEvenCorners - integDownAtOddCorners)- integUpAtEvenCorners = sumUp $ map (\pt -> chplEvalUp pt integUp) evenCorners- integUpAtOddCorners = sumUp $ map (\pt -> chplEvalUp pt integUp) oddCorners - integDownAtEvenCorners = sumDown $ map (\pt -> chplEvalDown pt integDown) evenCorners - integDownAtOddCorners = sumDown $ map (\pt -> chplEvalDown pt integDown) oddCorners- evenCorners = map (DBox.fromList) evenCornersL- oddCorners = map (DBox.fromList) oddCornersL- (evenCornersL, oddCornersL) =- allPairsCombinationsEvenOdd $ zip vars $ repeat (1,-1)- integUp = integrateByAllVars snd p vars- integDown = integrateByAllVars fst p vars- integrateByAllVars pick p [] = p- integrateByAllVars pick p (x : xs) =- integrateByAllVars pick ip xs- where- ip = pick $ chplIntegrate x p--- vars = chplGetVars p+--{-|+-- measure the volume between a polynomial and the zero axis on [-1,1]^n+---}+--chplVolumeAboveZero ::+-- (B.ERRealBase b, DomainBox box varid Int, Ord box, +-- DomainBoxMappable boxb boxbb varid b [ERInterval b]) =>+-- [varid] ->+-- ERChebPoly box b ->+-- (b,b)+--chplVolumeAboveZero vars p@(ERChebPoly coeffs) =+---- unsafePrint ("chplVolumeAboveZero: returning:" ++ show result) $+---- unsafePrint ("chplVolumeAboveZero: vars = " ++ show vars) $+-- result+-- where+-- result = +-- (- (integUpAtOddCorners - integDownAtEvenCorners), integUpAtEvenCorners - integDownAtOddCorners)+-- integUpAtEvenCorners = sumUp $ map (chplEvalUp integUp) evenCorners+-- integUpAtOddCorners = sumUp $ map (chplEvalUp integUp) oddCorners +-- integDownAtEvenCorners = sumDown $ map (chplEvalDown integDown) evenCorners +-- integDownAtOddCorners = sumDown $ map (chplEvalDown integDown) oddCorners+-- evenCorners = map (DBox.fromList) evenCornersL+-- oddCorners = map (DBox.fromList) oddCornersL+-- (evenCornersL, oddCornersL) =+-- allPairsCombinationsEvenOdd $ zip vars $ repeat (1,-1)+-- integUp = integrateByAllVars snd p vars+-- integDown = integrateByAllVars fst p vars+-- integrateByAllVars pick p [] = p+-- integrateByAllVars pick p (x : xs) =+-- integrateByAllVars pick ip xs+-- where+-- ip = pick $ chplIntegrate x p+---- vars = chplGetVars p - +-- --{-|--- Calculate approximations to the Chebyshev nodes.+-- Differentiate a polynomial using one of its variables. +-- +-- This is not implemented yet and will probably never be needed+-- because differentiation is not a computable operator+-- and thus we have to rely on automatic differentiation+-- when we need derivative enclosures. ---}---chebNodes ::--- (B.ERRealBase b) =>--- Granularity ->--- [[b]] -- ^ ith element is the ordered list of ith order Chebyshev nodes ---chebNodes gran =--- error "ERChebPoly: chebNodes: not implemented yet"- - -{-|- Differentiate a polynomial using one of its variables. --}-chplDifferentiate ::- (B.ERRealBase b, DomainBox box varid Int, Ord box) => - ERChebPoly box b ->- varid {-^ variable to differentiate over -} ->- ERChebPoly box b-chplDifferentiate (ERChebPoly coeffs) varName =- errorModule "chplDifferentiate: not implemented yet"+--chplDifferentiate ::+-- (B.ERRealBase b, DomainBox box varid Int, Ord box) => +-- ERChebPoly box b ->+-- varid {-^ variable to differentiate over -} ->+-- ERChebPoly box b+--chplDifferentiate (ERChebPoly coeffs) varName =+-- errorModule "chplDifferentiate: not implemented yet"
+ src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Reduce.hs view
@@ -0,0 +1,85 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE UndecidableInstances #-}+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Reduce+ Description : (internal) uniformly roudned polynomial size reductions + 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 pointwise rounding uniform over the whole unit domain.+-}++module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Reduce ++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.List as List+import qualified Data.Map as Map++chplReduceTermCount ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) =>+ Int -> + ERChebPoly box b -> + (ERChebPoly box b, ERChebPoly box b)+chplReduceTermCount maxTermCount p@(ERChebPoly coeffs) + | currentCount <= maxTermCount = (p,p)+ | otherwise =+ (ERChebPoly lessCoeffsDown, ERChebPoly lessCoeffsUp)+ where+ currentCount = chplCountTerms p+ lessCoeffsDown =+ Map.insertWith plusDown chplConstTermKey (- err) lessCoeffs+ lessCoeffsUp =+ Map.insertWith plusUp chplConstTermKey err lessCoeffs+ err = + sum $ map fst smallCoeffTerms+ lessCoeffs =+ Map.fromList $ map snd $ largeCoeffTerms+ (smallCoeffTerms, largeCoeffTerms) = + splitAt (Map.size coeffs - maxTermCount) $+ List.sort $ + map (\(t,c)->(abs c, (t,c))) $ Map.toList coeffs++chplReduceTermCountDown m = fst . chplReduceTermCount m+chplReduceTermCountUp m = snd . chplReduceTermCount m+++{-|+ 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)+ 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++
+ src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Ring.hs view
@@ -0,0 +1,218 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE UndecidableInstances #-}+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Ring+ Description : (internal) uniformly roudned pointwise ring operations + 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 addition and multiplication over polynomials + with pointwise rounding uniform over the whole unit domain.+-}+module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Ring++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++{-|+ Negate a polynomial exactly.+-}+chplNeg (ERChebPoly coeffs) =+ ERChebPoly $ Map.map negate coeffs++{-|+ Add a constant to a polynomial, rounding downwards and upwards. +-}+chplAddConst ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) => + b -> + ERChebPoly box b -> + (ERChebPoly box b, ERChebPoly box b, b)+ {-^ lower and upper bounds on the sum and an upper bound on their difference -}+chplAddConst c (ERChebPoly coeffs) =+ (ERChebPoly sumCoeffsDown, ERChebPoly sumCoeffsUp, err)+ where+ sumCoeffsUp =+ Map.insert chplConstTermKey newConstUp coeffs+ sumCoeffsDown =+ Map.insert chplConstTermKey newConstDown coeffs+ oldConst =+ case Map.lookup chplConstTermKey coeffs of+ Just c -> c+ Nothing -> 0+ newConstUp = oldConst `plusUp` c+ newConstDown = oldConst `plusDown` c+ err = newConstUp - newConstDown ++chplAddConstUp c p = (\(sumDown, sumUp, width) -> sumUp) $ chplAddConst c p+chplAddConstDown c p = (\(sumDown, sumUp, width) -> sumDown) $ chplAddConst c p++{-|+ Add two polynomials, rounding downwards and upwards. +-}+chplAdd ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) => + ERChebPoly box b -> + ERChebPoly box b -> + (ERChebPoly box b, ERChebPoly box b, b)+ {-^ lower and upper bounds on the sum and an upper bound on their difference -}+chplAdd (ERChebPoly coeffs1) (ERChebPoly coeffs2) =+ (ERChebPoly sumCoeffsDown, ERChebPoly sumCoeffsUp, 2 * maxError)+ where+ sumCoeffsUp =+ Map.insertWith plusUp chplConstTermKey maxError coeffsDown+ -- point-wise sum of polynomials with coeff rounding errors+ -- compensated for by enlarging the constant term+ sumCoeffsDown =+ Map.insertWith plusDown chplConstTermKey (- maxError) coeffsUp+ -- point-wise sum of polynomials with coeff rounding errors+ -- compensated for by enlarging the constant term+ coeffsUp =+ (Map.unionWith plusUp 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 plusUp 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 ++p1 +^ p2 = (\(sumDown, sumUp, width) -> sumUp) $ chplAdd p1 p2+p1 +. p2 = (\(sumDown, sumUp, width) -> sumDown) $ chplAdd p1 p2+p1 -^ p2 = p1 +^ (chplNeg p2)+p1 -. p2 = p1 +. (chplNeg p2)++{-|+ Multiply two polynomials, rounding downwards and upwards. +-}+chplMultiply ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) => + ERChebPoly box b -> + ERChebPoly box b -> + (ERChebPoly box b, ERChebPoly box b, b) + {-^ lower and upper bounds on the product and an upper bound on their difference -}+chplMultiply p1@(ERChebPoly coeffs1) p2@(ERChebPoly coeffs2) =+ case (chplGetConst p1, chplGetConst p2) of+ (Just c1, _) -> chplScale c1 p2+ (_, Just c2) -> chplScale c2 p1+ _ -> + (ERChebPoly prodCoeffsDown, ERChebPoly prodCoeffsUp, 2 * roundOffCompensation)+ where+ prodCoeffsUp =+ Map.insertWith plusUp chplConstTermKey roundOffCompensation $ + Map.map negate directProdCoeffsDownNeg+ prodCoeffsDown =+ Map.insertWith plusDown chplConstTermKey (- roundOffCompensation) $ + directProdCoeffsUp+ roundOffCompensation =+ Map.fold plusUp 0 $+ Map.unionWith plusUp directProdCoeffsUp directProdCoeffsDownNeg+ (directProdCoeffsUp, directProdCoeffsDownNeg) =+ foldl addCombiCoeff (Map.empty, Map.empty) combinedCoeffs+ where+ addCombiCoeff+ (prevCoeffsUp, prevCoeffsDownNeg) + (coeffUp, coeffDownNeg, (powersList, coeffCount)) =+ foldl addOnce (prevCoeffsUp, prevCoeffsDownNeg) powersList+ where+ addOnce (prevCoeffsUp, prevCoeffsDownNeg) powers =+ (Map.insertWith plusUp powers coeffUpFrac prevCoeffsUp, + Map.insertWith plusUp powers coeffDownNegFrac prevCoeffsDownNeg)+ coeffUpFrac = coeffUp / coeffCountB+ coeffDownNegFrac = coeffDownNeg / 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++p1 *^ p2 = (\(prodDown,prodUp,width) -> prodUp) $ chplMultiply p1 p2+p1 *. p2 = (\(prodDown,prodUp,width) -> prodDown) $ chplMultiply p1 p2++{-| 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, b)+ {-^ lower and upper bounds on the product and an upper bound on their difference -}+chplScale ratio p@(ERChebPoly coeffs) =+ case chplGetConst p of+ Just c -> + (chplConst cScaledDown, chplConst cScaledUp, cScaledUp - cScaledDown)+ where+ cScaledUp = ratio `timesUp` c+ cScaledDown = ratio `timesDown` c+ _ -> + (ERChebPoly coeffsDown, ERChebPoly coeffsUp, 2 * errBound)+ 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 p = (\(prodDown,prodUp,width) -> prodDown) $ chplScale r p+chplScaleUp r p = (\(prodDown,prodUp,width) -> prodUp) $ chplScale r p++{-|+ Multiply a polynomial by itself, rounding downwards and upwards.+-}+chplSquare ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) => + ERChebPoly box b ->+ (ERChebPoly box b, ERChebPoly box b)+chplSquare p =+ (p2Down, p2Up)+ where+ (p2Down, p2Up, wd) = chplMultiply p p
+ src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Bounds.hs view
@@ -0,0 +1,46 @@+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Bounds+ Description : (testing) properties of bounding operations+ Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable+ + Quickcheck properties of bounding operations, ie constant bounds and binary min/max.+-}+module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Bounds+where++import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Generate++import Data.Number.ER.BasicTypes++import Test.QuickCheck++prop_chplBounds_consistent (ixI, PSize30 (_,p)) =+ ixI >= 2 ==>+ ixI < 100 ==>+ chplAtKeyPointsCanBeLeq p pHigh+ &&+ chplAtKeyPointsCanBeLeq pLow p+ where+ pLow = chplConst cLow+ pHigh = chplConst cHigh+ (cLow, cHigh) = chplBounds ix p+ ix = int2effIx ixI++prop_chplMax_consistent + (Deg20Size20 maxDegree maxSize, PSize30 (_,p1), PSize30 (_, p2)) =+ chplAtKeyPointsPointwiseBinaryDownUpConsistent max p1 p2 (maxLow, maxHigh)+ where+ (maxLow, maxHigh) = chplMax maxDegree maxSize p1 p2++prop_chplMin_consistent (Deg20Size20 maxDegree maxSize, PSize30 (_,p1), PSize30 (_, p2)) =+ chplAtKeyPointsPointwiseBinaryDownUpConsistent min p1 p2 (minLow, minHigh)+ where+ (minLow, minHigh) = chplMin maxDegree maxSize p1 p2+
+ src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Compose.hs view
@@ -0,0 +1,114 @@+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Compose+ Description : (testing) properties of enclosure composition+ Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable+ + Quickcheck properties of polynomial enclosure composition.+-}+module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Compose+where++import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Compose+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Enclosure+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Ring+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Generate++import Data.Number.ER.Real.Approx.Interval+import qualified Data.Number.ER.Real.DomainBox as DBox+import Data.Number.ER.BasicTypes++import Data.Number.ER.Misc++import Test.QuickCheck++prop_enclCompose_ThinEncl_consistent+ reportFileName+ (Deg20Size20 maxDegree maxSize,+ varSelector,+ (PSize30 (n1,p1)),+ (PSize30 (n2,p2))) =+ compose_encl_consistent+ reportFileName + maxDegree maxSize+ varSelector+ n1 p1 n2 p2Encl+ where+ p2Encl = enclThin p2 ++prop_enclCompose_ThickEncl_consistent+ reportFileName+ (Deg20Size20 maxDegree maxSize,+ varSelector,+ (PSize30 (n1,p1)),+ (PSize30 (n21,p21), PSize30 (n22, p22))) =+ compose_encl_consistent+ reportFileName + maxDegree maxSize+ varSelector+ n1 p1 (n21, n22) p2Encl+ where+ p2Encl = makeThickEncl maxDegree maxSize p21 p22 ++prop_enclCompose_ParalEncl_consistent+ reportFileName+ (Deg20Size20 maxDegree maxSize,+ varSelector,+ (PSize30 (n1, p1)),+ (SmallRatio w2Num w2Denom, PSize30 (n2, p2))) =+ compose_encl_consistent + reportFileName+ maxDegree maxSize + varSelector+ n1 p1 ((w2Num, w2Denom), n2) p2Encl+ where+ p2Encl = makeParalEncl p2 w2Num w2Denom++compose_encl_consistent + reportFileName + maxDegree maxSize + varSelector+ p1Id p1 p2Id p2Encl@(p2LowNeg, p2High) =+-- unsafePrint+-- (+-- "compose_encl_consistent: "+-- ++ "\n p1 = " ++ show p1+-- ++ "\n substVar = " ++ show substVar+-- ++ "\n p2Low = " ++ show (chplNeg p2LowNeg)+-- ++ "\n p2High = " ++ show p2High+-- ++ "\n composition = " ++ show resEncl+-- ++ "\n**********************"+-- ) $+ enclAtKeyPointsConsistent+ reportFileName+ ((maxDegree, maxSize), varSelector, p1Id, p2Id)+ composeAtPointInner+ allVars+ resEncl+ where+ resEncl = enclCompose maxDegree maxSize p1 substVar p2Encl+ substVar = p1Vars !! (varSelector `mod` (length p1Vars))+ p1Vars = chplGetVars p1+ allVars = chplGetVars $ p1 +^ p2LowNeg +^ p2High+ p1Encl = (chplNeg p1, p1)+ composeAtPointInner point =+-- unsafePrintReturn+-- (+-- "\n point = " ++ show point+-- ++ "\n substVar = " ++ show substVar+-- ++ " substVal = " ++ show substVal+-- ++ "\n result = "+-- ) $+ enclRAEvalInner p1Encl pointWithSubst+ where+ pointWithSubst =+ DBox.insert substVar substVal $ DBox.map (\b -> ERInterval b b) point+ substVal =+ enclEvalInner p2Encl point+ +
+ src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Division.hs view
@@ -0,0 +1,78 @@+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Division+ Description : (testing) properties of basic enclosure operations+ Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable+ + Quickcheck properties of polynomial enclosure division.+-}+module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Division+where++import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Division+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Enclosure+--import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Generate++import Data.Number.ER.Real.Approx.Interval++import Data.Number.ER.BasicTypes++import Test.QuickCheck++prop_enclRecip_ThickEncl_consistent+ reportFileName+ (Deg20Size20 maxDegree maxSize,+ (Int20 ixInt, Int20 tauDegr),+ SmallRatio sepNum sepDenom,+ (isNegative, PSize30 (n1,p1), PSize30 (n2, p2))) =+ recip_encl_consistent+ reportFileName + maxDegree maxSize + ixInt tauDegr + sepNum sepDenom isNegative (n1, n2) preEncl+ where+ preEncl = makeThickEncl maxDegree maxSize p1 p2 ++prop_enclRecip_ParalEncl_consistent+ reportFileName+ (Deg20Size20 maxDegree maxSize,+ (Int20 ixInt, Int20 tauDegr),+ SmallRatio sepNum sepDenom,+ (isNegative, SmallRatio wNum wDenom, PSize30 (n, p))) =+ recip_encl_consistent + reportFileName+ maxDegree maxSize + ixInt tauDegr + sepNum sepDenom isNegative ((wNum, wDenom), n) preEncl+ where+ preEncl = makeParalEncl p wNum wDenom++recip_encl_consistent + reportFileName+ maxDegree maxSize + ixInt tauDegr + sepNum sepDenom isNegative pId preEncl =+ excludedZero ==>+ enclAtKeyPointsPointwiseUnaryDownUpConsistent+ reportFileName+ ((maxDegree, maxSize), (ixInt, tauDegr), (sepNum, sepDenom), (isNegative, pId)) + (intervalDivideInner 1) + pEncl resEncl+ where+ resEncl = enclRecip maxDegree maxSize ix tauDegr pEncl+ ix = int2effIx ixInt+ (excludedZero, pEncl) =+ enclRestrictRange ix rangeNoZero preEncl+ rangeNoZero+ | isNegative = (Nothing, Just (-sepB))+ | otherwise = (Just sepB, Nothing)+ sepB = abs sepNumB / sepDenomB+ sepNumB = fromInteger $ toInteger sepNum+ sepDenomB = fromInteger $ toInteger sepDenom+ +
+ src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Elementary.hs view
@@ -0,0 +1,120 @@+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Elementary+ Description : (testing) properties of basic enclosure operations+ Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable+ + Quickcheck properties of some elementary operations on primitive polynomial+ enclosures.+-}+module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Elementary+where++import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Elementary+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Enclosure+--import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Generate++import qualified Data.Number.ER.Real.Approx as RA+import Data.Number.ER.Real.Approx.Interval+import Data.Number.ER.Real.Arithmetic.Elementary++import Data.Number.ER.BasicTypes++import Test.QuickCheck++prop_enclExp_ThickEncl_consistent =+ encl_op_ThickEncl_consistent enclExp erExp_IR_Inner noDomainRestriction++prop_enclExp_ParalEncl_consistent =+ encl_op_ParalEncl_consistent enclExp erExp_IR_Inner noDomainRestriction+ +prop_enclExp_ThinEncl_consistent =+ encl_op_ThinEncl_consistent enclExp erExp_IR_Inner noDomainRestriction+ +prop_enclSine_ThickEncl_consistent =+ encl_op_ThickEncl_consistent enclSine erSine_IR_Inner sincosDomain++prop_enclSine_ParalEncl_consistent =+ encl_op_ParalEncl_consistent enclSine erSine_IR_Inner sincosDomain+ +prop_enclSine_ThinEncl_consistent =+ encl_op_ThinEncl_consistent enclSine erSine_IR_Inner sincosDomain+ +prop_enclCosine_ThickEncl_consistent =+ encl_op_ThickEncl_consistent enclCosine erCosine_IR_Inner sincosDomain++prop_enclCosine_ParalEncl_consistent =+ encl_op_ParalEncl_consistent enclCosine erCosine_IR_Inner sincosDomain+ +prop_enclCosine_ThinEncl_consistent =+ encl_op_ThinEncl_consistent enclCosine erCosine_IR_Inner sincosDomain+ +prop_enclAtan_ThickEncl_consistent =+ encl_op_ThickEncl_consistent enclAtan erATan_IR_Inner noDomainRestriction++prop_enclAtan_ParalEncl_consistent =+ encl_op_ParalEncl_consistent enclAtan erATan_IR_Inner noDomainRestriction+ +prop_enclAtan_ThinEncl_consistent =+ encl_op_ThinEncl_consistent enclAtan erATan_IR_Inner noDomainRestriction++sincosDomain = (Just (-1.57), Just 1.57) -- almost (-pi/2, pi/2)+noDomainRestriction = (Nothing, Nothing)+ +encl_op_ThickEncl_consistent+ opEncl opInner rangeRestriction+ reportFileName+ (Deg20Size20 maxDegree maxSize,+ (Int20 ixInt),+ (PSize30 (n1,p1), PSize30 (n2, p2))) =+ enclAtKeyPointsPointwiseUnaryDownUpConsistent+ reportFileName+ ((maxDegree, maxSize), ixInt, (n1, n2)) + (opInner ix) + pEncl resEncl+ where+ (succeeded, pEncl) = + enclRestrictRange ix rangeRestriction $ makeThickEncl maxDegree maxSize p1 p2 + resEncl = opEncl maxDegree maxSize ix pEncl+ ix = int2effIx ixInt+ +encl_op_ParalEncl_consistent+ opEncl opInner rangeRestriction+ reportFileName+ (Deg20Size20 maxDegree maxSize,+ (Int20 ixInt),+ (SmallRatio wNum wDenom, PSize30 (n, p))) =+ enclAtKeyPointsPointwiseUnaryDownUpConsistent + reportFileName+ ((maxDegree, maxSize), ixInt, ((wNum, wDenom), n)) + (opInner ix) + pEncl resEncl+ where+ (succeeded, pEncl) = + enclRestrictRange ix rangeRestriction $ makeParalEncl p wNum wDenom + resEncl = opEncl maxDegree maxSize ix pEncl+ ix = int2effIx ixInt+ +encl_op_ThinEncl_consistent+ opEncl opInner rangeRestriction+ reportFileName+ (Deg20Size20 maxDegree maxSize,+ (Int20 ixInt),+ (PSize30 (n, p))) =+ enclAtKeyPointsPointwiseUnaryDownUpConsistent + reportFileName+ ((maxDegree, maxSize), ixInt, n) + (opInner ix)+ pEncl resEncl+ where+ (succeeded, pEncl) = + enclRestrictRange ix rangeRestriction $ enclThin p + resEncl = opEncl maxDegree maxSize ix pEncl+ ix = int2effIx ixInt+ +
+ src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Enclosure.hs view
@@ -0,0 +1,106 @@+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Enclosure+ Description : (testing) properties of basic enclosure operations+ Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable+ + Quickcheck properties of basic enclosure operations, + mainly ring operations.+-}+module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Enclosure+where++import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Enclosure+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Generate++import Data.Number.ER.Real.Approx.Interval++prop_enclAdd_ThickEncls_consistent+ reportFileName+ (Deg20Size20 maxDegree maxSize,+ (PSize30 (n11,p11), PSize30 (n12, p12)),+ (PSize30 (n21,p21), PSize30 (n22, p22))) =+ enclAtKeyPointsPointwiseBinaryDownUpConsistent+ reportFileName+ ((maxDegree, maxSize), (n11, n12), (n21, n22))+ intervalPlusInner+ p1Encl p2Encl sumEncl+ where+ sumEncl = p1Encl +: p2Encl+ p1Encl = makeThickEncl maxDegree maxSize p11 p12 + p2Encl = makeThickEncl maxDegree maxSize p21 p22 + +prop_enclMultiply_ThickEncls_consistent+ reportFileName+ (Deg20Size20 maxDegree maxSize,+ (PSize30 (n11,p11), PSize30 (n12, p12)),+ (PSize30 (n21,p21), PSize30 (n22, p22))) =+ enclAtKeyPointsPointwiseBinaryDownUpConsistent+ reportFileName+ ((maxDegree, maxSize), (n11, n12), (n21, n22))+ intervalTimesInner+ p1Encl p2Encl prodEncl+ where+ prodEncl = enclMultiply maxDegree maxSize p1Encl p2Encl+ p1Encl = makeThickEncl maxDegree maxSize p11 p12 + p2Encl = makeThickEncl maxDegree maxSize p21 p22 + +prop_enclMultiply_ParalEncls_consistent+ reportFileName+ (Deg20Size20 maxDegree maxSize,+ (SmallRatio num1 denom1,+ PSize30 (n1,p1)),+ (SmallRatio num2 denom2,+ PSize30 (n2,p2))) =+ enclAtKeyPointsPointwiseBinaryDownUpConsistent + reportFileName+ ((maxDegree, maxSize), ((num1, denom1), n1), ((num2, denom2), n2))+ intervalTimesInner+ p1Encl p2Encl prodEncl+ where+ prodEncl = enclMultiply maxDegree maxSize p1Encl p2Encl+ p1Encl = makeParalEncl p1 num1 denom1+ p2Encl = makeParalEncl p2 num2 denom2+ +prop_enclScale_ThickEncl_consistent+ reportFileName+ (Deg20Size20 maxDegree maxSize,+ SmallRatio num denom,+ PSize30 (n1, p1), + PSize30 (n2, p2)) =+ enclAtKeyPointsPointwiseBinaryDownUpConsistent+ reportFileName + ((maxDegree, maxSize), (num, denom), (n1, n2))+ intervalTimesInner+ cEncl pEncl scaledEncl+ where+ scaledEncl = enclScale maxDegree maxSize cB pEncl+ pEncl = makeThickEncl maxDegree maxSize p1 p2 + cEncl = enclConst cB + cB = numB / denomB+ numB = fromInteger $ toInteger num+ denomB = fromInteger $ toInteger denom+ +prop_enclScale_ParalEncl_consistent+ reportFileName+ (Deg20Size20 maxDegree maxSize,+ SmallRatio cNum cDenom,+ (SmallRatio wNum wDenom, PSize30 (n, p))) =+ enclAtKeyPointsPointwiseBinaryDownUpConsistent+ reportFileName + ((maxDegree, maxSize), (cNum, cDenom), ((wNum, wDenom), n))+ intervalTimesInner + cEncl pEncl scaledEncl+ where+ scaledEncl = enclScale maxDegree maxSize cB pEncl+ pEncl = makeParalEncl p wNum wDenom + cEncl = enclConst cB + cB = cNumB / cDenomB+ cNumB = fromInteger $ toInteger cNum+ cDenomB = fromInteger $ toInteger cDenom+ +
+ src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Generate.hs view
@@ -0,0 +1,592 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Generate+ Description : (testing) generating polynomials for tests+ Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable+ + A collection of polynomials to pick from when testing.+-}+module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Generate+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.Reduce+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Eval+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Enclosure++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 Data.Number.ER.BasicTypes++import Data.Number.ER.Real.DefaultRepr+import Data.Number.ER.Real.DomainBox.IntMap+import Data.Number.ER.Real.Approx.Interval+import qualified Data.Number.ER.Real.Approx as RA+++import Test.QuickCheck hiding (two, three)++import qualified Data.Map as Map++{---------------------}+{----- Type synonyms for different polynomial generation distributions ----}+{---------------------}++type P = ERChebPoly (Box Int) BM++newtype PNoLimits = PNoLimits (Int, P) deriving (Show)+newtype PSize10Degree3 = PSize10Degree3 (Int, P) deriving (Show)+newtype PSize10Degree10 = PSize10Degree10 (Int, P) deriving (Show)+newtype PSize10 = PSize10 (Int, P) deriving (Show)+newtype PSize30 = PSize30 ((Int, Int), P) deriving (Show)++instance (Arbitrary PNoLimits)+ where+ arbitrary =+ elements $ map PNoLimits $ zip [0..] $ + polynomials1200ish id+ coarbitrary p =+ error "ERChebPoly: Generate: Arbitrary: coarbitrary not implemented for polynomials"++instance (Arbitrary PSize10Degree3) + where+ arbitrary =+ elements $ map PSize10Degree3 $ zip [0..] $ polynomials1200ishSize10Degree3 + coarbitrary (PSize10Degree3 p) =+ error "ERChebPoly: Generate: Arbitrary: coarbitrary not implemented for polynomials"++polynomials1200ishSize10Degree3 =+ polynomials1200ish $ chplReduceTermCountUp 10 . chplReduceDegreeUp 3++instance (Arbitrary PSize10Degree10) + where+ arbitrary =+ elements $ map PSize10Degree10 $ zip [0..] $ + polynomials1200ishSize10Degree10+ coarbitrary (PSize10Degree10 p) =+ error "ERChebPoly: Generate: Arbitrary: coarbitrary not implemented for polynomials"++polynomials1200ishSize10Degree10 =+ polynomials1200ish $ chplReduceTermCountUp 10 . chplReduceDegreeUp 10++instance (Arbitrary PSize10) + where+ arbitrary =+ elements $ map PSize10 $ zip [0..] $ polynomials1200ishSize10 + + coarbitrary (PSize10 p) =+ error "ERChebPoly: Generate: Arbitrary: coarbitrary not implemented for polynomials"++polynomials1200ishSize10 =+ polynomials1200ish $ chplReduceTermCountUp 10+ +instance (Arbitrary PSize30) + where+ arbitrary =+ sized arbitrarySized+ where+ arbitrarySized n + | n <= 28 =+ elements $ map PSize30 $ + zip (map (\n -> (0,n)) [0..]) $ + polynomials200ishSize30+ | otherwise =+ elements $ map PSize30 $ + zip (map (\n -> (1,n)) [0..]) $ + polynomials1200ishSize30+ coarbitrary (PSize30 p) =+ error "ERChebPoly: Generate: Arbitrary: coarbitrary not implemented for polynomials"++polynomials1200ishSize30 =+ polynomials1200ish $ chplReduceTermCountUp 30+ +polynomials200ishSize30 =+ polynomials200ishSmall $ chplReduceTermCountUp 30+ +data Int20 = Int20 Int deriving (Show)+ +instance (Arbitrary Int20)+ where+ arbitrary =+ do+ n <- choose (2,20)+ return $ Int20 n+ coarbitrary (Int20 n) =+ error "ERChebPoly: Generate: Arbitrary: coarbitrary not implemented for EffIx20"++data Deg20Size20 = Deg20Size20 Int Int deriving (Show)+ +instance (Arbitrary Deg20Size20)+ where+ arbitrary =+ do+ maxDegree <- choose (2,20)+ maxSize <- choose (10,20)+ return $ Deg20Size20 maxDegree maxSize+ coarbitrary (Deg20Size20 maxDegree maxSize) =+ error "ERChebPoly: Generate: Arbitrary: coarbitrary not implemented for Deg20Size20"++data SmallRatio = SmallRatio Int Int deriving (Show)+ +instance (Arbitrary SmallRatio)+ where+ arbitrary =+ do+ num <- choose (-1000000,1000000)+ denom <- choose (1,1000000)+ return $ SmallRatio num denom+ coarbitrary (SmallRatio num denom) =+ error "ERChebPoly: Generate: Arbitrary: coarbitrary not implemented for SmallRatio"+ + +{------------------}+{-------- Functions commonly used in tests: ----------}+{------------------}++chplAtKeyPointsCanBeLeq ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box, + DomainBoxMappable boxb boxbb varid b [ERInterval b], Show boxb) => + ERChebPoly box b ->+ ERChebPoly box b ->+ Bool+chplAtKeyPointsCanBeLeq p1 p2 =+ and $ map testPoint points+ where+ points = getKeyPoints (p1 +^ p2)+ testPoint point + | lower1 <= upper2 =+ True+ | otherwise =+ unsafePrint+ (+ "Failure at point = " ++ (show point)+ ) $+ False+ where+ lower1 = chplEvalDown p1 point+ upper2 = chplEvalUp p2 point + +getKeyPoints p =+ getKeyPointsForVars $ chplGetVars p+ +getKeyPointsForVars vars =+ points+ where+ points = map DBox.fromList $ allCombinations $ map getVarPoints varDoms+ varDoms = map (\v -> (v,unitInterval)) vars+ unitInterval = ERInterval (-1) 1+ getVarPoints (var, dom) =+ (var, [domLB, domMB, domRB])+ where+ ERInterval domLB domRB = dom+ domMB = (domLB + domRB)/2++chplAtKeyPointsPointwiseBinaryDownUpConsistent ::+ (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box, + DomainBoxMappable boxb boxbb varid b [ERInterval b], Show boxb) =>+ ((ERInterval b) -> (ERInterval b) -> (ERInterval b)) -> + ERChebPoly box b ->+ ERChebPoly box b ->+ (ERChebPoly box b, ERChebPoly box b) ->+ Bool+chplAtKeyPointsPointwiseBinaryDownUpConsistent raOp p1 p2 (resLow, resHigh) =+ and $ map testPoint points+ where+ points = getKeyPoints (p1 +^ p2)+ testPoint point + | ok = ok+ | otherwise =+ unsafePrint+ (+ "chplAtKeyPointsPointwiseBinaryDownUpConsistent failed:"+ ++ "\n point = " ++ show point+ ++ "\n raOpAtPointHigh = " ++ show raOpAtPointHigh+ ++ "\n raOpAtPointLow = " ++ show raOpAtPointLow+ ++ "\n resAtPointHigh = " ++ show resAtPointHigh+ ++ "\n resAtPointLow = " ++ show resAtPointLow+ )+ ok+ where+ ok = + raOpAtPointLow <= resAtPointHigh+ &&+ raOpAtPointHigh >= resAtPointLow+ resAtPointLow = chplEvalDown resLow point+ resAtPointHigh = chplEvalUp resHigh point+ raOpAtPoint@(ERInterval raOpAtPointLow raOpAtPointHigh) = + raOp p1AtPoint p2AtPoint + p1AtPoint = ERInterval p1AtPointLow p1AtPointHigh+ (p1AtPointLow, p1AtPointHigh) = chplEval p1 point+ p2AtPoint = ERInterval p2AtPointLow p2AtPointHigh+ (p2AtPointLow, p2AtPointHigh) = chplEval p2 point++makeThickEncl maxDegree maxSize p1 p2 =+ (pMax q1Neg q2Neg, pMax q1 q2)+ where+ q1Neg = chplNeg q1+ q2Neg = chplNeg q2+ q1 = p1 +^ p2Mp1ScaledDown+ q2 = p1 -^ p2Mp1ScaledDown+ p2Mp1ScaledDown =+ chplScaleUp (10/sizeB) p2Mp1+ where+ sizeB = max (abs upperB) (abs lowerB)+ (lowerB, upperB) = chplBounds 10 p2Mp1+ p2Mp1 = p2 -^ p1+ pMax = chplMaxUp maxDegree maxSize+ +makeParalEncl p num denom =+-- unsafePrintReturn+-- (+-- "makeThinEncl: result = "+-- )+ (pNeg, p +^ cP)+ where+ pNeg = chplNeg p+ cP = chplConst cB+ cB = abs $ numB / (1000 * denomB)+ numB = fromInteger $ toInteger num+ denomB = fromInteger $ toInteger denom+ +enclRestrictRange ix (Nothing, Nothing) pEncl = (True, pEncl)+enclRestrictRange ix (maybeLower, maybeUpper) preEncl =+ (succeeded, pEncl)+ where+ succeeded = lowerSucceeded && upperSucceeded+ lowerSucceeded =+ case maybeLower of+ Nothing -> True+ Just lower -> pLowerBound > lower + upperSucceeded =+ case maybeUpper of+ Nothing -> True+ Just upper -> pUpperBound < upper+ (pLowerBound, pUpperBound) = enclBounds ix pEncl+ pEncl =+ case (maybeLower, maybeUpper) of+ (Just lowerB, Nothing) ->+ case lowerB <= preLowerBoundB of+ True -> preEncl -- enclosure already in the range+ False -> -- a shift needed to get above the lower bound+ enclAddConst (lowerB - preLowerBoundB + sepB) preEncl+ (Nothing, Just upperB) ->+ case preUpperBoundB <= upperB of+ True -> preEncl -- enclosure already in the range+ False -> -- a shift needed to get below the upper bound+ enclAddConst (upperB - preUpperBoundB - sepB) preEncl+ (Just lowerB, Just upperB) ->+ case lowerB <= preLowerBoundB && preUpperBoundB <= upperB of+ True -> preEncl -- enclosure already in the range+ _ -> + case preWidthB + sepB <= widthB of+ True -> -- no scaling needed, only shifting by a constant to the centre of the range+ enclAddConst + (lowerB - preLowerBoundB + (preWidthB - widthB)/2) + preEncl+ _ -> -- full affine transformation needed+ enclAddConst+ (lowerB + sepB) $+ enclScaleNonneg -- scale preEncl so that it fits inside the range+ (widthB / saferPreWidthB) $+ enclAddConst -- shift preEncl so that it is non-negative and as close to 0 as safely possible+ (sepB - preLowerBoundB)+ preEncl+ where + widthB = upperB - lowerB+ saferPreWidthB = preWidthB + 2 * sepB+ sepB = preWidthB / 1000000+ preWidthB = preUpperBoundB - preLowerBoundB+ (preLowerBoundB, preUpperBoundB) = enclBounds ix preEncl+ + ++enclAtKeyPointsPointwiseBinaryDownUpConsistent ::+ (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box, + DomainBoxMappable boxb boxbb varid b [ERInterval b], Show boxb, Show testId) =>+ String {-^ report file name -} ->+ testId {-^ item to identify the random input given to the test -} ->+ ((ERInterval b) -> (ERInterval b) -> (ERInterval b)) ->+ {-^ this real approx operation has to return an inner approximation of the exact result set, + ie each number that the approximation supports is in the maximal extension -}+ (ERChebPoly box b, ERChebPoly box b) {-^ enclosure of argument 1 -} ->+ (ERChebPoly box b, ERChebPoly box b) {-^ enclosure of argument 2 -} ->+ (ERChebPoly box b, ERChebPoly box b) {-^ alleged enclosure of result -} ->+ Bool+enclAtKeyPointsPointwiseBinaryDownUpConsistent+ reportFileName testId+ raOpInner + p1Encl@(p1LowNeg, p1High) p2Encl@(p2LowNeg, p2High) resEncl =+ and $ map testPoint points+ where+ points = getKeyPoints (p1High +^ p2High +^ p1LowNeg +^ p2LowNeg)+ testPoint point + | result =+ unsafeReport reportFileName+ (+ show $ + (testId, point, p1OpInnerP2AtPoint, resAtPoint)+ ) + result+ | otherwise = + unsafePrint+ (+ "enclAtKeyPointsPointwiseBinaryDownUpConsistent failed"+ ++ "\n point = " ++ show point+ ++ "\n p1AtPoint = " ++ show p1AtPoint+ ++ "\n p2AtPoint = " ++ show p2AtPoint+ ++ "\n p1OpInnerP2AtPoint = " ++ show p1OpInnerP2AtPoint+ ++ "\n resAtPoint = " ++ show resAtPoint+ ) $+ result+ where+ result = p1OpInnerP2AtPoint `RA.refines` resAtPoint+ p1OpInnerP2AtPoint = p1AtPoint `raOpInner` p2AtPoint+ resAtPoint = enclEval resEncl point+-- resAtPoint = p1OpInnerP2AtPoint -- for dummy testing that never <<loop>>s+ p1AtPoint = normaliseERInterval $ enclEvalInner p1Encl point+ p2AtPoint = normaliseERInterval $ enclEvalInner p2Encl point++enclAtKeyPointsPointwiseUnaryDownUpConsistent ::+ (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box, + DomainBoxMappable boxb boxbb varid b [ERInterval b], Show boxb, Show testId) =>+ String {-^ report file name -} ->+ testId {-^ item to identify the random input given to the test -} ->+ ((ERInterval b) -> (ERInterval b)) ->+ {-^ this real approx operation has to return an inner approximation of the exact result set, + ie each number that the approximation supports is in the maximal extension -}+ (ERChebPoly box b, ERChebPoly box b) {-^ enclosure of argument -} ->+ (ERChebPoly box b, ERChebPoly box b) {-^ alleged enclosure of result -} ->+ Bool+enclAtKeyPointsPointwiseUnaryDownUpConsistent+ reportFileName testId+ raOpInner + pEncl@(pLowNeg, pHigh) resEncl =+ and $ map testPoint points+ where+ points = getKeyPoints (pHigh +^ pLowNeg)+ testPoint point + | result =+ unsafeReport reportFileName+ (+ show $ + (testId, point, opInnerPAtPoint, resAtPoint)+ )+ result + | otherwise = + unsafePrint+ (+ "enclAtKeyPointsPointwiseUnaryDownUpConsistent failed"+ ++ "\n point = " ++ show point+ ++ "\n pAtPoint = " ++ show pAtPoint+ ++ "\n opInnerPAtPoint = " ++ show opInnerPAtPoint+ ++ "\n resAtPoint = " ++ show resAtPoint+ ) $+ result+ where+ result = opInnerPAtPoint `RA.refines` resAtPoint+ opInnerPAtPoint = raOpInner pAtPoint+ resAtPoint = enclEval resEncl point+ pAtPoint = +-- normaliseERInterval $ + enclEvalInner pEncl point+++enclAtKeyPointsConsistent ::+ (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box, + DomainBoxMappable boxb boxbb varid b [ERInterval b], Show boxb, Show testId) =>+ String {-^ report file name -} ->+ testId {-^ item to identify the random input given to the test -} ->+ (boxb -> (ERInterval b)) ->+ {-^ this operation has to return an inner approximation of the exact result set, + ie each number that the approximation supports is a solution in the maximal extension -}+ [varid] {-^ variables to test over -} ->+ (ERChebPoly box b, ERChebPoly box b) {-^ alleged enclosure of result -} ->+ Bool+enclAtKeyPointsConsistent+ reportFileName testId+ opInner allVars resEncl@(resLowNeg, resHigh) =+ and $ map testPoint points+ where+ points = getKeyPointsForVars allVars+ testPoint point + | result =+ unsafeReport reportFileName+ (+ show $ + (testId, point, opInnerAtPoint, resAtPoint)+ )+ result + | otherwise = + unsafePrint+ (+ "enclAtKeyPointsConsistent failed"+ ++ "\n point = " ++ show point+ ++ "\n opInnerAtPoint = " ++ show opInnerAtPoint+ ++ "\n resAtPoint = " ++ show resAtPoint+ ) $+ result+ where+ result = opInnerAtPoint `RA.refines` resAtPoint+ opInnerAtPoint = opInner point+ resAtPoint = enclEval resEncl point+++{------------------}+{-------- A diverse collection of polynomials to pick from: ----------}+{------------------}++type E = (P,P)++vars :: [P]+vars = map chplVar [0..7]++varsE :: [E]+varsE = map (\p -> (chplNeg p, p)) vars++x0 = vars !! 0+x1 = vars !! 1+x2 = vars !! 2+x3 = vars !! 3+x4 = vars !! 4++x0E = varsE !! 0+x1E = varsE !! 1+x2E = varsE !! 2+x3E = varsE !! 3+x4E = varsE !! 4++one :: P+[mone, one, two, three, seven, thousand, million, tiny, huge] = + map chplConst + [-1,1,2,3,7,1000,1000000,10^^(-200),10^^200]++oneE :: E+[moneE, oneE, twoE, threeE, sevenE, thousandE, millionE, tinyE, hugeE] = + map (\ c -> (chplConst (-c), chplConst c))+ [-1,1,2,3,7,1000,1000000,10^^(-200),10^^200]++polynomials1200ish rdc =+ concat $ map (powers10 rdc) $+ concat $ map addConsts3 $+ concat $ map multConsts3 $+ polyBase13+ +polynomials200ish rdc =+ concat $ map (powers4 rdc) $+ concat $ map addConsts3 $+ concat $ map multConsts3 $+ polyBase5+ +polynomials40ish rdc =+ concat $ map (powers2 rdc) $+ concat $ map addConsts2 $+ concat $ map multConsts2 $+ polyBase5+ +polynomials200ishSmall rdc =+ concat $ map (powers4Small rdc) $+ concat $ map addConsts3 $+ concat $ map multConsts3 $+ polyBase5+ +polynomials40ishSmall rdc =+ concat $ map (powers2Small rdc) $+ concat $ map addConsts2 $+ concat $ map multConsts2 $+ polyBase5+ ++polyBase5 =+ [+ (two *^ x0) +^ x1+ ,+ (seven *^ x0) -^ x1+ ,+ (tiny *^ x0) +^ x1+ ,+ x0 -^ x1 *^ x2+ ,+ x0 -^ x1 +^ x2 -^ x3 +^ x4+ ]+ +polyBase13 =+ [+ x0+ ,+ x0 +^ x1+ ,+ x0 -^ x1+ ,+ (two *^ x0) +^ x1+ ,+ (two *^ x0) -^ x1+ ,+ (seven *^ x0) +^ x1+ ,+ (seven *^ x0) -^ x1+ ,+ (tiny *^ x0) +^ x1+ ,+ (tiny *^ x0) -^ x1+ ,+ x0 -^ x1 +^ x2+ ,+ x0 -^ x1 *^ x2+ ,+ x0 +^ x1 +^ x2 +^ x3 +^ x4+ ,+ x0 -^ x1 +^ x2 -^ x3 +^ x4+ ]+ +powersAll rdc p =+ powersAux [p, rdc $ p *^ p]+ where+ powersAux (pNHalfM1 : pNHalf : rest) = + pNHalfM1 : (powersAux $ (pNHalf : rest) ++ [pNM1, pN])+ where+ pNM1 = rdc $ pNHalf *^ pNHalfM1+ pN = rdc $ pNHalf *^ pNHalf++powersForExps rdc p exponents =+ map pw exponents+ where+ pw n = pws !! (n - 1)+ pws = powersAll rdc p++powers10 rdc p =+ powersForExps rdc p [1..10]++powers4 rdc p =+ powersForExps rdc p [1,3,5,7]+ +powers4Small rdc p =+ powersForExps rdc p [1,2,3,5]+ +powers2 rdc p =+ powersForExps rdc p [1,7]+ +powers2Small rdc p =+ powersForExps rdc p [1,3]+ +addConsts3 p =+ [p +^ one, p +^ three, p +^ seven]++multConsts3 p =+ [p *^ two, p *^ three, p *^ seven]+ +addConsts2 p =+ [p +^ one, p +^ three]++multConsts2 p =+ [p *^ two, p *^ seven]+
+ src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Reduce.hs view
@@ -0,0 +1,37 @@+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Reduce+ Description : (testing) properties of reduction operations+ Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable+ + Quickcheck properties of operations that reduce the size of polynomials.+-}+module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Reduce+where++import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Reduce+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Generate++import Test.QuickCheck++prop_chplReduceTermCount_consistent (PSize30 (_,p), Deg20Size20 _ maxSize) =+ maxSize < chplCountTerms p ==>+ chplAtKeyPointsCanBeLeq p pUp+ && + chplAtKeyPointsCanBeLeq pDown p+ where+ (pDown, pUp) = chplReduceTermCount maxSize p + ++prop_chplReduceDegree_consistent (PSize30 (_,p), Deg20Size20 maxDegree _) =+ maxDegree < chplGetDegree p ==>+ chplAtKeyPointsCanBeLeq p pUp+ && + chplAtKeyPointsCanBeLeq pDown p+ where+ (pDown, pUp) = chplReduceDegree maxDegree p
+ src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Ring.hs view
@@ -0,0 +1,47 @@+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Ring+ Description : (testing) properties of ring operations+ Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable+ + Quickcheck properties of ring operations, ie addition and multiplication.+-}+module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Ring+where++import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Ring+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Generate++prop_chplAdd_consistent (PSize30 (_,p1), PSize30 (_, p2)) =+ chplAtKeyPointsPointwiseBinaryDownUpConsistent (+) p1 p2 (sumLow, sumHigh)+ where+ (sumLow, sumHigh, _) = chplAdd p1 p2++prop_chplAddConst_consistent (SmallRatio num denom, PSize30 (_, p)) =+ chplAtKeyPointsPointwiseBinaryDownUpConsistent (+) cP p (sumLow, sumHigh)+ where+ (sumLow, sumHigh, _) = chplAddConst cB p+ cP = chplConst cB+ cB = numB / denomB+ numB = fromInteger $ toInteger num+ denomB = fromInteger $ toInteger denom++prop_chplMult_consistent (PSize30 (_,p1), PSize30 (_, p2)) =+ chplAtKeyPointsPointwiseBinaryDownUpConsistent (*) p1 p2 (prodLow, prodHigh)+ where+ (prodLow, prodHigh, _) = chplMultiply p1 p2++prop_chplScale_consistent (SmallRatio num denom, PSize30 (_, p)) =+ chplAtKeyPointsPointwiseBinaryDownUpConsistent (*) cP p (prodLow, prodHigh)+ where+ (prodLow, prodHigh, _) = chplScale cB p+ cP = chplConst cB+ cB = numB / denomB+ numB = fromInteger $ toInteger num+ denomB = fromInteger $ toInteger denom+
+ src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Run.hs view
@@ -0,0 +1,159 @@+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Run+ Description : (testing) running all polynomial tests in a batch+ Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable+ + Support for running all polynomial tests in a batch.+-}+module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Run+where++import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Generate+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Reduce+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Ring+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Bounds+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Enclosure+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Division+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Elementary+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Compose+--import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Integration++import qualified Data.Number.ER.RnToRm.UnitDom.Base as UFB+import qualified Data.Number.ER.Real.Base as B+import Data.Number.ER.Real.Approx.Interval+import Data.Number.ER.Real.Arithmetic.Elementary+import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox)++import Data.Number.ER.Real.DefaultRepr+import Data.Number.ER.Misc++import Test.QuickCheck+import Test.QuickCheck.Batch++import System.IO+import System.Directory+import qualified System.FilePath as FP+import Data.Time.Clock+import Data.Time.Calendar++initArith = B.initialiseBaseArithmetic (0::BM)++runPolynomTests =+ do+ (UTCTime (ModifiedJulianDay days) secs) <- getCurrentTime+ let folder = "tests-" ++ (show days) ++ "-" ++ (show $ floor $ toRational secs)+ createDirectory folder+-- mkRunTests "poly tests" chplTestOptions (chplTests folder)+ mkRunTests "poly tests" chplTestOptions (enclTests folder)+ +instance Show TestResult+ where+ show result =+ case result of+ TestOk msg ntest stamps ->+ msg ++ " " ++ show ntest ++ " " -- ++ show stamps+ TestExausted msg ntest stamps ->+ msg ++ " " ++ show ntest ++ " " -- ++ show stamps+ TestAborted exception ->+ "aborted: " ++ show exception+ TestFailed args ntest ->+ "failed after " ++ show ntest ++ " tests" + ++ "\n args = " ++ show args+ +mkRunTests testsetName options tests =+ do+ initArith+ mapM (mkRunTest $ length tests) $ zip [1..] tests+ return ()+ where+ mkRunTest testCount (n, (testName, test)) =+ do+ putStr testDescr+ result <- test options+ putStrLn $ " result: " ++ show result+-- runTests testDescr options [test]+ hFlush stdout+ where+ testDescr = + "(" ++ show n ++ "/" ++ show testCount ++ ") " ++ testsetName ++ ": " ++ testName ++ "\n" ++chplTestOptions = + TestOptions+ { +-- no_of_tests = 10+-- no_of_tests = 50+ no_of_tests = 100+-- no_of_tests = 200+ , + length_of_tests = 240 * 3600 -- ie 4h time limit+ ,+ debug_tests = False + }++chplTests folder =+ [+ ("reduce term count", run prop_chplReduceTermCount_consistent),+ ("reduce degree", run prop_chplReduceDegree_consistent),+ ("add two polys", run prop_chplAdd_consistent),+ ("add const to poly", run prop_chplAddConst_consistent),+ ("mult two polys", run prop_chplMult_consistent),+ ("scale poly", run prop_chplScale_consistent),+ ("bounds of poly", run prop_chplBounds_consistent),+ ("max of two polys", run prop_chplMax_consistent),+ ("min of two polys", run prop_chplMin_consistent)+ ]+enclTests folder =+ [+ ("add thick encls", run $ prop_enclAdd_ThickEncls_consistent $ addFolder "enclAdd_Thick"),+ ("mult paral encls", run $ prop_enclMultiply_ParalEncls_consistent $ addFolder "enclMultiply_Paral"),+ ("mult thick encls", run $ prop_enclMultiply_ThickEncls_consistent $ addFolder "enclMultiply_Thick"),+ ("scale paral encl", run $ prop_enclScale_ParalEncl_consistent $ addFolder "enclScale_Paral"),+ ("scale thick encl", run $ prop_enclScale_ThickEncl_consistent $ addFolder "enclScale_Thick"),+ ("recip paral encl", run $ prop_enclRecip_ParalEncl_consistent $ addFolder "enclRecip_Paral"),+ ("recip thick encl", run $ prop_enclRecip_ThickEncl_consistent $ addFolder "enclRecip_Thick"),+ ("compose thin encl", run $ prop_enclCompose_ThinEncl_consistent $ addFolder "enclCompose_Thin"),+ ("compose paral encl", run $ prop_enclCompose_ParalEncl_consistent $ addFolder "enclCompose_Paral"),+ ("compose thick encl", run $ prop_enclCompose_ThickEncl_consistent $ addFolder "enclCompose_Thick"),+ ("exp thin encl", run $ prop_enclExp_ThinEncl_consistent $ addFolder "enclExp_Thin"),+ ("exp paral encl", run $ prop_enclExp_ParalEncl_consistent $ addFolder "enclExp_Paral"),+ ("exp thick encl", run $ prop_enclExp_ThickEncl_consistent $ addFolder "enclExp_Thick"),+ ("sine thin encl", run $ prop_enclSine_ThinEncl_consistent $ addFolder "enclSine_Thin"),+ ("sine paral encl", run $ prop_enclSine_ParalEncl_consistent $ addFolder "enclSine_Paral"),+ ("sine thick encl", run $ prop_enclSine_ThickEncl_consistent $ addFolder "enclSine_Thick"),+ ("cosine thin encl", run $ prop_enclCosine_ThinEncl_consistent $ addFolder "enclCosine_Thin"),+ ("cosine paral encl", run $ prop_enclCosine_ParalEncl_consistent $ addFolder "enclCosine_Paral"),+ ("cosine thick encl", run $ prop_enclCosine_ThickEncl_consistent $ addFolder "enclCosine_Thick"),+ ("atan thin encl", run $ prop_enclAtan_ThinEncl_consistent $ addFolder "enclAtan_Thin"),+ ("atan paral encl", run $ prop_enclAtan_ParalEncl_consistent $ addFolder "enclAtan_Paral"),+ ("atan thick encl", run $ prop_enclAtan_ThickEncl_consistent $ addFolder "enclAtan_Thick")+ ]+ where+ addFolder name = FP.combine folder name+ ++-- failed tests:++--failed1 = +-- -- identified 19 Feb 9:33+-- -- fixed 19 Feb 16:50+-- prop_enclCompose_ThickEncl_consistent "a"+-- (Deg20Size20 4 18, 0,+-- PSize30 ((0,112), polynomials200ishSize30 !! 112),+-- (PSize30 ((0,57), polynomials200ishSize30 !! 57),+-- PSize30 ((0,18), polynomials200ishSize30 !! 18)+-- ) +-- )++failed2 = + -- identified 19 Feb 18:59 -- this one makes the automatic test abort with <<loop>>+ -- but runs ok when executed individually+ prop_enclMultiply_ParalEncls_consistent "a"+ (Deg20Size20 5 11,+ (SmallRatio 680377 535300, PSize30 ((1,1018), polynomials1200ishSize30 !! 1018)),+ (SmallRatio (-157647) 491208, PSize30 ((1,465), polynomials1200ishSize30 !! 465))+ )