diff --git a/AERN-RnToRm.cabal b/AERN-RnToRm.cabal --- a/AERN-RnToRm.cabal +++ b/AERN-RnToRm.cabal @@ -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 . . - 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 diff --git a/ChangeLog b/ChangeLog --- a/ChangeLog +++ b/ChangeLog @@ -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 diff --git a/src/Data/Number/ER/RnToRm/Approx.hs b/src/Data/Number/ER/RnToRm/Approx.hs --- a/src/Data/Number/ER/RnToRm/Approx.hs +++ b/src/Data/Number/ER/RnToRm/Approx.hs @@ -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 + diff --git a/src/Data/Number/ER/RnToRm/Approx/DomTransl.hs b/src/Data/Number/ER/RnToRm/Approx/DomTransl.hs --- a/src/Data/Number/ER/RnToRm/Approx/DomTransl.hs +++ b/src/Data/Number/ER/RnToRm/Approx/DomTransl.hs @@ -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 = diff --git a/src/Data/Number/ER/RnToRm/TestingDefs.hs b/src/Data/Number/ER/RnToRm/TestingDefs.hs --- a/src/Data/Number/ER/RnToRm/TestingDefs.hs +++ b/src/Data/Number/ER/RnToRm/TestingDefs.hs @@ -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 diff --git a/src/Data/Number/ER/RnToRm/UnitDom/Approx.hs b/src/Data/Number/ER/RnToRm/UnitDom/Approx.hs --- a/src/Data/Number/ER/RnToRm/UnitDom/Approx.hs +++ b/src/Data/Number/ER/RnToRm/UnitDom/Approx.hs @@ -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 + diff --git a/src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs b/src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs --- a/src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs +++ b/src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs @@ -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)) diff --git a/src/Data/Number/ER/RnToRm/UnitDom/Base.hs b/src/Data/Number/ER/RnToRm/UnitDom/Base.hs --- a/src/Data/Number/ER/RnToRm/UnitDom/Base.hs +++ b/src/Data/Number/ER/RnToRm/UnitDom/Base.hs @@ -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 + diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs --- a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs @@ -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 + diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Basic.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Basic.hs --- a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Basic.hs +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Basic.hs @@ -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 $ diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Bounds.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Bounds.hs --- a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Bounds.hs +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Bounds.hs @@ -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 diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Compose.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Compose.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Compose.hs @@ -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 + diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Division.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Division.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Division.hs @@ -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 diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs --- a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs @@ -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 diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Enclosure.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Enclosure.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Enclosure.hs @@ -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 + + diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Eval.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Eval.hs --- a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Eval.hs +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Eval.hs @@ -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 diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Field.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Field.hs deleted file mode 100644 --- a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Field.hs +++ /dev/null @@ -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" - diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Integration.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Integration.hs --- a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Integration.hs +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Integration.hs @@ -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" diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Reduce.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Reduce.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Reduce.hs @@ -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 + + diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Ring.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Ring.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Ring.hs @@ -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 diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Bounds.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Bounds.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Bounds.hs @@ -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 + diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Compose.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Compose.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Compose.hs @@ -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 + + diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Division.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Division.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Division.hs @@ -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 + + diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Elementary.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Elementary.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Elementary.hs @@ -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 + + diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Enclosure.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Enclosure.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Enclosure.hs @@ -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 + + diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Generate.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Generate.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Generate.hs @@ -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 <>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] + diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Reduce.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Reduce.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Reduce.hs @@ -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 diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Ring.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Ring.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Ring.hs @@ -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 + diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Run.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Run.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Run.hs @@ -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 <> + -- 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)) + )