{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE UndecidableInstances #-}
{-|
    Module      :  Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Field
    Description :  (internal) field operations 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 field arithmetic over polynomials 
    with rounding consistent over the whole domain.
-}
module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Field 

where

import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic

import qualified Data.Number.ER.Real.Base as B
import qualified Data.Number.ER.Real.DomainBox as DBox
import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainIntBox)
import Data.Number.ER.Misc

import qualified Data.Map as Map

chplAffine ::
    (B.ERRealBase b, DomainBox box varid Int, Ord box) =>
    b -> 
    Map.Map varid b ->
    ERChebPoly box b
chplAffine at0 varCoeffs =
    ERChebPoly $ 
        Map.insert chplConstTermKey at0 $
            Map.mapKeys (\ i -> DBox.singleton i 1) varCoeffs
    
{-|
    Convert a polynomial to a lower-order one that is dominated by (resp. dominates)
    it closely on the domain [-1,1].
-}
chplReduceDegree ::
    (B.ERRealBase b, DomainBox box varid Int, Ord box) => 
    Int {-^ new maximal order -} ->
    ERChebPoly box b -> 
    (ERChebPoly box b, ERChebPoly box b) {-^ lower and upper bounds with limited degree -}
chplReduceDegree maxOrder (ERChebPoly coeffs) =
    (ERChebPoly newCoeffsDown, ERChebPoly newCoeffsUp)
--    errorModule "chplSetMaxOrder: not implemented yet"
    where
    newCoeffsUp =
        Map.insertWith plusUp chplConstTermKey highOrderCompensation coeffsLowOrder
    newCoeffsDown =
        Map.insertWith plusDown chplConstTermKey (-highOrderCompensation) coeffsLowOrder
    highOrderCompensation =
        Map.fold (\ new prev -> prev + (abs new)) 0 coeffsHighOrder
    (coeffsHighOrder, coeffsLowOrder) =        
        Map.partitionWithKey (\ k v -> chplTermOrder k > maxOrder) coeffs

chplReduceDegreeDown m = fst . chplReduceDegree m
chplReduceDegreeUp m = snd . chplReduceDegree m

instance (B.ERRealBase b, DomainBox box varid Int, Ord box) => Num (ERChebPoly box b)
    where
    fromInteger n =
        ERChebPoly $ Map.singleton chplConstTermKey (fromInteger n)
    abs (ERChebPoly coeffs) =
        errorModule "abs of a polynomial not implemented, use UFB.max instead"
    signum (ERChebPoly coeffs) =
        errorModule "signum of a polynomial not implemented, use RA.leqReals instead"
    --------- negation ----------
    negate (ERChebPoly coeffs) =
        ERChebPoly $ Map.map negate coeffs
    --------- addition ----------
    (ERChebPoly coeffs1) + (ERChebPoly coeffs2) =
        ERChebPoly sumCoeffs
        where
        sumCoeffs =
            Map.insertWith (+) chplConstTermKey maxError coeffsDown
            -- point-wise sum of polynomials with coeff rounding errors
            -- compensated for by enlarging the constant term
        coeffsUp =
            (Map.unionWith (+) coeffs1 coeffs2)
            -- point-wise sum of polynomials with coeffs rounded upwards
        coeffsDown =
            (Map.unionWith (\c1 c2 -> - ((- c1) + (- c2))) coeffs1 coeffs2)
            -- point-wise sum of polynomials with coeffs rounded upwards
        maxError =
            Map.fold (+) 0 $ 
                Map.intersectionWith (-) coeffsUp coeffsDown
            -- addition must round upwards on interval [-1,1]
                    -- non-constant terms are multiplied by quantities in [-1,1] 
                    -- and thus can make the result drop below the exact result
                    -- -> to compensate add the rounding difference to the constant term 
    --------- multiplication ----------
    (ERChebPoly coeffs1) * (ERChebPoly coeffs2) =
        ERChebPoly prodCoeffs
        where        
        prodCoeffs =
            Map.insertWith (+) chplConstTermKey roundOffCompensation $ 
                Map.map negate directProdCoeffsDown
        roundOffCompensation =
            Map.fold (+) 0 $
                Map.unionWith (+) directProdCoeffsDown directProdCoeffsUp
        (directProdCoeffsUp, directProdCoeffsDown) =
            foldl addCombiCoeff (Map.empty, Map.empty) combinedCoeffs
            where
            addCombiCoeff
                    (prevCoeffsUp, prevCoeffsDown) 
                    (coeffUp, coeffDown, (powersList, coeffCount)) =
                foldl addOnce (prevCoeffsUp, prevCoeffsDown) powersList
                where
                addOnce (prevCoeffsUp, prevCoeffsDown) powers =
                    (Map.insertWith (+) powers coeffUpFrac prevCoeffsUp, 
                     Map.insertWith (+) powers coeffDownFrac prevCoeffsDown)
                coeffUpFrac = coeffUp / coeffCountB
                coeffDownFrac = coeffDown / coeffCountB
                coeffCountB = fromInteger coeffCount
        combinedCoeffs =
            [   -- (list of triples)
                (
                    (c1 * c2) -- upwards rounded product
                ,
                    ((- c1) * c2) -- downwards rounded negated product
                ,
                    combinePowers powers1 powers2
                )
            |
                (powers1, c1) <- coeffs1List,
                (powers2, c2) <- coeffs2List
            ]
        combinePowers powers1 powers2 =
            (combinedPowers, 2 ^ (length sumsDiffs)) 
            where
            combinedPowers =
                map (DBox.fromAscList . (filter $ \ (k,v) -> v > 0)) $
                    allPairsCombinations $ 
                        sumsDiffs
            sumsDiffs = 
                -- associative list with the sum and difference of powers for each variable
                zipWith (\(k,s) (_,d) -> (k,(s,d)))
                    (DBox.toAscList $ DBox.unionWith (\a b -> (a + b)) powers1 powers2)
                    (DBox.toAscList $ DBox.unionWith (\a b -> abs (a - b)) powers1 powers2)
        coeffs1List =
            Map.toList coeffs1
        coeffs2List =
            Map.toList coeffs2


-- | multiply a polynomial by a scalar rounding downwards and upwards 
chplScale ::
    (B.ERRealBase b, DomainBox box varid Int, Ord box) =>
    b -> 
    (ERChebPoly box b) -> 
    (ERChebPoly box b, ERChebPoly box b)
chplScale ratio (ERChebPoly coeffs) =
    (ERChebPoly coeffsDown, ERChebPoly coeffsUp)
    where
    coeffsDown = 
        Map.insertWith plusDown chplConstTermKey (- errBound) coeffsScaled 
    coeffsUp = 
        Map.insertWith plusUp chplConstTermKey errBound coeffsScaled
    (errBound, coeffsScaled) =
        Map.mapAccum processTerm 0 coeffs
    processTerm errBoundPrev coeff =
        (errBoundPrev + errBoundHere, coeffScaledUp)
        where
        errBoundHere = coeffScaledUp - coeffScaledDown
        coeffScaledDown = ratio `timesDown` coeff
        coeffScaledUp = ratio `timesUp` coeff    

chplScaleDown r = fst . chplScale r
chplScaleUp r = snd . chplScale r

-- | multiply a polynomial by a scalar interval rounding downwards and upwards 
chplScaleApprox :: 
    (B.ERRealBase b, DomainBox box varid Int, Ord box) =>
    (b, b) -> 
    (ERChebPoly box b) -> 
    (ERChebPoly box b, ERChebPoly box b)
chplScaleApprox (ratioDown, ratioUp) (ERChebPoly coeffs) =
    (ERChebPoly coeffsDown, ERChebPoly coeffsUp)
    where
    coeffsDown =
        Map.insertWith plusDown chplConstTermKey (- errBound) coeffsScaled 
    coeffsUp = 
        Map.insertWith plusUp chplConstTermKey errBound coeffsScaled
    (errBound, coeffsScaled) =
        Map.mapAccum processTerm 0 coeffs
    processTerm errBoundPrev coeff =
        (errBoundPrev + errBoundHere, coeffScaledUp)
        where
        errBoundHere = coeffScaledUp - coeffScaledDown
        (coeffScaledDown, coeffScaledUp)
            | coeff >= 0 = 
                (ratioDown `timesDown` coeff, ratioUp `timesUp` coeff)
            | coeff < 0 = 
                (ratioUp `timesDown` coeff, ratioDown `timesUp` coeff)


instance (B.ERRealBase b, DomainBox box varid Int, Ord box) => Fractional (ERChebPoly box b)
    where
    fromRational r =
        ERChebPoly $ Map.singleton chplConstTermKey (fromRational r)
    --------- division ----------
    _ / _ =
        errorModule "for division use chplRecip from module Elementary"    
    
instance (B.ERRealBase b, DomainBox box varid Int, Ord box) => Ord (ERChebPoly box b)
    where
    compare _ _ =
        errorModule "cannot compare polynomials, consider using RA.leqReals instead"
    
--instance (B.ERRealBase b, DomainBox box varid Int, Ord box) => Real (ERChebPoly box b)
--    where
--    toRational _ =
--        errorModule "toRational: cannot convert polynomial to rational"    
--    
--instance (B.ERRealBase b, DomainBox box varid Int, Ord box) => RealFrac (ERChebPoly box b)
--    where
--    properFraction _ =
--        errorModule "properFraction: rounding of polynomials not implemented"