summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorMichalKonecny <>2009-02-23 16:40:58 (GMT)
committerLuite Stegeman <luite@luite.com>2009-02-23 16:40:58 (GMT)
commit2bbd45fcf507f0956e0bfc3845af20aa1f230705 (patch)
tree0bfa88efd2e04186f62d96250929621c5e755ec2
parent207aa45f6ada4a47828546ac9580b67555382326 (diff)
version 0.4.90.4.9
-rw-r--r--AERN-RnToRm.cabal50
-rw-r--r--ChangeLog14
-rw-r--r--src/Data/Number/ER/RnToRm/Approx.hs129
-rw-r--r--src/Data/Number/ER/RnToRm/Approx/DomTransl.hs89
-rw-r--r--src/Data/Number/ER/RnToRm/TestingDefs.hs31
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/Approx.hs85
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs563
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/Base.hs539
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs94
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Basic.hs126
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Bounds.hs348
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Compose.hs138
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Division.hs168
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs667
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Enclosure.hs311
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Eval.hs99
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Field.hs228
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Integration.hs120
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Reduce.hs85
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Ring.hs218
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Bounds.hs46
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Compose.hs114
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Division.hs78
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Elementary.hs120
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Enclosure.hs106
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Generate.hs592
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Reduce.hs37
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Ring.hs47
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Tests/Run.hs159
29 files changed, 3895 insertions, 1506 deletions
diff --git a/AERN-RnToRm.cabal b/AERN-RnToRm.cabal
index 726e7d3..ce95d1f 100644
--- 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 @@ Maintainer: mik@konecny.aow.cz
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 @@ Description:
with Taylor Models is included in the
paper <http://www-users.aston.ac.uk/~konecnym/papers/cfv08.html>.
.
- Simple examples of usage can be found in module @Demo.hs@ in folder @tests@.
+ Simple examples of usage can be found in folder @tests@.
Extra-source-files:
tests/Demo.hs tests/ISin3.hs
Data-files:
ChangeLog
-Flag containers-in-base
- Default: False
-
Library
hs-source-dirs: src
- if flag(containers-in-base)
- Build-Depends:
- base < 3, binary >= 0.4, html >= 1.0, AERN-Real >= 0.9.7
- else
- Build-Depends:
- base >= 3, containers, binary >= 0.4, html >= 1.0, AERN-Real >= 0.9.7
+ Build-Depends:
+ AERN-Real >= 0.9.9, base >= 3, base < 4, containers, binary >= 0.4, html >= 1.0, QuickCheck >= 1.2, QuickCheck < 2, time, filepath, directory
Exposed-modules:
Data.Number.ER.RnToRm,
- Data.Number.ER.RnToRm.BisectionTree.Path,
- Data.Number.ER.RnToRm.BisectionTree.Integration,
- Data.Number.ER.RnToRm.BisectionTree,
- Data.Number.ER.RnToRm.DefaultRepr,
- Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic,
- Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Elementary,
- Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Field,
Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Integration,
+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Compose,
Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Eval,
Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds,
+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic,
+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Ring,
+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Enclosure,
+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Elementary,
+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Reduce,
+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Division,
+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Compose,
+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Bounds,
+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Ring,
+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Enclosure,
+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Elementary,
+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Generate,
+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Reduce,
+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Division,
+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Run,
Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom,
Data.Number.ER.RnToRm.UnitDom.Base,
Data.Number.ER.RnToRm.UnitDom.Approx.Interval,
Data.Number.ER.RnToRm.UnitDom.Approx,
+ Data.Number.ER.RnToRm.TestingDefs,
+ Data.Number.ER.RnToRm.DefaultRepr,
+ Data.Number.ER.RnToRm.BisectionTree.Integration,
+ Data.Number.ER.RnToRm.BisectionTree.Path,
+ Data.Number.ER.RnToRm.BisectionTree,
+ Data.Number.ER.RnToRm.Approx.DomEdges,
Data.Number.ER.RnToRm.Approx.DomTransl,
Data.Number.ER.RnToRm.Approx.PieceWise,
- Data.Number.ER.RnToRm.Approx.DomEdges,
Data.Number.ER.RnToRm.Approx.Tuple,
- Data.Number.ER.RnToRm.Approx,
- Data.Number.ER.RnToRm.TestingDefs
+ Data.Number.ER.RnToRm.Approx
diff --git a/ChangeLog b/ChangeLog
index 4798a59..584c5aa 100644
--- 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
index b5dc932..86be60c 100644
--- a/src/Data/Number/ER/RnToRm/Approx.hs
+++ b/src/Data/Number/ER/RnToRm/Approx.hs
@@ -18,7 +18,9 @@ module Data.Number.ER.RnToRm.Approx
(
ERFnApprox(..),
ERFnDomApprox(..),
- bisectUnbisectDepth
+ bisectUnbisectDepth,
+ keyPointsConsistencyCheck,
+ keyPointsPointwiseConsistencyCheck
)
where
@@ -29,6 +31,8 @@ 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
{-|
@@ -48,7 +52,7 @@ import qualified Data.Map as Map
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 @@ class
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 @@ class
-}
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 @@ class
-}
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 @@ class
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 @@ bisectUnbisectDepth depth f =
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
+ \ No newline at end of file
diff --git a/src/Data/Number/ER/RnToRm/Approx/DomTransl.hs b/src/Data/Number/ER/RnToRm/Approx/DomTransl.hs
index c01795e..d2c8282 100644
--- a/src/Data/Number/ER/RnToRm/Approx/DomTransl.hs
+++ b/src/Data/Number/ER/RnToRm/Approx/DomTransl.hs
@@ -25,7 +25,6 @@ module Data.Number.ER.RnToRm.Approx.DomTransl
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.Real.DomainBox (VariableID(..), DomainIntBox, DomainBoxMap
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 @@ instance
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 @@ instance
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 @@ instance
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 @@ instance
(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 @@ instance
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 @@ instance
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
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 @@ instance
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 @@ instance
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
index b9b8469..378e17c 100644
--- a/src/Data/Number/ER/RnToRm/TestingDefs.hs
+++ b/src/Data/Number/ER/RnToRm/TestingDefs.hs
@@ -31,8 +31,13 @@ fapdConst01 = (FA.const DBox.noinfo [0 RA.\/ 1]) :: (FAPD B)
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)
@@ -70,16 +75,30 @@ testIntegrE =
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
index 1e0b989..a529fd9 100644
--- 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 @@ import qualified Data.Map as Map
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 @@ class (ERFnApprox box varid domra ranra fa) =>
[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 @@ class (ERFnApprox box varid domra ranra fa) =>
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
+ \ No newline at end of file
diff --git a/src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs b/src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs
index 48f9f69..5f409fd 100644
--- 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 Data.Number.ER.Real.Arithmetic.Elementary
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 @@ data ERFnInterval fb ra =
|
ERFnInterval
{
- erfnUpper :: fb,
erfnLowerNeg :: fb,
+ erfnUpper :: fb,
erfnContext :: ERFnContext,
erfnGlobal :: ra
}
@@ -89,24 +90,26 @@ data ERFnContext =
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 @@ instance
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 @@ instance
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)
- 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
+ ERFnInterval ln h ctxt (gl1 * gl2)
+ where
+ (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 @@ erfnintLeq left right
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 @@ instance
-- 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 @@ instance
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 @@ instance
-- ++ "\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 @@ sincos isSine maxDegree ix u l
(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 @@ sincos isSine maxDegree ix u l
(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 @@ sincos isSine maxDegree ix u l
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
- where
- (lCosDown, lCosUp) = UFB.cos maxDegree ix (l `plusUp` shiftPoly)
- (uCosDown, uCosUp) = UFB.cos maxDegree ix (u `plusDown` shiftPoly)
- shiftPoly = UFB.const shiftLOB
+ 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
+ (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 @@ instance
normalise $
ERFnInterval
{
- erfnUpper = fbH,
erfnLowerNeg = fbLNeg,
+ erfnUpper = fbH,
erfnContext = context,
erfnGlobal = val
}
@@ -694,8 +777,8 @@ instance
normalise $
ERFnInterval
{
- erfnUpper = fbH,
erfnLowerNeg = fbLNeg,
+ erfnUpper = fbH,
erfnContext = context,
erfnGlobal =
UFB.raFromEndpoints fbH
@@ -723,6 +806,42 @@ instance
{
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 @@ instance
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 @@ instance
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 @@ instance
-- ++ "\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
index 46d4752..23011c7 100644
--- 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 Data.Number.ER.BasicTypes
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
index a100bab..d236527 100644
--- 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 @@ module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom
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 qualified Data.Number.ER.Real.Base as B
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 @@ x1 = chplVar 1 :: P
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
index e720865..6135ee4 100644
--- 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 @@ data ERChebPoly box b =
{
chplCoeffs :: (Map.Map (TermKey box) b)
}
- deriving (Eq, Typeable, Data)
+ deriving (Eq, Ord, Typeable, Data)
type TermKey box = box
@@ -138,7 +138,7 @@ chplConst val =
(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 ::
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 @@ chplNormalise (ERChebPoly coeffs) =
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 @@ instance (B.ERRealBase b, DomainBox box varid Int, Ord box) => Show (ERChebPoly
-}
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 @@ chplShow showChebyshevBasis (ERChebPoly coeffs)
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 @@ chebToXBasis ::
(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
index f10847a..3e9b949 100644
--- 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 @@ module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds
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 @@ import Data.List
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 @@ chplUpperBoundQuadr ::
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 @@ chplUpperBoundQuadr ix p@(ERChebPoly coeffs) =
(and $ map maybeInUnit $ DBox.elems peak)
,
erintv_right $
- chplEvalApprox makeInterval peak p
+ chplRAEval makeInterval p peak
)
Nothing -> (False, undefined)
where
@@ -167,7 +180,7 @@ chplUpperBoundQuadr ix p@(ERChebPoly coeffs) =
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 @@ chplUpperBoundQuadr ix p@(ERChebPoly coeffs) =
_ ->
[(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 @@ chplNonneg = chplNonnegCubic
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) `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
- 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 *)
+ (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
index 0000000..9dd7353
--- /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
+ \ No newline at end of file
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
index 0000000..e9e5700
--- /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
index cbebc67..d1157bd 100644
--- 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 @@ module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Elementary
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 Data.Number.ER.Misc
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
- )
+ enclAddConst 1 $
+ (chplConst 0, pSquareHigh)
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)
- 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
+ coeff = coeffBase / (coeffBase + 1)
- -- define such translated T_i's for all i >= 0:
- trTis =
- map (mapPair (chplReduceDegreeDown maxDegree, chplReduceDegreeUp maxDegree)) $
- chebyEvalTsRoundDownUp trT1Dn
-
- -- 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
index 0000000..d6b9d3a
--- /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
index 7907aac..6e728e1 100644
--- 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 @@ module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Eval
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 @@ import qualified Data.Map as Map
-}
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 @@ chebyEvalTsExact val =
{-|
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 @@ chplEvalApprox b2ra vals (ERChebPoly coeffs) =
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 @@ chplPartialEvalApprox ra2endpts substitutions (ERChebPoly coeffs) =
(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
index 0603f9b..0000000
--- 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
index 461103b..400cf5a 100644
--- 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 @@ where
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 @@ chplIntegrate ::
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 @@ chplIntegrate x (ERChebPoly coeffs) =
| 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 @@ chplIntegrate x (ERChebPoly coeffs) =
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 @@ chplIntegrate x (ERChebPoly coeffs) =
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
index 0000000..61b0e12
--- /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
index 0000000..7bf93c5
--- /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
index 0000000..31a448a
--- /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
index 0000000..983ebed
--- /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
+
+ \ No newline at end of file
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
index 0000000..9d2198f
--- /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
+
+ \ No newline at end of file
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
index 0000000..9e8f470
--- /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
+
+ \ No newline at end of file
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
index 0000000..2fe024c
--- /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
+
+ \ No newline at end of file
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
index 0000000..24a503e
--- /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 <<loop>>s
+ p1AtPoint = normaliseERInterval $ enclEvalInner p1Encl point
+ p2AtPoint = normaliseERInterval $ enclEvalInner p2Encl point
+
+enclAtKeyPointsPointwiseUnaryDownUpConsistent ::
+ (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box,
+ DomainBoxMappable boxb boxbb varid b [ERInterval b], Show boxb, Show testId) =>
+ String {-^ report file name -} ->
+ testId {-^ item to identify the random input given to the test -} ->
+ ((ERInterval b) -> (ERInterval b)) ->
+ {-^ this real approx operation has to return an inner approximation of the exact result set,
+ ie each number that the approximation supports is in the maximal extension -}
+ (ERChebPoly box b, ERChebPoly box b) {-^ enclosure of argument -} ->
+ (ERChebPoly box b, ERChebPoly box b) {-^ alleged enclosure of result -} ->
+ Bool
+enclAtKeyPointsPointwiseUnaryDownUpConsistent
+ reportFileName testId
+ raOpInner
+ pEncl@(pLowNeg, pHigh) resEncl =
+ and $ map testPoint points
+ where
+ points = getKeyPoints (pHigh +^ pLowNeg)
+ testPoint point
+ | result =
+ unsafeReport reportFileName
+ (
+ show $
+ (testId, point, opInnerPAtPoint, resAtPoint)
+ )
+ result
+ | otherwise =
+ unsafePrint
+ (
+ "enclAtKeyPointsPointwiseUnaryDownUpConsistent failed"
+ ++ "\n point = " ++ show point
+ ++ "\n pAtPoint = " ++ show pAtPoint
+ ++ "\n opInnerPAtPoint = " ++ show opInnerPAtPoint
+ ++ "\n resAtPoint = " ++ show resAtPoint
+ ) $
+ result
+ where
+ result = opInnerPAtPoint `RA.refines` resAtPoint
+ opInnerPAtPoint = raOpInner pAtPoint
+ resAtPoint = enclEval resEncl point
+ pAtPoint =
+-- normaliseERInterval $
+ enclEvalInner pEncl point
+
+
+enclAtKeyPointsConsistent ::
+ (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box,
+ DomainBoxMappable boxb boxbb varid b [ERInterval b], Show boxb, Show testId) =>
+ String {-^ report file name -} ->
+ testId {-^ item to identify the random input given to the test -} ->
+ (boxb -> (ERInterval b)) ->
+ {-^ this operation has to return an inner approximation of the exact result set,
+ ie each number that the approximation supports is a solution in the maximal extension -}
+ [varid] {-^ variables to test over -} ->
+ (ERChebPoly box b, ERChebPoly box b) {-^ alleged enclosure of result -} ->
+ Bool
+enclAtKeyPointsConsistent
+ reportFileName testId
+ opInner allVars resEncl@(resLowNeg, resHigh) =
+ and $ map testPoint points
+ where
+ points = getKeyPointsForVars allVars
+ testPoint point
+ | result =
+ unsafeReport reportFileName
+ (
+ show $
+ (testId, point, opInnerAtPoint, resAtPoint)
+ )
+ result
+ | otherwise =
+ unsafePrint
+ (
+ "enclAtKeyPointsConsistent failed"
+ ++ "\n point = " ++ show point
+ ++ "\n opInnerAtPoint = " ++ show opInnerAtPoint
+ ++ "\n resAtPoint = " ++ show resAtPoint
+ ) $
+ result
+ where
+ result = opInnerAtPoint `RA.refines` resAtPoint
+ opInnerAtPoint = opInner point
+ resAtPoint = enclEval resEncl point
+
+
+{------------------}
+{-------- A diverse collection of polynomials to pick from: ----------}
+{------------------}
+
+type E = (P,P)
+
+vars :: [P]
+vars = map chplVar [0..7]
+
+varsE :: [E]
+varsE = map (\p -> (chplNeg p, p)) vars
+
+x0 = vars !! 0
+x1 = vars !! 1
+x2 = vars !! 2
+x3 = vars !! 3
+x4 = vars !! 4
+
+x0E = varsE !! 0
+x1E = varsE !! 1
+x2E = varsE !! 2
+x3E = varsE !! 3
+x4E = varsE !! 4
+
+one :: P
+[mone, one, two, three, seven, thousand, million, tiny, huge] =
+ map chplConst
+ [-1,1,2,3,7,1000,1000000,10^^(-200),10^^200]
+
+oneE :: E
+[moneE, oneE, twoE, threeE, sevenE, thousandE, millionE, tinyE, hugeE] =
+ map (\ c -> (chplConst (-c), chplConst c))
+ [-1,1,2,3,7,1000,1000000,10^^(-200),10^^200]
+
+polynomials1200ish rdc =
+ concat $ map (powers10 rdc) $
+ concat $ map addConsts3 $
+ concat $ map multConsts3 $
+ polyBase13
+
+polynomials200ish rdc =
+ concat $ map (powers4 rdc) $
+ concat $ map addConsts3 $
+ concat $ map multConsts3 $
+ polyBase5
+
+polynomials40ish rdc =
+ concat $ map (powers2 rdc) $
+ concat $ map addConsts2 $
+ concat $ map multConsts2 $
+ polyBase5
+
+polynomials200ishSmall rdc =
+ concat $ map (powers4Small rdc) $
+ concat $ map addConsts3 $
+ concat $ map multConsts3 $
+ polyBase5
+
+polynomials40ishSmall rdc =
+ concat $ map (powers2Small rdc) $
+ concat $ map addConsts2 $
+ concat $ map multConsts2 $
+ polyBase5
+
+
+polyBase5 =
+ [
+ (two *^ x0) +^ x1
+ ,
+ (seven *^ x0) -^ x1
+ ,
+ (tiny *^ x0) +^ x1
+ ,
+ x0 -^ x1 *^ x2
+ ,
+ x0 -^ x1 +^ x2 -^ x3 +^ x4
+ ]
+
+polyBase13 =
+ [
+ x0
+ ,
+ x0 +^ x1
+ ,
+ x0 -^ x1
+ ,
+ (two *^ x0) +^ x1
+ ,
+ (two *^ x0) -^ x1
+ ,
+ (seven *^ x0) +^ x1
+ ,
+ (seven *^ x0) -^ x1
+ ,
+ (tiny *^ x0) +^ x1
+ ,
+ (tiny *^ x0) -^ x1
+ ,
+ x0 -^ x1 +^ x2
+ ,
+ x0 -^ x1 *^ x2
+ ,
+ x0 +^ x1 +^ x2 +^ x3 +^ x4
+ ,
+ x0 -^ x1 +^ x2 -^ x3 +^ x4
+ ]
+
+powersAll rdc p =
+ powersAux [p, rdc $ p *^ p]
+ where
+ powersAux (pNHalfM1 : pNHalf : rest) =
+ pNHalfM1 : (powersAux $ (pNHalf : rest) ++ [pNM1, pN])
+ where
+ pNM1 = rdc $ pNHalf *^ pNHalfM1
+ pN = rdc $ pNHalf *^ pNHalf
+
+powersForExps rdc p exponents =
+ map pw exponents
+ where
+ pw n = pws !! (n - 1)
+ pws = powersAll rdc p
+
+powers10 rdc p =
+ powersForExps rdc p [1..10]
+
+powers4 rdc p =
+ powersForExps rdc p [1,3,5,7]
+
+powers4Small rdc p =
+ powersForExps rdc p [1,2,3,5]
+
+powers2 rdc p =
+ powersForExps rdc p [1,7]
+
+powers2Small rdc p =
+ powersForExps rdc p [1,3]
+
+addConsts3 p =
+ [p +^ one, p +^ three, p +^ seven]
+
+multConsts3 p =
+ [p *^ two, p *^ three, p *^ seven]
+
+addConsts2 p =
+ [p +^ one, p +^ three]
+
+multConsts2 p =
+ [p *^ two, p *^ seven]
+
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
index 0000000..4d954b6
--- /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
index 0000000..fae2a10
--- /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
index 0000000..edca3f6
--- /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 <<loop>>
+ -- but runs ok when executed individually
+ prop_enclMultiply_ParalEncls_consistent "a"
+ (Deg20Size20 5 11,
+ (SmallRatio 680377 535300, PSize30 ((1,1018), polynomials1200ishSize30 !! 1018)),
+ (SmallRatio (-157647) 491208, PSize30 ((1,465), polynomials1200ishSize30 !! 465))
+ )