{-| 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) => EffortIndex -> ira -> ira erSqr_IR = erSqr_R erSqr_R :: (RA.ERIntApprox 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) => EffortIndex -> Integer -> ira -> ira erPow_IR = erPow_R erPow_R :: (RA.ERIntApprox 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) => 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) sqrtExtrema ix x | 0 `RA.refines` x = [0] | otherwise = [] erSqrtContFr_R :: (RA.ERIntApprox ira) => EffortIndex -> ira -> ira erSqrtContFr_R ix a | aR == 0 = 0 | aL == 1/0 = 1/0 | aR < 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) => 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) => EffortIndex -> Integer -> ira -> ira erRoot_R = erRootNewton_R erRoot_IR :: (RA.ERIntApprox 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) => 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 ) 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))