{-# LANGUAGE FlexibleContexts #-}
{-|
    Module      :  Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds
    Description :  (internal) bounds of single and multiple polynomials
    Copyright   :  (c) 2007-2008 Michal Konecny
    License     :  BSD3

    Maintainer  :  mik@konecny.aow.cz
    Stability   :  experimental
    Portability :  portable

    Internal module for "Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom".
    
    Implementation of various functions related to the bounds of polynomials.    
-}
module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds 
where

import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Ring
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Reduce
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Eval
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Derivative

import qualified Data.Number.ER.Real.Approx as RA
import qualified Data.Number.ER.Real.Base as B
import Data.Number.ER.Real.Approx.Interval
import Data.Number.ER.Real.Arithmetic.LinearSolver
import qualified Data.Number.ER.BasicTypes.DomainBox as DBox
import Data.Number.ER.BasicTypes.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox)
import Data.Number.ER.BasicTypes
import Data.Number.ER.Misc

import qualified Data.Map as Map

import Data.List

{-|
    Find an upper bound on a polynomial over the 
    unit domain [-1,1]^n.  

    Quick method that does not converge to exact result with increasing 
    effort index.
-}
chplUpperBound ::
    (B.ERRealBase b, 
     DomainBox box varid Int, Ord box, Show varid,
     DomainIntBox boxra varid (ERInterval b),
     DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) => 
    EffortIndex {-^ how hard to try -} ->
    ERChebPoly box b ->
    b
chplUpperBound ix p = snd $ chplBounds ix p

{-|
    Find a lower bound on a polynomial over the 
    unit domain [-1,1]^n.  

    Quick method that does not converge to exact result with increasing 
    effort index.
-}
chplLowerBound ::
    (B.ERRealBase b, 
     DomainBox box varid Int, Ord box, Show varid,
     DomainIntBox boxra varid (ERInterval b),
     DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) => 
    EffortIndex {-^ how hard to try -} ->
    ERChebPoly box b ->
    b
chplLowerBound ix p = fst $ chplBounds ix p

{-|
    Find both lower and upper bounds on a polynomial over the 
    unit domain [-1,1]^n.  

    Quick method that does not converge to exact result with increasing 
    effort index.
-}
chplBounds ::
    (B.ERRealBase b, 
     DomainBox box varid Int, Ord box, Show varid,
     DomainIntBox boxra varid (ERInterval b),
     DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) => 
    EffortIndex {-^ how hard to try -} ->
    ERChebPoly box b ->
    (b,b)
chplBounds = 
    chplBoundsAffine

{-|
    Find an upper bound on a polynomial over the 
    unit domain [-1,1]^n.  
-}
chplUpperBoundExpensive ::
    (B.ERRealBase b, 
     DomainBox box varid Int, Ord box, Show varid,
     DomainIntBox boxra varid (ERInterval b),
     DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) => 
    EffortIndex {-^ how hard to try -} ->
    ERChebPoly box b ->
    b
chplUpperBoundExpensive ix p = snd $ chplBoundsExpensive ix p

{-|
    Find a lower bound on a polynomial over the 
    unit domain [-1,1]^n.  
    
    Quick method that does not converge to exact result with increasing 
    effort index.
-}
chplLowerBoundExpensive ::
    (B.ERRealBase b, 
     DomainBox box varid Int, Ord box, Show varid,
     DomainIntBox boxra varid (ERInterval b),
     DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) => 
    EffortIndex {-^ how hard to try -} ->
    ERChebPoly box b ->
    b
chplLowerBoundExpensive ix p = fst $ chplBoundsExpensive ix p

{-|
    Find both lower and upper bounds on a polynomial over the 
    unit domain [-1,1]^n.
    
    Quick method that does not converge to exact result with increasing 
    effort index.
-}
chplBoundsExpensive ::
    (B.ERRealBase b, 
     DomainBox box varid Int, Ord box, Show varid,
     DomainIntBox boxra varid (ERInterval b),
     DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) => 
    EffortIndex {-^ how hard to try -} ->
    ERChebPoly box b ->
    (b,b)
chplBoundsExpensive = chplBoundsByDerivative

{-|
    Find bounds on a polynomial over the unit domain [-1,1]^n.
    
    Fast but inaccurate method, in essence
    taking the maximum of the upper affine reduction.
-}
chplBoundsAffine ::
    (B.ERRealBase b, DomainBox box varid Int, Ord box) => 
    EffortIndex {-^ how hard to try -} ->
    ERChebPoly box b ->
    (b,b)
chplBoundsAffine ix p@(ERChebPoly coeffs) =
--    unsafePrintReturn
--    (
--        "chplBoundsAffine:"
--        ++ "\n p = " ++ show p
--        ++ "\n noConstCoeffAbsSum = " ++ show noConstCoeffAbsSum
--        ++ "\n result = "
--    )    
    result
    where
    result =
        (constTerm `plusDown` (- noConstCoeffAbsSum),
         constTerm `plusUp` noConstCoeffAbsSum)
    noConstCoeffAbsSum = Map.fold plusUp 0 absCoeffs
    absCoeffs = Map.map abs $ Map.delete chplConstTermKey coeffs
    constTerm = Map.findWithDefault 0 chplConstTermKey coeffs


{-|
    Find a close upper bound of a polynomial over the 
    unit domain [-1,1]^n.
    
    Approximates all local extrema and computes their maximum.  
-}
chplBoundsByDerivative ::
    (B.ERRealBase b, 
     DomainBox box varid Int, Ord box, Show varid,
     DomainIntBox boxra varid (ERInterval b),
     DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) => 
    EffortIndex {-^ how hard to try looking for peaks -} ->
    ERChebPoly box b ->
    (b,b)
chplBoundsByDerivative ix p =
--    unsafePrint
--    (
--        "chplBoundsByDerivative: "
--        ++ "\n extremaValues = " ++ show extremaValues
--    ) $
    (lowerBound, upperBound)
    where
    lowerBound = foldl1 min $ map fst extremaValues 
    upperBound = foldl1 max $ map snd extremaValues
    ra2bb (ERInterval l r) = (l,r)
    b2ra b = ERInterval b b
    extremaValues = 
        collectValuesOnFaces vars varDerivatives (p,p)
        where
        vars = chplGetVars p
        varDerivatives = -- var |-> (lower, upper) bounds on partial derivative
            Map.fromList $
                map getDerivatives vars
        getDerivatives var =
            (var,
                chplBall2DownUp $
                    ballDifferentiate p var)
    collectValuesOnFaces varsSpecialise varDerivatives (pDown, pUp) =
--        unsafePrint
--        (
--            "chplBoundsByDerivative: collectValuesOnFaces: "
--            ++ "\n vars = " ++ (show $ Map.keys varDerivatives)
--            ++ "\n valuesThisFace = " ++ show valuesThisFace
--        ) $
        valuesThisFace ++ (valuesSubFaces varsSpecialise)
        where
        valuesThisFace =
            collectExtremeValues varDerivatives (pDown, pUp)
        valuesSubFaces [] = []
        valuesSubFaces (var : vars) =
            (collectValuesOnFaces vars varDerivativesNoVarL (pDownNoVarL, pUpNoVarL))
            ++
            (collectValuesOnFaces vars varDerivativesNoVarR (pDownNoVarR, pUpNoVarR))
            ++
            (valuesSubFaces vars)
            where
            (pDownNoVarR, pUpNoVarR) = substVarR (pDown, pUp)
            (pDownNoVarL, pUpNoVarL) = substVarL (pDown, pUp)
            substVarL = substVar (-1)
            substVarR = substVar 1
            substVar val (pDown, pUp) =
                (fst $ chplPartialRAEval ra2bb pDown $ DBox.singleton var val,
                 snd $ chplPartialRAEval ra2bb pUp $ DBox.singleton var val)
            varDerivativesNoVarL =
                Map.map substVarL varDerivativesNoVar 
            varDerivativesNoVarR =
                Map.map substVarR varDerivativesNoVar 
            varDerivativesNoVar = 
                Map.delete var varDerivatives
    collectExtremeValues varDerivatives (pDown, pUp)
        | null varsNoConst =
--            unsafePrint
--            (
--                "chplBoundsByDerivative: collectExtremeValues:" 
--                ++ "\n null varsNoConst"
--                ++ "\n varDerivatives = " ++ show varDerivatives
--            )
            -- corner or near constant function 
            [pEvalAt unitDomBox] 
        | otherwise =
--            unsafePrint
--            (
--                "chplBoundsByDerivative: collectExtremeValues:" 
--                ++ "\n varDerivatives = " ++ show varDerivatives
--                ++ "\n boxesWithPotentialExtrema = " ++ show boxesWithPotentialExtrema
--            ) $
            map pEvalAt boxesWithPotentialExtrema
        where
        boxesWithPotentialExtrema = 
            paveFindBoxes [(unitDomBox,0)] 
        varDerivativesNoZeros =
            Map.filter (not . isConstWithZero) varDerivatives
            where
            isConstWithZero (pDown, pUp) =
                (snd $ chplBoundsAffine ix pDown) <= 0
                &&
                (fst $ chplBoundsAffine ix pUp) >= 0
--                case (chplGetConst pDown, chplGetConst pUp) of
--                    (Just cDown, Just cUp) ->
--                        cDown <= 0 && cUp >= 0 
--                    _ -> False 
        vars = Map.keys varDerivatives
        varsNoConst = Map.keys varDerivativesNoZeros
        varsNoConstLength = length varsNoConst
        pEvalAt = evalAt (pDown, pUp)
        evalAt (pDown,pUp) box =
            (fst $ ra2bb $ chplRAEval b2ra pDown box,
             snd $ ra2bb $ chplRAEval b2ra pUp box)
        unitDomBox =
            DBox.fromList $ zip vars (repeat unitInterval)
        unitInterval = ((-1) RA.\/ 1)
        maxDepth = fromInteger $ toInteger $ max 3 ix
        keepBox box =
            and $ map evalDeriv $ Map.elems varDerivativesNoZeros
            where
            evalDeriv derivBounds = hasZero $ evalAt derivBounds box  
            hasZero (l,h) = l <= 0 && h >= 0
        paveFindBoxes [] = [] 
        paveFindBoxes boxes@((box, depth) : boxesRest) 
            | keepBox box =
                case depth < maxDepth of
                    True ->
                        paveFindBoxes ((boxL, newDepth) : (boxR, newDepth) : boxesRest)
                    False ->
                        box : (paveFindBoxes boxesRest)
            | otherwise =
                paveFindBoxes boxesRest
            where
            var = varsNoConst !! (depth `mod` varsNoConstLength)
            (boxL, boxR) = DBox.split box var Nothing
            newDepth = depth + 1


--{-|
--    Find a close upper bound on a quadratic polynomial over the 
--    unit domain [-1,1]^n.  
--
--    Much slower and somewhat more accurate method, in essence
--    taking the maximum of the upper quadratic reduction.
--    
--    !!! Not properly tested !!!
---}
--chplUpperBoundQuadr ::
--    (B.ERRealBase b, DomainBox box varid Int, Ord box,
--     DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b],
--     DomainBoxMappable boxra boxra varid (ERInterval b) (ERInterval b), 
--     DomainIntBox boxra varid (ERInterval b), Num varid, Enum varid) => 
--    EffortIndex {-^ how hard to try looking for peaks -} ->
--    ERChebPoly box b ->
--    b
--chplUpperBoundQuadr ix p@(ERChebPoly coeffs) =
--    quadBound (coeffsQ, vars)
--    where
--    pQ@(ERChebPoly coeffsQ) = chplReduceDegreeUp 2 p
--    vars = chplGetVars pQ
--    quadBound (coeffs, vars)
--        | null vars =
--            Map.findWithDefault 0 chplConstTermKey coeffs
--        | hasInteriorPeak =
--            foldl max peakValue edgeBounds
--        | otherwise =
--            foldl1 max edgeBounds
--        where
--        edgeBounds =
--            map quadBound $ concat $ map removeVar vars
--        (hasInteriorPeak, peakValue) =
--            case maybePeak of
--                Just peak ->
--                    (noPositiveSquare -- if any term x^2 has a positive coeff, there is no peak  
--                     &&
--                     (and $ map maybeInUnit $ DBox.elems peak)
--                    ,
--                     erintv_right $
--                     chplRAEval makeInterval p peak
--                    )
--                Nothing -> (False, undefined)
--            where
--            noPositiveSquare =
--                and $ map (<= 0) $ map getQuadCoeff vars
--            getQuadCoeff var = 
--                Map.findWithDefault 0 (DBox.singleton var 2) coeffs
--            maybeInUnit r =
--                case (RA.compareReals r (-1), RA.compareReals (1) r) of
--                    (Just LT, _) -> False -- ie r < -1
--                    (_, Just LT) -> False -- ie r > 1
--                    _ -> True
--        maybePeak =
--            linearSolver
--                (map derivZeroLinearEq vars)
--                (DBox.fromList $ map (\v -> (v,(-1) RA.\/ 1)) vars)
--                (2^^(-ix))
--            where
--            derivZeroLinearEq var =
--                (linCoeffs, - constCoeff)
--                where
--                constCoeff =
--                    makeInterval $
--                    Map.findWithDefault 0 (DBox.singleton var 1) coeffs
--                      -- recall T_1(x) = x, T_1'(x) = 1
--                linCoeffs =
--                    DBox.fromList $
--                        (var, 4 * quadCoeff) -- T_2(x) = 2*x^2 - 1; T_2'(x) = 4*x
--                        : (map getVarVarCoeff $ var `delete` vars)
--                quadCoeff =
--                    makeInterval $
--                    Map.findWithDefault 0 (DBox.singleton var 2) coeffs
--                getVarVarCoeff var2 =
--                    (var2,
--                      makeInterval $
--                      Map.findWithDefault 0 (DBox.fromList [(var,1), (var2,1)]) coeffs)
--        makeInterval b = ERInterval b b
--        removeVar var =
--            [(substVar True, newVars), 
--             (substVar False, newVars)]
--            where
--            newVars = var `delete` vars
--            substVar isOne =
--                chplCoeffs $
--                    foldl (+^) (chplConst 0) $ 
--                        map (makeMonomial isOne) $ 
--                            Map.toList coeffs
--            makeMonomial isOne (term, coeff) =
--                ERChebPoly $ Map.fromList $
--                case (DBox.toList term) of
--                    [(v,2)] | v == var ->
--                        [(chplConstTermKey, coeff)]
--                    [(v,1)] | v == var ->
--                        [(chplConstTermKey, 
--                          case isOne of True -> coeff; False -> - coeff)]
--                    [(v1,1), (v2,1)] | v1 == var ->
--                        [(DBox.fromList [(v2,1)], 
--                          case isOne of True -> coeff; False -> - coeff)]
--                    [(v1,1), (v2,1)] | v2 == var ->
--                        [(DBox.fromList [(v1,1)], 
--                          case isOne of True -> coeff; False -> - coeff)]
--                    _ ->
--                        [(term, coeff)]

{-|
     Approximate from below and  from above the pointwise maximum of two polynomials
-}
chplMax ::
    (B.ERRealBase b, 
     DomainBox box varid Int, Ord box, Show varid,
     DomainIntBox boxra varid (ERInterval b),
     DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) => 
    Int {-^ maximum polynomial degree -} -> 
    Int {-^ maximum term count -} -> 
    ERChebPoly box b ->
    ERChebPoly box b ->
    (ERChebPoly box b, ERChebPoly box b)
chplMax maxDegree maxSize p1 p2 =
    (p1 +. differenceDown, p1 +^ differenceUp)
    where
    (differenceDown, _) = chplNonneg maxDegree maxSize p2MinusP1Down
    (_, differenceUp) = chplNonneg maxDegree maxSize $ p2MinusP1Up
    (p2MinusP1Down, p2MinusP1Up) = chplBall2DownUp $ ballAdd p2 (chplNeg p1)

chplMaxDn m s a b = fst $ chplMax m s a b
chplMaxUp m s a b = snd $ chplMax m s a b
chplMinDn m s a b = fst $ chplMin m s a b
chplMinUp m s a b = snd $ chplMin m s a b

{-|
     Approximate from below and  from above the pointwise minimum of two polynomials
-}
chplMin ::
    (B.ERRealBase b, 
     DomainBox box varid Int, Ord box, Show varid,
     DomainIntBox boxra varid (ERInterval b),
     DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) => 
    Int {-^ maximum polynomial degree -} -> 
    Int {-^ maximum term count -} -> 
    ERChebPoly box b ->
    ERChebPoly box b ->
    (ERChebPoly box b, ERChebPoly box b)
chplMin m s a b =
    (chplNeg u,chplNeg l)
    where
    (l,u) = chplMax m s (chplNeg a) (chplNeg b)

chplNonnegDown m s p = fst $ chplNonneg m s p
chplNonnegUp m s p = snd $ chplNonneg m s p 
chplNonposDown m s p = fst $ chplNonpos m s p
chplNonposUp m s p = snd $ chplNonpos m s p 

chplNonpos m s p =
    (chplNeg h, chplNeg l)
    where
    (l,h) = chplNonneg m s (chplNeg p)

{-|
     Approximate the function max(0,p(x)) by a polynomial from below
     and from above. 
-}
chplNonneg ::
    (B.ERRealBase b, 
     DomainBox box varid Int, Ord box, Show varid,
     DomainIntBox boxra varid (ERInterval b),
     DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) => 
    Int {-^ maximum polynomial degree -} -> 
    Int {-^ maximum term count -} -> 
    ERChebPoly box b ->
    (ERChebPoly box b, ERChebPoly box b)
chplNonneg = chplNonnegCubic

{-|
    A version of 'chplNonneg' using a cubic approximation. 
-}
chplNonnegCubic ::
    (B.ERRealBase b, 
     DomainBox box varid Int, Ord box, Show varid,
     DomainIntBox boxra varid (ERInterval b),
     DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) => 
    Int {-^ maximum polynomial degree -} -> 
    Int {-^ maximum term count -} -> 
    ERChebPoly box b ->
    (ERChebPoly box b, ERChebPoly box b)
chplNonnegCubic maxDegree maxSize p
    | upperB <= 0 = (chplConst 0, chplConst 0)
    | lowerB >= 0 = (p, p)
    | not allInterimsBounded = (chplConst (B.plusInfinity), chplConst (B.plusInfinity))
    | otherwise = -- ie lowerB < 0 < upperB: polynomial may be crossing 0...
--        unsafePrintReturn
--        (
--            "chplNonnegCubic:"
--            ++ "\n p = " ++ show p
--            ++ "\n maxDegree = " ++ show maxDegree
--            ++ "\n maxSize = " ++ show maxSize
--            ++ "\n upperB = " ++ show upperB
--            ++ "\n lowerB = " ++ show lowerB
--            ++ "\n a0 = " ++ show a0
--            ++ "\n a1 = " ++ show a1
--            ++ "\n a2 = " ++ show a2
--            ++ "\n a3 = " ++ show a3
--            ++ "\n b = " ++ show b
--            ++ "\n rb = " ++ show rb
--            ++ "\n correctionB = " ++ show correctionB
--            ++ "\n valueAt0B = " ++ show valueAt0B
--            ++ "\n result = "
--        )
        -- work out the cubic polynomial (a3*x^3 + a2*x^2 + a1*x + a0) / b 
        -- that hits 0 at lowerB with derivative 0 
        -- and hits upperB at upperB with derivative 1 
        (chplAddConstDown (- valueAt0B) cubicAppliedOnPDown, 
         chplAddConstUp correctionB cubicAppliedOnPUp)
    where    
    (lowerB, upperB) = chplBounds 10 p
    (cubicAppliedOnPDown, cubicAppliedOnPUp, width) =
        p0 `scaleByPositiveConsts` (rbLo, rbHi)
        where
        p0 = (multiplyByP p1) `addConsts` (a0Lo, a0Hi) -- ie p*(p*(p * a3 + a2) + a1) + a0 enclosure
        p1 = (multiplyByP p2) `addConsts` (a1Lo, a1Hi) -- ie p*(p * a3 + a2) + a1 enclosure
        p2 = (multiplyByP p3) `addConsts` (a2Lo, a2Hi) -- ie p * a3 + a2 enclosure
        p3 = (chplConst a3Lo, chplConst a3Hi, a3Hi - a3Lo) -- ie a3 enclosure
    multiplyByP (lo,hi,wd) =
        (ploRed, phiRed, pwd)
        where
        ploRed = reduceDgSzDown plo
        phiRed = reduceDgSzUp phi 
        pwd = chplUpperBound 10 $ phiRed -^ ploRed 
        (plo, phi, _) = chplTimesLoHi p (lo,hi,wd)
    reduceDgSzUp =
        chplReduceTermCountUp maxSize . chplReduceDegreeUp maxDegree
    reduceDgSzDown =
        chplReduceTermCountDown maxSize . chplReduceDegreeDown maxDegree
    addConsts (lo, hi, wd) (cLo, cHi) =
        (alo, ahi, wd + wdlo + wdhi)
        where
        (alo, _, wdlo) = chplBall2DownUpWd $ ballAddConst cLo lo 
        (_, ahi, wdhi) = chplBall2DownUpWd $ ballAddConst cHi hi 
    scaleByPositiveConsts (lo, hi, wd) (cLo, cHi) =
        (slo, shi, wd + wdlo + wdhi)
        where
        (slo, _, wdlo) = chplBall2DownUpWd $ ballScale cLo lo 
        (_, shi, wdhi) = chplBall2DownUpWd $ ballScale cHi hi 
    
    -- convert interval coefficients to pairs of bounds:
    ERInterval rbLo rbHi = rb
    ERInterval a3Lo a3Hi = a3
    ERInterval a2Lo a2Hi = a2
    ERInterval a1Lo a1Hi = a1
    ERInterval a0Lo a0Hi = a0
    allInterimsBounded = 
        and $ map RA.isBounded [w, s, rb, a0, a1, a2, a3, correction]
    {-
      The cubic polynomial's coefficients are calculated by solving a system of 4 linear eqs.
      The generic solution is as follows:
         b = (r - l)^3   always positive
         a3 = -(r + l)
         a2 = 2*(r^2 + r*l + l^2)
         a1 = -l*(4*r^2 + r*l + l^2)
         a0 = 2*r^2*l^2
    -}
    rb = recip b
    b = w3 -- = w^3 -- see below
    w = r - l
    r = ERInterval upperB upperB
    l = ERInterval lowerB lowerB
    --
    a3 = - s
    s = r + l
    --
    a2 = 2 * (r2PrlPl2)
    r2PrlPl2 = s2 - rl
    rl = r * l
    --
    a1 = (- l) * (r2PrlPl2 + 3*r2)
    a0 = 2*r2*l2
    -- interval arithmetic tricks to speed it up and reduce dependency errors:
    w3 = ERInterval (wLo * wLo * wLo) (wHi * wHi * wHi) -- x^3 is monotone 
    ERInterval wLo wHi = w
    s2 = ERInterval (max 0 s2Lo) s2Hi
    s2Lo = min sLo2 sHi2 
    s2Hi = max sLo2 sHi2
    sLo2 = sLo * sLo
    sHi2 = sHi * sHi 
    ERInterval sLo sHi = s    
    r2 = ERInterval (upperB `timesDown` upperB) (upperB `timesUp` upperB)    
    l2 = ERInterval (lowerB `timesDown` lowerB) (lowerB `timesUp` lowerB)
    {- 
        The cubic polynomial may sometimes fail to dominate
        x or sometimes it dips below 0.
        Work out the amount by which it has to be lifted up
        to fix these problems. 
    -}
    ERInterval _ correctionB = correction
    correction =
        case (RA.compareReals (2 * r2) (l*s), RA.compareReals (2 * l2) (r*s)) of
            (Just LT, _) ->
                (peak0 * (peak0 * (peak0 * (-a3) - a2) - a1) - a0) / b
            (_, Just LT) ->
                ((peakP * (peakP * (peakP * (-a3) - a2) - a1) - a0) / b) + peakP
            _ -> 0
        where
        peak0 = (l + 4*r2/s) / 3 
        peakP = (r + 4*l2/s) / 3
    {-
        The same cubic polynomial can be used as a lower bound when
        we subtract its value at 0 rounded upwards.
    -}
    valueAt0B = 
        case a0 / b of
            ERInterval lo hi -> hi

{-|
    Multiply a polynomial by an enclosure (with non-negated lower bound).
-}
chplTimesLoHi ::
    (B.ERRealBase b, 
     DomainBox box varid Int, Ord box, Show varid,
     DomainIntBox boxra varid (ERInterval b),
     DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b]) => 
    ERChebPoly box b ->
    (ERChebPoly box b, ERChebPoly box b, b) ->
    (ERChebPoly box b, ERChebPoly box b, b)
chplTimesLoHi p1 (p2Low, p2High, p2Width) =
    (prodMid -. (chplConst width), 
     prodMid +^ (chplConst width),
     2 * width)
    where
    prodMid = prodLowUp
    (prodLowDown, prodLowUp, prodLowWidth) = 
        chplBall2DownUpWd $ ballMultiply p1 p2Low
    (prodHighDown, prodHighUp, prodHighWidth) = 
        chplBall2DownUpWd $ ballMultiply p1 p2High
    width = 
        p1Norm `timesUp` p2Width `plusUp` prodLowWidth `plusUp` prodHighWidth
    p1Norm = 
        max (abs $ p1LowerBound) (abs $ p1UpperBound)
    (p1LowerBound, p1UpperBound) = 
        chplBounds ix p1
    ix = 10