{-|
    Module      :  Data.Number.ER.Real.Arithmetic.Elementary
    Description :  some elementary functions
    Copyright   :  (c) Michal Konecny, Amin Farjudian, Jan Duracz
    License     :  BSD3

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

    Some important elementary functions for real approximations
    and their maximal extensions for interval approximations.
-}
module Data.Number.ER.Real.Arithmetic.Elementary
(   
    -- * specialised exponentiation
    erSqr_R,
    erSqr_IR,
    erPow_R,
    erPow_IR,
    erSqrt_R,
    erSqrt_IR,
    erRoot_R,
    erRoot_IR,
    -- * exponentiation and logarithm 
    erExp_R,
    erExp_IR,
    erLog_R,
    erLog_IR,
    -- * trigonometrics
    erSine_R,
    erSine_IR,
    erCosine_R,
    erCosine_IR,
    erATan_R,
    erATan_IR,
    erPi_R
)
where

import qualified Data.Number.ER.Real.Approx as RA
import Data.Number.ER.BasicTypes

import Data.Number.ER.Real.Arithmetic.Taylor
-- import Data.Number.ER.Real.Arithmetic.Newton

import Data.Number.ER.Misc

{-
    sqr
-}

erSqr_IR ::
    (RA.ERIntApprox ira, Ord ira) =>
    EffortIndex -> 
    ira -> ira
erSqr_IR = erSqr_R

erSqr_R ::
    (RA.ERIntApprox ira, Ord ira) =>
    EffortIndex -> 
    ira -> ira
erSqr_R ix a
    | RA.isEmpty a =
        RA.emptyApprox
    | otherwise = 
        max 0 $ a' * a'
    where
    a' = RA.setMinGranularity gran a
    gran = effIx2gran ix
    
{-
    integer exponentiation x ^ p
-}

erPow_IR ::
    (RA.ERIntApprox ira, Ord ira) =>
    EffortIndex -> 
    Integer ->
    ira -> ira
erPow_IR = erPow_R

erPow_R ::
    (RA.ERIntApprox ira, Ord ira) =>
    EffortIndex ->
    Integer ->
    ira -> ira
erPow_R ix p a
    | RA.isEmpty a =
        RA.emptyApprox
    | p < 0 =
        1 / erPow_R ix (-p) a
    | p == 0 = 
        1
    | even p =
        erPow_R ix (div p 2) (erSqr_R ix a')
    | otherwise =
        a' * (erPow_R ix (div (p - 1) 2) (erSqr_R ix a'))
    where
    a' = RA.setMinGranularity gran a
    gran = effIx2gran ix

{-
    sqrt
-}

erSqrt_R ::
    (RA.ERIntApprox ira, Ord ira) => 
    EffortIndex -> ira -> ira
erSqrt_R = erSqrtNewton_R  
    
erSqrt_IR ::
    (RA.ERIntApprox ira, Ord ira) => 
    EffortIndex -> ira -> ira
erSqrt_IR =
    RA.maxExtensionR2R 
        sqrtExtrema
        (\ ix x -> erSqrt_R ix x)

sqrtExtrema ix x
    | 0 `RA.refines` x = [0]
    | otherwise = []
        


erSqrtContFr_R ::
    (RA.ERIntApprox ira, Ord ira) => 
    EffortIndex -> ira -> ira
erSqrtContFr_R ix a
    | aR == 0 = 0
    | aL == 1/0 = 1/0
    | aR `RA.ltSingletons` 0 = RA.emptyApprox
    | otherwise =
        contFrIter (ix + 3) $
            RA.setMinGranularity gran $ max 0 (0 RA.\/ a) 
    where
    gran = effIx2gran ix
    (aL, aR) = RA.bounds a
    aM1 = a - 1
    
    contFrIter i x_i
        | i == 0 =
            x_i
        | otherwise =
            1 + (aM1 / (x_iPlus1 + 1))
        where
        x_iPlus1 = contFrIter (i - 1) x_i
            
erSqrtNewton_R ::
    (RA.ERIntApprox ira, Ord ira) => 
    EffortIndex -> ira -> ira
erSqrtNewton_R ix a
    | RA.isEmpty a = RA.emptyApprox
    | aR == 0 = 0
    | aL == 1/0 = 1/0
    | aR < 0 = RA.emptyApprox
    | otherwise =        
        x_i RA.\/ (a/x_i)
    where
    gran = effIx2gran ix
    (aL, aR) = RA.bounds a
    aM1 = a - 1
    
    x_i = 
        newtonIter ((ix `div` 100) + 5) $
                RA.setMinGranularity gran $ max 0 aR 
    newtonIter i x_i
        | i == 0 = x_i
        | otherwise =
                snd $ RA.bounds $
                    (x_iMinus1 + a / (x_iMinus1)) / 2
        where
        x_iMinus1 = newtonIter (i - 1) x_i

{-
    pth root x ^ (1/p)
-}

erRoot_R ::
    (RA.ERIntApprox ira, Ord ira) => 
    EffortIndex -> Integer -> ira -> ira
erRoot_R = erRootNewton_R    
    
erRoot_IR ::
    (RA.ERIntApprox ira, Ord ira) => 
    EffortIndex -> Integer -> ira -> ira
erRoot_IR ix p =
    RA.maxExtensionR2R 
        (rootExtrema p)
        (\ ix x -> erRoot_R ix p x) $
            ix

rootExtrema p ix x
    | 0 `RA.refines` x && even p = [0]
    | otherwise = []

erRootNewton_R ::
    (RA.ERIntApprox ira, Ord ira) => 
    EffortIndex -> Integer -> ira -> ira
erRootNewton_R ix p a
    | RA.isEmpty a = RA.emptyApprox
    | aR == 0 = 0
    | aL == 1/0 = 1/0
    | aR < 0 && even p = RA.emptyApprox
    | aR < 0 = - erRootNewton_R ix p (-a)
    | p > 0 =
        x_i RA.\/ (a/x_i_pow_p_minus_1)
    | otherwise =   
        1 / (erRootNewton_R ix (-p) a) -- TODO: check extremes
    where
    gran = effIx2gran ix
    (aL, aR) = RA.bounds a
    aM1 = a - 1
    pIRA = fromInteger p
    pIRA_minus_1 = pIRA - 1
    
    (x_i, x_i_pow_p_minus_1) = 
        newtonIter (ix + 5) $
                RA.setMinGranularity gran $ max 0 aR
 
    newtonIter i x_0
        | i == 0 = 
            (x_0, x_0_pow_p_minus_1)
        | otherwise =
            (x_i, x_i_pow_p_minus_1)
            
        where
        (x_iMinus1, x_iMinus1_pow_p_minus_1) = 
            newtonIter (i - 1) x_0
        x_i =
                snd $ RA.bounds $
                    (pIRA_minus_1 * x_iMinus1 + a / x_iMinus1_pow_p_minus_1) / pIRA
        x_i_pow_p_minus_1 =
                erPow_R ix (p - 1) x_i
        x_0_pow_p_minus_1 =
                erPow_R ix (p - 1) x_0

{-
    e^x and log
-}

erExp_R :: 
    (RA.ERIntApprox ira) => 
    EffortIndex -> ira -> ira
    
erExp_R = erExp_Tay_Opt_R

{- 
    exp as derived from Taylor series is already a maximal extension
    because it does not suffer from the wrapping effect - all
    functions used are monotone - all Taylor coeffs are positive
-}
erExp_IR :: 
    (RA.ERIntApprox ira) => 
    EffortIndex -> ira -> ira
    
erExp_IR ix x
    | 0 `RA.refines` x || (-1/0) `RA.refines` x=
        RA.maxExtensionR2R
            (\ ix x -> [])
            (\ ix x -> erExp_R ix x)
            ix x
    | otherwise =
        erExp_R ix x


{- Log using Newton -}

erLog_R :: 
    (RA.ERIntApprox ira) => 
    EffortIndex -> ira -> ira
    
erLog_R =
    logDivSeries_R 
--    erLog_IR -- intervals are more efficient for log than singletons 

erLog_IR ::
    (RA.ERIntApprox ira) => 
    EffortIndex -> ira -> ira
    
erLog_IR =
    RA.maxExtensionR2R
        logExtrema
        (\ ix x -> logDivSeries_R ix x)
        
logExtrema ix x
    | 0 `RA.refines` x = [-1/0]
    | otherwise = []
        
{-| log using a fast converging series, designed to be used with singletons -}
logDivSeries_R ::
    (RA.ERIntApprox ira) => EffortIndex -> ira -> ira 
logDivSeries_R ix x 
    | RA.isEmpty posx = RA.emptyApprox
    | posx `RA.refines` 0 = -1/0 
    | posx `RA.refines` (1/0) = 1/0
    | otherwise =
        case RA.compareReals posx 1 of
            Just LT ->
                - (logDivSeries_R ix (recip posx))
            _ ->
                nearLogx + 2 * t * (series ix (RA.setMinGranularity gran 1))
    where
    gran = effIx2gran ix
    posx = (RA.setMinGranularity gran x) RA./\ (0 RA.\/ (1/0))
    nearLogx =
        0.69314718055994530941 * (fromInteger $ intLogUp 2 $ xCeiling)
    remNearLogx =
        posx / (erExp_R ix nearLogx) -- should be very close to 1
    xCeiling = 
        snd $ RA.integerBounds posx
    t = 
        ((remNearLogx - 1) / (remNearLogx + 1)) -- the range of this expression is [-1,1] 
            RA./\ ((-1) RA.\/ 1) -- correction of wrapping 
    tsqare = abs $ t * t -- the range is [0,1]
    series termsCount currentDenominator 
        | termsCount > 0 =
            (recip currentDenominator) + tsqare * (series (termsCount - 1) (currentDenominator + 2))
        | otherwise =
            (recip currentDenominator)
            * (1 RA.\/ (recip $ 1 - tsqare)) -- [1,1/(1-t^2)] is a valid error bound
        
--{- log using Newton -}
--    
--logNewton_RA
--    :: (RA.ERIntApprox ira)
--    => EffortIndex
--    -> ra -- must not be below 1
--    -> ra
--    
--logNewton_RA i x = 
--    case compareReals posx 1 of
--        Just LT ->
--            - (logNewton_RA i (recip posx))
--        _ ->    
--            erNewton_FullArgs 
--                ( \ i y -> (erExp_RA i y) - posx, erExp_RA) 
--                (RA.setMinGranularity gran nearLogx) 
--                (RA.setMinGranularity gran 1) 
--                (fromInteger $ toInteger i)
--                i
--    where
--    gran = effIx2gran i
--    posx = 
--        RA.setMinGranularity gran x /\ (ira2ra $ 0 RA.\/ (1/0))
--    nearLogx =                    
--        0.69314718055994530941 * (fromInteger $ intLog 2 $ xCeiling)
--    xCeiling 
--        | RA.isEmpty posx = 1 -- choice of constant irrelevant
--        | otherwise =
--            snd $ RA.iraIntegerBounds $ ra2ira posx


{-
    sin(x) and cos(x)
-}

erSine_R :: 
    (RA.ERIntApprox ira) => 
    EffortIndex -> ira -> ira

erSine_R ix x
    | (1/0) `RA.refines` x || (-1/0) `RA.refines` x = 
        (-1) RA.\/ 1  
    | otherwise =
        erSine_Tay_Opt_R ix x

erCosine_R :: 
    (RA.ERIntApprox ira) => 
    EffortIndex -> ira -> ira
     
erCosine_R ix x
    | (1/0) `RA.refines` x || (-1/0) `RA.refines` x =   
        (-1) RA.\/ 1  
    | otherwise =
        erCosine_Tay_Opt_R ix x 


{- Sine using generic Taylor (see Taylor for an optimised version) -}

erSine_Tay_R :: 
    (RA.ERIntApprox ira) =>
    EffortIndex -> ira -> ira

erSine_Tay_R ix x
    | (1/0) `RA.refines` x || (-1/0) `RA.refines` x = 
        (-1) RA.\/ 1  
    | otherwise =
        erTaylor_R ix sine_coefSeq sine_error 0 x

sine_coefSeq :: 
    (RA.ERIntApprox ira) => 
    Int -> ira

sine_coefSeq n
  | n `mod` 4 == 0 = 0
  | n `mod` 4 == 1 = 1
  | n `mod` 4 == 2 = 0
  | n `mod` 4 == 3 = -1
  
sine_error n = (-1) RA.\/ 1  

{- maximal extensions -}

erSine_IR ::
    (RA.ERIntApprox ira) =>
    EffortIndex -> ira -> ira 
    
erSine_IR = 
    RA.maxExtensionR2R sineExtremes erSine_R
    
erCosine_IR ::
    (RA.ERIntApprox ira) =>
    EffortIndex -> ira -> ira 
    
erCosine_IR = 
    RA.maxExtensionR2R cosineExtremes erCosine_R
        
sineExtremes ix x 
    | RA.isBounded x =
        alternatingExtremes 1 (-1) ix scaledX
    | otherwise = [-1,1]
    where
    scaledX = (x / (erPi_R ix)) - 0.5
    
cosineExtremes ix x
    | RA.isBounded x =
        alternatingExtremes 1 (-1) ix scaledX
    | otherwise = [-1,1]
    where
    scaledX = (x / (erPi_R ix))
    
alternatingExtremes extr0 extr1 ix scaledX
    | extremesCount >= 2 = [extr0,extr1]  
    | extremesCount == 1 && even minExtremeN = [extr0]
    | extremesCount == 1 = [extr1]
    | otherwise = []
    where
    extremesCount = 1 + maxExtremeN - minExtremeN
    (xFloor, xCeiling) = RA.integerBounds scaledX
    minExtremeN = 
        case RA.compareReals (fromInteger $ toInteger xFloor) scaledX of
            Just LT -> (xFloor + 1)
            _ -> xFloor
    maxExtremeN =
        case RA.compareReals scaledX (fromInteger $ toInteger xCeiling) of
            Just LT -> xCeiling - 1
            _ -> xCeiling
        

{-
    tan(x), atan(x) and pi
-}

erATan_R :: 
    (RA.ERIntApprox ira) => 
    EffortIndex -> ira -> ira
    
erATan_R = atanEuler_R

erATan_IR ::
    (RA.ERIntApprox ira) =>
    EffortIndex -> ira -> ira 
    
erATan_IR =
    RA.maxExtensionR2R atanExtremes erATan_R

atanExtremes ix x = []

{- 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]
-}

atanEuler_R ::
    (RA.ERIntApprox ira) => 
    EffortIndex -> ira -> ira

atanEuler_R ix x 
    | RA.isEmpty x = RA.emptyApprox
    | otherwise =
        (x / xSquarePlus1) * (series ix (RA.setMinGranularity gran 2))
    where
    gran = effIx2gran ix
    series termsCount coeffBase 
        | termsCount > 0 =
            1 + xSquareOverXSquarePlus1 * coeff * (series (termsCount - 1) (coeffBase + 2))
        | otherwise =
            1 + xSquare * (coeff RA.\/ 1)
        where
        coeff = coeffBase / (coeffBase + 1)
    xSquare = abs $ x * x
    xSquarePlus1 = xSquare + 1
    xSquareOverXSquarePlus1 = xSquare / xSquarePlus1
    
--{- atan using Newton -}
--
--atanNewton_RA :: 
--    (RA.ERIntApprox ira) => 
--    EffortIndex -> ra -> ra
--    
--atanNewton_RA i x = 
--    erNewton_FullArgs 
--        ( \ i y -> (erTan_RA i y) - x, erTanDeriv_RA) 
--        (RA.setMinGranularity (effIx2gran i) (x))
--        (RA.setMinGranularity (effIx2gran i) 1) 
--        (fromInteger $ toInteger i)
--        i

{- tan -}

erTan_R :: 
    (RA.ERIntApprox ira) => 
    EffortIndex -> ira -> ira
    
erTan_R ix x =
    (erSine_R ix x) / (erCosine_R ix x)

erTanDeriv_R ix x = 
    recip $ abs $ cosx * cosx
    where
    cosx = erCosine_R ix x


{- pi -}

{-|
    pi using Bellard's formula
    
    Convergence properties:
    
    * shrinking sequence
     
    * rate at least 2^(-i).
    
-}
erPi_R :: 
    (RA.ERIntApprox ira) => 
    EffortIndex -> ira
erPi_R = piBellard_R

-- | pi using atan 
piAtan_R ::
    (RA.ERIntApprox ira) => 
    EffortIndex -> ira
piAtan_R ix =
    (*) 4 $ atanEuler_R ix 1

{-|
    pi using Bellard's formula
    (see <http://en.wikipedia.org/wiki/Computing_π>)
    
    Convergence properties:
    
    * shrinking sequence
     
    * rate at least 2^(-i).
    
-}
piBellard_R ::
    (RA.ERIntApprox ira) => 
    EffortIndex -> ira
piBellard_R ix =
    r1over64 * (sum $ reverse $ bellardTerms 0 (10 + (ix `div` 10)) (1,z,z))
    {- 
      sum from the smallest to the largest 
      (got this trick from Martin Escardo who said he got it from Andrej Bauer)
      
      the rounding error dominates the truncation error to such
      a degree that the truncation error can be safely left out
      
      each bellard term contributes 10 binary digits that the following terms
      do not influence
    -} 
    where
    gran = max 0 (effIx2gran ix) + 10
    r1over64 = (RA.setMinGranularity gran 1) / 64
    r1over1024 = (RA.setMinGranularity gran 1) / 1024
    z = RA.setMinGranularity gran 0
    bellardTerms n nMax (mult, r4n, r10n)
        | n >= nMax = []
        | otherwise =
             termN : rest
        where
        rest = 
            bellardTerms (n + 1) nMax (- mult * r1over1024, r4n + 4, r10n + 10)
        termN = 
            mult * bellardSum
        bellardSum =
            -- sum from the smallest to the largest
            (recip $ r10n + 9)
            - (recip $ r4n + 3)
            - 4 * ((recip $ r10n + 7) + (recip $ r10n + 5))
            - (64 / (r10n + 3))
            - (32 / (r4n + 1))
            + (256 / (r10n + 1))