module Data.Number.ER.Real.Arithmetic.Elementary
(
erSqr_R,
erSqr_IR,
erPow_R,
erPow_IR,
erSqrt_R,
erSqrt_IR,
erSqrt_IR_Inner,
erRoot_R,
erRoot_IR,
erRoot_IR_Inner,
erExp_R,
erExp_IR,
erExp_IR_Inner,
erLog_R,
erLog_IR,
erLog_IR_Inner,
erSine_R,
erSine_IR,
erSine_IR_Inner,
erCosine_R,
erCosine_IR,
erCosine_IR_Inner,
erATan_R,
erATan_IR,
erATan_IR_Inner,
erPi_R
)
where
import qualified Data.Number.ER.Real.Approx as RA
import Data.Number.ER.BasicTypes
import qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI
import Data.Number.ER.Real.Arithmetic.Taylor
import Data.Number.ER.Misc
erSqr_IR ::
(RA.ERIntApprox ira, Ord ira) =>
EffortIndex ->
ira -> ira
erSqr_IR =
RA.maxExtensionR2R
sqrExtrema
erSqr_R
where
sqrExtrema ix x
| 0 `RA.refines` x = [0]
| otherwise = []
erSqr_R ::
(RA.ERIntApprox ira, Ord ira) =>
EffortIndex ->
ira -> ira
erSqr_R ix a =
max 0 $ a' * a'
where
a' = RA.setMinGranularityOuter gran a
gran = effIx2gran ix
erPow_IR ::
(RA.ERIntApprox ira, Ord ira) =>
EffortIndex ->
Integer ->
ira -> ira
erPow_IR ix n x =
RA.maxExtensionR2R
powExtrema
(\ ix x -> erPow_R ix n x)
ix x
where
powExtrema ix x
| even n && 0 `RA.refines` x = [0]
| otherwise = []
erPow_R ::
(RA.ERIntApprox ira, Ord ira) =>
EffortIndex ->
Integer ->
ira -> ira
erPow_R ix p a
| 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.setMinGranularityOuter gran a
gran = effIx2gran ix
erSqrt_R ::
(RA.ERIntApprox ira) =>
EffortIndex -> ira -> ira
erSqrt_R = erSqrtNewton_R
erSqrt_IR ::
(RA.ERIntApprox ira) =>
EffortIndex -> ira -> ira
erSqrt_IR =
RA.maxExtensionR2R
sqrtExtrema
(\ ix x -> erSqrt_R ix x)
erSqrt_IR_Inner ::
(RA.ERIntApprox ira) =>
EffortIndex -> ira -> ira
erSqrt_IR_Inner =
RA.maxExtensionInnerR2R
sqrtExtrema
(\ ix x -> erSqrt_R ix x)
sqrtExtrema ix x = fst $ sqrtExtremaAndDirections ix x
sqrtExtremaAndDirections ix x =
case RA.compareReals 0 x of
Just LT -> ([], (Just True, Just True))
Just GT -> ([], (Nothing, Nothing))
_ -> ([0], (Nothing, Just True))
erSqrtContFr_R ::
(RA.ERIntApprox ira) =>
EffortIndex -> ira -> ira
erSqrtContFr_R ix a
| aR == 0 = 0
| aL == RA.plusInfinity = RA.plusInfinity
| aR `RA.ltSingletons` 0 = RA.topApprox
| otherwise =
contFrIter (ix + 3) $
RA.setMinGranularityOuter gran $ 0 RA.\/ aR
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) =>
EffortIndex -> ira -> ira
erSqrtNewton_R ix a
| aR == 0 = 0
| aL == RA.plusInfinity = RA.plusInfinity
| aR `RA.ltSingletons` 0 = RA.topApprox
| otherwise =
x_i RA.\/ (a/x_i)
where
gran = effIx2gran ix
(aL, aR) = RA.bounds a
aM1 = a 1
x_i =
newtonIter ((ix `div` 10) + 5) $
RA.setMinGranularityOuter gran 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
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
erRoot_IR_Inner ::
(RA.ERIntApprox ira, Ord ira) =>
EffortIndex -> Integer -> ira -> ira
erRoot_IR_Inner ix p =
RA.maxExtensionInnerR2R
(rootExtrema p)
(\ ix x -> erRoot_R ix p x) $
ix
rootExtrema p ix x = fst $ rootExtremaAndDirections p ix x
rootExtremaAndDirections p ix x
| odd p = ([], (Just True, Just True))
| otherwise =
case RA.compareReals 0 x of
Just LT -> ([], (Just True, Just True))
Just GT -> ([], (Nothing, Nothing))
_ -> ([0], (Nothing, Just True))
erRootNewton_R ::
(RA.ERIntApprox ira, Ord ira) =>
EffortIndex -> Integer -> ira -> ira
erRootNewton_R ix p a
| aR == 0 = 0
| aL == RA.plusInfinity = RA.plusInfinity
| aR < 0 && even p = RA.topApprox
| 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)
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.setMinGranularityOuter 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
erExp_R ::
(RA.ERIntApprox ira, Ord ira) =>
EffortIndex -> ira -> ira
erExp_R ix x
| RA.isBounded x =
erPow_IR ix n $
erExp_Tay_Opt_R ix xNear0
| x `RA.refines` (RA.plusInfinity) = 0
| (RA.plusInfinity) `RA.refines` x =
0 RA.\/ (erExp_R ix (snd $ RA.bounds x))
| otherwise = RA.bottomApprox
where
(xNear0, n) = scaleNear0 (x,1)
scaleNear0 (xPrev, nPrev) =
case xPrev `RA.refines` ((1) RA.\/ 1) of
True -> (xPrev, nPrev)
False -> scaleNear0 (xNext, nNext)
where
xNext = xPrev / 2
nNext = 2 * nPrev
erExp_IR ::
(RA.ERIntApprox ira, Ord ira) =>
EffortIndex -> ira -> ira
erExp_IR =
RA.maxExtensionR2R
noExtrema
erExp_R
erExp_IR_Inner ::
(RA.ERIntApprox ira, Ord ira) =>
EffortIndex -> ira -> ira
erExp_IR_Inner =
RA.maxExtensionInnerR2R
noExtrema
erExp_R
noExtrema ix x = []
erLog_R ::
(RA.ERIntApprox ira, Ord ira) =>
EffortIndex -> ira -> ira
erLog_R =
logDivSeries_R
erLog_IR ::
(RA.ERIntApprox ira, Ord ira) =>
EffortIndex -> ira -> ira
erLog_IR =
RA.maxExtensionR2R
logExtrema
(\ ix x -> logDivSeries_R ix x)
erLog_IR_Inner ::
(RA.ERIntApprox ira, Ord ira) =>
EffortIndex -> ira -> ira
erLog_IR_Inner =
RA.maxExtensionInnerR2R
logExtrema
(\ ix x -> logDivSeries_R ix x)
logExtrema ix x = fst $ logExtremaAndDirections ix x
logExtremaAndDirections ix x =
case RA.compareReals 0 x of
Just LT -> ([], (Just True, Just True))
Just GT -> ([], (Nothing, Nothing))
_ -> ([RA.plusInfinity], (Nothing, Just True))
logDivSeries_R ::
(RA.ERIntApprox ira, Ord ira) => EffortIndex -> ira -> ira
logDivSeries_R ix x
| posx `RA.refines` 0 = RA.plusInfinity
| 0 `RA.refines` posx = RA.bottomApprox
| posx `RA.refines` (RA.plusInfinity) = RA.plusInfinity
| otherwise =
case RA.compareReals posx 1 of
Just LT ->
negate $
(logDivSeries_R ix posxRecipL)
RA.\/
(logDivSeries_R ix posxRecipR)
_ ->
nearLogx + 2 * t * (series ix (RA.setMinGranularityOuter gran 1))
where
gran = effIx2gran ix
posx = (RA.setMinGranularityOuter gran x) RA./\ (0 RA.\/ (RA.plusInfinity))
(posxRecipL, posxRecipR) = RA.bounds $ recip posx
nearLogx =
0.69314718055994530941 * (fromInteger $ intLogUp 2 $ xCeiling)
remNearLogx =
posx / (erExp_R ix nearLogx)
xCeiling =
snd $ RA.integerBounds posx
t =
((remNearLogx 1) / (remNearLogx + 1))
RA./\ ((1) RA.\/ 1)
tsquare = abs $ t * t
series termsCount currentDenominator
| termsCount > 0 =
(recip currentDenominator) + tsquare * (series (termsCount 1) (currentDenominator + 2))
| otherwise =
(recip currentDenominator)
* (1 RA.\/ (recip $ 1 tsquare))
--logNewton_RA
erSine_R ::
(RA.ERIntApprox ira) =>
EffortIndex -> ira -> ira
erSine_R ix x =
case (RA.isBounded x) of
True | xNear0 `RA.refines` plusMinusPiHalf ->
erSine_Tay_Opt_R ix xNear0
True | (xNear0 piHalf) `RA.refines` plusMinusPiHalf ->
erCosine_Tay_Opt_R ix (xNear0 piHalf)
True | (xNear0 + piHalf) `RA.refines` plusMinusPiHalf ->
negate $ erCosine_Tay_Opt_R ix (xNear0 + piHalf)
True | (xNear0 pi) `RA.refines` plusMinusPiHalf ->
negate $ erSine_Tay_Opt_R ix (xNear0 pi)
_ ->
(1) RA.\/ 1
where
xNear0 = x k * 2 * pi
k = fromInteger $ toInteger kEI
(kEI,_) = RA.integerBounds $ 0.5 + (x / (2*pi))
plusMinusPiHalf = ( piHalf) RA.\/ piHalf
piHalf = pi / 2
pi = erPi_R ix
erCosine_R ::
(RA.ERIntApprox ira) =>
EffortIndex -> ira -> ira
erCosine_R ix x =
case (RA.isBounded x) of
True | xNear0 `RA.refines` plusMinusPiHalf ->
erCosine_Tay_Opt_R ix xNear0
True | (xNear0 piHalf) `RA.refines` plusMinusPiHalf ->
negate $ erSine_Tay_Opt_R ix (xNear0 piHalf)
True | (xNear0 + piHalf) `RA.refines` plusMinusPiHalf ->
erSine_Tay_Opt_R ix (xNear0 + piHalf)
True | (xNear0 pi) `RA.refines` plusMinusPiHalf ->
negate $ erCosine_Tay_Opt_R ix (xNear0 pi)
_ ->
(1) RA.\/ 1
where
xNear0 = x k * 2 * pi
k = fromInteger $ toInteger kEI
(kEI,_) = RA.integerBounds $ 0.5 + (x / (2*pi))
plusMinusPiHalf = ( piHalf) RA.\/ piHalf
piHalf = pi / 2
pi = erPi_R ix
erSine_Tay_R ::
(RA.ERIntApprox ira) =>
EffortIndex -> ira -> ira
erSine_Tay_R ix x
| (RA.plusInfinity) `RA.refines` x || (RA.plusInfinity) `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
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
erSine_IR_Inner ::
(RA.ERIntApprox ira) =>
EffortIndex -> ira -> ira
erSine_IR_Inner =
RA.maxExtensionInnerR2R sineExtremes erSine_R
erCosine_IR_Inner ::
(RA.ERIntApprox ira) =>
EffortIndex -> ira -> ira
erCosine_IR_Inner =
RA.maxExtensionInnerR2R cosineExtremes erCosine_R
sineExtremes ix x = fst $ sineExtremesAndDirections ix x
cosineExtremes ix x = fst $ cosineExtremesAndDirections ix x
sineExtremesAndDirections ix x
| RA.isBounded x =
alternatingExtremes 1 (1) ix scaledX
| otherwise = ([1,1], (Nothing, Nothing))
where
scaledX = (x / (erPi_R ix)) 0.5
cosineExtremesAndDirections ix x
| RA.isBounded x =
alternatingExtremes 1 (1) ix scaledX
| otherwise = ([1,1], (Nothing, Nothing))
where
scaledX = (x / (erPi_R ix))
alternatingExtremes extrHigh extrLow ix scaledX
| extremesCount == 1 && even minExtremeN =
([extrHigh], (Just True, Just False))
| extremesCount == 1 =
([extrLow], (Just False, Just True))
| extremesCount >= 2 =
([extrHigh,extrLow], (Just $ even minExtremeN, Just $ odd maxExtremeN))
| otherwise =
([], (Just isIncreasing, Just isIncreasing))
where
extremesCount = 1 + maxExtremeN minExtremeN
isIncreasing = even maxExtremeN
(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
erATan_R ::
(RA.ERIntApprox ira) =>
EffortIndex -> ira -> ira
erATan_R = atanEuler_R
erATan_IR ::
(RA.ERIntApprox ira) =>
EffortIndex -> ira -> ira
erATan_IR =
RA.maxExtensionR2R noExtrema erATan_R
erATan_IR_Inner ::
(RA.ERIntApprox ira) =>
EffortIndex -> ira -> ira
erATan_IR_Inner =
RA.maxExtensionInnerR2R noExtrema erATan_R
atanEuler_R ::
(RA.ERIntApprox ira) =>
EffortIndex -> ira -> ira
atanEuler_R ix x
| x `RA.refines` RA.plusInfinity = RA.plusInfinity
| x `RA.refines` ( RA.plusInfinity) = RA.plusInfinity
| not $ RA.isBounded x = RA.bottomApprox
| x `RA.refines` ((1.5) RA.\/ 1.5) =
(x / xSquarePlus1) * (series ix (RA.setMinGranularityOuter gran 2))
| otherwise =
2 * (atanEuler_R ix $ x / (1 + sqrtXQuarePlus1))
where
gran = effIx2gran ix
series termsCount coeffBase
| termsCount > 0 =
1 + xSquareOverXSquarePlus1 * coeff * (series (termsCount 1) (coeffBase + 2))
| otherwise =
1 + xSquare * (0 RA.\/ 1)
where
coeff = coeffBase / (coeffBase + 1)
xSquare = abs $ x * x
xSquarePlus1 = xSquare + 1
xSquareOverXSquarePlus1 = xSquare / xSquarePlus1
sqrtXQuarePlus1 =
iterateIx 10 EI.MinusInfinity
where
iterateIx ix prevPrec
| prevPrec == currentPrec = result
| otherwise =
iterateIx (ix * 2) currentPrec
where
result = erSqrt_R ix xSquarePlus1
currentPrec = RA.getPrecision result
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
erPi_R ::
(RA.ERIntApprox ira) =>
EffortIndex -> ira
erPi_R = piBellard_R
piAtan_R ::
(RA.ERIntApprox ira) =>
EffortIndex -> ira
piAtan_R ix =
(*) 4 $ atanEuler_R ix 1
piBellard_R ::
(RA.ERIntApprox ira) =>
EffortIndex -> ira
piBellard_R ix =
r1over64 * (sum $ reverse $ bellardTerms 0 (10 + (ix `div` 10)) (1,z,z))
where
gran = max 0 (effIx2gran ix) + 10
r1over64 = (RA.setMinGranularityOuter gran 1) / 64
r1over1024 = (RA.setMinGranularityOuter gran 1) / 1024
z = RA.setMinGranularityOuter 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 =
(recip $ r10n + 9)
(recip $ r4n + 3)
4 * ((recip $ r10n + 7) + (recip $ r10n + 5))
(64 / (r10n + 3))
(32 / (r4n + 1))
+ (256 / (r10n + 1))