diff --git a/AERN-RnToRm.cabal b/AERN-RnToRm.cabal --- a/AERN-RnToRm.cabal +++ b/AERN-RnToRm.cabal @@ -1,5 +1,5 @@ Name: AERN-RnToRm -Version: 0.3.0.3 +Version: 0.4 Cabal-Version: >= 1.2 Build-Type: Simple License: BSD3 diff --git a/ChangeLog b/ChangeLog --- a/ChangeLog +++ b/ChangeLog @@ -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 diff --git a/src/Data/Number/ER/RnToRm/Approx/DomEdges.hs b/src/Data/Number/ER/RnToRm/Approx/DomEdges.hs --- a/src/Data/Number/ER/RnToRm/Approx/DomEdges.hs +++ b/src/Data/Number/ER/RnToRm/Approx/DomEdges.hs @@ -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) => diff --git a/src/Data/Number/ER/RnToRm/Approx/DomTransl.hs b/src/Data/Number/ER/RnToRm/Approx/DomTransl.hs --- a/src/Data/Number/ER/RnToRm/Approx/DomTransl.hs +++ b/src/Data/Number/ER/RnToRm/Approx/DomTransl.hs @@ -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, diff --git a/src/Data/Number/ER/RnToRm/Approx/PieceWise.hs b/src/Data/Number/ER/RnToRm/Approx/PieceWise.hs --- a/src/Data/Number/ER/RnToRm/Approx/PieceWise.hs +++ b/src/Data/Number/ER/RnToRm/Approx/PieceWise.hs @@ -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) = diff --git a/src/Data/Number/ER/RnToRm/Approx/Tuple.hs b/src/Data/Number/ER/RnToRm/Approx/Tuple.hs --- a/src/Data/Number/ER/RnToRm/Approx/Tuple.hs +++ b/src/Data/Number/ER/RnToRm/Approx/Tuple.hs @@ -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) => diff --git a/src/Data/Number/ER/RnToRm/TestingDefs.hs b/src/Data/Number/ER/RnToRm/TestingDefs.hs --- a/src/Data/Number/ER/RnToRm/TestingDefs.hs +++ b/src/Data/Number/ER/RnToRm/TestingDefs.hs @@ -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! diff --git a/src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs b/src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs --- a/src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs +++ b/src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs @@ -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) => diff --git a/src/Data/Number/ER/RnToRm/UnitDom/Base.hs b/src/Data/Number/ER/RnToRm/UnitDom/Base.hs --- a/src/Data/Number/ER/RnToRm/UnitDom/Base.hs +++ b/src/Data/Number/ER/RnToRm/UnitDom/Base.hs @@ -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) -} -> diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs --- a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs @@ -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 = diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Bounds.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Bounds.hs --- a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Bounds.hs +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Bounds.hs @@ -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) diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs --- a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs +++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs @@ -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: