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

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))

```