module Data.Number.ER.Real.Arithmetic.Elementary
(
erSqr_R,
erSqr_IR,
erPow_R,
erPow_IR,
erSqrt_R,
erSqrt_IR,
erRoot_R,
erRoot_IR,
erExp_R,
erExp_IR,
erLog_R,
erLog_IR,
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.Misc
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
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
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
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)
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
erExp_R ::
(RA.ERIntApprox ira) =>
EffortIndex -> ira -> ira
erExp_R = erExp_Tay_Opt_R
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
erLog_R ::
(RA.ERIntApprox ira) =>
EffortIndex -> ira -> ira
erLog_R =
logDivSeries_R
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 = []
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)
xCeiling =
snd $ RA.integerBounds posx
t =
((remNearLogx 1) / (remNearLogx + 1))
RA./\ ((1) RA.\/ 1)
tsqare = abs $ t * t
series termsCount currentDenominator
| termsCount > 0 =
(recip currentDenominator) + tsqare * (series (termsCount 1) (currentDenominator + 2))
| otherwise =
(recip currentDenominator)
* (1 RA.\/ (recip $ 1 tsqare))
--logNewton_RA
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
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
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
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 = []
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
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.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 =
(recip $ r10n + 9)
(recip $ r4n + 3)
4 * ((recip $ r10n + 7) + (recip $ r10n + 5))
(64 / (r10n + 3))
(32 / (r4n + 1))
+ (256 / (r10n + 1))