diff --git a/AERN-RnToRm.cabal b/AERN-RnToRm.cabal
--- a/AERN-RnToRm.cabal
+++ b/AERN-RnToRm.cabal
@@ -1,5 +1,5 @@
 Name:           AERN-RnToRm
-Version:        0.4.1
+Version:        0.4.2
 Cabal-Version:  >= 1.2
 Build-Type:     Simple
 License:        BSD3
@@ -10,7 +10,7 @@
 Stability:      experimental
 Category:       Data, Math
 Synopsis:       polynomial function enclosures (PFEs) approximating exact real functions
-Tested-with:    GHC ==6.8.2
+Tested-with:    GHC ==6.8.3
 Description:
     AERN-RnToRm provides
     datatypes and abstractions for approximating functions
@@ -34,7 +34,9 @@
     .
     Simple examples of usage can be found in module @Demo.hs@ in folder @tests@.
 Extra-source-files:
-    ChangeLog tests/Demo.hs tests/ISin3.hs
+    tests/Demo.hs tests/ISin3.hs
+Data-files:
+    ChangeLog
 
 Flag containers-in-base
     Default: False
@@ -43,10 +45,10 @@
   hs-source-dirs:  src
   if flag(containers-in-base)
     Build-Depends:
-      base < 3, binary >= 0.4, AERN-Real >= 0.9.7
+      base < 3, binary >= 0.4, html >= 1.0, AERN-Real >= 0.9.7
   else
     Build-Depends:
-      base >= 3, containers, binary >= 0.4, AERN-Real >= 0.9.7
+      base >= 3, containers, binary >= 0.4, html >= 1.0, AERN-Real >= 0.9.7
   Exposed-modules:
     Data.Number.ER.RnToRm,
     Data.Number.ER.RnToRm.BisectionTree.Path,
@@ -69,13 +71,4 @@
     Data.Number.ER.RnToRm.Approx.Tuple,
     Data.Number.ER.RnToRm.Approx,
     Data.Number.ER.RnToRm.TestingDefs  
-    
-  Extensions: 
-    CPP,
-    DeriveDataTypeable,
-    FlexibleContexts,
-    FlexibleInstances,
-    FunctionalDependencies,
-    MultiParamTypeClasses,
-    UndecidableInstances
-  
+
diff --git a/ChangeLog b/ChangeLog
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,6 @@
+0.4.2: 1 December 2008
+    * substantially improved division by a constant PFE
+    * added proper handling of overflown coefficients
 0.4.1: 30 September 2008
     * updated to work with AERN-Real 0.9.7    
 0.4: 20 August 2008
diff --git a/src/Data/Number/ER/RnToRm/Approx.hs b/src/Data/Number/ER/RnToRm/Approx.hs
--- a/src/Data/Number/ER/RnToRm/Approx.hs
+++ b/src/Data/Number/ER/RnToRm/Approx.hs
@@ -120,8 +120,8 @@
     partialIntersect ::
         EffortIndex -> 
         box {-^ the subdomain; defined by clipping the range of some variables -} ->
-        fa {-^ function to improve by intersecting its subdomain -} -> 
         fa {-^ the enclosure to be used on the subdomain (but defined on the whole domain) -} ->
+        fa {-^ function to improve by intersecting its subdomain -} -> 
         fa
     {-|
         Intersect two enclosures and measure the global improvement as one number.
diff --git a/src/Data/Number/ER/RnToRm/Approx/DomEdges.hs b/src/Data/Number/ER/RnToRm/Approx/DomEdges.hs
--- a/src/Data/Number/ER/RnToRm/Approx/DomEdges.hs
+++ b/src/Data/Number/ER/RnToRm/Approx/DomEdges.hs
@@ -31,6 +31,9 @@
 import Data.Number.ER.Misc
 import Data.Number.ER.PlusMinus
 
+import Data.Number.ER.ShowHTML
+import qualified Text.Html as H
+
 import Data.Typeable
 import Data.Generics.Basics
 import Data.Binary
@@ -152,6 +155,39 @@
         domB = FA.dom fa
 
 instance
+    (FA.ERFnDomApprox box varid domra ranra fa, H.HTML fa) =>
+    H.HTML (ERFnDomEdgesApprox varid fa)
+    where
+    toHtml f@(ERFnDomEdgesApprox mainEncl edges) 
+        | Map.null edges =
+            H.toHtml mainEncl
+        | otherwise =
+            H.toHtml $
+                abovesTable [H.border 2] [leftEdgesHtml, H.toHtml mainEncl, rightEdgesHtml]
+        where
+        leftEdgesHtml = 
+            H.toHtml $ H.simpleTable [H.border 1] [] leftRows
+        rightEdgesHtml = 
+            H.toHtml $ H.simpleTable [H.border 1] [] rightRows
+        leftRows =
+            map (makeRow fst) leftEdges
+        rightRows =
+            map (makeRow snd) rightEdges
+        (leftEdges, rightEdges) =
+             partition (\((_, pm),_) -> pm == Minus) $ Map.toList edges
+        makeRow pickEndpoint ((varId, _), faEdge) =
+            [
+                H.stringToHtml $ showVar varId ++ "=" ++ show varValueLeft
+            ,
+                H.toHtml faEdge
+            ]
+            where
+            varValueLeft = 
+                pickEndpoint $ 
+                RA.bounds $ DBox.lookup "RnToRm.Approx.DomEdges: " varId domB
+        domB = FA.dom mainEncl
+        
+instance
     (FA.ERFnApprox box varid domra ranra fa) =>
     Eq (ERFnDomEdgesApprox varid fa)
     where
@@ -183,7 +219,7 @@
 
 
 instance 
-    (FA.ERFnDomApprox box varid domra ranra fa, VariableID varid) =>
+    (FA.ERFnDomApprox box varid domra ranra fa, VariableID varid, Show box) =>
     RA.ERApprox (ERFnDomEdgesApprox varid fa)
     where
     initialiseBaseArithmetic _ =
@@ -211,9 +247,44 @@
             Map.intersectionWith (RA.intersectMeasureImprovement ix) edges1 edges2 
     leqReals fa1 fa2 =
         RA.leqReals (erfnMainVolume fa1) (erfnMainVolume fa2)
-
+    refines
+            f1@(ERFnDomEdgesApprox mainEncl1 edges1) 
+            f2@(ERFnDomEdgesApprox mainEncl2 edges2) =
+--        unsafePrint 
+--        (
+--            "ERFnDomEdgesApprox: refines: "
+--            ++ "\n domB = " ++ show (FA.dom mainEncl1) 
+--            ++ "\n mainRefines = " ++ show mainRefines 
+--            ++ "\n mainEncl1 = " ++ show mainEncl1 
+--            ++ "\n mainEncl2 = " ++ show mainEncl2 
+--        )$
+        mainRefines && edgesRefine
+        where
+        mainRefines =
+            mainEncl1 `RA.refines` mainEncl2
+        edgesRefine =
+            Map.fold (&&) True $
+                Map.intersectionWith RA.refines edges1 edges2  
+    compareApprox
+            f1@(ERFnDomEdgesApprox mainEncl1 edges1) 
+            f2@(ERFnDomEdgesApprox mainEncl2 edges2) =
+        compareComposeMany
+        [
+            RA.compareApprox mainEncl1 mainEncl2
+        ,
+            compareListsWith compareEdges (Map.toAscList edges1) (Map.toAscList edges2)
+        ]
+        where
+        compareEdges ((v1,pm1),fa1) ((v2,pm2),fa2) =
+            compareComposeMany
+            [
+                compare v1 v2,
+                compare pm1 pm2,
+                RA.compareApprox fa1 fa2
+            ]
 instance 
-    (FA.ERFnDomApprox box varid domra ranra fa, RA.ERIntApprox fa, VariableID varid) =>
+    (FA.ERFnDomApprox box varid domra ranra fa, RA.ERIntApprox fa, 
+     VariableID varid, Show box) =>
     RA.ERIntApprox (ERFnDomEdgesApprox varid fa)
     where
 --    doubleBounds = :: ira -> (Double, Double) 
@@ -246,7 +317,8 @@
     f1 \/ f2 = edgesLift2 (RA.\/) f1 f2
 
 instance 
-    (FA.ERFnDomApprox box varid domra ranra fa, RAEL.ERApproxElementary fa, VariableID varid) =>
+    (FA.ERFnDomApprox box varid domra ranra fa, RAEL.ERApproxElementary fa, 
+     VariableID varid, Show box) =>
     RAEL.ERApproxElementary (ERFnDomEdgesApprox varid fa)
     where
     abs ix = edgesLift1 $ RAEL.abs ix
@@ -257,7 +329,8 @@
     atan ix = edgesLift1 $ RAEL.atan ix
         
 instance 
-    (FA.ERFnDomApprox box varid domra ranra fa, VariableID varid) =>
+    (FA.ERFnDomApprox box varid domra ranra fa, 
+     VariableID varid, Show box) =>
     FA.ERFnApprox box varid domra ranra (ERFnDomEdgesApprox varid fa)
     where
     check prgLocation (ERFnDomEdgesApprox mainEncl edges) =
@@ -312,7 +385,7 @@
         partialIntersectEdge (var, pm) edge1 edge2 
             | withinSubstitutions =
                 FA.partialIntersect ix substitutions edge1 edge2
-            | otherwise = edge1
+            | otherwise = edge2
             where
             withinSubstitutions =
                 (varDomEndpoint pm) `RA.refines` varVal
@@ -322,7 +395,7 @@
             varDomEndpoint Minus = varDomLO
             varDomEndpoint Plus = varDomHI
             (varDomLO, varDomHI) = RA.bounds varDom
-            varDom = DBox.lookup "DomEdges: partialIntersect: " var $ FA.dom mainEncl1 
+            varDom = DBox.lookup "DomEdges: partialIntersect: " var $ FA.dom mainEncl2 
     eval ptBox (ERFnDomEdgesApprox mainEncl edges) 
         | null edgeVals =
             mainVal
@@ -364,7 +437,8 @@
             substitutionsNew = DBox.delete varID substitutionsPrev
             
 instance 
-    (FA.ERFnDomApprox box varid domra ranra fa, VariableID varid) =>
+    (FA.ERFnDomApprox box varid domra ranra fa, 
+     VariableID varid, Show box) =>
     FA.ERFnDomApprox box varid domra ranra (ERFnDomEdgesApprox varid fa)
     where
     dom (ERFnDomEdgesApprox mainEncl edges) = FA.dom mainEncl
diff --git a/src/Data/Number/ER/RnToRm/Approx/DomTransl.hs b/src/Data/Number/ER/RnToRm/Approx/DomTransl.hs
--- a/src/Data/Number/ER/RnToRm/Approx/DomTransl.hs
+++ b/src/Data/Number/ER/RnToRm/Approx/DomTransl.hs
@@ -35,6 +35,8 @@
 import Data.Number.ER.BasicTypes
 import Data.Number.ER.Misc
 
+import qualified Text.Html as H
+
 import Data.Typeable
 import Data.Generics.Basics
 import Data.Binary
@@ -90,6 +92,13 @@
     
 instance
     (RA.ERIntApprox domra) =>
+    Ord (DomTransl domra)
+    where
+    compare (DomTransl _ _ _ _ dom1) (DomTransl _ _ _ _ dom2) =
+        RA.compareApprox dom1 dom2
+    
+instance
+    (RA.ERIntApprox domra) =>
     Show (DomTransl domra)
     where
     show (DomTransl fromA fromB toA toB dom) =
@@ -98,6 +107,15 @@
         "  fromUnit = " ++ show fromA ++ " * x + " ++ show fromB ++ "\n" ++
         "  toUnit = " ++ show toA ++ " * x + " ++ show toB ++ "\n"
 
+instance
+    (RA.ERIntApprox domra, H.HTML domra) =>
+    H.HTML (DomTransl domra)
+    where
+    toHtml (DomTransl fromA fromB toA toB dom) =
+        error "DomTransl: toHtml not implemented yet"
+
+
+
 dtrIdentity ::
     (RA.ERIntApprox ira) =>
     DomTransl ira
@@ -177,9 +195,9 @@
 
 #endif
 
-dtrBShow dtrs =
+dtrBShow dtrB =
     concatWith "," $
-    map showOne $ DBox.toList dtrs
+    map showOne $ DBox.toList dtrB
     where
     showOne (var, dtr) =
         showVar var ++ " in " ++ show (dtrDom dtr)
@@ -194,24 +212,42 @@
         "\nERFnDomTranslApprox" ++
         show ufaDom ++
 --        show ufa ++
-        "\n dom = [" ++
-        (concatWith ", " $ map showVarDom $ DBox.toList $ dtrBToDomB dtrB)
-        ++ "]"
+        "\n dom = [" ++ dtrBShow dtrB ++ "]"
         where
         ufaDom =
-            FA.composeThin ufa $ 
-                Map.fromAscList $ 
-                    map mkToUnitUFA $ 
-                        DBox.toAscList dtrB
+            translateUfaToDom ufa dtrB
 --        gr = 20 + (RA.getGranularity ufa)
-        mkToUnitUFA (var, tr) =
-            (var, UFA.affine [co] (Map.singleton var [sl]))
-            where
-            sl = FA.domra2ranra ufa $ dtrToUnitSlope tr
-            co = FA.domra2ranra ufa $ dtrToUnitConst tr
-        showVarDom (varID, varDom) =
-            showVar varID ++ " -> " ++ show varDom 
+
+translateUfaToDom ufa dtrB =
+    FA.composeThin ufa $ 
+        Map.fromAscList $ 
+            map mkToUnitUFA $ 
+                 DBox.toAscList dtrB
+    where
+    mkToUnitUFA (var, tr) =
+        (var, UFA.affine [co] (Map.singleton var [sl]))
+        where
+        sl = FA.domra2ranra ufa $ dtrToUnitSlope tr
+        co = FA.domra2ranra ufa $ dtrToUnitConst tr
+
+
 instance 
+    (UFA.ERUnitFnApprox box varid domra ranra ufa,
+     DomainBoxMappable dtrbox box varid (DomTransl domra) domra, 
+     DomainBoxMappable box dtrbox varid domra (DomTransl domra),
+     H.HTML ufa) =>
+    H.HTML (ERFnDomTranslApprox dtrbox varid ufa domra)
+    where
+    toHtml (ERFnDomTranslApprox ufa dtrB) =
+        H.toHtml $ translateUfaToDom ufa dtrB
+--        H.toHtml $
+--            abovesTable []
+--                [
+--                    H.toHtml $ "enclosure on dom = [" ++ dtrBShow dtrB ++ "]",
+--                    H.toHtml $ translateUfaToDom ufa dtrB
+--                ]
+
+instance 
     (UFA.ERUnitFnApprox box varid domra ranra ufa, 
      Eq dtrbox) =>
     Eq (ERFnDomTranslApprox dtrbox varid ufa domra)
@@ -258,7 +294,7 @@
 
 instance 
     (UFA.ERUnitFnApprox box varid domra ranra ufa
-    , DomainBoxMappable dtrbox box varid (DomTransl domra) domra, Eq dtrbox) =>
+    , DomainBoxMappable dtrbox box varid (DomTransl domra) domra, Eq dtrbox, Ord dtrbox) =>
     RA.ERApprox (ERFnDomTranslApprox dtrbox varid ufa domra)
     where
     initialiseBaseArithmetic _ =
@@ -286,11 +322,18 @@
         vars = DBox.keys dtrB
     leqReals fa1 fa2 =
         RA.leqReals (erfnUnitApprox fa1) (erfnUnitApprox fa2)
-
+    refines fa1 fa2 =
+        RA.refines (erfnUnitApprox fa1) (erfnUnitApprox fa2)
+    compareApprox (ERFnDomTranslApprox ufa1 dtrB1) (ERFnDomTranslApprox ufa2 dtrB2) =
+        compareComposeMany
+            [ 
+                RA.compareApprox ufa1 ufa2,
+                compare dtrB1 dtrB2
+            ]
 
 instance
     (UFA.ERUnitFnApprox box varid domra ranra ufa, RA.ERIntApprox ufa 
-    , DomainBoxMappable dtrbox box varid (DomTransl domra) domra, Eq dtrbox) =>
+    , DomainBoxMappable dtrbox box varid (DomTransl domra) domra, Eq dtrbox, Ord dtrbox) =>
     RA.ERIntApprox (ERFnDomTranslApprox dtrbox varid ufa domra)
     where
 --    doubleBounds = :: ira -> (Double, Double) 
@@ -313,7 +356,7 @@
 
 instance
     (UFA.ERUnitFnApprox box varid domra ranra ufa, RAEL.ERApproxElementary ufa
-    , DomainBoxMappable dtrbox box varid (DomTransl domra) domra, Eq dtrbox) =>
+    , DomainBoxMappable dtrbox box varid (DomTransl domra) domra, Eq dtrbox, Ord dtrbox) =>
     RAEL.ERApproxElementary (ERFnDomTranslApprox dtrbox varid ufa domra)
     where
     abs ix (ERFnDomTranslApprox ufa dtrB) =
@@ -331,10 +374,10 @@
 
 instance 
     (UFA.ERUnitFnApprox box varid domra ranra ufa, 
-     DomainIntBox box varid domra, 
      DomainBoxMappable dtrbox box varid (DomTransl domra) domra, 
+     DomainIntBox box varid domra, 
      DomainBoxMappable box dtrbox varid domra (DomTransl domra), 
-     Eq dtrbox) =>
+     Eq dtrbox, Ord dtrbox) =>
     FA.ERFnApprox box varid domra ranra (ERFnDomTranslApprox dtrbox varid ufa domra)
     where
     check prgLocation (ERFnDomTranslApprox ufa dtrB) =
@@ -345,6 +388,10 @@
         FA.ranra2domra (erfnUnitApprox fa) r
     setMaxDegree maxDegree (ERFnDomTranslApprox ufa dtrB) =
         ERFnDomTranslApprox (FA.setMaxDegree maxDegree ufa) dtrB
+    getMaxDegree (ERFnDomTranslApprox ufa _) =
+        FA.getMaxDegree ufa
+    getRangeApprox (ERFnDomTranslApprox ufa dtrB) =
+        FA.getRangeApprox ufa
     volume (ERFnDomTranslApprox ufa dtrB) =
         DBox.fold 
             (\tr vol -> vol * (FA.domra2ranra ufa $ dtrFromUnitSlope tr)) 
@@ -355,12 +402,12 @@
         (ERFnDomTranslApprox (FA.scale ratio ufa) dtrB)
     partialIntersect ix substitutions f1 f2 
         | insideSubstitutions = f1 RA./\ f2
-        | otherwise = f1
+        | otherwise = f2
         where
         insideSubstitutions =
             and $ map snd $
                 DBox.zipWith (RA.refines) dom1 substitutions
-        dom1 = FA.dom f1
+        dom1 = FA.dom f2
     eval ptBox (ERFnDomTranslApprox ufa dtrB) =
         FA.eval (domToUnit dtrB ptBox) ufa
     partialEval substitutions (ERFnDomTranslApprox ufa dtrB) =
@@ -368,6 +415,7 @@
         where
         dtrBNoVars =
             DBox.difference dtrB substitutions
+    
 
 --instance 
 --    (UFA.ERUnitFnApprox box varid domra ranra ufa, 
@@ -393,7 +441,7 @@
      DomainIntBox box varid domra,
      DomainBoxMappable dtrbox box varid (DomTransl domra) domra, 
      DomainBoxMappable box dtrbox varid domra (DomTransl domra), 
-     Eq dtrbox) =>
+     Eq dtrbox, Ord dtrbox) =>
     FA.ERFnDomApprox box varid domra ranra (ERFnDomTranslApprox dtrbox varid ufa domra)
     where
     dom (ERFnDomTranslApprox ufa dtrB) = dtrBToDomB dtrB
diff --git a/src/Data/Number/ER/RnToRm/Approx/PieceWise.hs b/src/Data/Number/ER/RnToRm/Approx/PieceWise.hs
--- a/src/Data/Number/ER/RnToRm/Approx/PieceWise.hs
+++ b/src/Data/Number/ER/RnToRm/Approx/PieceWise.hs
@@ -41,6 +41,8 @@
 import Data.Number.ER.BasicTypes
 import Data.Number.ER.Misc
 
+import qualified Text.Html as H
+
 import Data.Typeable
 import Data.Generics.Basics
 import Data.Binary
@@ -133,6 +135,13 @@
     show f@(ERFnPiecewise bistr) =
         "\nERFnPiecewise:" ++ show bistr
 
+instance 
+    (FA.ERFnDomApprox box varid domra ranra fa, RA.ERIntApprox fa, H.HTML fa) =>
+    H.HTML (ERFnPiecewise box varid domra fa)
+    where
+    toHtml (ERFnPiecewise bistr) =
+        H.toHtml bistr
+
 instance
     (FA.ERFnDomApprox box varid domra ranra fa) =>
     Eq (ERFnPiecewise box varid domra fa)
@@ -202,7 +211,10 @@
         leqTuple [] = Just True
         leqTuple _ = 
             error $ "ERFnTuple: leqReals not implemented" 
-
+    refines f1@(ERFnPiecewise bistr1) f2@(ERFnPiecewise bistr2) =
+        and $ BISTR.collectValues $ pwbistrZipWith (RA.refines) 10 bistr1 bistr2
+    compareApprox f1@(ERFnPiecewise bistr1) f2@(ERFnPiecewise bistr2) =
+        BISTR.compare RA.compareApprox RA.compareApprox bistr1 bistr2
 instance 
     (FA.ERFnDomApprox box varid domra ranra fa, RA.ERIntApprox fa, VariableID varid) =>
     RA.ERIntApprox (ERFnPiecewise box varid domra fa)
@@ -287,17 +299,31 @@
             [FA.partialIntersect ix substitutions val1 val2]
         updateTouch = updateInside
         updateAway dom [val1, val2] =
-            [val1]
+            [val2]
     eval ptBox (ERFnPiecewise bistr) =
-        foldl1 (zipWith (RA.\/)) $ 
+        foldl1 (zipWith (RA./\)) $ 
+--            unsafePrintReturn ("ERFnPiecewise: eval: vals = ")$
             map (\fa -> FA.eval ptBox fa) $
-                BISTR.collectValues $ BISTR.lookupSubtreeDom bistr ptBox 
+                concat $ map BISTR.collectValues $ 
+                    BISTR.lookupSubtreeDoms bistr ptBox 
     partialEval substitutions f@(ERFnPiecewise bistr) = 
         pwLift1 (FA.partialEval substitutions) (ERFnPiecewise bistrNoVars)
         where
         bistrNoVars =
             BISTR.removeVars substitutions bistr
-
+    intersectMeasureImprovement ix f1@(ERFnPiecewise bistr1) f2@(ERFnPiecewise bistr2) =
+        (intersection, improvementRA)
+        where
+        (intersection, _) = RA.intersectMeasureImprovement ix f1 f2
+        improvementRA 
+            | 0 `RA.refines` intersectionVolume && 0 `RA.refines` f1Volume = 1
+--                error $ 
+--                    "ERFnInterval: intersectMeasureImprovement: inconsistent result: " 
+--                    ++ show intersection
+            | otherwise =
+                 f1Volume / intersectionVolume
+        intersectionVolume = FA.volume intersection
+        f1Volume = FA.volume f1
 instance 
     (FA.ERFnDomApprox box varid domra ranra fa, RA.ERIntApprox fa, Show box, 
      DomainBoxMappable box box varid domra domra) =>
@@ -323,9 +349,19 @@
                 Just pt -> pt
     unBisect var (ERFnPiecewise bistr1, ERFnPiecewise bistr2) =
         ERFnPiecewise $ 
-            BISTR.Node (depth1 - 1) dom var domVarMid bistr1 bistr2
+            BISTR.Node (depth - 1) dom var domVarMid bistr1Dp bistr2Dp
         where
+        depth = max depth1 depth2
         depth1 = BISTR.bistrDepth bistr1
+        depth2 = BISTR.bistrDepth bistr2
+        bistr1Dp 
+            | depth1 == depth = bistr1
+            | otherwise =
+                BISTR.setDepth depth bistr1
+        bistr2Dp 
+            | depth2 == depth = bistr2
+            | otherwise =
+                BISTR.setDepth depth bistr2
         dom1 = BISTR.bistrDom bistr1
         dom2 = BISTR.bistrDom bistr2
         dom = DBox.unionWith (RA.\/) dom1 dom2
@@ -360,13 +396,13 @@
                     (Just faLO, Nothing) ->
                         FA.partialIntersect ix 
                             (DBox.singleton x domLO) 
-                            bottomFn 
                             (ERFnPiecewise (BISTR.Leaf depth domB faLO)) 
+                            bottomFn 
                     (Nothing, Just faHI) -> 
                         FA.partialIntersect ix 
                             (DBox.singleton x domHI) 
-                            bottomFn 
                             (ERFnPiecewise (BISTR.Leaf depth domB faHI))
+                            bottomFn 
             bottomFn =
                 ERFnPiecewise $ BISTR.Leaf depth domB $ FA.bottomApprox domB (FA.getTupleSize fD)
             (domLO, domHI) = 
@@ -442,13 +478,13 @@
                     (Just faLO, Nothing) ->
                         FA.partialIntersect ix 
                             (DBox.singleton x domLO) 
-                            (ERFnPiecewise bistrP) 
                             (ERFnPiecewise (BISTR.Leaf depth domB faLO)) 
+                            (ERFnPiecewise bistrP) 
                     (Nothing, Just faHI) -> 
                         FA.partialIntersect ix 
                             (DBox.singleton x domHI) 
-                            (ERFnPiecewise bistrP) 
                             (ERFnPiecewise (BISTR.Leaf depth domB faHI)) 
+                            (ERFnPiecewise bistrP) 
             (domLO, domHI) = 
                 RA.bounds $ 
                     DBox.lookup "ERFnPieceWise: integrateMeasureImprovement: zipOutsideRange: " x domB
@@ -497,7 +533,7 @@
             (faIsect, faImpr) = 
                 FA.integrateMeasureImprovement ix faD x integdomBox domLO faPadj
             faPadj =
-                FA.partialIntersect ix (DBox.singleton x domLO) faP faLO
+                FA.partialIntersect ix (DBox.singleton x domLO) faLO faP
             faHI =
                 FA.partialEval (DBox.singleton x domHI) faIsect
             (domLO, domHI) = 
@@ -516,7 +552,7 @@
             (faIsect, faImpr) = 
                 FA.integrateMeasureImprovement ix faD x integdomBox domHI faPadj
             faPadj =
-                FA.partialIntersect ix (DBox.singleton x domHI) faP faHI
+                FA.partialIntersect ix (DBox.singleton x domHI) faHI faP
             faLO =
                 FA.partialEval (DBox.singleton x domLO) faIsect
             (domLO, domHI) = 
diff --git a/src/Data/Number/ER/RnToRm/Approx/Tuple.hs b/src/Data/Number/ER/RnToRm/Approx/Tuple.hs
--- a/src/Data/Number/ER/RnToRm/Approx/Tuple.hs
+++ b/src/Data/Number/ER/RnToRm/Approx/Tuple.hs
@@ -27,7 +27,11 @@
 import qualified Data.Number.ER.Real.Approx.Elementary as RAEL
 import qualified Data.Number.ER.Real.DomainBox as DBox
 import Data.Number.ER.BasicTypes
+import Data.Number.ER.Misc
 
+import Data.Number.ER.ShowHTML
+import qualified Text.Html as H
+
 import Data.Typeable
 import Data.Generics.Basics
 import Data.Binary
@@ -117,6 +121,18 @@
         showFA (fnname, fa) =
             "\n>>> Function " ++ show fnname ++ ":" ++ show fa
 
+instance 
+    (FA.ERFnDomApprox box varid domra ranra fa, H.HTML fa) =>
+    H.HTML (ERFnTuple fa)
+    where
+    toHtml (ERFnTuple []) = H.toHtml "[]"
+--    toHtml (ERFnTuple [f]) = H.toHtml f
+    toHtml (ERFnTuple fs) =
+        H.toHtml $
+            abovesTable [H.border 2] $
+--                [H.toHtml "Function tuple"] ++ 
+                map H.toHtml fs
+
 instance
     (FA.ERFnApprox box varid domra ranra fa) =>
     Eq (ERFnTuple fa)
@@ -159,6 +175,8 @@
     setGranularity gran = tuplesLift1 (RA.setGranularity gran) 
     setMinGranularity gran = tuplesLift1 (RA.setMinGranularity gran)
     f1 /\ f2 = tuplesLift2 "ERFnTuple: /\\: " (RA./\) f1 f2
+    refines f1@(ERFnTuple fas1) f2@(ERFnTuple fas2) =
+        and $ zipWith RA.refines fas1 fas2
     intersectMeasureImprovement ix f1@(ERFnTuple fas1) f2@(ERFnTuple fas2)
         | length fas1 == length fas2 =
             (ERFnTuple fasIsect, ERFnTuple fasImpr)
@@ -174,7 +192,9 @@
         where
         leqTuple [] = Just True
         leqTuple _ = 
-            error $ "ERFnTuple: leqReals not implemented" 
+            error $ "ERFnTuple: leqReals not implemented"
+    compareApprox f1@(ERFnTuple fas1) f2@(ERFnTuple fas2) =
+        compareListsWith RA.compareApprox fas1 fas2
 
 instance 
     (FA.ERFnDomApprox box varid domra ranra fa, RA.ERIntApprox fa) =>
diff --git a/src/Data/Number/ER/RnToRm/BisectionTree.hs b/src/Data/Number/ER/RnToRm/BisectionTree.hs
--- a/src/Data/Number/ER/RnToRm/BisectionTree.hs
+++ b/src/Data/Number/ER/RnToRm/BisectionTree.hs
@@ -27,6 +27,7 @@
     removeVars,
     sync2,
     syncMany,
+    setDepth,
     split,
     mapWithDom,
     mapLeaves,
@@ -36,7 +37,8 @@
     combineWith,
     collectValues,
     collectDomValues,
-    lookupSubtreeDom
+    compare,
+    lookupSubtreeDoms
 )
 where
 
@@ -50,6 +52,9 @@
 
 import Data.Number.ER.Misc
 
+import Data.Number.ER.ShowHTML
+import qualified Text.Html as H
+
 import Data.Typeable
 import Data.Generics.Basics
 import Data.Binary
@@ -127,6 +132,31 @@
     showVD (v,d) =
         showVar v ++ "->" ++ show d
 
+instance (Show d, H.HTML v, DomainBox box varid d) => 
+    H.HTML (BisectionTree box varid d v)
+    where
+    toHtml (Leaf depth dom val) =
+        H.toHtmlFromList $
+            [
+                H.toHtml $ concatWith "," (Prelude.map showVD $ DBox.toList dom)
+            ,
+                H.primHtml " &rarr; "
+            , 
+                H.toHtml val
+            ]
+        where
+        showVD (v,d) =
+            showVar v ++ " in " ++ show d
+    toHtml (Node depth dom dir pt lo hi) =
+        H.toHtml $
+        besidesTable [H.border 2]
+            [
+             abovesTable [] [H.toHtml $ "(" ++ showVar dir ++ ")"]
+            ,
+             abovesTable [] [H.toHtml lo, H.toHtml hi]
+            ]
+
+
 isLeaf ::
     BisectionTree box varid d v ->
     Bool     
@@ -153,6 +183,23 @@
     (EffortIndex -> Depth -> (BisectionTree box varid d v) -> v) 
     
     
+setDepth ::
+    Depth ->
+    BisectionTree box varid d v ->
+    BisectionTree box varid d v
+setDepth depth bistr 
+    | isLeaf bistr =
+        bistr { bistrDepth = depth }
+    | otherwise =
+        bistr
+        {
+            bistrLO = setDepth depthInc $ bistrLO bistr,
+            bistrHI = setDepth depthInc $ bistrHI bistr
+        }
+    where
+    depthInc = depth + 1
+
+    
 split ::
     (RA.ERIntApprox d, DomainBox box varid d) =>
     ValueSplitter box varid d v ->
@@ -559,41 +606,55 @@
     linear ordering on bisection trees
 -}        
 compare ::
-    (Ord d, Ord varid) =>
+    (Ord varid, DomainBox box varid d) =>
+    (d -> d -> Ordering) ->
     (v -> v -> Ordering) ->
     (BisectionTree box varid d v) ->
     (BisectionTree box varid d v) ->
     Ordering
-compare compValues (Leaf _ _ _) (Node _ _ _ _ _ _) = LT
-compare compValues (Node _ _ _ _ _ _) (Leaf _ _ _) = GT
-compare compValues (Leaf _ _ val1) (Leaf _ _ val2) =
-    compValues val1 val2
-compare compValues 
-        (Node _ _ dir1 pt1 lo1 hi1) 
-        (Node _ _ dir2 pt2 lo2 hi2) =
-    compareComposeMany $
-    [Prelude.compare dir1 dir2,
-     Prelude.compare pt1 pt2,
-     compare compValues lo1 lo2,
-     compare compValues hi1 hi2]
+compare compDoms compValues (Leaf _ _ _) (Node _ _ _ _ _ _) = LT
+compare compDoms compValues (Node _ _ _ _ _ _) (Leaf _ _ _) = GT
+compare compDoms compValues (Leaf depth1 dom1 val1) (Leaf depth2 dom2 val2) =
+    compareComposeMany 
+        [
+            Prelude.compare depth1 depth2,
+            DBox.compare compDoms dom1 dom2,
+            compValues val1 val2
+        ]
+compare compDoms compValues 
+        (Node depth1 dom1 dir1 pt1 lo1 hi1) 
+        (Node depth2 dom2 dir2 pt2 lo2 hi2) =
+    compareComposeMany 
+        [
+            Prelude.compare dir1 dir2,
+            compDoms pt1 pt2,
+            compare compDoms compValues lo1 lo2,
+            compare compDoms compValues hi1 hi2
+        ]
 
 {-|
-    lookup the smallest subtree whose domain covers the given rectangle
+    lookup all maximal subtrees whose domain intersect the given rectangle
 -}
-lookupSubtreeDom ::
+lookupSubtreeDoms ::
     (RA.ERIntApprox d, DomainBox box varid d) =>
     (BisectionTree box varid d v) ->
     box {-^ domain to look up within the tree -} ->
-    (BisectionTree box varid d v)
-lookupSubtreeDom origBistr dom = 
+    [BisectionTree box varid d v]
+lookupSubtreeDoms origBistr dom = 
     lk origBistr
     where
-    lk bistr@(Leaf _ _ _) = bistr
+    lk bistr@(Leaf _ _ _) = [bistr]
     lk bistr@(Node _ _ _ _ lo hi)
-        | and $ Prelude.map snd $ DBox.zipWithDefault RA.bottomApprox (RA.refines) dom domHI = lk hi
-        | and $ Prelude.map snd $ DBox.zipWithDefault RA.bottomApprox (RA.refines) dom domLO = lk lo
-        | otherwise = bistr
+        | loDisjoint = lk hi
+        | hiDisjoint = lk lo
+        | otherwise = (lk lo) ++ (lk hi)
         where
+        loDisjoint =
+            and $ Prelude.map snd $ 
+                DBox.zipWithDefault RA.bottomApprox (RA.isDisjoint) dom domLO  
+        hiDisjoint =
+            and $ Prelude.map snd $ 
+                DBox.zipWithDefault RA.bottomApprox (RA.isDisjoint) dom domHI  
         domLO = bistrDom lo
         domHI = bistrDom hi
 
diff --git a/src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs b/src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs
--- a/src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs
+++ b/src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs
@@ -41,6 +41,9 @@
 
 import Data.Number.ER.Misc
 
+import Data.Number.ER.ShowHTML
+import qualified Text.Html as H
+
 import qualified Data.Map as Map
 
 import Data.Typeable
@@ -120,6 +123,24 @@
 
 instance
     (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>
+    H.HTML (ERFnInterval fb ra) 
+    where
+    toHtml (ERFnIntervalAny ctxt) =
+        H.toHtml "ERFnIntervalAny"
+    toHtml (ERFnInterval h ln 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)]
+                        ]
+--                ]
+
+instance
+    (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>
     Eq (ERFnInterval fb ra)
     where
     (ERFnInterval h1 ln1 ctxt1 gl1) 
@@ -135,9 +156,9 @@
     compare 
             (ERFnInterval h1 ln1 ctxt1 gl1) 
             (ERFnInterval h2 ln2 ctxt2 gl2) =
-        error "ERFnInterval: comparison not implemented; consider leqReals from class ERApprox instead"
+        error "ERFnInterval: comparison not implemented; consider leqReals or compareApprox from class ERApprox instead"
     compare _ _ =
-        error "ERFnInterval: comparison not implemented; consider leqReals from class ERApprox instead"
+        error "ERFnInterval: comparison not implemented; consider leqReals or compareApprox from class ERApprox instead"
     
     
 instance 
@@ -149,6 +170,7 @@
     negate (ERFnInterval h ln ctxt gl) =
         (ERFnInterval ln h ctxt (negate gl))
     (ERFnInterval h1 ln1 ctxt1 gl1) + (ERFnInterval h2 ln2 ctxt2 gl2) =
+        normalise $
         ERFnInterval (h1 + h2) (ln1 + ln2) ctxt (gl1 + gl2)
         where
         ctxt = erfnContextUnify ctxt1 ctxt2
@@ -156,6 +178,7 @@
         where
         ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2)
     (ERFnInterval h1 ln1 ctxt1 gl1) * (ERFnInterval h2 ln2 ctxt2 gl2) =
+        normalise $
         ERFnInterval h ln ctxt (gl1 * gl2)
         where
         (h, ln) =
@@ -196,6 +219,7 @@
     recip f@(ERFnIntervalAny _) = f
     recip (ERFnInterval h ln ctxt gl) 
         | certainNoZero =
+            normalise $
             ERFnInterval lRecipUp hnRecipUp ctxt (recip gl)
         | otherwise = ERFnIntervalAny ctxt
         where
@@ -237,7 +261,11 @@
             (ctxt { erfnCoeffGranularity = max gran (erfnCoeffGranularity ctxt) }) gl
 --    getPrecision (ERFnIntervalAny _) = 0
 --    getPrecision f = intLog 2 (1 + (fst $ RA.integerBounds (FA.volume f))) -- wrong! 
-    (ERFnInterval h1 ln1 ctxt1 gl1) /\ (ERFnInterval h2 ln2 ctxt2 gl2) =
+    f1@(ERFnInterval h1 ln1 ctxt1 gl1) /\ f2@(ERFnInterval h2 ln2 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)
         where
         ctxt = erfnContextUnify ctxt1 ctxt2
@@ -254,6 +282,21 @@
         where
         ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2)
     leqReals = erfnintLeq
+    refines _ (ERFnIntervalAny _) = True
+    refines (ERFnIntervalAny _) _ = False
+    refines (ERFnInterval h1 ln1 _ _) (ERFnInterval h2 ln2 _ _) = 
+        (UFB.upperBound 10 (h2 - h1) >= 0)
+        &&
+        (UFB.upperBound 10 (ln2 - ln1) >= 0)
+    compareApprox (ERFnIntervalAny _) (ERFnIntervalAny _) = EQ
+    compareApprox (ERFnIntervalAny _) _ = LT
+    compareApprox _ (ERFnIntervalAny _) = GT
+    compareApprox (ERFnInterval h1 ln1 _ _) (ERFnInterval h2 ln2 _ _) =
+        compareComposeMany
+        [
+            UFB.compareApprox h1 h2,
+            UFB.compareApprox ln1 ln2
+        ]
 
 erfnintLeq left right
     | left `isClearlyBelow` right = Just True
@@ -295,7 +338,11 @@
     bounds (ERFnInterval u ln c g) =
         (ERFnInterval (negate ln) ln c g,
          ERFnInterval u (negate u) c g) 
-    (ERFnInterval u1 ln1 c1 g1) \/ (ERFnInterval u2 ln2 c2 g2) =
+    f1@(ERFnInterval u1 ln1 c1 g1) \/ f2@(ERFnInterval u2 ln2 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)
         where
         u = UFB.maxUp maxDegree u1 u2
@@ -320,7 +367,11 @@
     where
     -- default abs does not work because we do not have Prelude.abs
     abs _ f@(ERFnIntervalAny _) = f
-    abs _ (ERFnInterval u ln c g) =
+    abs _ f@(ERFnInterval u ln c g) =
+---- #ifdef RUNTIME_CHECKS
+----         check ("ERFnInterval: abs:\n f:\n" ++ show f ++ "\n result:\n") $
+---- #endif
+        normalise $
         ERFnInterval maxulnUp maxunl0Dn c (abs g)
         where
         maxDegree = erfnMaxDegree c
@@ -329,17 +380,17 @@
             fst $ UFB.max maxDegree 0 $
                 fst $ UFB.max maxDegree (- u) (- ln)
     exp ix f@(ERFnIntervalAny _) = f
-    exp ix (ERFnInterval u ln c g) =
+    exp ix f@(ERFnInterval u ln c g) = 
+        normalise $
         ERFnInterval uExp lExpNeg c (RAEL.exp ix g)
         where
         maxDegree = erfnMaxDegree c
---        ix = int2effIx maxDegree
         uExp = snd $ UFB.exp maxDegree ix u
         lExpNeg = 
             negate $ fst $ UFB.exp maxDegree ix (negate ln) 
     sin ix f@(ERFnIntervalAny c) = 
         ERFnInterval 1 1 c ((-1) RA.\/ 1)
-    sin ix (ERFnInterval u ln c g) =
+    sin ix f@(ERFnInterval u ln c g) =
 --        unsafePrint
 --        (
 --            "ERFnInterval: RAEL.sin: "
@@ -348,13 +399,17 @@
 --            ++ "\n uSin = " ++ show uSin
 --            ++ "\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)
         where
         (lSin, uSin) = sincos True maxDegree ix u (-ln)  
         maxDegree = erfnMaxDegree c
     cos ix f@(ERFnIntervalAny c) =
         ERFnInterval 1 1 c ((-1) RA.\/ 1)
-    cos ix (ERFnInterval u ln c g) =
+    cos ix f@(ERFnInterval u ln c g) =
 --        unsafePrint
 --        (
 --            "ERFnInterval: RAEL.cos: "
@@ -363,13 +418,14 @@
 --            ++ "\n uCos = " ++ show uCos
 --            ++ "\n lCosNeg = " ++ show lCosNeg
 --        ) $
+        normalise $
         ERFnInterval uCos (- lCos) c (RAEL.cos ix g)
         where
         (lCos, uCos) = sincos False maxDegree ix u (-ln) 
         maxDegree = erfnMaxDegree c
     atan ix f@(ERFnIntervalAny c) =
         ERFnInterval 1 1 c ((-1) RA.\/ 1)
-    atan ix (ERFnInterval u ln c g) =
+    atan ix f@(ERFnInterval u ln c g) =
 --        unsafePrint
 --        (
 --            "ERFnInterval: RAEL.atan: "
@@ -378,6 +434,7 @@
 --            ++ "\n uAtan = " ++ show uAtan
 --            ++ "\n lAtanNeg = " ++ show lAtanNeg
 --        ) $
+        normalise $
         ERFnInterval uAtan lAtanNeg c (RAEL.atan ix g)
         where
         maxDegree = erfnMaxDegree c
@@ -514,18 +571,25 @@
         (pLO `plusDown` (- 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
+    | otherwise = ERFnIntervalAny c 
     
+check callerLocation f@(ERFnIntervalAny c) = f
+check callerLocation f@(ERFnInterval u ln c g) =
+    ERFnInterval 
+        (UFB.check (callerLocation ++ "upper: ") u) 
+        (UFB.check (callerLocation ++ "neg lower: ") ln) 
+        c g 
 
+
 instance 
     (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>
     FA.ERFnApprox boxra varid ra ra (ERFnInterval fb ra)
     where
-    check callerLocation f@(ERFnIntervalAny c) = f
-    check callerLocation (ERFnInterval u ln c g) =
-        ERFnInterval 
-            (UFB.check (callerLocation ++ "upper: ") u) 
-            (UFB.check (callerLocation ++ "neg lower: ") ln) 
-            c g 
+    check = check
     domra2ranra _ = id
     ranra2domra _ = id
     setMaxDegree maxDegr (ERFnIntervalAny c) =
@@ -536,16 +600,26 @@
             (UFB.reduceDegreeUp maxDegr ln)
             (c { erfnMaxDegree = maxDegr } )
             g
-    getRangeApprox (ERFnIntervalAny _) = RA.bottomApprox 
+    getMaxDegree (ERFnIntervalAny c) =
+        erfnMaxDegree c
+    getMaxDegree (ERFnInterval _ _ c _) =
+        erfnMaxDegree c
+    getRangeApprox (ERFnIntervalAny _) = 
+        RA.bottomApprox 
     getRangeApprox (ERFnInterval u ln c g) =
         UFB.raFromEndpoints u
         (
-         (UFB.upperBound 10 u)
-        ,
          (- (UFB.upperBound 10 ln))
+        ,
+         (UFB.upperBound 10 u)
         )
-    scale ratio f@(ERFnIntervalAny c) = f
+    scale ratio f@(ERFnIntervalAny c) = 
+        f
     scale ratio f@(ERFnInterval u ln 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
@@ -560,7 +634,8 @@
         up = UFB.evalApprox ptBox u
         lo = negate $ UFB.evalApprox ptBox ln
     partialEval substitutions f@(ERFnIntervalAny c) = f
-    partialEval substitutions (ERFnInterval u ln c g) =
+    partialEval substitutions f@(ERFnInterval u ln c g) =
+        normalise $
         (ERFnInterval uP lnP c g)
         where
         uP = UFB.partialEvalApproxUp substitutions u
@@ -589,6 +664,10 @@
         ERFnIntervalAny erfnContextDefault
     const [val] 
         | RA.isBounded val =
+---- #ifdef RUNTIME_CHECKS
+----             check ("ERFnInterval: const:\n") $
+---- #endif
+            normalise $
             ERFnInterval
             {
                 erfnUpper = fbH,
@@ -609,6 +688,10 @@
             }
     affine [val] coeffsSingletons
         | RA.isBounded val && (and $ map (RA.isBounded . head) $ Map.elems coeffsSingletons) =
+---- #ifdef RUNTIME_CHECKS
+----             check ("ERFnInterval: affine:\n") $
+---- #endif
+            normalise $
             ERFnInterval
             {
                 erfnUpper = fbH,
@@ -665,7 +748,12 @@
             Just LT -> (f1, 1) -- intersection made it worse, keep original
             _ ->  (intersection, improvementRA)
         where
-        intersection = f1 RA./\ f2
+        intersection = 
+---- #ifdef RUNTIME_CHECKS
+----             check ("ERFnInterval: intersectMeasureImprovement:\n f1:\n" ++ show f1 ++ "\n f2:\n" ++ show f2 ++ "\n intersection:\n") $
+---- #endif
+            normalise $
+            f1 RA./\ f2
         improvementRA 
             | 0 `RA.refines` intersectionVolume && 0 `RA.refines` f1Volume = 1
 --                error $ 
@@ -682,9 +770,10 @@
 --        where
 --        result =
             UFB.raFromEndpoints u $ UFB.volumeAboveZero vars (u + ln)
+    integrate _ f@(ERFnIntervalAny c) _ _ _ = f 
     integrate 
-            ix (ERFnInterval u ln c g) x 
-            origin (ERFnInterval uInit lnInit cInit gInit) =
+            ix fD@(ERFnInterval u ln c g) x 
+            origin fI@(ERFnInterval uInit lnInit cInit gInit) =
 --        unsafePrint
 --        (
 --            "ERFnInterval: integrate: " 
@@ -704,6 +793,10 @@
 --            ++ "\n uIov = " ++ show uIov
 --            ++ "\n lnIov = " ++ show lnIov
 --        )
+---- #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)
         where
         -- perform raw integration of both bounds:
diff --git a/src/Data/Number/ER/RnToRm/UnitDom/Base.hs b/src/Data/Number/ER/RnToRm/UnitDom/Base.hs
--- a/src/Data/Number/ER/RnToRm/UnitDom/Base.hs
+++ b/src/Data/Number/ER/RnToRm/UnitDom/Base.hs
@@ -39,7 +39,15 @@
     where
     initialiseBaseArithmetic :: ufb -> IO ()
     initialiseBaseArithmetic _ =
-    	B.initialiseBaseArithmetic (0 :: b)  
+    	B.initialiseBaseArithmetic (0 :: b)
+    {-|
+        Check internal consistency, typically absence of NaN.
+    -}
+    isValid :: ufb -> Bool
+    {-|
+        A linear ordering, which can be syntactic and rather arbitrary. 
+    -}
+    compareApprox :: ufb -> ufb -> Ordering
     {-| 
         Check internal consistency of the function and report problem if any.
     -}
diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs
--- a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs
+++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs
@@ -59,7 +59,9 @@
      DomainIntBox boxra varid (ERInterval rb)) =>
     (UFB.ERUnitFnBase boxb boxra varid rb (ERInterval rb) (ERChebPoly box rb))
     where
+    isValid = chplHasNoNaNOrInfty
     check = chplCheck
+    compareApprox = chplCompareApprox
     getGranularity = chplGetGranularity
     setMinGranularity = chplSetMinGranularity
     setGranularity = chplSetGranularity
diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Basic.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Basic.hs
--- a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Basic.hs
+++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Basic.hs
@@ -47,23 +47,36 @@
         chplCoeffs :: (Map.Map (TermKey box) b)
     }
     deriving (Eq, Typeable, Data)
+
+type TermKey box = box
     
 instance (Ord a, Binary a, Binary b) => Binary (ERChebPoly a b) where
   put (ERChebPoly a) = put a
   get = get >>= \a -> return (ERChebPoly a)
     
-    
-chplCheck prgLocation p@(ERChebPoly coeffs)
-    | ok = p
-    | otherwise = 
-        unsafePrint (prgLocation ++ " problem with p = \n" ++ show p) p
+chplHasNoNaN p@(ERChebPoly coeffs) =
+    Map.fold (&&) True $ Map.map coeffOK coeffs
     where
-    ok = 
-        Map.fold (&&) True $ Map.map coeffOK coeffs
     coeffOK c =
         not $ B.isERNaN c
-    
-type TermKey box = box
+
+chplHasNoNaNOrInfty p@(ERChebPoly coeffs) =
+    Map.fold (&&) True $ Map.map coeffOK coeffs
+    where
+    coeffOK = isBounded
+    isBounded c 
+        | B.isERNaN c = False
+        | B.isPlusInfinity c = False
+        | B.isPlusInfinity (- c) = False
+        | otherwise = True
+
+chplCheck prgLocation p
+    | chplHasNoNaNOrInfty p = p
+    | otherwise =
+        unsafePrint (prgLocation ++ " problem with p = \n" ++ show p) p
+
+chplCompareApprox (ERChebPoly coeffs1) (ERChebPoly coeffs2) =
+    compare coeffs1 coeffs2 
     
 chplConstTermKey :: (DomainBox box varid d) => box
 chplConstTermKey = DBox.noinfo
@@ -87,6 +100,20 @@
     d
 chplGetDegree (ERChebPoly coeffs) =
     foldl max 0 $ map chplTermOrder $ Map.keys coeffs
+
+{-|
+    If the polynomial is constant, return the constant,
+    otherwise return Nothing.
+-}
+chplGetConst ::
+    (B.ERRealBase b, DomainBox box varid d, Num d, Ord d) =>
+    (ERChebPoly box b) ->
+    Maybe b
+chplGetConst (ERChebPoly coeffs) =
+    case Map.keys coeffs of
+        [key] | chplIsConstTermKey key ->
+            Just $ head $ Map.elems coeffs
+        _ -> Nothing
 
     
 -- chplGetArity = length . chplGetVars  
diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs
--- a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs
+++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs
@@ -57,6 +57,17 @@
     ERChebPoly box b ->
     (ERChebPoly box b, ERChebPoly box b)
 chplExp maxDegree ix p =
+--    unsafePrint
+--    ( 
+--        "chplExp:" ++ 
+--        "\n expM = " ++ show expM ++
+--        "\n upperB = " ++ show upperB ++
+--        "\n lowerB = " ++ show lowerB ++
+--        "\n a_int = " ++ show a_int ++
+--        "\n expNear0Dn pNear0Dn = " ++ show (expNear0Dn pNear0Dn) ++
+--        "\n chplPow maxDegree (expNear0Up pNear0Up) 2000 = " ++ show (chplPow maxDegree (expNear0Up pNear0Up) 2000)
+--    ) 
+--    $ 
     (expDownwards, expUpwards + valueAtRDnNeg + (chplConst expRUp))
     where
     expUpwards =
@@ -69,7 +80,8 @@
     r = (upperB - lowerB) / 2 
     expMUp = erintv_right expM 
     expMDn = erintv_left expM
-    expM = erExp_R ix (ERInterval m m)
+    expM =
+        erExp_R ix (ERInterval m m)
     pNear0Up = (p - (chplConst m)) * (chplConst $ recip a_base)
     pNear0Dn = - (((-p) + (chplConst m)) * (chplConst $ recip a_base))
     a_base = fromInteger a_int
@@ -449,7 +461,9 @@
     EffortIndex {-^ minimum approx degree -} -> 
     ERChebPoly box b ->
     (ERChebPoly box b, ERChebPoly box b)
-chplRecip maxDegree ix p
+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
@@ -470,7 +484,12 @@
     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
     
diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Field.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Field.hs
--- a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Field.hs
+++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Field.hs
@@ -86,7 +86,7 @@
             (Map.unionWith (+) coeffs1 coeffs2)
             -- point-wise sum of polynomials with coeffs rounded upwards
         coeffsDown =
-            (Map.unionWith (\c1 c2 -> - ((- c1) + (- c2))) coeffs1 coeffs2)
+            (Map.unionWith plusDown coeffs1 coeffs2)
             -- point-wise sum of polynomials with coeffs rounded upwards
         maxError =
             Map.fold (+) 0 $ 
@@ -199,6 +199,8 @@
                 (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)
@@ -212,7 +214,7 @@
 instance (B.ERRealBase b, DomainBox box varid Int, Ord box) => Ord (ERChebPoly box b)
     where
     compare _ _ =
-        errorModule "cannot compare polynomials, consider using RA.leqReals instead"
+        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
diff --git a/tests/ISin3.hs b/tests/ISin3.hs
--- a/tests/ISin3.hs
+++ b/tests/ISin3.hs
@@ -35,7 +35,7 @@
         head $ 
             AERNFunc.eval (AERNFunc.unary 1) $ 
                 AERNFunc.integrateUnary 0 
-                    (sin3 60 0 21 60) 
+                    (sin3 120 0 42 160) 
                     0 (0 AERN.\/ 1) [0]
     sin3 ix depth deg gran = 
         AERN.sin ix $