exact-real 0.9.0 → 0.10.0
raw patch · 3 files changed
+56/−45 lines, 3 filesdep +memoizePVP ok
version bump matches the API change (PVP)
Dependencies added: memoize
API changes (from Hackage documentation)
+ Data.CReal.Internal: crMemoize :: (Int -> Integer) -> CReal n
Files
- exact-real.cabal +3/−2
- readme.md +2/−0
- src/Data/CReal/Internal.hs +51/−43
exact-real.cabal view
@@ -1,5 +1,5 @@ name: exact-real-version: 0.9.0+version: 0.10.0 synopsis: Exact real arithmetic description: A type to represent exact real number using a fast binary Cauchy sequence@@ -28,7 +28,8 @@ Data.CReal.Internal build-depends: base >= 4.8 && < 4.9,- integer-gmp < 1.1.0.0+ integer-gmp < 1.1.0.0,+ memoize >= 0.7 && < 0.8 hs-source-dirs: src default-language:
readme.md view
@@ -1,6 +1,8 @@ exact-real ========== +[](https://hackage.haskell.org/package/exact-real)+ Exact real arithmetic implemented by fast binary Cauchy sequences.
src/Data/CReal/Internal.hs view
@@ -41,6 +41,7 @@ , cosBounded -- * Utilities for operating inside CReals+ , crMemoize , powerSeries , alternateSign @@ -62,6 +63,7 @@ import GHC.Integer.Logarithms (integerLog2#, integerLogBase#) import GHC.TypeLits import Numeric (readSigned, readFloat)+import Data.Function.Memoize (memoize) {-# ANN module "HLint: ignore Reduce duplication" #-} @@ -77,6 +79,12 @@ -- a function from Ints to Integers. newtype CReal (n :: Nat) = CR (Int -> Integer) +-- | 'crMemoize' takes a fast binary Cauchy sequence and returns a CReal+-- represented by that sequence which will memoize the values at each+-- precision. This is essential for getting good performance.+crMemoize :: (Int -> Integer) -> CReal n+crMemoize = CR . memoize+ -- | crealPrecision x returns the type level parameter representing x's default -- precision. --@@ -124,27 +132,27 @@ -- 1.0 instance Num (CReal n) where {-# INLINE fromInteger #-}- fromInteger i = CR (\p -> i * 2 ^ p)+ fromInteger i = crMemoize (\p -> i * 2 ^ p) {-# INLINE negate #-}- negate (CR x) = CR (negate . x)+ negate (CR x) = crMemoize (negate . x) {-# INLINE abs #-}- abs (CR x) = CR (abs . x)+ abs (CR x) = crMemoize (abs . x) {-# INLINE (+) #-}- CR x1 + CR x2 = CR (\p -> let n1 = x1 (p + 2)- n2 = x2 (p + 2)- in (n1 + n2) /. 4)+ CR x1 + CR x2 = crMemoize (\p -> let n1 = x1 (p + 2)+ n2 = x2 (p + 2)+ in (n1 + n2) /. 4) {-# INLINE (*) #-}- CR x1 * CR x2 = CR (\p -> let s1 = log2 (abs (x1 0) + 2) + 3- s2 = log2 (abs (x2 0) + 2) + 3- n1 = x1 (p + s2)- n2 = x2 (p + s1)- in (n1 * n2) /. 2^(p + s1 + s2) )+ CR x1 * CR x2 = crMemoize (\p -> let s1 = log2 (abs (x1 0) + 2) + 3+ s2 = log2 (abs (x2 0) + 2) + 3+ n1 = x1 (p + s2)+ n2 = x2 (p + s1)+ in (n1 * n2) /. 2^(p + s1 + s2) ) - signum x = CR (\p -> signum (x `atPrecision` p) * 2^p)+ signum x = crMemoize (\p -> signum (x `atPrecision` p) * 2^p) -- | Taking the reciprocal of zero will not terminate instance Fractional (CReal n) where@@ -153,9 +161,9 @@ {-# INLINE recip #-} -- TODO: Make recip 0 throw an error (if, for example, it would take more -- than 4GB of memory to represent the result)- recip (CR x) = CR (\p -> let s = findFirstMonotonic ((3 <=) . abs . x)- n = x (p + 2 * s + 2)- in 2^(2 * p + 2 * s + 2) /. n)+ recip (CR x) = crMemoize (\p -> let s = findFirstMonotonic ((3 <=) . abs . x)+ n = x (p + 2 * s + 2)+ in 2^(2 * p + 2 * s + 2) /. n) instance Floating (CReal n) where -- TODO: Could we use something faster such as Ramanujan's formula@@ -179,8 +187,8 @@ | l > 0 -> let a = x `shiftR` l in logBounded a + fromIntegral l *. ln2 - sqrt (CR x) = CR (\p -> let n = x (2 * p)- in isqrt n)+ sqrt (CR x) = crMemoize (\p -> let n = x (2 * p)+ in isqrt n) -- | This will diverge when the base is not positive x ** y = exp (log x * y)@@ -266,7 +274,7 @@ isDenormalized _ = False isNegativeZero _ = False isIEEE _ = False- atan2 y x = CR (\p ->+ atan2 y x = crMemoize (\p -> let y' = y `atPrecision` p x' = x `atPrecision` p θ = if | x' > 0 -> atan (y/x)@@ -292,8 +300,8 @@ instance KnownNat n => Ord (CReal n) where compare x y = let p = crealPrecision x in compare ((x - y) `atPrecision` p) 0- max (CR x) (CR y) = CR (\p -> max (x p) (y p))- min (CR x) (CR y) = CR (\p -> min (x p) (y p))+ max (CR x) (CR y) = crMemoize (\p -> max (x p) (y p))+ min (CR x) (CR y) = crMemoize (\p -> min (x p) (y p)) -------------------------------------------------------------------------------- -- Some utility functions@@ -330,33 +338,33 @@ -- | A more efficient multiply with the restriction that the first argument -- must be in the closed range [-1..1] mulBoundedL :: CReal n -> CReal n -> CReal n-mulBoundedL (CR x1) (CR x2) = CR (\p -> let s1 = 4- s2 = log2 (abs (x2 0) + 2) + 3- n1 = x1 (p + s2)- n2 = x2 (p + s1)- in (n1 * n2) /. 2^(p + s1 + s2))+mulBoundedL (CR x1) (CR x2) = crMemoize (\p -> let s1 = 4+ s2 = log2 (abs (x2 0) + 2) + 3+ n1 = x1 (p + s2)+ n2 = x2 (p + s1)+ in (n1 * n2) /. 2^(p + s1 + s2)) -- | A more efficient multiply with the restriction that both values must be -- in the closed range [-1..1] mulBounded :: CReal n -> CReal n -> CReal n-mulBounded (CR x1) (CR x2) = CR (\p -> let s1 = 4- s2 = 4- n1 = x1 (p + s2)- n2 = x2 (p + s1)- in (n1 * n2) /. 2^(p + s1 + s2))+mulBounded (CR x1) (CR x2) = crMemoize (\p -> let s1 = 4+ s2 = 4+ n1 = x1 (p + s2)+ n2 = x2 (p + s1)+ in (n1 * n2) /. 2^(p + s1 + s2)) -- | A more efficient 'recip' with the restriction that the input must have -- absolute value greater than or equal to 1 recipBounded :: CReal n -> CReal n-recipBounded (CR x) = CR (\p -> let s = 2- n = x (p + 2 * s + 2)- in 2^(2 * p + 2 * s + 2) /. n)+recipBounded (CR x) = crMemoize (\p -> let s = 2+ n = x (p + 2 * s + 2)+ in 2^(2 * p + 2 * s + 2) /. n) -- | Return the square of the input, more efficient than @('*')@ square :: CReal n -> CReal n-square (CR x) = CR (\p -> let s = log2 (abs (x 0) + 2) + 3- n = x (p + s)- in (n * n) /. 2^(p + 2 * s))+square (CR x) = crMemoize (\p -> let s = log2 (abs (x 0) + 2) + 3+ n = x (p + s)+ in (n * n) /. 2^(p + 2 * s)) -- -- Bounded exponential functions and expPosNeg@@ -407,7 +415,7 @@ atanBounded x = let q = scanl' (*) 1 [n % (n + 1) | n <- [2,4..]] d = 1 + x .*. x rd = recipBounded d- in CR (\p -> ((x .*. rd) .* powerSeries q (+1) (x .*. x .*. rd)) `atPrecision` p)+ in (x .*. rd) .* powerSeries q (+1) (x .*. x .*. rd) -- -- Multiplication with powers of two@@ -421,10 +429,10 @@ -- -- This can be faster than doing the division shiftR :: CReal n -> Int -> CReal n-shiftR (CR x) n = CR (\p -> let p' = p - n- in if p' >= 0- then x p'- else x 0 /. 2^(-p'))+shiftR (CR x) n = crMemoize (\p -> let p' = p - n+ in if p' >= 0+ then x p'+ else x 0 /. 2^(-p')) -- | @x \`shiftL\` n@ is equal to @x@ multiplied by 2^@n@ --@@ -540,8 +548,8 @@ -- >>> powerSeries [1 % (n!) | n <- [0..]] (max 5) 1 :: CReal 218 -- 2.718281828459045235360287471352662497757247093699959574966967627724 powerSeries :: [Rational] -> (Int -> Int) -> CReal n -> CReal n-powerSeries q termsAtPrecision (CR x) =- CR (\p -> let t = termsAtPrecision p+powerSeries q termsAtPrecision (CR x) = crMemoize+ (\p -> let t = termsAtPrecision p d = log2 (toInteger t) + 2 p' = p + d p'' = p' + d