{-|
    Module      :  Data.Number.ER.Real.Arithmetic.Taylor
    Description :  implementation of Taylor series
    Copyright   :  (c) Amin Farjudian, Michal Konecny
    License     :  BSD3

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

    Taylor series related functions.
-}
module Data.Number.ER.Real.Arithmetic.Taylor where

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


erTaylor_R
    :: (RA.ERIntApprox ira)
    => EffortIndex
    -> (Int -> ira) -- ^ coefficients of the Taylor series
    -> (Int -> ira) -- ^ function to estimate the n'th derivative between a and x
    -> ira -- ^ centre of the Taylor Expansion
    -> ira 
    -> ira
erTaylor_R ix coefSeq derivBounds a x =
    erTaylor_R_FullArgs coefSeq derivBounds n a gran x
    where
    n = fromInteger ix
    gran = fromInteger $ toInteger $ ix

erTaylor_R_FullArgs
    :: (RA.ERIntApprox ira)
    => (Int -> ira)  -- ^ coefficients of the Taylor series
    -> (Int -> ira) -- ^ function to estimate the n'th derivative between a and x
    -> Int -- ^ use this many elements of the series (+ account for error appropriately)
    -> ira -- ^ centre of the Taylor Expansion
    -> Granularity -- ^ make all constants have this granularity, thus influencing rounding errors
    -> ira 
    -> ira
erTaylor_R_FullArgs coefSeq derivBounds n a gran x = 
    rec_apTaylor (RA.setMinGranularity gran 0) 0
    where
    rec_apTaylor i j
        | n > j = (coefSeq(j)) + 
                        ((x - a)/(i+1)) * (rec_apTaylor (i+1) (j+1))
        | n == j = derivBounds n
        | otherwise = 
            error "Data.Number.ER.Real.Arithmetic.Taylor.hs: erTaylor_RA_FullArgs: The index n cannot be negative"

{-|
    A Taylor series for exponentiation.    
-}
erExp_Tay_Opt_R
    :: (RA.ERIntApprox ira)
    => EffortIndex
    -> ira
    -> ira
erExp_Tay_Opt_R ix x 
    | RA.isEmpty x = RA.emptyApprox
    | x `RA.refines` (-1/0) = 0 -- -infty is not handled well by the Taylor formula
    | otherwise = 1 + (te ix x (RA.setMinGranularity gran 1))
    where
    gran = effIx2gran ix
    te steps x i
        | steps > 0 =
            (x/i) * (1 + (te (steps - 1) x (i + 1)))
        | steps == 0 = 
            errorBound
            where
            errorBound = 
                (x/i) * ithDerivBound
            ithDerivBound 
                | xCeiling == EI.MinusInfinity = -- certainly -infty:
                    0
                | xCeiling < 0 = -- certainly negative:
                    pow26xFloor RA.\/ 1
                | xFloor > 0 = -- certainly positive:
                    1 RA.\/ pow28xCeiling
                | otherwise = -- could contain 0:
                    pow26xFloor RA.\/ pow28xCeiling
                where
                (xFloor, xCeiling) = RA.integerBounds x
                pow26xFloor 
                    | xFloor == EI.MinusInfinity =
                        0
                    | otherwise = 
                        ((RA.setMinGranularity gran 26)/10) ^^ xFloor 
                            -- lower estimate of e^x
                pow28xCeiling 
                    | xCeiling == EI.PlusInfinity =
                        (1/0)
                    | otherwise = 
                        ((RA.setMinGranularity gran 28)/10) ^^ xCeiling 
                            -- upper estimate of e^x

{-
 The sine and cosine are implemented in almost exactly the same way 
-}

{-|
    A Taylor series for sine.    
-}
erSine_Tay_Opt_R
    :: (RA.ERIntApprox ira)
    => EffortIndex
    -> ira
    -> ira
erSine_Tay_Opt_R ix x = taylor_seg ix x (RA.setMinGranularity gran 1)
	where
		gran = effIx2gran ix
		taylor_seg i x n -- 'i' for iterator
			| i > 0  = x - ((x*x)/((n+1)*(n+2))) * (taylor_seg (i-2) x (n+2))
			| otherwise = errorRegion
				where 
					errorRegion = (-1) RA.\/ (1)
		
{-|
    A Taylor series for cosine.    
-}
erCosine_Tay_Opt_R 
    :: (RA.ERIntApprox ira) 
    => EffortIndex 
    -> ira
    -> ira
erCosine_Tay_Opt_R ix x = taylor_seg ix x (RA.setMinGranularity gran 1)
	where
		gran = effIx2gran ix
		taylor_seg i x n -- 'i' for iterator
			| i > 0  = 1 - ((x*x)/(n*(n+1))) * (taylor_seg (i-2) x (n+2))
			| otherwise = errorRegion
				where 
					errorRegion = (-1) RA.\/ (1)

   
				
{-| Natural Logarithm: The following is a code for obtaining natural
 	logarithm using taylor series. Note that it only works for 
 	x in [ 1, 2]. For other values, a scaling by factors of e^q is
 	best to do, i.e. if x is not in [1,2], then find some ratioal number				
 	q where exp(q) * x is in [1,2]. Then you have:
 	log ( exp(q) * m)  = q + log(m)
-}

{-| Coefficients of the taylor series around a=1 -}
--logTayCoefs
--	:: (RA.ERIntApprox ira)
--	=>	Int -- up to how many terms of the Taylor series is desired
--	-> Int
--    -> ra
--	
--logTayCoefs n i 
----	| i < 0 = error "ERTaylor.logTayCoefs: Negative n for the n-th term of Taylor series for logarithm"
--	| i == 0 = 0
--	| i == n = fromInteger $ toInteger $ amFact(n-1)	
--	| otherwise = fromInteger $ toInteger $ ((-1)^(i-1) * amFact(i-1))
--	where
--		amFact (m) = product [1..m]
		
--logTay_RA
--	:: (RA.ERIntApprox ira)
--	=> EffortIndex
--    -> ra
--    -> ra
--	
--logTay_RA i = erTaylor_RA_FullArgs (logTayCoefs $fromInteger $ toInteger $ i) 
--                (100000) (setMinGranularity (effIx2gran i) 1) (effIx2gran i)
--	
--	
--logTay 
--	:: (RA.ERIntApprox ira)	
--	=> (ConvergRealSeq ra)
--    -> (ConvergRealSeq ra)
--logTay = convertFuncRA2Seq logTay_RA