{-# LANGUAGE FlexibleContexts #-}
{-|
    Module      :  Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Elementary
    Description :  (internal) elementary functions applied to 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 elementary functions applied to polynomials.
-}
module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Elementary 
where

import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Eval
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Field
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds

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

import qualified Data.Map as Map

{-|
    Approximate the pointwise square root of a polynomial 
    by another polynomial from below and from above. 
-}
chplSqrt ::
    (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)
chplSqrt maxDegree ix p =
    error "ERChebPoly: chplSqrt: not implemented yet"

{-|
    Approximate the pointwise exponential of a polynomial 
    by another polynomial from below and from above. 
-}
chplExp ::
    (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) => 
    Int {-^ maximum polynomial degree -} -> 
    EffortIndex {-^ minimum approx Taylor degree -} -> 
    ERChebPoly box b ->
    (ERChebPoly box b, ERChebPoly box b)
chplExp maxDegree ix p =
    (expDownwards, expUpwards + valueAtRDnNeg + (chplConst expRUp))
    where
    expUpwards =
        (chplConst expMUp) * (chplPow maxDegree (expNear0Up pNear0Up) a_int) 
    expDownwards =
        (chplConst expMDn) * (chplPow maxDegree (expNear0Dn pNear0Dn) a_int) 
    upperB = chplUpperBoundAffine ix p 
    lowerB = - (chplUpperBoundAffine ix (- p))
    m = (upperB + lowerB) / 2
    r = (upperB - lowerB) / 2 
    expMUp = erintv_right expM 
    expMDn = erintv_left expM
    expM = erExp_R ix (ERInterval m m)
    pNear0Up = (p - (chplConst m)) * (chplConst $ recip a_base)
    pNear0Dn = - (((-p) + (chplConst m)) * (chplConst $ recip a_base))
    a_base = fromInteger a_int
    a_int = max 1 $ floor r -- could this be too high?
    expNear0Up p0 =
        expAux p0 1 (B.setGranularity coeffGr 1)
    expNear0Dn p0 =
        negate $ expAux p0 1 (B.setGranularity coeffGr (-1))
    expAux p0 nextDegree thisCoeff
            | nextDegree > taylorDegree =
                chplConst thisCoeff
            | otherwise =
                snd $ chplReduceDegree maxDegree $
                (chplConst thisCoeff) + p0 * (expAux p0 (nextDegree + 1) nextCoeff)
            where
            nextCoeff = 
                thisCoeff / (fromInteger nextDegree)
    taylorDegree = 1 + 2 * (ix `div` 6)
    coeffGr = effIx2gran $ 10 + 3 * taylorDegree
    expRUp = erintv_right expR
    expR = erExp_R ix (ERInterval r r)
    valueAtRDnNeg = 
        expAux (chplConst r) 1 (B.setGranularity coeffGr (-1))

    
{-|
    Approximate the pointwise integer power of a polynomial by another polynomial from above. 
-}
chplPow ::
    (B.ERRealBase b, Integral i, DomainBox box varid Int, Ord box) => 
    Int {-^ maximum polynomial degree -} -> 
    ERChebPoly box b ->
    i ->
    ERChebPoly box b
chplPow maxDegree p n
    | n == 0 =
        chplConst 1
    | n == 1 =
        p 
    | even n =
        snd $ chplReduceDegree maxDegree $ powHalfN * powHalfN
    | odd n =
        snd $ chplReduceDegree maxDegree $ 
            p * 
            (snd $ chplReduceDegree maxDegree $
             powHalfN * powHalfN)
    where
    powHalfN =
        chplPow maxDegree p halfN
    halfN = n `div` 2
    
{-|
    Approximate the pointwise natural logarithm of a polynomial 
    by another polynomial from below and from above. 
-}
chplLog ::
    (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)
chplLog maxDegree ix p =
    error "ERChebPoly: chplLog: not implemented yet"

{-|
    Approximate the pointwise sine of a polynomial 
    by another polynomial from below and from above.
    
    Assuming the polynomial range is [-pi/2, pi/2]. 
-}
chplSine ::
    (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)
chplSine maxDegree ix p =
--        unsafePrint
--        (
--            "ERChebPoly: sineTaylor: "
--            ++ "\n p = " ++ show p
--            ++ "\n ranLargerEndpoint = " ++ show ranLargerEndpoint
--            ++ "\n sineUp = " ++ show sineUp
--            ++ "\n sineDown = " ++ show sineDown
--        ) $
        (sineDown, sineUp)
        where
        (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
            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.
    
    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: chplCosine: "
--            ++ "\n p = " ++ show p
--            ++ "\n ranLargerEndpoint = " ++ show ranLargerEndpoint
--            ++ "\n cosineUp = " ++ show cosineUp
--            ++ "\n cosineDown = " ++ show cosineDown
--        ) $
        (cosineDown, cosineUp)
        where
        (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
    
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
--                )
        ((thisCoeffDownP, thisCoeffUpP), nextDegree, (nextCoeffDown, nextCoeffUp))
    | otherwise =
--                unsafePrint
--                (
--                    "ERChebPoly: chplSine: taylorAux: "
--                    ++ "\n thisCoeff = " ++ show thisCoeff
--                    ++ "\n nextDegree = " ++ show nextDegree
--                    ++ "\n errorTermCoeff = " ++ show errorTermCoeff
--                    ++ "\n errorTermDegree = " ++ show errorTermDegree
--                )
        ((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    
    as described in [Mason & Handscomb 2003, p 62]. 
-}
chplRecip ::
    (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) => 
    Int {-^ maximum polynomial degree -} -> 
    EffortIndex {-^ minimum approx degree -} -> 
    ERChebPoly box b ->
    (ERChebPoly box b, ERChebPoly box b)
chplRecip maxDegree ix p
    | upperB < 0 = -- range negative
        (\(lo, hi) -> (-hi, -lo)) $ chplRecip maxDegree ix (negate p)
    | lowerB > 0 = -- range positive
--        unsafePrint
--        (
--            "ERChebPoly: chplRecip: "
--            ++ "\n k = " ++ show k
--            ++ "\n lowerB = " ++ show lowerB
--            ++ "\n tau = " ++ (show $ recip tauInv)
--        )
        (resDn, resUp)
    | otherwise = -- cannot establish 0 freedom
        error $
             "ERChebPoly: chplRecip: "
             ++ "cannot deal with estimated range " ++ show ranp
             ++ "of polynomial: \n" ++ show p 
    where
    ranp = ERInterval lowerB upperB
    lowerB = - (chplUpperBoundAffine ix (- p))
    upperB = chplUpperBoundAffine ix p
     
    tauDegree = effIx2int (max 2 $ ix `div` 3)
    coeffGr = effIx2gran $ ix
    
    -- translate p to have range above 1:
    k = intLogUp 2 $ ceiling (1/lowerB) -- 2^k * lowerB >= 1
    upperBtr = upperB * 2^k -- upper bound of translated poly
    (pAbove1Dn, pAbove1Up) = -- p multiplied by 2^k; range in [1,upperBtr]    
        chplScale (2^k) p
        
    -- translate T_1 to domain [0, upperBtr] and apply it to (pAbove1 - 1):
    -- T'_1 = nu * (p - 1) - 1
    trT1Dn = 
        (chplScaleDown nuLOB (pAbove1Dn - 1)) - 1
    trT1Up =
        (chplScaleUp nuHIB (pAbove1Up - 1)) - 1
    nu = recip nuInv -- auxiliary constant
    ERInterval nuLOB nuHIB = nu
    nuInv = (RA.setMinGranularity coeffGr (ERInterval upperBtr upperBtr)) / 2
    nuPlus1 = nu + 1
    nuInvPlus1 = nuInv + 1
    nuInvDiv2 = nuInv / 2
        
    -- define such translated T_i's for all i >= 0:
    trTis =
        map (mapPair (chplReduceDegreeDown maxDegree, chplReduceDegreeUp maxDegree)) $ 
            chebyEvalTsRoundDownUp trT1Dn 
        
    -- construct the result from interval coefficients:
    resDn =
        chplScaleDown (2^k) $
            (-tauAbsUpPoly) `plusDown` 
                (chplScaleUp tauAbsDnB $
                    sumDown $
                        (- errPoly) : (zipWith scaleDn cis trTis))
    resUp =
        chplScaleUp (2^k) $
            (tauAbsUpPoly) `plusUp` 
                (chplScaleUp tauAbsUpB $
                    sumUp $
                        (errPoly) : (zipWith scaleUp cis trTis))
                        
    scaleDn c (trTDn, trTUp) 
        | r >= 0 = chplScaleDown r trTDn
        | otherwise = chplScaleDown r trTUp
        where
        r = c * tauSign
    scaleUp c (trTDn, trTUp) 
        | r >= 0 = chplScaleUp r trTUp
        | otherwise = chplScaleUp r trTDn
        where
        r = c * tauSign
         
    tauAbsUpPoly = chplConst $ tauAbsUpB
    tauSign = 
        case RA.compareReals tauInv 0 of
            Just GT -> 1
            Just LT -> -1
    ERInterval tauAbsDnB tauAbsUpB = abs $ recip tauInv
    cis =
        map (\(ERInterval lo hi) -> hi) c0n 
    errPoly = chplConst err
    err =
        foldl1 plusUp $
            map (\(ERInterval lo hi) -> hi - lo) c0n
                
    -- work out the coefficients in interval arithmetic using the tau method:
    c0n = c0 : c1n
    tauInv = c0 * nuInvPlus1 + c1 * nuInvDiv2
    c0 = - c1 * nuPlus1 - c2/2
    (c1 : c2 : _) = c1n
    c1n = reverse $ take n $ csRev
    n = tauDegree
    csRev =
        cn : cnM1 : (csAux cn cnM1)
        where
        cn = 1
        cnM1 = - 2 * nuPlus1
    csAux cn cnM1 =
        cnM2 : (csAux cnM1 cnM2)
        where
        cnM2 = - cn - 2 * nuPlus1 * cnM1