AERN-Real-0.9.5: src/Data/Number/ER/Real/Arithmetic/Elementary.hs
{-|
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 $ intLog 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))