exact-real 0.7.1.0 → 0.8.0.0

raw patch · 7 files changed

+346/−53 lines, 7 filesdep +criteriondep ~basePVP ok

version bump matches the API change (PVP)

Dependencies added: criterion

Dependency ranges changed: base

API changes (from Hackage documentation)

Files

• bench/Bench.hs +36/−0
• exact-real.cabal +35/−19
• readme.md +78/−0
• src/Data/CReal/Internal.hs +76/−22
• test/BoundedFunctions.hs +46/−0
• test/RealFrac.hs +66/−3
• test/Test.hs +9/−9

@@ -0,0 +1,36 @@+module Main where++import Criterion.Main+import Data.CReal.Internal++main :: IO ()+main = defaultMain [ bgroup "pi" [ bench "0"    $ whnf (pi `atPrecision`) 0+                                 , bench "4"    $ whnf (pi `atPrecision`) 4+                                 , bench "16"   $ whnf (pi `atPrecision`) 16+                                 , bench "64"   $ whnf (pi `atPrecision`) 64+                                 , bench "256"  $ whnf (pi `atPrecision`) 256+                                 , bench "1024" $ whnf (pi `atPrecision`) 1024+                                 ]+                   , bgroup "sin 1" [ bench "0"    $ whnf (sin 1 `atPrecision`) 0+                                    , bench "4"    $ whnf (sin 1 `atPrecision`) 4+                                    , bench "16"   $ whnf (sin 1 `atPrecision`) 16+                                    , bench "64"   $ whnf (sin 1 `atPrecision`) 64+                                    , bench "256"  $ whnf (sin 1 `atPrecision`) 256+                                    , bench "1024" $ whnf (sin 1 `atPrecision`) 1024+                                 ]+                   , bgroup "sin (π/4)" [ bench "0"    $ whnf (sin (pi/4) `atPrecision`) 0+                                        , bench "4"    $ whnf (sin (pi/4) `atPrecision`) 4+                                        , bench "16"   $ whnf (sin (pi/4) `atPrecision`) 16+                                        , bench "64"   $ whnf (sin (pi/4) `atPrecision`) 64+                                        , bench "256"  $ whnf (sin (pi/4) `atPrecision`) 256+                                        , bench "1024" $ whnf (sin (pi/4) `atPrecision`) 1024+                                 ]+                   , bgroup "asin (π/4)" [ bench "0"    $ whnf (asin (pi/4) `atPrecision`) 0+                                         , bench "4"    $ whnf (asin (pi/4) `atPrecision`) 4+                                         , bench "16"   $ whnf (asin (pi/4) `atPrecision`) 16+                                         , bench "64"   $ whnf (asin (pi/4) `atPrecision`) 64+                                         , bench "256"  $ whnf (asin (pi/4) `atPrecision`) 256+                                         , bench "1024" $ whnf (asin (pi/4) `atPrecision`) 1024+                                 ]+                                 ]+

@@ -1,30 +1,34 @@-name:                exact-real-version:             0.7.1.0-synopsis:            Exact real arithmetic+name:         exact-real+version:      0.8.0.0+synopsis:     Exact real arithmetic description:   A type to represent exact real number using a fast binary Cauchy sequence-license:             MIT-license-file:        LICENSE-author:              Joe Hermaszewski-maintainer:          Joe Hermaszewski <keep.it.real@monoid.al>-homepage:            http://github.com/expipiplus1/exact-real-bug-reports:         http://github.com/expipiplus1/exact-real/issues-copyright:           2015 Joe Hermaszewski-category:            Math-build-type:          Simple+license:      MIT+license-file: LICENSE+author:       Joe Hermaszewski+maintainer:   Joe Hermaszewski <keep.it.real@monoid.al>+homepage:     http://github.com/expipiplus1/exact-real+bug-reports:  http://github.com/expipiplus1/exact-real/issues+copyright:    2015 Joe Hermaszewski+category:     Math+build-type:   Simple extra-source-files:   .gitignore   readme.md   stack.yaml cabal-version:       >=1.10 +source-repository head+  type: git+  location: https://github.com/expipiplus1/exact-real+ library   exposed-modules:     Data.CReal     Data.CReal.Internal   build-depends:-    base >=4.8 && <4.9,-    integer-gmp < 1.1.0.0+    base        >= 4.8 && < 4.9,+    integer-gmp           < 1.1.0.0   hs-source-dirs:     src   default-language:@@ -32,6 +36,22 @@   ghc-options:     -Wall +benchmark bench+  default-language:+    Haskell2010+  type:+    exitcode-stdio-1.0+  ghc-options:+    -Wall -threaded -O2+  hs-source-dirs:+    bench+  main-is:+    Bench.hs+  build-depends:+    base        >= 4 && < 5,+    criterion   >= 1.1 && < 1.2,+    exact-real+ test-suite test   default-language:     Haskell2010@@ -44,6 +64,7 @@   main-is:     Test.hs   other-modules:+    BoundedFunctions,     Data.CReal.Extra,     Data.Monoid.Extra,     Data.Ratio.Extra,@@ -86,9 +107,4 @@     directory >= 1.0 && < 1.3,     doctest   >= 0.8 && < 0.11,     filepath  >= 1.3 && < 1.5--source-repository head-  type: git-  location: https://github.com/expipiplus1/exact-real- 

@@ -2,3 +2,81 @@ ==========  Exact real arithmetic implemented by fast binary Cauchy sequences.+++Motivating Example+-------------------++Compare evaluating Euler's identity with a `Float`:++``` haskell+λ> let i = 0 :+ 1+λ> exp (i * pi) + 1 :: Complex Float+0.0 :+ (-8.742278e-8)+```++... and with a `CReal`++``` haskell+λ> import Data.CReal+λ> let i = 0 :+ 1+λ> exp (i * pi) + 1 :: Complex (CReal 0)+0 :+ 0+```++Implementation+--------------++The basic operations have explanations and proofs of correctness+[here][correctness].++`CReal`'s phantom type parameter `n :: Nat` represents the precision at which+values should be evaluated at when converting to a less precise representation.+For instance the definition of `x == y` in the instance for `Eq` evaluates `x -+y` at precision `n` and compares the resulting `Integer` to zero. I think that+this is the most reasonable solution to the fact that lots of of operations+(such as equality) are not computable on the reals but we want to pretend that+they are for the sake of writing useful programs. Please see the+[Caveats](#caveats) section for more information.++The `CReal` type is an instance of `Num`, `Fractional`, `Floating`, `Real`,+`RealFrac`, `RealFloat`, `Eq`, `Ord`, `Show` and `Read`. The only functions not+implemented are a handful from `RealFloat` which assume the number is+implemented with a mantissa and exponent.++There is a comprehensive test suite to test the properties of these classes.++The performance isn't terrible on most operations but it's obviously not nearly+as speedy as performing the operations on `Float` or `Double`. The only two+super slow functions are `asinh` and `atanh` at the moment.+++Caveats+-------++The implementation is not without its caveats however. The big gotcha is that+although internally the `CReal n`s are represented exactly, whenever a value is+extracted to another type such as a `Rational` or `Float` it is evaluated to+within `2^-p` of the true value.++For example when using the `CReal 0` type (numbers within 1 of the true value)+one can produce the following:++``` haskell+λ> 0.5 == (1 :: CReal 0)+True+λ> 0.5 * 2 == (1 :: CReal 0) * 2+False+```++Contributing+------------++Contributions and bug reports are welcome!++Please feel free to contact me on GitHub or as "jophish" on freenode.++-Joe++[goldberg]: http://www.validlab.com/goldberg/paper.pdf "What Every Computer Scientist Should Know About Floating-Point Arithmetic"+[correctness]: http://www.imada.sdu.dk/~kornerup/RNC4/papers/p07.ps "The Correctness of an Implementation of Exact Arithmetic"

@@ -14,6 +14,13 @@   , atPrecision   , crealPrecision +  , (.*)+  , (*.)+  , (.*.)+  , mulBounded+  , mulBoundedL+  , recipBounded+   , expBounded   , logBounded @@ -44,6 +51,8 @@ import GHC.TypeLits import Numeric (readSigned, readFloat) +{-# ANN module "HLint: ignore Reduce duplication" #-}+ -- $setup -- >>> :set -XDataKinds @@ -98,10 +107,13 @@ -- >>> signum (0.1 :: CReal 3) -- 1.0 instance Num (CReal n) where+  {-# INLINE fromInteger #-}   fromInteger i = CR (\p -> i * 2 ^ p) +  {-# INLINE negate #-}   negate (CR x) = CR (negate . x) +  {-# INLINE abs #-}   abs (CR x) = CR (abs . x)    {-# INLINE (+) #-}@@ -120,8 +132,7 @@  -- | Taking the reciprocal of zero will not terminate instance Fractional (CReal n) where-  -- This should be in base-  fromRational n = fromInteger (numerator n) / fromInteger (denominator n)+  fromRational n = fromInteger (numerator n) *. recipBounded (fromInteger (denominator n))    {-# INLINE recip #-}   -- TODO: Make recip 0 throw an error (if, for example, it would take more@@ -144,10 +155,13 @@   -- | Range reduction on the principle that ln (a * b) = ln a + ln b   log x = let CR o = x               l = log2 (o 2) - 2-              a = x `shiftR` l-          in if | l < 0  -> - log (recip x)+          in if   -- x <= 0.75+                | l < 0  -> - log (recip x)+                  -- 0.75 <= x <= 2                 | l == 0 -> logBounded x-                | l > 0  -> logBounded a + fromIntegral l * ln2+                  -- x >= 2+                | 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)@@ -173,19 +187,23 @@                    ,          cosBounded . (piBy4-)]           in (fs !! octant) offset -  -- TODO: use multiplyBounded here-  tan x = sin x / cos x+  tan x = sin x .* recip (cos x) -  asin x = 2 * atan (x / (1 + sqrt (1 - x*x)))+  asin x = 2 * atan (x .*. recipBounded (1 + sqrt (1 - x.*.x)))    acos x = pi/2 - asin x -  atan x = let -- q is x to the nearest 1/4+  atan x = let -- q is 4 times x to within 1/4                q = x `atPrecision` 2-           in if | q <  -4 -> atanBounded (negate (recip x)) - pi / 2-                 | q == -4 -> -pi / 4 - atanBounded ((x + 1) / (x - 1))-                 | q ==  4 -> pi / 4 + atanBounded ((x - 1) / (x + 1))-                 | q >   4 -> pi / 2 - atanBounded (recip x)+           in if   -- x <= -1+                 | q <  -4 -> atanBounded (negate (recipBounded x)) - pi / 2+                   -- -1.25 <= x <= -0.75+                 | q == -4 -> -pi / 4 - atanBounded ((x + 1) .*. recipBounded (x - 1))+                   -- 0.75 <= x <= 1.25+                 | q ==  4 -> pi / 4 + atanBounded ((x - 1) .*. recipBounded (x + 1))+                   -- x >= 1+                 | q >   4 -> pi / 2 - atanBounded (recipBounded x)+                   -- -0.75 <= x <= 0.75                  | otherwise -> atanBounded x    -- TODO: benchmark replacing these with their series expansion@@ -204,8 +222,9 @@                  in x `atPrecision` p % 2^p  instance KnownNat n => RealFrac (CReal n) where-  properFraction x = let n = x `atPrecision` 0-                         f = x - fromIntegral n+  properFraction x = let p = crealPrecision x+                         n = (x `atPrecision` p) `quot` 2^p+                         f =  x - fromInteger n                      in (fromInteger n, f)  -- | Several of the functions in this class ('floatDigits', 'floatRange',@@ -273,6 +292,40 @@ ln2 = logBounded 2  --+-- Bounded multiplication+--++infixl 7 `mulBounded`, `mulBoundedL`, .*, *., .*.++(.*), (*.), (.*.) :: CReal n -> CReal n -> CReal n+(.*) = mulBoundedL+(*.) = flip mulBoundedL+(.*.) = mulBounded++-- | The first argument to @mulBoundedL@ must be in the 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))++-- | Both arguments to @mulBounded@ must be in the 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))++-- | The absolute value of the argument to @recipBounded@ must be 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)++-- -- Bounded exponential functions -- @@ -281,11 +334,11 @@ expBounded x = let q = [1 % (n!) | n <- [0..]]                in powerSeries q (max 5) x --- | The input must be in [1..2]+-- | The input must be in [2/3..2] logBounded :: CReal n -> CReal n logBounded x = let q = [1 % n | n <- [1..]]-                   y = (x - 1) / x-               in y * powerSeries q (*2) y+                   y = (x - 1) .* recip x+               in y .* powerSeries q id y  -- -- Bounded trigonometric functions@@ -294,18 +347,19 @@ -- | The input to sinBounded must be in (-1..1) sinBounded :: CReal n -> CReal n sinBounded x = let q = alternateSign (scanl' (*) 1 [ 1 % (n*(n+1)) | n <- [2,4..]])-               in x * powerSeries q (max 1) (x*x)+               in x * powerSeries q (max 1) (x .*. x)  -- | The input to cosBounded must be in (-1..1) cosBounded :: CReal n -> CReal n cosBounded x = let q = alternateSign (scanl' (*) 1 [1 % (n*(n+1)) | n <- [1,3..]])-               in powerSeries q (max 1) (x*x)+               in powerSeries q (max 1) (x .*. x)  -- | The input to atanBounded must be in [-1..1] atanBounded :: CReal n -> CReal n atanBounded x = let q = scanl' (*) 1 [n % (n + 1) | n <- [2,4..]]-                    d = 1 + x * x-                in CR (\p -> ((x/d) * powerSeries q (+1) (x*x/d)) `atPrecision` p)+                    d = 1 + x .*. x+                    rd = recipBounded d+                in CR (\p -> ((x .*. rd) .* powerSeries q (+1) (x .*. x .*. rd)) `atPrecision` p)  -- -- Multiplication with powers of two

@@ -0,0 +1,46 @@+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeFamilies #-}+module BoundedFunctions+  ( boundedFunctions+  ) where++import Test.Tasty (testGroup, TestTree)+import Test.Tasty.QuickCheck ((==>), testProperty)+import Test.QuickCheck.Extra (BiunitInterval(..), UnitInterval(..))+import Test.QuickCheck.Checkers ((=-=))+import GHC.TypeLits (KnownNat)++import Data.CReal.Internal+import Data.CReal.Extra ()++boundedFunctions :: forall a n. (KnownNat n, a ~ CReal n) =>  a -> TestTree+boundedFunctions _ = testGroup "bounded functions" ts+  where ts = [ testProperty "mulBounded"+                  (\(BiunitInterval x) (BiunitInterval y) ->+                    (x :: a) * y =-= x .*. y)+             , testProperty "mulBoundedL"+                  (\(BiunitInterval x) y ->+                    (x :: a) * y =-= x .* y)+             , testProperty "mulBoundedR"+                  (\x (BiunitInterval y) ->+                    (x :: a) * y =-= x *. y)+             , testProperty "recipBounded"+                  (\x -> (abs x >= 1) ==> recip (x::a) =-= recipBounded x)+             , testProperty "expBounded"+                  (\(BiunitInterval x) ->+                    exp (x :: a) =-= expBounded x)+             , testProperty "logBounded"+                  (\(UnitInterval x) ->+                    let x' = (x :: a) * (2 - 2/3) + 2/3+                    in log x' =-= logBounded x')+             , testProperty "sinBounded"+                  (\(BiunitInterval x) ->+                    sin (x :: a) =-= sinBounded x)+             , testProperty "cosBounded"+                  (\(BiunitInterval x) ->+                    cos (x :: a) =-= cosBounded x)+             , testProperty "atanBounded"+                  (\(BiunitInterval x) ->+                    atan (x :: a) =-= atanBounded x)+             ]+

@@ -4,6 +4,7 @@   ( realFrac   ) where +import Data.Function (on) import Data.Ratio.Extra () import Test.QuickCheck.Checkers (EqProp, (=-=)) import Test.Tasty (testGroup, TestTree)@@ -13,8 +14,20 @@ realFrac :: forall a. (Arbitrary a, EqProp a, Show a, RealFrac a) =>             a -> TestTree realFrac x = testGroup "Test RealFrac instance" ts-  where ts = [ properFractionLaws "properFraction laws" x ]+  where ts = [ properFractionLaws "properFraction laws" x+             , truncateLaws "truncate laws" x+             , roundLaws "round laws" x+             , ceilingLaws "ceiling laws" x+             , floorLaws "floor laws" x+             ] +-- Thses are used to cope with CReal 0 being a little weird with comparisons+-- between non-integers+infix 4 <., >.+(<.), (>.) :: (Ord a, Num a) => a -> a -> Bool+(<.) = (<) `on` (*2)+(>.) = (>) `on` (*2)+ -- | This tests a slightly different law for n having the same sign as x properFractionLaws :: forall a. (Arbitrary a, EqProp a, Show a, RealFrac a) =>                       String -> a -> TestTree@@ -24,10 +37,60 @@                                    in x =-= fromInteger n + f)              , testProperty "n has same sign or is zero"                             (\x -> let (n, _) = properFraction (x :: a)-                                   in n == 0 || sign x == sign (n::Int))+                                   in n == 0 || sign x == sign (n::Integer))              , testProperty "abs f < 1"                             (\x -> let (_::Int, f) = properFraction (x :: a)-                                   in abs f < 1)+                                   in abs f <. 1)+             , testProperty "f has same sign or is zero"+                            (\x -> let (_::Int, f) = properFraction (x :: a)+                                   in f == 0 || sign x == sign f)+             ]++truncateLaws :: forall a. (Arbitrary a, EqProp a, Show a, RealFrac a) =>+                           String -> a -> TestTree+truncateLaws s _ = testGroup s ts+  where ts = [ testProperty "abs (truncate x) <= abs x"+                            (\x -> let t = truncate (x :: a)+                                   in fromInteger (abs t) <= abs x)+             , testProperty "abs (truncate x) + 1 > abs x"+                            (\x -> let t = truncate (x :: a)+                                   in fromInteger (abs t + 1) >. abs x)+             , testProperty "truncate x has same sign or is zero"+                            (\x -> let t = truncate (x :: a)+                                   in t == 0 || sign x == sign (t::Integer))+             ]++roundLaws :: forall a. (Arbitrary a, EqProp a, Show a, RealFrac a) =>+                           String -> a -> TestTree+roundLaws s _ = testGroup s ts+  where ts = [ testProperty "abs (round x - x) <= 0.5"+                            (\x -> let r = round (x :: a)+                                       in abs (fromInteger r - x) <= 0.5)+             , testProperty "round to even if eqiudistant"+                            (\i -> let x = fromInteger i + 0.5 :: a+                                   in even (round x :: Integer))+             ]++ceilingLaws :: forall a. (Arbitrary a, EqProp a, Show a, RealFrac a) =>+                           String -> a -> TestTree+ceilingLaws s _ = testGroup s ts+  where ts = [ testProperty "ceiling x - 1 < x"+                            (\x -> let c = ceiling (x :: a)+                                   in fromInteger c - 1 <. x)+             , testProperty "ceiling x >= x"+                            (\x -> let c = ceiling (x :: a)+                                   in fromInteger c >= x)+             ]++floorLaws :: forall a. (Arbitrary a, EqProp a, Show a, RealFrac a) =>+                           String -> a -> TestTree+floorLaws s _ = testGroup s ts+  where ts = [ testProperty "floor x + 1 > x"+                            (\x -> let f = floor (x :: a)+                                   in fromInteger f + 1 >. x)+             , testProperty "floor x <= x"+                            (\x -> let f = floor (x :: a)+                                   in fromInteger f <= x)              ]  data Sign = Positive

@@ -6,12 +6,13 @@  import Data.Ratio ((%)) import Test.Tasty (testGroup, TestTree)-import Test.Tasty.QuickCheck (Positive(..), testProperty, (===), Property)+import Test.Tasty.QuickCheck (Positive(..), testProperty, (===), Property, (==>), (.&&.), testProperty) import Test.Tasty.TH (defaultMainGenerator)  import Data.CReal.Internal import Data.CReal.Extra () +import BoundedFunctions (boundedFunctions) import Floating (floating) import Ord (ord) import Read (read')@@ -23,11 +24,6 @@ -- precisions type Precision = 10 -infixr 1 ==>-(==>) :: Bool -> Bool -> Bool-False ==> _ = True-True ==> b = b- {-# ANN test_floating "HLint: ignore Use camelCase" #-} test_floating :: [TestTree] test_floating = [floating (undefined :: CReal Precision)]@@ -38,7 +34,7 @@  {-# ANN test_real "HLint: ignore Use camelCase" #-} test_real :: [TestTree]-test_real = [ real (\x -> 1 % toInteger (crealPrecision (x::CReal Precision))) ]+test_real = [ real (\x -> 1 % toInteger (max 1 (crealPrecision (x::CReal Precision)))) ]  {-# ANN test_realFrac "HLint: ignore Use camelCase" #-} test_realFrac :: [TestTree]@@ -52,9 +48,9 @@ test_read :: [TestTree] test_read = [ read' (undefined :: CReal Precision) ] -prop_decimalDigits :: Positive Int -> Bool+prop_decimalDigits :: Positive Int -> Property prop_decimalDigits (Positive p) = let d = decimalDigitsAtPrecision p-                                  in 10^d >= (2^p :: Integer) &&+                                  in 10^d >= (2^p :: Integer) .&&.                                      (d > 0 ==> 10^(d-1) < (2^p :: Integer))  prop_showIntegral :: Integer -> Property@@ -70,6 +66,10 @@ prop_showNumDigits (Positive places) x =   let s = rationalToDecimal places x   in length (dropWhile (/= '.') s) === places + 1++{-# ANN test_boundedFunctions "HLint: ignore Use camelCase" #-}+test_boundedFunctions :: [TestTree]+test_boundedFunctions = [ boundedFunctions (undefined :: CReal Precision) ]  main :: IO () main = $(defaultMainGenerator)