packages feed

AERN-RnToRm 0.3.0.3 → 0.4

raw patch · 12 files changed

+517/−204 lines, 12 files

Files

AERN-RnToRm.cabal view
@@ -1,5 +1,5 @@ Name:           AERN-RnToRm-Version:        0.3.0.3+Version:        0.4 Cabal-Version:  >= 1.2 Build-Type:     Simple License:        BSD3
ChangeLog view
@@ -1,3 +1,6 @@+0.4: 20 August 2008+    * fixed several serious bugs in sin and cos+    * added arctan 0.3.0.3: 7 August 2008     * revamped package description to make it much shorter and linked it       to the main module
src/Data/Number/ER/RnToRm/Approx/DomEdges.hs view
@@ -251,6 +251,7 @@     log ix = edgesLift1 $ RAEL.log ix     sin ix = edgesLift1 $ RAEL.sin ix     cos ix = edgesLift1 $ RAEL.cos ix+    atan ix = edgesLift1 $ RAEL.atan ix          instance      (FA.ERFnDomApprox box varid domra ranra fa, VariableID varid) =>
src/Data/Number/ER/RnToRm/Approx/DomTransl.hs view
@@ -323,6 +323,8 @@         ERFnDomTranslApprox (RAEL.sin ix ufa) dtrB     cos ix (ERFnDomTranslApprox ufa dtrB) =         ERFnDomTranslApprox (RAEL.cos ix ufa) dtrB+    atan ix (ERFnDomTranslApprox ufa dtrB) =+        ERFnDomTranslApprox (RAEL.atan ix ufa) dtrB  instance      (UFA.ERUnitFnApprox box varid domra ranra ufa,
src/Data/Number/ER/RnToRm/Approx/PieceWise.hs view
@@ -232,6 +232,7 @@     log ix = pwLift1 $ RAEL.log ix     sin ix = pwLift1 $ RAEL.sin ix     cos ix = pwLift1 $ RAEL.cos ix+    atan ix = pwLift1 $ RAEL.atan ix      instance      (FA.ERFnDomApprox box varid domra ranra fa, @@ -343,12 +344,12 @@                 integrateOriginHigher                 [bistrD,  bistrInit]         zipOutsideRange maybeFromL maybeFromR [bistrD, bistrInit] =-            unsafePrint-            (-                "ERFnPiecewise: integrateMeasureImprovement: zipOutsideRange: "-                ++ "\n domB = " ++ show domB-                ++ "\n bottomFn = " ++ show bottomFn-            )+--            unsafePrint+--            (+--                "ERFnPiecewise: integrateMeasureImprovement: zipOutsideRange: "+--                ++ "\n domB = " ++ show domB+--                ++ "\n bottomFn = " ++ show bottomFn+--            )             [bistrPadj]             where             (ERFnPiecewise bistrPadj) =
src/Data/Number/ER/RnToRm/Approx/Tuple.hs view
@@ -203,6 +203,7 @@     log ix = tuplesLift1 $ RAEL.log ix     sin ix = tuplesLift1 $ RAEL.sin ix     cos ix = tuplesLift1 $ RAEL.cos ix+    atan ix = tuplesLift1 $ RAEL.atan ix          instance      (FA.ERFnDomApprox box varid domra ranra fa) =>
src/Data/Number/ER/RnToRm/TestingDefs.hs view
@@ -13,6 +13,7 @@ module Data.Number.ER.RnToRm.TestingDefs where  import Data.Number.ER.RnToRm.DefaultRepr+import Data.Number.ER.Real.DefaultRepr  import qualified Data.Number.ER.RnToRm.Approx as FA import qualified Data.Number.ER.RnToRm.UnitDom.Approx as UFA@@ -50,7 +51,7 @@ fapeUConst13 = (FA.const (DBox.unary $ (0)RA.\/1) [1 RA.\/ 3]) :: FAPE fapeUConst13InitPt = FA.partialIntersect 1 (DBox.unary 0) fapeUConst13 fapeUConst1  -fapwUUX0 = (FA.proj (DBox.fromAscList [(0,(1) RA.\/ 1)]) 0) :: FAPWP+fapwUUX0 = (FA.proj (DBox.fromAscList [(0,(-1) RA.\/ 1)]) 0) :: FAPWP fapwUUX1 = (FA.proj (DBox.fromAscList [(1,(-1) RA.\/ 1)]) 1) :: FAPWP  fapwUX0 = (FA.proj (DBox.fromAscList [(0,(0) RA.\/ 1)]) 0) :: FAPWP@@ -66,7 +67,14 @@ testIntegrP =      FA.integrateMeasureImprovement 1 (FA.setMaxDegree 0 fapwUConst13InitPt) 0 (DBox.unary $ 0 RA.\/ 0.5) 0 fapwUConst13InitPt -x = FA.setMaxDegree 4 fapwUX0+x = +--    FA.bisectUnbisectDepth 1 $+    FA.setMaxDegree 4 +    fapwUUX0+    +xLR = +    snd $ FA.bisect 0 Nothing $ fst $ FA.bisect 0 Nothing $ x+     fn1 = (1 + x) RA.\/ (1 + 3*x) fn2 = FA.integrateUnary 0 fn1 0 (0 RA.\/ 1) [1] fn3 = FA.integrateUnary 0 fn2 0 (0 RA.\/ 1) [1] -- this seems wrong!
src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs view
@@ -124,9 +124,9 @@     where     (ERFnInterval h1 ln1 ctxt1 gl1)              == (ERFnInterval h2 ln2 ctxt2 gl2) =-        error "ERFnInterval: equality not implemented yet"+        error "ERFnInterval: equality not implemented"     _ == _ =-        error "ERFnInterval: equality not implemented yet"+        error "ERFnInterval: equality not implemented"  instance      (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>@@ -135,10 +135,11 @@     compare              (ERFnInterval h1 ln1 ctxt1 gl1)              (ERFnInterval h2 ln2 ctxt2 gl2) =-        error "ERFnInterval: comparison not implemented yet"+        error "ERFnInterval: comparison not implemented; consider leqReals from class ERApprox instead"     compare _ _ =-        error "ERFnInterval: comparison not implemented yet"+        error "ERFnInterval: comparison not implemented; consider leqReals from class ERApprox instead"     +     instance      (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>     Num (ERFnInterval fb ra)@@ -312,7 +313,7 @@         ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2)  instance-    (UFB.ERUnitFnBase boxb boxra varid b ra fb, RAEL.ERApproxElementary ra) =>+    (UFB.ERUnitFnBase boxb boxra varid b ra fb, RAEL.ERApproxElementary ra, RealFrac b) =>     RAEL.ERApproxElementary (ERFnInterval fb ra)      where     -- default abs does not work because we do not have Prelude.abs@@ -345,31 +346,173 @@ --            ++ "\n uSin = " ++ show uSin --            ++ "\n lSinNeg = " ++ show lSinNeg --        ) $-        ERFnInterval uSin lSinNeg c (RAEL.sin ix g)+        ERFnInterval uSin (- lSin) c (RAEL.sin ix g)         where+        (lSin, uSin) = sincos True maxDegree ix u (-ln)           maxDegree = erfnMaxDegree c---        ix = int2effIx maxDegree-        uSin = snd $ UFB.sin maxDegree ix u-        lSinNeg = -            negate $ fst $ UFB.sin maxDegree ix (negate ln)      cos ix f@(ERFnIntervalAny c) =         ERFnInterval 1 1 c ((-1) RA.\/ 1)     cos ix (ERFnInterval u ln c g) = --        unsafePrint --        (---            "ERFnInterval: RAEL.sin: "+--            "ERFnInterval: RAEL.cos: " --            ++ "\n u = " ++ show u --            ++ "\n ln = " ++ show ln---            ++ "\n uSin = " ++ show uSin---            ++ "\n lSinNeg = " ++ show lSinNeg+--            ++ "\n uCos = " ++ show uCos+--            ++ "\n lCosNeg = " ++ show lCosNeg --        ) $-        ERFnInterval uCos lCosNeg c (RAEL.cos ix g)+        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) =+--        unsafePrint+--        (+--            "ERFnInterval: RAEL.atan: "+--            ++ "\n u = " ++ show u+--            ++ "\n ln = " ++ show ln+--            ++ "\n uAtan = " ++ show uAtan+--            ++ "\n lAtanNeg = " ++ show lAtanNeg+--        ) $+        ERFnInterval uAtan lAtanNeg c (RAEL.atan ix g)+        where+        maxDegree = erfnMaxDegree c --        ix = int2effIx maxDegree-        uCos = snd $ UFB.cos maxDegree ix u-        lCosNeg = -            negate $ fst $ UFB.cos maxDegree ix (negate ln) +        uAtan = snd $ UFB.atan maxDegree ix u+        lAtanNeg = +            negate $ fst $ UFB.atan maxDegree ix (negate ln) ++sincos ::+    (UFB.ERUnitFnBase boxb boxra varid b ra fb, RAEL.ERApproxElementary ra, RealFrac b) =>+    Bool {-^ True iff sine, False iff cosine -} -> +    Int {-^ maximum representation degree -} -> +    EffortIndex {-^ how hard to try to eliminate truncation errors -} -> +    fb ->+    fb ->+    (fb, fb)+sincos isSine maxDegree ix u l+    -- p - 2k*pi range within [-pi/2, pi/2]: +    | ranfNear0 `RA.refines` plusMinusPiHalf =+--        unsafePrint+--        (+--            "ERFnInterval: sincos: [-pi/2, pi/2]: "+--            ++ "\n u = " ++ show u+--            ++ "\n l = " ++ show l+--            ++ "\n ranf = " ++ show ranf+--            ++ "\n k = " ++ show k+--            ++ "\n ranfNear0 = " ++ show ranfNear0+--        ) $+        case isSine of+            True -> sineShifted (- k2pi)+            False -> cosineShifted (- k2pi)+    -- p - 2k*pi range within [0, pi]: +    | (ranfNear0 - piHalf) `RA.refines` plusMinusPiHalf =+--        unsafePrint+--        (+--            "ERFnInterval: sincos: [0, pi]: "+--            ++ "\n u = " ++ show u+--            ++ "\n l = " ++ show l+--            ++ "\n ranf = " ++ show ranf+--            ++ "\n k = " ++ show k+--            ++ "\n ranfNear0 = " ++ show ranfNear0+--        ) $+        case isSine of+            -- use sin(x) = cos(x - pi/2) and cos(x) = - sin(x - pi/2):+            True -> cosineShifted (- k2pi - piHalf)+            False -> sineShiftedNegated (- k2pi - piHalf)+    -- p - 2k*pi range within [-pi, 0]: +    | (ranfNear0 + piHalf) `RA.refines` plusMinusPiHalf =+--        unsafePrint+--        (+--            "ERFnInterval: sincos: [-pi, 0]: "+--            ++ "\n u = " ++ show u+--            ++ "\n l = " ++ show l+--            ++ "\n ranf = " ++ show ranf+--            ++ "\n k = " ++ show k+--            ++ "\n ranfNear0 = " ++ show ranfNear0+--        ) $+        case isSine of+            -- use sin(x) = - cos(x + pi/2) and cos(x) = sin(x + pi/2):+            True -> cosineShiftedNegated (-k2pi + piHalf)+            False -> sineShifted (-k2pi + piHalf)+    -- p - 2k*pi range within [pi/2, 3pi/2]: +    | (ranfNear0 - pi) `RA.refines` plusMinusPiHalf =+--        unsafePrint+--        (+--            "ERFnInterval: sincos: [pi/2, 3pi/2]: "+--            ++ "\n u = " ++ show u+--            ++ "\n l = " ++ show l+--            ++ "\n ranf = " ++ show ranf+--            ++ "\n k = " ++ show k+--            ++ "\n ranfNear0 = " ++ show ranfNear0+--        ) $+        -- use sin(x) = - sin(x - pi) and cos(x) = - cos(x - pi)+        case isSine of+            True -> sineShiftedNegated (- k2pi - pi)+            False -> cosineShiftedNegated (- k2pi - pi)+    | otherwise = +--        unsafePrint+--        (+--            "ERFnInterval: sincos: big range: "+--            ++ "\n u = " ++ show u+--            ++ "\n l = " ++ show l+--            ++ "\n ranf = " ++ show ranf+--            ++ "\n k = " ++ show k+--            ++ "\n ranfNear0 = " ++ show ranfNear0+--        ) $+        (UFB.const (-1), UFB.const 1)+--    (expDownwards, expUpwards + valueAtRDnNeg + (UFB.const expRUp))+    where+    ranfNear0 = ranf - k2pi+    k2pi = k * 2 * pi+    plusMinusPiHalf = (-piHalfLO) RA.\/ piHalfLO+    pi = RAEL.pi ix  +    piHalf = pi / 2+    (piHalfLO, piHalfHI) = RA.bounds piHalf+    ranf = +        ERInterval +            (UFB.lowerBound 10 l) +            (UFB.upperBound 10 u)+    k = +        fromInteger $ floor $ +            case (pi + ranf) / (2 * pi) of ERInterval lo hi -> lo++    sineShiftedNegated shift =+        boundsNegate $ sineShifted shift+        +    cosineShiftedNegated shift =+        boundsNegate $ cosineShifted shift++    boundsNegate (pLO, pHI) = (- pHI, - pLO)+        +    sineShifted shift =+        boundsAddErr shiftWidthB (lSinDown, uSinUp)+        where+        lSinDown = fst $ UFB.sin maxDegree ix (l `plusUp` shiftPoly)+        uSinUp = snd $ UFB.sin maxDegree ix (u `plusDown` shiftPoly)  +        shiftPoly = UFB.const shiftLOB+        ERInterval shiftLOB shiftHIB = shift+        shiftWidthB = shiftHIB - shiftLOB+    +    cosineShifted shift =+        boundsAddErr shiftWidthB $ +            (UFB.minDown maxDegree lCosDown uCosDown,+             UFB.maxUp maxDegree lCosUp uCosUp +                + (snd $ UFB.scale 0.5 (u-l))) -- important near 0+        where+        (lCosDown, lCosUp) = UFB.cos maxDegree ix (l `plusUp` shiftPoly)+        (uCosDown, uCosUp) = UFB.cos maxDegree ix (u `plusDown` shiftPoly)  +        shiftPoly = UFB.const shiftLOB+        ERInterval shiftLOB shiftHIB = shift+        shiftWidthB = shiftHIB - shiftLOB+    +    boundsAddErr errB (pLO, pHI) =+        (pLO `plusDown` (- errPoly), pHI + errPoly)+        where+        errPoly = UFB.const errB+      instance      (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>
src/Data/Number/ER/RnToRm/UnitDom/Base.hs view
@@ -237,7 +237,8 @@         ufb {-^ p(x) -} ->          (ufb, ufb)     {-| -        Approximate @sin(p(x))@ from below and from above.+        Approximate @sin(p(x))@ from below and from above,+        assuming the range of p is within [-pi/2,pi/2].     -}     sin ::          Int {-^ max degree for result -} -> @@ -245,9 +246,18 @@         ufb {-^ p(x) -} ->          (ufb, ufb)     {-|-        Approximate @cos(p(x))@ from below and from above.+        Approximate @cos(p(x))@ from below and from above,+        assuming the range of p is within [-pi/2,pi/2].     -}     cos :: +        Int {-^ max degree for result -} -> +        EffortIndex {-^ how hard to try when approximating cos as a polynomial -} -> +        ufb {-^ p(x) -} -> +        (ufb, ufb)+    {-|+        Approximate @atan(p(x))@ from below and from above.+    -}+    atan ::          Int {-^ max degree for result -} ->          EffortIndex {-^ how hard to try when approximating cos as a polynomial -} ->          ufb {-^ p(x) -} ->
src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs view
@@ -23,7 +23,7 @@ module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom (     ERChebPoly(..), TermKey-) +) where  import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic@@ -80,8 +80,9 @@     sqrt = chplSqrt     exp = chplExp     log = chplLog-    sin = chplSineCosine True-    cos = chplSineCosine False+    sin = chplSine+    cos = chplCosine+    atan = chplAtan     eval = chplEval     evalApprox ufb x = chplEvalApprox (\ b -> ERInterval b b) ufb x     partialEvalApprox substitutions ufb =
src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Bounds.hs view
@@ -187,6 +187,16 @@                     _ ->                         [(term, coeff)] +chplMaxDn m a b = fst $ chplMax m a b+chplMaxUp m a b = snd $ chplMax m a b+chplMinDn m a b = fst $ chplMin m a b+chplMinUp m a b = snd $ chplMin m a b++chplMin m a b =+    (-u,-l)+    where+    (l,u) = chplMax m (-a) (-b)+ {-|      Approximate from below and  from above the pointwise maximum of two polynomials -}@@ -197,7 +207,7 @@     ERChebPoly box b ->     (ERChebPoly box b, ERChebPoly box b) chplMax maxDegree p1 p2 =-    (- (-p1 - differenceDown), p1 + differenceUp)+    (p1 `plusDown` differenceDown, p1 `plusUp` differenceUp)     where     (differenceDown, differenceUp) = chplNonneg maxDegree $ p2 - p1 @@ -241,9 +251,9 @@         p2 = multiplyByP p3 + (chplConst a2) -- ie p * a3 + a2         p3 = chplConst a3     multiplyByPUp =-        snd . chplReduceDegree maxDegree . (p *)+        chplReduceDegreeUp maxDegree . (p *)     multiplyByPDown =-        fst . chplReduceDegree maxDegree . (p *)+        chplReduceDegreeDown maxDegree . (p *)     {-       The cubic polynomial's coefficients are calculated by solving a system of 4 linear eqs.       The generic solution is as follows:@@ -291,3 +301,37 @@         we subtract its value at 0 rounded upwards.     -}     valueAt0 = chplConst $ a0 / b++{-|+    Multiply a thin enclosure by a non-thin enclosure+-}+chplThinTimesEncl ::+    (B.ERRealBase b, DomainBox box varid Int, Ord box) => +    Int {-^ maximum polynomial degree -} -> +    ERChebPoly box b ->+    (ERChebPoly box b, ERChebPoly box b) ->+    (ERChebPoly box b, ERChebPoly box b)+chplThinTimesEncl maxDegree p1 (p2LO, p2HI) =+    (prodLO, prodHI)+    where+    prodHI =+        chplMaxUp maxDegree +            (p1 `timesUp` p2HI)+            (p1 `timesUp` p2LO) -- beware: p1 can be negative+    prodLO =+        negate $+        chplMaxUp maxDegree +            (p1n `timesUp` p2HI)+            (p1n `timesUp` p2LO)+    p1n = negate p1++{-|+    Safely multiply a polynomial by itself.+-}+chplSquare ::+    (B.ERRealBase b, DomainBox box varid Int, Ord box) => +    Int {-^ maximum polynomial degree -} -> +    ERChebPoly box b ->+    (ERChebPoly box b, ERChebPoly box b)+chplSquare maxDegree p =+    (p `timesDown` p, p `timesUp` p)
src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs view
@@ -136,188 +136,132 @@  {-|     Approximate the pointwise sine of a polynomial -    by another polynomial from below and from above. +    by another polynomial from below and from above.+    +    Assuming the polynomial range is [-pi/2, pi/2].  -}-chplSineCosine ::+chplSine ::     (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>-    Bool {-^ True iff sine, False iff cosine -} ->      Int {-^ maximum polynomial degree -} -> -    EffortIndex {-^ minimum approx Taylor degree -} -> +    EffortIndex {-^ how hard to try (determines Taylor degree and granularity) -} ->      ERChebPoly box b ->     (ERChebPoly box b, ERChebPoly box b)-chplSineCosine isSine maxDegree ix p-    -- p - 2k*pi range within [-pi/2, pi/2]: -    | ranfNear0 `RA.refines` plusMinusPiHalf =---        unsafePrint---        (---            "ERChebPoly: chplSineCosine: [-pi/2, pi/2]: "---            ++ "\n p = " ++ show p---            ++ "\n ranf = " ++ show ranf---            ++ "\n k = " ++ show k---            ++ "\n ranfNear0 = " ++ show ranfNear0---        ) $-        case isSine of-            True -> sineShifted (- k2pi)-            False -> cosineShifted (- k2pi)-    -- p - 2k*pi range within [0, pi]: -    | (ranfNear0 - piHalf) `RA.refines` plusMinusPiHalf =+chplSine maxDegree ix p = --        unsafePrint --        (---            "ERChebPoly: chplSineCosine: [0, pi]: "+--            "ERChebPoly: sineTaylor: " --            ++ "\n p = " ++ show p---            ++ "\n ranf = " ++ show ranf---            ++ "\n k = " ++ show k---            ++ "\n ranfNear0 = " ++ show ranfNear0+--            ++ "\n ranLargerEndpoint = " ++ show ranLargerEndpoint+--            ++ "\n sineUp = " ++ show sineUp+--            ++ "\n sineDown = " ++ show sineDown --        ) $-        case isSine of-            -- use sin(x) = cos(x - pi/2) and cos(x) = - sin(x - pi/2):-            True -> cosineShifted (- k2pi - piHalf)-            False -> sineShiftedNegated (- k2pi - piHalf)-    -- p - 2k*pi range within [-pi, 0]: -    | (ranfNear0 + piHalf) `RA.refines` plusMinusPiHalf =-        case isSine of-            -- use sin(x) = - cos(x + pi/2) and cos(x) = sin(x + pi/2):-            True -> cosineShiftedNegated (-k2pi + piHalf)-            False -> sineShifted (-k2pi + piHalf)-    -- p - 2k*pi range within [pi/2, 3pi/2]: -    | (ranfNear0 - pi) `RA.refines` plusMinusPiHalf =-        -- use sin(x) = - sin(x - pi) and cos(x) = - cos(x - pi)-        case isSine of-            True -> sineShiftedNegated (- k2pi - pi)-            False -> cosineShiftedNegated (- k2pi - pi)-    | otherwise = (chplConst (-1), chplConst 1)---    (expDownwards, expUpwards + valueAtRDnNeg + (chplConst expRUp))-    where-    ranfNear0 = ranf - k2pi -    k2pi = k * 2 * pi-    plusMinusPiHalf = (-piHalfLO) RA.\/ piHalfLO-    pi = RAEL.pi ix  -    piHalf = pi / 2-    (piHalfLO, piHalfHI) = RA.bounds piHalf-    ranf = -        ERInterval -            (negate $ chplUpperBoundAffine 10 (-p)) -            (chplUpperBoundAffine 10 p)-    k = -        fromInteger $ floor $ -            case (pi + ranf) / (2 * pi) of ERInterval lo hi -> lo-            -    sineShiftedNegated shift =-        boundsNegate $ sineShifted shift-        -    cosineShiftedNegated shift =-        boundsNegate $ cosineShifted shift--    boundsNegate (pLO, pHI) = (- pHI, - pLO)-        -    sineShifted shift =-        boundsAddErr shiftWidthB $ sineTaylor (p + shiftPoly) (ranf + shift)-        where-        shiftPoly = chplConst shiftLOB-        ERInterval shiftLOB shiftHIB = shift-        shiftWidthB = shiftHIB - shiftLOB-    -    cosineShifted shift =-        boundsAddErr shiftWidthB $ cosineTaylor (p + shiftPoly) (ranf + shift)-        where-        shiftPoly = chplConst shiftLOB-        ERInterval shiftLOB shiftHIB = shift-        shiftWidthB = shiftHIB - shiftLOB-    -    boundsAddErr errB (pLO, pHI) =-        (pLO `plusDown` (- errPoly), pHI + errPoly)-        where-        errPoly = chplConst errB-    -    sineTaylor x xran =         (sineDown, sineUp)         where-        sineUp =-            chplReduceDegreeUp maxDegree $ -                x * sineUpTaylor + (chplConst sineUpErrorBound)-        (sineUpTaylor, sineUpErrorTermDegree, sineUpErrorTermCoeff) =-            taylorAux x 1 (B.setGranularity coeffGr 1)-        sineUpErrorBound =-            case sineUpErrorBoundRA of ERInterval lo hi -> hi-            where-            sineUpErrorBoundRA =        -                (xranLargerEndpoint ^ (1 + sineUpErrorTermDegree)) * sineUpErrorTermCoeffRA-            sineUpErrorTermCoeffRA =-                abs $-                ERInterval sineUpErrorTermCoeff sineUpErrorTermCoeff-        sineDown = -            negate $ chplReduceDegreeUp maxDegree $ -                x * sineDownTaylorNeg + (chplConst $ sineDownErrorBound)-        (sineDownTaylorNeg, sineDownErrorTermDegree, sineDownErrorTermCoeff) =-            taylorAux x 1 (B.setGranularity coeffGr (-1))-        sineDownErrorBound =-            case sineDownErrorBoundRA of ERInterval lo hi -> hi+        (sineDown, sineUp) =+            boundsAddErr sineErrorBound $+            chplThinTimesEncl maxDegree p (sineDownTaylor, sineUpTaylor) +        ((sineDownTaylor, sineUpTaylor), +         sineErrorTermDegree, +         (sineErrorTermCoeffDown, sineErrorTermCoeffUp)) =+            sincosTaylorAux True (chplSquare maxDegree p) taylorDegree 1 (one, one)+        one = B.setGranularity coeffGr 1+        sineErrorBound =+            case sineErrorBoundRA of ERInterval lo hi -> hi             where-            sineDownErrorBoundRA =-                (xranLargerEndpoint ^ (1 + sineDownErrorTermDegree)) * sineDownErrorTermCoeffRA-            sineDownErrorTermCoeffRA =-                abs $-                ERInterval sineDownErrorTermCoeff sineDownErrorTermCoeff-        xranLargerEndpoint =        -            max (abs xranLO) (abs xranHI)-        (xranLO, xranHI) = RA.bounds xran+            sineErrorBoundRA =        +                (ranLargerEndpointRA ^ (sineErrorTermDegree)) * sineErrorTermCoeffRA+            sineErrorTermCoeffRA =+                ERInterval sineErrorTermCoeff sineErrorTermCoeff+            sineErrorTermCoeff =+                max (abs sineErrorTermCoeffDown) (abs sineErrorTermCoeffUp)+        ranLargerEndpointRA =+            ERInterval ranLargerEndpoint ranLargerEndpoint+        ranLargerEndpoint =+            max (abs ranLO) (abs ranHI)+        ranLO = negate $ chplUpperBoundAffine ix (-p)+        ranHI = chplUpperBoundAffine ix p+        taylorDegree = effIx2int $ ix `div` 3+        coeffGr = effIx2gran $ ix+        +boundsAddErr errB (pLO, pHI) =+    (pLO `plusDown` (- errPoly), pHI + errPoly)+    where+    errPoly = chplConst errB     -    cosineTaylor x xran =+{-|+    Approximate the pointwise sine of a polynomial +    by another polynomial from below and from above.+    +    Assuming the polynomial range is [-pi/2, pi/2]. +-}+chplCosine ::+    (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>+    Int {-^ maximum polynomial degree -} -> +    EffortIndex {-^ how hard to try (determines Taylor degree and granularity) -} -> +    ERChebPoly box b ->+    (ERChebPoly box b, ERChebPoly box b)+chplCosine maxDegree ix p = --        unsafePrint --        (---            "ERChebPoly.Elementary: chplSineCosine: cosineTaylor: "---            ++ "\n xran = " ++ show xran---            ++ "\n cosineUpErrorBound = " ++ show cosineUpErrorBound---            ++ "\n cosineUpErrorTermDegree = " ++ show cosineUpErrorTermDegree---            ++ "\n cosineUpErrorTermCoeff = " ++ show cosineUpErrorTermCoeff---            ++ "\n xranLargerEndpoint = " ++ show xranLargerEndpoint---        )+--            "ERChebPoly: chplCosine: "+--            ++ "\n p = " ++ show p+--            ++ "\n ranLargerEndpoint = " ++ show ranLargerEndpoint+--            ++ "\n cosineUp = " ++ show cosineUp+--            ++ "\n cosineDown = " ++ show cosineDown+--        ) $         (cosineDown, cosineUp)         where-        cosineUp =-            chplReduceDegreeUp maxDegree $ -                cosineUpTaylor + (chplConst cosineUpErrorBound)-        (cosineUpTaylor, cosineUpErrorTermDegree, cosineUpErrorTermCoeff) =-            taylorAux x 0 (B.setGranularity coeffGr 1)-        cosineUpErrorBound-            | odd (cosineUpErrorTermDegree `div` 2) = 0-            | otherwise =-                case cosineUpErrorBoundRA of ERInterval lo hi -> hi-                where-                cosineUpErrorBoundRA =        -                    (xranLargerEndpoint ^ (cosineUpErrorTermDegree)) * cosineUpErrorTermCoeffRA-                cosineUpErrorTermCoeffRA =-                    abs $-                    ERInterval cosineUpErrorTermCoeff cosineUpErrorTermCoeff-        cosineDown = -            negate $ chplReduceDegreeUp maxDegree $ -                cosineDownTaylorNeg + (chplConst $ cosineDownErrorBound)-        (cosineDownTaylorNeg, cosineDownErrorTermDegree, cosineDownErrorTermCoeff) =-            taylorAux x 0 (B.setGranularity coeffGr (-1))-        cosineDownErrorBound -            | even (cosineDownErrorTermDegree `div` 2) = 0-            | otherwise =-                case cosineDownErrorBoundRA of ERInterval lo hi -> hi-                where-                cosineDownErrorBoundRA =-                    (xranLargerEndpoint ^ (cosineDownErrorTermDegree)) * cosineDownErrorTermCoeffRA-                cosineDownErrorTermCoeffRA =-                    abs $-                    ERInterval cosineDownErrorTermCoeff cosineDownErrorTermCoeff-        xranLargerEndpoint =        -            max (abs xranLO) (abs xranHI)-        (xranLO, xranHI) = RA.bounds xran+        (cosineDown, cosineUp) =+            boundsAddErr cosineErrorBound $+            (cosineDownTaylor, cosineUpTaylor) +        ((cosineDownTaylor, cosineUpTaylor), +         cosineErrorTermDegree, +         (cosineErrorTermCoeffDown, cosineErrorTermCoeffUp)) =+            sincosTaylorAux True (chplSquare maxDegree p) taylorDegree 0 (one, one)+        one = B.setGranularity coeffGr 1+        cosineErrorBound =+            case cosineErrorBoundRA of ERInterval lo hi -> hi+            where+            cosineErrorBoundRA =+                (ranLargerEndpointRA ^ (cosineErrorTermDegree)) * cosineErrorTermCoeffRA+            cosineErrorTermCoeffRA =+                ERInterval cosineErrorTermCoeff cosineErrorTermCoeff+            cosineErrorTermCoeff =+                max (abs cosineErrorTermCoeffDown) (abs cosineErrorTermCoeffUp)+        ranLargerEndpointRA =+            ERInterval ranLargerEndpoint ranLargerEndpoint+        ranLargerEndpoint =+            max (abs ranLO) (abs ranHI)+        ranLO = negate $ chplUpperBoundAffine ix (-p)+        ranHI = chplUpperBoundAffine ix p+        taylorDegree = effIx2int $ ix `div` 3+        coeffGr = effIx2gran $ ix     -    taylorAux p0 thisDegree thisCoeff-            | nextDegree > taylorDegree =+sincosTaylorAux ::+    (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>+    Bool -> +    (ERChebPoly box b, ERChebPoly box b) ->+    Int {-^ how far to go in the Taylor series -} ->+    Int {-^ degree of the term now being constructed -} ->+    (b,b) -> +    ((ERChebPoly box b, ERChebPoly box b),+     Int,+     (b,b))+    {-^ +        Bounds for the series result and information about the first discarded term,+        from which some bound on the uniform error can be deduced.+    -} +sincosTaylorAux resultPositive pSquares@(pSquareDown, pSquareUp) +        maxDegree thisDegree (thisCoeffDown, thisCoeffUp)+    | nextDegree > maxDegree = --                unsafePrint --                ( --                    "ERChebPoly: chplSine: taylorAux: " --                    ++ "\n thisCoeff = " ++ show thisCoeff --                    ++ "\n nextDegree = " ++ show nextDegree --                )-                (chplConst thisCoeff, nextDegree, nextCoeff)-            | otherwise =+        ((thisCoeffDownP, thisCoeffUpP), nextDegree, (nextCoeffDown, nextCoeffUp))+    | otherwise = --                unsafePrint --                ( --                    "ERChebPoly: chplSine: taylorAux: "@@ -326,19 +270,174 @@ --                    ++ "\n errorTermCoeff = " ++ show errorTermCoeff --                    ++ "\n errorTermDegree = " ++ show errorTermDegree --                )-                (chplReduceDegreeUp maxDegree $-                    (chplConst thisCoeff) + p0 * p0 * rest,-                 errorTermDegree, errorTermCoeff) -            where-            (rest, errorTermDegree, errorTermCoeff) =-                taylorAux p0 nextDegree nextCoeff-            nextDegree = thisDegree + 2-            nextCoeff = -                thisCoeff / (fromInteger $ negate $ nextDegree * (nextDegree - 1))-    taylorDegree = ix `div` 3-    coeffGr = effIx2gran $ ix+        ((resultDown, resultUp), errorTermDegree, errorTermCoeffs) +    where+    thisCoeffDownP = chplConst thisCoeffDown+    thisCoeffUpP = chplConst thisCoeffUp+    resultDown+                | resultPositive = +                -- ie rest's ideal value is negative and thisCoeff is positive+                    chplReduceDegreeDown maxDegree $+                        thisCoeffDownP `plusDown` (pSquareUp `timesDown` restDown)+                | otherwise =+                -- ie rest's ideal value is positive and thisCoeff is negative+                    chplReduceDegreeDown maxDegree $+                        thisCoeffDownP `plusDown` (pSquareDown `timesDown` restDown)+    resultUp+                | resultPositive = +                -- ie rest's ideal value is negative and thisCoeff is positive+                    chplReduceDegreeUp maxDegree $+                        thisCoeffUpP `plusUp` (pSquareDown `timesUp` restUp)+                | otherwise =+                -- ie rest's ideal value is positive and thisCoeff is negative+                    chplReduceDegreeUp maxDegree $+                        thisCoeffUpP `plusUp` (pSquareUp `timesUp` restUp)+    ((restDown, restUp), errorTermDegree, errorTermCoeffs) =+        sincosTaylorAux (not resultPositive) pSquares maxDegree nextDegree (nextCoeffDown, nextCoeffUp)+    nextDegree = thisDegree + 2+    nextCoeffUp+                | resultPositive = +                    thisCoeffDown / nextCoeffDenominator -- positive / negative+                | otherwise = +                    thisCoeffUp / nextCoeffDenominator -- negative / negative+    nextCoeffDown +                | resultPositive = +                    thisCoeffUp `divDown` nextCoeffDenominator -- positive / negative+                | otherwise = +                    thisCoeffDown `divDown` nextCoeffDenominator -- negative / negative+    nextCoeffDenominator =+        fromInteger $ toInteger $ negate $ nextDegree * (nextDegree - 1)+    divDown a b = negate $ a / (negate b)   {-|+    Approximate the pointwise natural logarithm of a polynomial +    by another polynomial from below and from above. +-}+chplAtan ::+    (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) => +    Int {-^ maximum polynomial degree -} -> +    EffortIndex {-^  ?? -} -> +    ERChebPoly box b ->+    (ERChebPoly box b, ERChebPoly box b)+{- atan using Euler's series: +    x / (1 + x^2) * (1 + t*2*1/(2*1 + 1)*(1 + t*2*2/(2*2 + 1)*(1 + ... (1 + t*2*n/(2*n+1)*(1 + ...)))))+    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]+-}+chplAtan maxDegree ix p +    | avoidingDivBy0 = +--        unsafePrint+--        (+--            "ERChebPoly.Elementary: chplAtan: "+--             ++ "\n maxDegree = " ++ show maxDegree+--             ++ "\n p = " ++ show p+--             ++ "\n pSquareDn = " ++ show pSquareDn+--             ++ "\n pSquareUp = " ++ show pSquareUp+--             ++ "\n pOverPSquarePlus1Dn = " ++ show pOverPSquarePlus1Dn+--             ++ "\n pOverPSquarePlus1Up = " ++ show pOverPSquarePlus1Up+--             ++ "\n preresDn = " ++ show preresDn+--             ++ "\n preresUp = " ++ show preresUp+--             ++ "\n resDn = " ++ show resDn+--             ++ "\n resUp = " ++ show resUp+--        )+        (resDn, resUp)+    | otherwise =+        (chplConst (-2), chplConst 2) -- this is always safe...    +    where+    avoidingDivBy0 = +        (chplUpperBoundAffine ix (- pSquarePlus1Dn) < 0)+        &&+        (chplUpperBoundAffine ix (- pSquarePlus1Up) < 0)+    resDn = +        negate $+        chplMaxUp maxDegree +            (chplReduceDegreeUp maxDegree $ +                pOverPSquarePlus1DnNeg `timesUp` preresDn) -- beware: pOverPSquarePlus1Dn can be negative+            (chplReduceDegreeUp maxDegree $+                pOverPSquarePlus1DnNeg `timesUp` preresUp)+        where+        pOverPSquarePlus1DnNeg = negate pOverPSquarePlus1Dn+    resUp = +        chplMaxUp maxDegree +            (chplReduceDegreeUp maxDegree $+                pOverPSquarePlus1Up `timesUp` preresDn) -- beware: pOverPSquarePlus1Up can be negative+            (chplReduceDegreeUp maxDegree $+                pOverPSquarePlus1Up `timesUp` preresUp)+    (preresDn, preresUp) = seriesDnUp termsCount 2+    termsCount = max 0 $ ix `div` 3+    gran = effIx2gran ix+    seriesDnUp termsCount coeffBase +        | termsCount > 0 =+            (+             chplReduceDegreeDown maxDegree $+             1 `plusDown` +                (pSquareOverPSquarePlus1Dn -- >=0 +                    `timesDown` (chplConst coeffDn) -- >=0 +                    `timesDown` restDn -- >=0+                )+            ,+             chplReduceDegreeUp maxDegree $+             1 `plusUp`+                (pSquareOverPSquarePlus1Up -- >=0 +                    `timesUp` (chplConst coeffUp) -- >=0 +                    `timesUp` restUp -- >=0+                )+            )+        | otherwise =+            (+             1 `plusDown` (pSquareDn `timesDown` (chplConst coeffDn)) -- both >=0+            ,+             1 `plusUp` pSquareUp+            )+        where+        (restDn, restUp) = seriesDnUp (termsCount - 1) (coeffBase + 2)+        coeffUp = coeffBaseB / (coeffBaseB `plusDown` 1)+        coeffDn = negate $ coeffBaseB / (negate $ coeffBaseB `plusUp` 1)+        coeffBaseB = B.setMinGranularity gran $ fromInteger coeffBase+    (pSquareDn, pSquareUp) = chplSquare maxDegree p+    pSquarePlus1Dn = pSquareDn `plusDown` 1+    pSquarePlus1Up = pSquareUp `plusUp` 1+    recipPSquarePlus1Dn = chplRecipDn maxDegree ix pSquarePlus1Up+    recipPSquarePlus1Up = chplRecipUp maxDegree ix pSquarePlus1Dn+--        -- safely compute the square of an enclosure:+--        pSquareDn = chplMinDn m pUpTDnpUp (chplMinDn m pDnTDnpUp pDnTDnpDn)+--        pSquareUp = chplMaxUp m pUpTUppUp (chplMaxUp m pDnTUppUp pDnTUppDn) +--        pUpTDnpUp = pUp `timesDown` pUp+--        pDnTDnpUp = pDn `timesDown` pUp+--        pDnTDnpDn = pDn `timesDown` pDn+--        pUpTUppUp = pUp `timesUp` pUp+--        pDnTUppUp = pDn `timesUp` pUp+--        pDnTUppDn = pDn `timesUp` pDn+--        mMinus1 = m - 1+    pSquareOverPSquarePlus1Up = +        pSquareUp `timesUp` recipPSquarePlus1Up -- both >=0+    pSquareOverPSquarePlus1Dn = +        pSquareDn `timesDown` recipPSquarePlus1Dn -- both >=0 (one enclosure may dip below 0, not a problem)+--        negate $+--        chplMaxUp maxDegree+--            (pSquareDnNeg `timesUp` recipPSquarePlus1Up) -- beware: pSquareDn may dip below 0+--            (pSquareDnNeg `timesUp` recipPSquarePlus1Dn)+--        where+--        pSquareDnNeg = negate pSquareDn+    pOverPSquarePlus1Up =+        chplMaxUp maxDegree +            (p `timesUp` recipPSquarePlus1Up)+            (p `timesUp` recipPSquarePlus1Dn) -- beware: p can be negative+    pOverPSquarePlus1Dn =+        negate $+        chplMaxUp maxDegree+            (pn `timesUp` recipPSquarePlus1Up) -- beware: pn can be positive+            (pn `timesUp` recipPSquarePlus1Dn)+        where+        pn = negate p++chplRecipDn m i = fst . chplRecip m i+chplRecipUp m i = snd . chplRecip m i++{-|     Approximate the pointwise cosine of a polynomial      by another polynomial from below and from above     using the tau method    @@ -372,7 +471,7 @@     lowerB = - (chplUpperBoundAffine ix (- p))     upperB = chplUpperBoundAffine ix p      -    tauDegree = effIx2int (ix `div` 3)+    tauDegree = effIx2int (max 2 $ ix `div` 3)     coeffGr = effIx2gran $ ix          -- translate p to have range above 1: