AERN-Real 0.9.8 → 0.9.9

raw patch · 10 files changed

+421/−151 lines, 10 filesdep +timedep ~basePVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: time

Dependency ranges changed: base

API changes (from Hackage documentation)

Files

• AERN-Real.cabal +5/−17
• ChangeLog +9/−0
• src/Data/Number/ER/Misc.hs +18/−0
• src/Data/Number/ER/Real/Approx.hs +52/−12
• src/Data/Number/ER/Real/Approx/Interval.hs +112/−41
• src/Data/Number/ER/Real/Arithmetic/Elementary.hs +206/−66
• src/Data/Number/ER/Real/Arithmetic/Integration.hs +1/−1
• src/Data/Number/ER/Real/Arithmetic/Taylor.hs +16/−14
• src/Data/Number/ER/Real/DomainBox.hs +1/−0
• src/Data/Number/ER/Real/DomainBox/IntMap.hs +1/−0

@@ -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,

@@ -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

@@ -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

@@ -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+                                   

@@ -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

@@ -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

@@ -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 

@@ -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)      				

@@ -59,6 +59,7 @@     where     noinfo :: box     isNoinfo :: box -> Bool+    size :: box -> Int     {-| constructor using 'defaultVar' -}     unary :: val -> box     singleton :: varid -> val -> box

@@ -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