diff --git a/AERN-Real.cabal b/AERN-Real.cabal
--- a/AERN-Real.cabal
+++ b/AERN-Real.cabal
@@ -1,5 +1,5 @@
 Name:           AERN-Real
-Version:        0.9.8
+Version:        0.9.9
 Cabal-Version:  >= 1.2
 Build-Type:     Simple
 License:        BSD3
@@ -10,7 +10,7 @@
 Stability:      experimental
 Category:       Data, Math
 Synopsis:       arbitrary precision interval arithmetic for approximating exact real numbers
-Tested-with:    GHC ==6.8.3
+Tested-with:    GHC ==6.10.1
 Description:
     Datatypes and abstractions for approximating exact real numbers
     and a basic arithmetic over such approximations.
@@ -31,30 +31,18 @@
 Data-files:
     ChangeLog
 
-Flag containers-in-base
-    Default: False
-
 Flag use-hmpfr
     Default: False
 
 Library
   hs-source-dirs:  src
-  if flag(containers-in-base)
-    if flag(use-hmpfr)
+  if flag(use-hmpfr)
       Build-Depends:
-        base < 3, binary >= 0.4, html >= 1.0, regex-compat >= 0.71, stm, hmpfr == 0.1.3
+        base >= 3, base < 4, containers, binary >= 0.4, html >= 1.0, regex-compat >= 0.71, stm, time, hmpfr == 0.1.3
       cpp-options: -DUSE_MPFR
-    else
-      Build-Depends:
-        base < 3, binary >= 0.4, html >= 1.0, regex-compat >= 0.71, stm
   else
-    if flag(use-hmpfr)
       Build-Depends:
-        base >= 3, containers, binary >= 0.4, html >= 1.0, regex-compat >= 0.71, stm, hmpfr == 0.1.3
-      cpp-options: -DUSE_MPFR
-    else
-      Build-Depends:
-        base >= 3, containers, binary >= 0.4, html >= 1.0, regex-compat >= 0.71, stm
+        base >= 3, base < 4, containers, binary >= 0.4, html >= 1.0, regex-compat >= 0.71, stm, time
   Exposed-modules:
     Data.Number.ER,
     Data.Number.ER.Real,
diff --git a/ChangeLog b/ChangeLog
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,12 @@
+0.9.9: 23 February 2009
+    * Small changes needed in other AERN packages:
+        * New operation for domain boxes: get its dimension.
+    * Exponentiation, sine, cosine and arctan signinificantly improved for arguments further away from 0.
+    * Fixed a bug in sine Taylor series error term.
+    * Some interval arithmetic operations now have also "inner" versions
+      that approximate the maximal extension of the operation from inside
+      (useful for testing the normal "outer" versions).
+
 0.9.8: 1 December 2008
     * added instance of the HTML class for intervals
     * added syntactic comparison of variable-indexed domain boxes
diff --git a/src/Data/Number/ER/Misc.hs b/src/Data/Number/ER/Misc.hs
--- a/src/Data/Number/ER/Misc.hs
+++ b/src/Data/Number/ER/Misc.hs
@@ -14,6 +14,7 @@
 
 import Data.List
 import System.IO.Unsafe
+import Data.Time.Clock.POSIX
 
 unsafePrint msg val =
     unsafePerformIO $
@@ -24,6 +25,18 @@
 unsafePrintReturn msg a =
     unsafePrint (msg ++ show a) a
 
+unsafeReport fileName msg val =
+    unsafePerformIO $
+        do
+        stamp <- getPOSIXTime
+        appendFile fileName $ showStamp stamp ++ ":"
+        appendFile fileName $ msg ++ "\n"
+        return val
+    where
+    showStamp stamp =
+        padTo18 $ show stamp
+    padTo18 s = s ++ (replicate (18 - (length s)) ' ')
+
 {-|
     Compose as when defining the lexicographical ordering.
 -}
@@ -246,6 +259,9 @@
 plusUp, plusDown, timesUp, timesDown :: 
     (Num t) =>
     t -> t -> t
+divideUp, divideDown :: 
+    (Fractional t) =>
+    t -> t -> t
 sumUp, sumDown, productDown, productUp :: 
     (Num t) =>
     [t] -> t
@@ -257,6 +273,8 @@
 timesDown c1 c2 = - ((- c1) * c2)
 productUp = foldl timesUp 1
 productDown = foldl timesDown 1
+divideUp c1 c2 = c1 / c2
+divideDown c1 c2 = - ((- c1) / c2)
 
 {- parsing -}
 readMaybe :: (Read a) => String -> Maybe a
diff --git a/src/Data/Number/ER/Real/Approx.hs b/src/Data/Number/ER/Real/Approx.hs
--- a/src/Data/Number/ER/Real/Approx.hs
+++ b/src/Data/Number/ER/Real/Approx.hs
@@ -29,7 +29,6 @@
 (
     ERApprox(..),
     ERIntApprox(..),
-    bounds2ira,
     effIx2ra,
     splitIRA,
 --    checkShrinking,
@@ -38,7 +37,8 @@
     ltSingletons,
     equalIntervals,
     exactMiddle,
-    maxExtensionR2R
+    maxExtensionR2R,
+    maxExtensionInnerR2R
 )
 where
 
@@ -162,6 +162,8 @@
     defaultBisectPt :: ira -> ira
     -- | returns thin approximations of endpoints, in natural order 
     bounds :: ira -> (ira, ira)
+    -- | make an interval from thin approximations of endpoints 
+    fromBounds :: (ira, ira) -> ira
     {-|
          meet, usually constructing interval from approximations of its endpoints
          
@@ -171,14 +173,6 @@
     -}
     (\/) :: ira -> ira -> ira
     
-{-| 
-    Inverse of 'bounds'.
--}
-bounds2ira ::
-    (ERIntApprox ira) =>
-    (ira, ira) -> ira
-bounds2ira (a,b) = a \/ b
-    
 {-|
     Assuming the arguments are singletons, equality is decidable.
 -}
@@ -295,13 +289,59 @@
         {- ^ a function behaving well on sequences that intersect to a non-empty interval -}
 maxExtensionR2R getExtremes f ix x
     | getPrecision x < effIx2prec ix =
-        (f ix xL) \/ (f ix xR) \/ 
-        (foldl (\/) emptyApprox $ getExtremes ix x)
+        foldl1 (\/) $ [f ix xL, f ix xR] ++ (getExtremes ix x)
     -- x is thin enough (?), don't bother evaluating by endpoints and extrema:
     | otherwise =
         f ix x
     where
     (xL, xR) = bounds x
+        
+{-| 
+    This produces a function that computes the maximal extension of the
+    given function.  A maximal extension function has the property:
+    f(I) = { f(x) | x in I }.  Here we get this property only for the
+    limit function for its 'EffortIndex' tending to infinity.
+-}
+maxExtensionInnerR2R ::
+    (ERIntApprox ira) =>
+    (EffortIndex -> ira -> ([ira], (Maybe Bool, Maybe Bool)))
+        {-^ returns a safe approximation of all extrema within the interval
+            and an indication whether the function is increasing or decreasing 
+            at the endpoints of the queried real approximation -} ->
+    (EffortIndex -> ira -> ira) 
+        {-^ a function behaving well on sequences that intersect to a point -} ->
+    (EffortIndex -> ira -> ira)
+        {- ^ a function behaving well on sequences that intersect to a non-empty interval -}
+maxExtensionInnerR2R getExtremesAndDirections f ix x =
+    case (isIncreasing, isDecreasing, compareReals leftVal rightVal) of
+--        (True, _, Just GT) -> emptyApprox
+--        (True, _, Nothing) -> emptyApprox
+--        (_, True, Just LT) -> emptyApprox
+--        (_, True, Nothing) -> emptyApprox
+        (True, _, _) -> fromBounds (leftVal, rightVal)
+        (_, True, _) -> fromBounds (rightVal, leftVal)
+        _ -> 
+            (/\) ((-1) \/ 1) $ 
+                foldl1 (\/) $ [leftVal, rightVal] ++ extremes
+    where
+    (extremes, (maybeLowIncreasing, maybeHighIncreasing)) =
+        getExtremesAndDirections ix x 
+    (Just lowIsIncreasing) = maybeLowIncreasing 
+    (isIncreasing, isDecreasing)
+        | null extremes = (lowIsIncreasing, not lowIsIncreasing)
+        | otherwise = (False, False)
+    leftVal =
+        case maybeLowIncreasing of
+            Just True -> snd $ bounds $ f ix xL
+            Just False -> fst $ bounds $ f ix xL
+            Nothing -> emptyApprox
+    rightVal =
+        case maybeHighIncreasing of
+            Just True -> fst $ bounds $ f ix xR
+            Just False -> snd $ bounds $ f ix xR
+            Nothing -> emptyApprox
+    (xL, xR) = bounds x
+        
         
         
         
diff --git a/src/Data/Number/ER/Real/Approx/Interval.hs b/src/Data/Number/ER/Real/Approx/Interval.hs
--- a/src/Data/Number/ER/Real/Approx/Interval.hs
+++ b/src/Data/Number/ER/Real/Approx/Interval.hs
@@ -16,7 +16,10 @@
 module Data.Number.ER.Real.Approx.Interval 
 (
     ERInterval(..),
-    normaliseERInterval
+    normaliseERInterval,
+    intervalTimesInner,
+    intervalPlusInner,
+    intervalDivideInner
 )
 where
 
@@ -355,59 +358,118 @@
     negate ERIntervalEmpty = ERIntervalEmpty
     negate ERIntervalAny = ERIntervalAny
     {- addition -}
-    (ERInterval l1 r1) + (ERInterval l2 r2) =
-        normaliseERInterval $
-        ERInterval 
-            (-((-l1) + (-l2))) -- reverse the rounding mode
-            (r1 + r2)
-    ERIntervalAny + i2 = ERIntervalAny
-    i1 + ERIntervalAny = ERIntervalAny
-    ERIntervalEmpty + i2 = ERIntervalEmpty
-    i1 + ERIntervalEmpty = ERIntervalEmpty
+    i1 + i2 = fst $ intervalPlusOuterInner i1 i2
     {- multiplication -}
-    (ERInterval l1 r1) * (ERInterval l2 r2)
-        | haveNan = ERIntervalAny
-        | otherwise =
-            normaliseERInterval $
-            ERInterval minProd maxProd
-        where
-        haveNan = or $ map B.isERNaN (prodsL ++ prodsR)
-        minProd = foldl1 min prodsL
-        maxProd = foldl1 max prodsR
-        prodsL = [-((-l1) * l2), -((-l1) * r2), -((-r1) * l2), -((-r1) * r2)]
-        prodsR = [l1 * l2, l1 * r2, r1 * l2, r1 * r2]
-    ERIntervalAny * i2 = ERIntervalAny
-    i1 * ERIntervalAny = ERIntervalAny
-    ERIntervalEmpty * i2 = ERIntervalEmpty
-    i1 * ERIntervalEmpty = ERIntervalEmpty
+    i1 * i2 = fst $ intervalTimesOuterInner i1 i2 
 
+{-|
+    Add two real approximations, assuming the approximations are `inner'
+    as opposed to `outer':
+    
+    * `outer': the approximation contains all the number(s) of interest
+    * `inner': all numbers eligible for the approximation are numbers of interest
+-}
+intervalPlusInner ::
+    (B.ERRealBase b) =>
+    (ERInterval b) -> 
+    (ERInterval b) -> 
+    (ERInterval b) 
+intervalPlusInner i1 i2 = snd $ intervalPlusOuterInner i1 i2
+
+{-|
+    Multiply two real approximations, assuming the approximations are `inner'
+    as opposed to `outer':
+    
+    * `outer': the approximation contains all the number(s) of interest
+    * `inner': all numbers eligible for the approximation are numbers of interest
+-}
+intervalTimesInner ::
+    (B.ERRealBase b) =>
+    (ERInterval b) -> 
+    (ERInterval b) -> 
+    (ERInterval b) 
+intervalTimesInner i1 i2 = snd $ intervalTimesOuterInner i1 i2
+
+intervalPlusOuterInner (ERInterval l1 r1) (ERInterval l2 r2) =
+        (normaliseERInterval $
+         ERInterval (l1 `plusDown` l2) (r1 `plusUp` r2),
+         ERInterval (l1 `plusUp` l2) (r1 `plusDown` r2))
+intervalPlusOuterInner ERIntervalAny i2 = (ERIntervalAny, ERIntervalAny)
+intervalPlusOuterInner l1 ERIntervalAny = (ERIntervalAny, ERIntervalAny)
+intervalPlusOuterInner ERIntervalEmpty i2 = (ERIntervalEmpty, ERIntervalEmpty)
+intervalPlusOuterInner l1 ERIntervalEmpty = (ERIntervalEmpty, ERIntervalEmpty)
+
+intervalTimesOuterInner (ERInterval l1 r1) (ERInterval l2 r2)
+    | haveNan = (ERIntervalAny, ERIntervalAny)
+    | otherwise =
+        (normaliseERInterval $
+         ERInterval minProdOuter maxProdOuter, 
+         ERInterval minProdInner maxProdInner)
+    where
+    haveNan = or $ map B.isERNaN (prodsUp ++ prodsDown)
+    minProdOuter = foldl1 min prodsDown
+    maxProdOuter = foldl1 max prodsUp
+    minProdInner = foldl1 min prodsUp
+    maxProdInner = foldl1 max prodsDown
+    prodsDown = [l1 `timesDown` l2, l1 `timesDown` r2, r1 `timesDown` l2, r1 `timesDown` r2]
+    prodsUp = [l1 `timesUp` l2, l1 `timesUp` r2, r1 `timesUp` l2, r1 `timesUp` r2]
+intervalTimesOuterInner ERIntervalAny i2 = (ERIntervalAny, ERIntervalAny)
+intervalTimesOuterInner l1 ERIntervalAny = (ERIntervalAny, ERIntervalAny)
+intervalTimesOuterInner ERIntervalEmpty i2 = (ERIntervalEmpty, ERIntervalEmpty)
+intervalTimesOuterInner l1 ERIntervalEmpty = (ERIntervalEmpty, ERIntervalEmpty)
+
 instance (B.ERRealBase b) => Fractional (ERInterval b) where
     fromRational rat =
         (fromInteger $ numerator rat)
         / (fromInteger $ denominator rat)
     {- division -}
-    (ERInterval l1 r1) / (ERInterval l2 r2)
-        | l2 < 0 && r2 > 0 = ERIntervalAny
+    i1 / i2 =
+        fst $ intervalDivideOuterInner i1 i2
+        
+intervalDivideInner ::
+    (B.ERRealBase b) =>
+    (ERInterval b) -> 
+    (ERInterval b) -> 
+    (ERInterval b) 
+intervalDivideInner i1 i2 = snd $ intervalDivideOuterInner i1 i2
+
+intervalDivideOuterInner (ERInterval l1 r1) (ERInterval l2 r2)
+        | l2 < 0 && r2 > 0 = (ERIntervalAny, ERIntervalAny)
         | haveNan = 
 --            unsafePrint "ERInterval: /: haveNan" $ 
-            ERIntervalAny
+            (ERIntervalAny, ERIntervalAny)
         | l2 == 0 && r2 > 0 && 1/l2 < 0 = -- minus 0
-            (ERInterval l1 r1) / (ERInterval (-l2) r2) -- correct it to +0
+            intervalDivideOuterInner (ERInterval l1 r1) (ERInterval (-l2) r2) -- correct it to +0
         | r2 == 0 && l2 < 0 && 1/r2 > 0 = -- plus 0
-            (ERInterval l1 r1) / (ERInterval l2 (-r2)) -- correct it to -0
+            intervalDivideOuterInner (ERInterval l1 r1) (ERInterval l2 (-r2)) -- correct it to -0
         | otherwise =
-            normaliseERInterval $
-            ERInterval minDiv maxDiv
+--            unsafePrintReturn
+--            (
+--                "intervalDivideOuterInner: "
+--                ++ "\n divsUp = " ++ show divsUp
+--                ++ "\n divsDown = " ++ show divsDown
+--                ++ "\n result = "
+--            )
+            (
+             normaliseERInterval $
+             ERInterval minDivOuter maxDivOuter
+            ,
+             ERInterval minDivInner maxDivInner
+            )
         where
-        haveNan = or $ map B.isERNaN (divsL ++ divsR)
-        minDiv = foldl1 min divsL
-        maxDiv = foldl1 max divsR
-        divsL = [-(l1 / (-l2)), -(l1 / (-r2)), -(r1 / (-l2)), -(r1 / (-r2))]
-        divsR = [l1 / l2, l1 / r2, r1 / l2, r1 / r2]
-    ERIntervalAny / i2 = ERIntervalAny
-    i1 / ERIntervalAny = ERIntervalAny
-    ERIntervalEmpty / i2 = ERIntervalEmpty
-    i1 / ERIntervalEmpty = ERIntervalEmpty
+        haveNan = or $ map B.isERNaN (divsUp ++ divsDown)
+        minDivOuter = foldl1 min divsDown
+        maxDivOuter = foldl1 max divsUp
+        minDivInner = foldl1 min divsUp
+        maxDivInner = foldl1 max divsDown
+        divsDown = [l1 `divideDown` l2, l1 `divideDown` r2, r1 `divideDown` l2, r1 `divideDown` r2]
+        divsUp = [l1 `divideUp` l2, l1 `divideUp` r2, r1 `divideUp` l2, r1 `divideUp` r2]
+intervalDivideOuterInner ERIntervalAny i2 = (ERIntervalAny, ERIntervalAny)
+intervalDivideOuterInner i1 ERIntervalAny = (ERIntervalAny, ERIntervalAny)
+intervalDivideOuterInner ERIntervalEmpty i2 = (ERIntervalEmpty, ERIntervalEmpty)
+intervalDivideOuterInner i1 ERIntervalEmpty = (ERIntervalEmpty, ERIntervalEmpty)
+
+     
             
 instance (B.ERRealBase b, RealFrac b) => RA.ERApprox (ERInterval b) where
     initialiseBaseArithmetic _ =
@@ -473,6 +535,8 @@
     {- refines -}
     refines _ ERIntervalAny = True
     refines ERIntervalEmpty _ = True
+    refines ERIntervalAny (ERInterval l r) 
+        | B.isPlusInfinity r && B.isPlusInfinity (-l) = True
     refines ERIntervalAny _ = False
     refines _ ERIntervalEmpty = False
     refines (ERInterval l1 r1) (ERInterval l2 r2) =
@@ -538,6 +602,13 @@
     bounds ERIntervalEmpty = (ERIntervalEmpty, ERIntervalEmpty)
     bounds (ERInterval l r) = 
         (ERInterval l l, ERInterval r r)
+    {- RA.fromBounds -}
+    fromBounds (ERInterval l1 r1, ERInterval l2 r2) 
+        | l1 == r1 && l2 == r2 = ERInterval l1 l2
+    fromBounds i1i2 =
+        error $
+            "ER.Real.Approx.Interval: fromBounds: bounds not thin: "
+            ++ show i1i2
 
 instance (B.ERRealBase b, RealFrac b) => RAEL.ERApproxElementary (ERInterval b)
 -- all operations here have appropriate default implementations
diff --git a/src/Data/Number/ER/Real/Arithmetic/Elementary.hs b/src/Data/Number/ER/Real/Arithmetic/Elementary.hs
--- a/src/Data/Number/ER/Real/Arithmetic/Elementary.hs
+++ b/src/Data/Number/ER/Real/Arithmetic/Elementary.hs
@@ -20,20 +20,27 @@
     erPow_IR,
     erSqrt_R,
     erSqrt_IR,
+    erSqrt_IR_Inner,
     erRoot_R,
     erRoot_IR,
+    erRoot_IR_Inner,
     -- * exponentiation and logarithm 
     erExp_R,
     erExp_IR,
+    erExp_IR_Inner,
     erLog_R,
     erLog_IR,
+    erLog_IR_Inner,
     -- * trigonometrics
     erSine_R,
     erSine_IR,
+    erSine_IR_Inner,
     erCosine_R,
     erCosine_IR,
+    erCosine_IR_Inner,
     erATan_R,
     erATan_IR,
+    erATan_IR_Inner,
     erPi_R
 )
 where
@@ -54,7 +61,14 @@
     (RA.ERIntApprox ira, Ord ira) =>
     EffortIndex -> 
     ira -> ira
-erSqr_IR = erSqr_R
+erSqr_IR =
+    RA.maxExtensionR2R
+        sqrExtrema
+        erSqr_R
+    where
+    sqrExtrema ix x 
+        | 0 `RA.refines` x = [0]
+        | otherwise = [] 
 
 erSqr_R ::
     (RA.ERIntApprox ira, Ord ira) =>
@@ -78,8 +92,17 @@
     EffortIndex -> 
     Integer ->
     ira -> ira
-erPow_IR = erPow_R
+erPow_IR ix n x = 
+    RA.maxExtensionR2R
+        powExtrema
+        (\ ix x -> erPow_R ix n x)
+        ix x
+    where
+    powExtrema ix x 
+        | even n && 0 `RA.refines` x = [0]
+        | otherwise = [] 
 
+
 erPow_R ::
     (RA.ERIntApprox ira, Ord ira) =>
     EffortIndex ->
@@ -105,26 +128,38 @@
 -}
 
 erSqrt_R ::
-    (RA.ERIntApprox ira, Ord ira) => 
+    (RA.ERIntApprox ira) => 
     EffortIndex -> ira -> ira
 erSqrt_R = erSqrtNewton_R  
     
 erSqrt_IR ::
-    (RA.ERIntApprox ira, Ord ira) => 
+    (RA.ERIntApprox ira) => 
     EffortIndex -> ira -> ira
 erSqrt_IR =
     RA.maxExtensionR2R 
         sqrtExtrema
         (\ ix x -> erSqrt_R ix x)
 
-sqrtExtrema ix x
-    | 0 `RA.refines` x = [0]
-    | otherwise = []
+erSqrt_IR_Inner ::
+    (RA.ERIntApprox ira) => 
+    EffortIndex -> ira -> ira
+erSqrt_IR_Inner =
+    RA.maxExtensionInnerR2R 
+        sqrtExtremaAndDirections
+        (\ ix x -> erSqrt_R ix x)
+
+sqrtExtrema ix x = fst $ sqrtExtremaAndDirections ix x
         
+sqrtExtremaAndDirections ix x =
+    case RA.compareReals 0 x of
+        Just LT -> ([], (Just True, Just True))
+        Just GT -> ([], (Nothing, Nothing))
+        _ -> ([0], (Nothing, Just True))
+        
 
 
 erSqrtContFr_R ::
-    (RA.ERIntApprox ira, Ord ira) => 
+    (RA.ERIntApprox ira) => 
     EffortIndex -> ira -> ira
 erSqrtContFr_R ix a
     | aR == 0 = 0
@@ -132,7 +167,7 @@
     | aR `RA.ltSingletons` 0 = RA.emptyApprox
     | otherwise =
         contFrIter (ix + 3) $
-            RA.setMinGranularity gran $ max 0 (0 RA.\/ a) 
+            RA.setMinGranularity gran $ 0 RA.\/ aR -- assuming aR >= 0 
     where
     gran = effIx2gran ix
     (aL, aR) = RA.bounds a
@@ -147,14 +182,14 @@
         x_iPlus1 = contFrIter (i - 1) x_i
             
 erSqrtNewton_R ::
-    (RA.ERIntApprox ira, Ord ira) => 
+    (RA.ERIntApprox ira) => 
     EffortIndex -> ira -> ira
 erSqrtNewton_R ix a
     | RA.isEmpty a = RA.emptyApprox
     | aR == 0 = 0
     | aL == 1/0 = 1/0
-    | aR < 0 = RA.emptyApprox
-    | otherwise =        
+    | aR `RA.ltSingletons` 0 = RA.emptyApprox
+    | otherwise =
         x_i RA.\/ (a/x_i)
     where
     gran = effIx2gran ix
@@ -162,8 +197,8 @@
     aM1 = a - 1
     
     x_i = 
-        newtonIter ((ix `div` 100) + 5) $
-                RA.setMinGranularity gran $ max 0 aR 
+        newtonIter ((ix `div` 10) + 5) $
+                RA.setMinGranularity gran aR -- assuming aR >= 0 
     newtonIter i x_i
         | i == 0 = x_i
         | otherwise =
@@ -190,10 +225,24 @@
         (\ ix x -> erRoot_R ix p x) $
             ix
 
-rootExtrema p ix x
-    | 0 `RA.refines` x && even p = [0]
-    | otherwise = []
+erRoot_IR_Inner ::
+    (RA.ERIntApprox ira, Ord ira) => 
+    EffortIndex -> Integer -> ira -> ira
+erRoot_IR_Inner ix p =
+    RA.maxExtensionInnerR2R 
+        (rootExtremaAndDirections p)
+        (\ ix x -> erRoot_R ix p x) $
+            ix
+rootExtrema p ix x = fst $ rootExtremaAndDirections p ix x
 
+rootExtremaAndDirections p ix x
+    | odd p = ([], (Just True, Just True))
+    | otherwise =
+        case RA.compareReals 0 x of
+            Just LT -> ([], (Just True, Just True))
+            Just GT -> ([], (Nothing, Nothing))
+            _ -> ([0], (Nothing, Just True))
+
 erRootNewton_R ::
     (RA.ERIntApprox ira, Ord ira) => 
     EffortIndex -> Integer -> ira -> ira
@@ -240,34 +289,48 @@
 -}
 
 erExp_R :: 
-    (RA.ERIntApprox ira) => 
+    (RA.ERIntApprox ira, Ord ira) => 
     EffortIndex -> ira -> ira
     
-erExp_R = erExp_Tay_Opt_R
+erExp_R ix x 
+    | RA.isBounded x =
+        erPow_IR ix n $ 
+        erExp_Tay_Opt_R ix xNear0
+    | x `RA.refines` (-1/0) = 0
+    | (-1/0) `RA.refines` x =
+        0 RA.\/ (erExp_R ix (snd $ RA.bounds x))
+    | otherwise = RA.bottomApprox
+    where
+    (xNear0, n) = scaleNear0 (x,1)
+    scaleNear0 (xPrev, nPrev) =
+        case xPrev `RA.refines` ((-1) RA.\/ 1) of
+            True -> (xPrev, nPrev)
+            False -> scaleNear0 (xNext, nNext)
+        where
+        xNext = xPrev / 2
+        nNext = 2 * nPrev
 
-{- 
-    exp as derived from Taylor series is already a maximal extension
-    because it does not suffer from the wrapping effect - all
-    functions used are monotone - all Taylor coeffs are positive
--}
 erExp_IR :: 
-    (RA.ERIntApprox ira) => 
+    (RA.ERIntApprox ira, Ord ira) => 
     EffortIndex -> ira -> ira
     
-erExp_IR ix x
-    | 0 `RA.refines` x || (-1/0) `RA.refines` x=
-        RA.maxExtensionR2R
-            (\ ix x -> [])
-            (\ ix x -> erExp_R ix x)
-            ix x
-    | otherwise =
-        erExp_R ix x
+erExp_IR =
+    RA.maxExtensionR2R
+        (\ ix x -> [])
+        erExp_R
 
+erExp_IR_Inner :: 
+    (RA.ERIntApprox ira, Ord ira) => 
+    EffortIndex -> ira -> ira
+erExp_IR_Inner =
+    RA.maxExtensionInnerR2R
+        (\ ix x -> ([], (Just True, Just True)))
+        erExp_R
 
 {- Log using Newton -}
 
 erLog_R :: 
-    (RA.ERIntApprox ira) => 
+    (RA.ERIntApprox ira, Ord ira) => 
     EffortIndex -> ira -> ira
     
 erLog_R =
@@ -275,7 +338,7 @@
 --    erLog_IR -- intervals are more efficient for log than singletons 
 
 erLog_IR ::
-    (RA.ERIntApprox ira) => 
+    (RA.ERIntApprox ira, Ord ira) => 
     EffortIndex -> ira -> ira
     
 erLog_IR =
@@ -283,13 +346,26 @@
         logExtrema
         (\ ix x -> logDivSeries_R ix x)
         
-logExtrema ix x
-    | 0 `RA.refines` x = [-1/0]
-    | otherwise = []
+erLog_IR_Inner ::
+    (RA.ERIntApprox ira, Ord ira) => 
+    EffortIndex -> ira -> ira
+    
+erLog_IR_Inner =
+    RA.maxExtensionInnerR2R
+        logExtremaAndDirections
+        (\ ix x -> logDivSeries_R ix x)
         
+logExtrema ix x = fst $ logExtremaAndDirections ix x
+        
+logExtremaAndDirections ix x =
+    case RA.compareReals 0 x of
+        Just LT -> ([], (Just True, Just True))
+        Just GT -> ([], (Nothing, Nothing))
+        _ -> ([-1/0], (Nothing, Just True))
+        
 {-| log using a fast converging series, designed to be used with singletons -}
 logDivSeries_R ::
-    (RA.ERIntApprox ira) => EffortIndex -> ira -> ira 
+    (RA.ERIntApprox ira, Ord ira) => EffortIndex -> ira -> ira 
 logDivSeries_R ix x 
     | RA.isEmpty posx = RA.emptyApprox
     | posx `RA.refines` 0 = -1/0 
@@ -312,13 +388,13 @@
     t = 
         ((remNearLogx - 1) / (remNearLogx + 1)) -- the range of this expression is [-1,1] 
             RA./\ ((-1) RA.\/ 1) -- correction of wrapping 
-    tsqare = abs $ t * t -- the range is [0,1]
+    tsquare = abs $ t * t -- the range is [0,1]
     series termsCount currentDenominator 
         | termsCount > 0 =
-            (recip currentDenominator) + tsqare * (series (termsCount - 1) (currentDenominator + 2))
+            (recip currentDenominator) + tsquare * (series (termsCount - 1) (currentDenominator + 2))
         | otherwise =
             (recip currentDenominator)
-            * (1 RA.\/ (recip $ 1 - tsqare)) -- [1,1/(1-t^2)] is a valid error bound
+            * (1 RA.\/ (recip $ 1 - tsquare)) -- [1,1/(1-t^2)] is a valid error bound
         
 --{- log using Newton -}
 --    
@@ -355,27 +431,58 @@
     sin(x) and cos(x)
 -}
 
-erSine_R :: 
+erSine_R ::
     (RA.ERIntApprox ira) => 
     EffortIndex -> ira -> ira
 
-erSine_R ix x
-    | (1/0) `RA.refines` x || (-1/0) `RA.refines` x = 
-        (-1) RA.\/ 1  
-    | otherwise =
-        erSine_Tay_Opt_R ix x
+erSine_R ix x =
+    case (RA.isBounded x) of
+        True | xNear0 `RA.refines` plusMinusPiHalf ->
+            erSine_Tay_Opt_R ix xNear0
+        True | (xNear0 - piHalf) `RA.refines` plusMinusPiHalf ->
+            erCosine_Tay_Opt_R ix (xNear0 - piHalf)
+        True | (xNear0 + piHalf) `RA.refines` plusMinusPiHalf ->
+            negate $ erCosine_Tay_Opt_R ix (xNear0 + piHalf)
+        True | (xNear0 - pi) `RA.refines` plusMinusPiHalf ->
+            negate $ erSine_Tay_Opt_R ix (xNear0 - pi)
+        _ ->
+            (-1) RA.\/ 1
+    where
+    xNear0 = x - k * 2 * pi -- should be in [-pi,3*pi/2]
+    k = fromInteger $ toInteger kEI
+    (kEI,_) = RA.integerBounds $ 0.5 + (x / (2*pi))
 
+    plusMinusPiHalf = (- piHalf) RA.\/ piHalf
+    piHalf = pi / 2
+    pi = erPi_R ix
+    
+
 erCosine_R :: 
     (RA.ERIntApprox ira) => 
     EffortIndex -> ira -> ira
      
-erCosine_R ix x
-    | (1/0) `RA.refines` x || (-1/0) `RA.refines` x =   
-        (-1) RA.\/ 1  
-    | otherwise =
-        erCosine_Tay_Opt_R ix x 
+erCosine_R ix x =
+    case (RA.isBounded x) of
+        True | xNear0 `RA.refines` plusMinusPiHalf ->
+            erCosine_Tay_Opt_R ix xNear0
+        True | (xNear0 - piHalf) `RA.refines` plusMinusPiHalf ->
+            negate $ erSine_Tay_Opt_R ix (xNear0 - piHalf)
+        True | (xNear0 + piHalf) `RA.refines` plusMinusPiHalf ->
+            erSine_Tay_Opt_R ix (xNear0 + piHalf)
+        True | (xNear0 - pi) `RA.refines` plusMinusPiHalf ->
+            negate $ erCosine_Tay_Opt_R ix (xNear0 - pi)
+        _ ->
+            (-1) RA.\/ 1
+    where
+    xNear0 = x - k * 2 * pi -- should be in [-pi,3*pi/2]
+    k = fromInteger $ toInteger kEI
+    (kEI,_) = RA.integerBounds $ 0.5 + (x / (2*pi))
 
+    plusMinusPiHalf = (- piHalf) RA.\/ piHalf
+    piHalf = pi / 2
+    pi = erPi_R ix
 
+
 {- Sine using generic Taylor (see Taylor for an optimised version) -}
 
 erSine_Tay_R :: 
@@ -416,27 +523,49 @@
 erCosine_IR = 
     RA.maxExtensionR2R cosineExtremes erCosine_R
         
-sineExtremes ix x 
+erSine_IR_Inner ::
+    (RA.ERIntApprox ira) =>
+    EffortIndex -> ira -> ira 
+    
+erSine_IR_Inner = 
+    RA.maxExtensionInnerR2R sineExtremesAndDirections erSine_R
+    
+erCosine_IR_Inner ::
+    (RA.ERIntApprox ira) =>
+    EffortIndex -> ira -> ira 
+    
+erCosine_IR_Inner = 
+    RA.maxExtensionInnerR2R cosineExtremesAndDirections erCosine_R
+        
+sineExtremes ix x = fst $ sineExtremesAndDirections ix x
+cosineExtremes ix x = fst $ cosineExtremesAndDirections ix x
+        
+sineExtremesAndDirections ix x 
     | RA.isBounded x =
         alternatingExtremes 1 (-1) ix scaledX
-    | otherwise = [-1,1]
+    | otherwise = ([-1,1], (Nothing, Nothing))
     where
     scaledX = (x / (erPi_R ix)) - 0.5
     
-cosineExtremes ix x
+cosineExtremesAndDirections ix x
     | RA.isBounded x =
         alternatingExtremes 1 (-1) ix scaledX
-    | otherwise = [-1,1]
+    | otherwise = ([-1,1], (Nothing, Nothing))
     where
     scaledX = (x / (erPi_R ix))
     
-alternatingExtremes extr0 extr1 ix scaledX
-    | extremesCount >= 2 = [extr0,extr1]  
-    | extremesCount == 1 && even minExtremeN = [extr0]
-    | extremesCount == 1 = [extr1]
-    | otherwise = []
+alternatingExtremes extrHigh extrLow ix scaledX
+    | extremesCount == 1 && even minExtremeN = 
+        ([extrHigh], (Just True, Just False)) -- increasing, decreasing
+    | extremesCount == 1 =
+        ([extrLow], (Just False, Just True)) -- decreasing, increasing
+    | extremesCount >= 2 = 
+        ([extrHigh,extrLow], (Just $ even minExtremeN, Just $ odd maxExtremeN))  
+    | otherwise = 
+        ([], (Just isIncreasing, Just isIncreasing))
     where
     extremesCount = 1 + maxExtremeN - minExtremeN
+    isIncreasing = even maxExtremeN
     (xFloor, xCeiling) = RA.integerBounds scaledX
     minExtremeN = 
         case RA.compareReals (fromInteger $ toInteger xFloor) scaledX of
@@ -467,13 +596,22 @@
 
 atanExtremes ix x = []
 
+erATan_IR_Inner ::
+    (RA.ERIntApprox ira) =>
+    EffortIndex -> ira -> ira 
+
+erATan_IR_Inner =
+    RA.maxExtensionInnerR2R atanExtremesAndDirections erATan_R
+
+atanExtremesAndDirections ix x = ([], (Just True, Just True))
+
 {- 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 + ...)))))
+    (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]
 -}
 
 atanEuler_R ::
@@ -482,15 +620,17 @@
 
 atanEuler_R ix x 
     | RA.isEmpty x = RA.emptyApprox
-    | otherwise =
+    | x `RA.refines` ((-1.5) RA.\/ 1.5) =
         (x / xSquarePlus1) * (series ix (RA.setMinGranularity gran 2))
+    | otherwise = -- too far from 0, needs atan(x) = 2*atan(x/(1+sqrt(1+x^2)))
+        2 * (atanEuler_R ix $ x / (1 + (erSqrt_R ix $ 1 + x * x)))
     where
     gran = effIx2gran ix
     series termsCount coeffBase 
         | termsCount > 0 =
             1 + xSquareOverXSquarePlus1 * coeff * (series (termsCount - 1) (coeffBase + 2))
         | otherwise =
-            1 + xSquare * (coeff RA.\/ 1)
+            1 + xSquare * (0 RA.\/ 1)
         where
         coeff = coeffBase / (coeffBase + 1)
     xSquare = abs $ x * x
diff --git a/src/Data/Number/ER/Real/Arithmetic/Integration.hs b/src/Data/Number/ER/Real/Arithmetic/Integration.hs
--- a/src/Data/Number/ER/Real/Arithmetic/Integration.hs
+++ b/src/Data/Number/ER/Real/Arithmetic/Integration.hs
@@ -26,7 +26,7 @@
 import Data.Number.ER.Real.Arithmetic.Elementary
 
 testIntegr1 :: 
-    (RA.ERIntApprox ira) => 
+    (RA.ERIntApprox ira, Ord ira) => 
     (ConvergRealSeq ira)
 testIntegr1 = integrateCont erExp_IR 0 1
 
diff --git a/src/Data/Number/ER/Real/Arithmetic/Taylor.hs b/src/Data/Number/ER/Real/Arithmetic/Taylor.hs
--- a/src/Data/Number/ER/Real/Arithmetic/Taylor.hs
+++ b/src/Data/Number/ER/Real/Arithmetic/Taylor.hs
@@ -108,14 +108,16 @@
     => EffortIndex
     -> ira
     -> ira
-erSine_Tay_Opt_R ix x = taylor_seg ix x (RA.setMinGranularity gran 1)
-	where
-		gran = effIx2gran ix
-		taylor_seg i x n -- 'i' for iterator
-			| i > 0  = x - ((x*x)/((n+1)*(n+2))) * (taylor_seg (i-2) x (n+2))
-			| otherwise = errorRegion
-				where 
-					errorRegion = (-1) RA.\/ (1)
+erSine_Tay_Opt_R ix x = 
+    taylor_seg ix x (RA.setMinGranularity gran 1)
+    where
+    gran = effIx2gran ix
+    taylor_seg i x n -- 'i' for iterator
+        | i > 0  = x - (x*x)/((n+1)*(n+2)) * (taylor_seg (i-2) x (n+2))
+        | otherwise = errorRegion
+        where 
+        errorRegion = (- x) RA.\/ x
+                    
 		
 {-|
     A Taylor series for cosine.    
@@ -127,12 +129,12 @@
     -> ira
 erCosine_Tay_Opt_R ix x = taylor_seg ix x (RA.setMinGranularity gran 1)
 	where
-		gran = effIx2gran ix
-		taylor_seg i x n -- 'i' for iterator
-			| i > 0  = 1 - ((x*x)/(n*(n+1))) * (taylor_seg (i-2) x (n+2))
-			| otherwise = errorRegion
-				where 
-					errorRegion = (-1) RA.\/ (1)
+    gran = effIx2gran ix
+    taylor_seg i x n -- 'i' for iterator
+        | i > 0  = 1 - ((x*x)/(n*(n+1))) * (taylor_seg (i-2) x (n+2))
+        | otherwise = errorRegion
+        where 
+        errorRegion = (-1) RA.\/ (1)
 
    
 				
diff --git a/src/Data/Number/ER/Real/DomainBox.hs b/src/Data/Number/ER/Real/DomainBox.hs
--- a/src/Data/Number/ER/Real/DomainBox.hs
+++ b/src/Data/Number/ER/Real/DomainBox.hs
@@ -59,6 +59,7 @@
     where
     noinfo :: box
     isNoinfo :: box -> Bool
+    size :: box -> Int
     {-| constructor using 'defaultVar' -}
     unary :: val -> box
     singleton :: varid -> val -> box
diff --git a/src/Data/Number/ER/Real/DomainBox/IntMap.hs b/src/Data/Number/ER/Real/DomainBox/IntMap.hs
--- a/src/Data/Number/ER/Real/DomainBox/IntMap.hs
+++ b/src/Data/Number/ER/Real/DomainBox/IntMap.hs
@@ -45,6 +45,7 @@
     where
     noinfo = IMap.empty
     isNoinfo = IMap.null
+    size = IMap.size
     unary r = IMap.singleton defaultVar r
     singleton = IMap.singleton
     toList = IMap.toList