summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorMichalKonecny <>2008-08-20 21:25:16 (GMT)
committerLuite Stegeman <luite@luite.com>2008-08-20 21:25:16 (GMT)
commitc8332b5c28876c7de09677eae9feaa3ef559dbd7 (patch)
treedfbcd859b2e49cdb6d29be016b69c3a32960ba24
parent9533a4bcb1f771e4b61b43cfb8b1cedb781d7a15 (diff)
version 0.40.4
-rw-r--r--AERN-RnToRm.cabal2
-rw-r--r--ChangeLog3
-rw-r--r--src/Data/Number/ER/RnToRm/Approx/DomEdges.hs1
-rw-r--r--src/Data/Number/ER/RnToRm/Approx/DomTransl.hs2
-rw-r--r--src/Data/Number/ER/RnToRm/Approx/PieceWise.hs13
-rw-r--r--src/Data/Number/ER/RnToRm/Approx/Tuple.hs1
-rw-r--r--src/Data/Number/ER/RnToRm/TestingDefs.hs12
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs177
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/Base.hs14
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs7
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Bounds.hs50
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs439
12 files changed, 517 insertions, 204 deletions
diff --git a/AERN-RnToRm.cabal b/AERN-RnToRm.cabal
index d695337..f58e4ea 100644
--- 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
index bd4f9ce..0280297 100644
--- 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
index b5176fe..15a5fb5 100644
--- a/src/Data/Number/ER/RnToRm/Approx/DomEdges.hs
+++ b/src/Data/Number/ER/RnToRm/Approx/DomEdges.hs
@@ -251,6 +251,7 @@ instance
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
index 780770f..02ea605 100644
--- a/src/Data/Number/ER/RnToRm/Approx/DomTransl.hs
+++ b/src/Data/Number/ER/RnToRm/Approx/DomTransl.hs
@@ -323,6 +323,8 @@ instance
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
index 968b506..dce8881 100644
--- a/src/Data/Number/ER/RnToRm/Approx/PieceWise.hs
+++ b/src/Data/Number/ER/RnToRm/Approx/PieceWise.hs
@@ -232,6 +232,7 @@ instance
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 @@ instance
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
index 4b37c15..be2e88d 100644
--- a/src/Data/Number/ER/RnToRm/Approx/Tuple.hs
+++ b/src/Data/Number/ER/RnToRm/Approx/Tuple.hs
@@ -203,6 +203,7 @@ instance
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
index 3b1b57f..58b0784 100644
--- 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 @@ fapeUConst1 = (FA.const (DBox.unary $ (0)RA.\/1) [1]) :: FAPE
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 @@ testIntegrE =
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
index 86bdc3f..4e62932 100644
--- a/src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs
+++ b/src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs
@@ -124,9 +124,9 @@ instance
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,9 +135,10 @@ instance
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) =>
@@ -312,7 +313,7 @@ instance
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 @@ instance
-- ++ "\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
index 0a041b3..d3dfd40 100644
--- a/src/Data/Number/ER/RnToRm/UnitDom/Base.hs
+++ b/src/Data/Number/ER/RnToRm/UnitDom/Base.hs
@@ -237,7 +237,8 @@ class
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,13 +246,22 @@ class
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) -} ->
+ (ufb, ufb)
{-|
Evaluate at a point, rounding upwards and downwards.
-}
diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs
index 809f92b..9304cf4 100644
--- 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 @@ instance
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
index 68b2699..f10847a 100644
--- 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 @@ chplUpperBoundQuadr ix p@(ERChebPoly coeffs) =
_ ->
[(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 @@ chplMax ::
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 @@ chplNonnegCubic maxDegree p
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 @@ chplNonnegCubic maxDegree p
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
index 20b80fd..b3ad053 100644
--- 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 @@ chplLog maxDegree ix p =
{-|
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 =
+chplSine maxDegree ix p =
-- unsafePrint
-- (
--- "ERChebPoly: chplSineCosine: [-pi/2, pi/2]: "
+-- "ERChebPoly: sineTaylor: "
-- ++ "\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 =
--- unsafePrint
--- (
--- "ERChebPoly: chplSineCosine: [0, pi]: "
--- ++ "\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
+ (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
- 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
- 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
+
+{-|
+ Approximate the pointwise sine of a polynomial
+ by another polynomial from below and from above.
- cosineTaylor x xran =
+ 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,17 +270,172 @@ chplSineCosine isSine maxDegree ix p
-- ++ "\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
@@ -372,7 +471,7 @@ chplRecip maxDegree ix p
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: