diff --git a/bench/Bench.hs b/bench/Bench.hs new file mode 100644 --- /dev/null +++ b/bench/Bench.hs @@ -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 + ] + ] + diff --git a/exact-real.cabal b/exact-real.cabal --- a/exact-real.cabal +++ b/exact-real.cabal @@ -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 -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 +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 - diff --git a/readme.md b/readme.md --- a/readme.md +++ b/readme.md @@ -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" diff --git a/src/Data/CReal/Internal.hs b/src/Data/CReal/Internal.hs --- a/src/Data/CReal/Internal.hs +++ b/src/Data/CReal/Internal.hs @@ -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 diff --git a/test/BoundedFunctions.hs b/test/BoundedFunctions.hs new file mode 100644 --- /dev/null +++ b/test/BoundedFunctions.hs @@ -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) + ] + diff --git a/test/RealFrac.hs b/test/RealFrac.hs --- a/test/RealFrac.hs +++ b/test/RealFrac.hs @@ -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 diff --git a/test/Test.hs b/test/Test.hs --- a/test/Test.hs +++ b/test/Test.hs @@ -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)