exact-real 0.2.1.0 → 0.3.0.0

raw patch · 9 files changed

+313/−23 lines, 9 filesdep +tasty-hunitPVP ok

version bump matches the API change (PVP)

Dependencies added: tasty-hunit

API changes (from Hackage documentation)

Files

• exact-real.cabal +3/−1
• src/Data/CReal.hs +4/−0
• src/Data/CReal/Internal.hs +189/−3
• test/Data/CReal/Extra.hs +20/−1
• test/Floating.hs +66/−0
• test/Fractional.hs +5/−1
• test/Test.hs +10/−13
• test/Test/QuickCheck/Classes/Extra.hs +3/−1
• test/Test/QuickCheck/Extra.hs +13/−3

@@ -1,5 +1,5 @@ name:                exact-real-version:             0.2.1.0+version:             0.3.0.0 synopsis:            Exact real arithmetic description:         please see readme.md license:             MIT@@ -40,6 +40,7 @@   main-is:     Test.hs   other-modules:+    Floating,     Fractional,     Num,     Ord,@@ -54,6 +55,7 @@     tasty            >= 0.10 && < 0.12,     tasty-th         >= 0.1  && < 0.2,     tasty-quickcheck >= 0.8  && < 0.9,+    tasty-hunit      >= 0.9  && < 0.10,     QuickCheck       >= 2.8  && < 2.9,     checkers         >= 0.4  && < 0.5,     random           >= 1.0  && < 1.2,

@@ -1,3 +1,7 @@+-----------------------------------------------------------------------------+-- | This module exports everything you need to use exact real numbers+----------------------------------------------------------------------------+ module Data.CReal   ( CReal   , atPrecision

@@ -1,23 +1,43 @@+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE KindSignatures #-} {-# LANGUAGE MagicHash #-} {-# LANGUAGE MultiWayIf #-}-{-# LANGUAGE KindSignatures #-}-{-# LANGUAGE DataKinds #-}+{-# LANGUAGE PostfixOperators #-} +-----------------------------------------------------------------------------+-- | This module exports a bunch of utilities for working inside the CReal+-- datatype. One should be careful to maintain the CReal invariant when using+-- these functions+---------------------------------------------------------------------------- module Data.CReal.Internal   ( CReal(..)   , atPrecision   , crealPrecision +  , expBounded+  , logBounded++  , atanBounded+  , sinBounded+  , cosBounded+   , shiftL   , shiftR +  , powerSeries+  , alternateSign+   , (/.)   , log2   , log10+  , isqrt++  , showAtPrecision   , decimalDigitsAtPrecision   , rationalToDecimal   ) where +import Data.List (scanl') import Data.Ratio (numerator,denominator,(%)) import GHC.Base (Int(..)) import GHC.Integer.Logarithms (integerLog2#, integerLogBase#)@@ -31,7 +51,7 @@ default ()  -- | The type CReal represents a fast binary Cauchy sequence. This is--- a Cauchy sequence with the property that the pth element will be within+-- a Cauchy sequence with the invariant that the pth element will be within -- 2^-p of the true value. Internally this sequence is represented as -- a function from Ints to Integers. newtype CReal (n :: Nat) = CR (Int -> Integer)@@ -106,6 +126,74 @@                                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+  pi = 4 * piBy4++  exp x = let CR o = x / ln2+              l = o 0+              y = x - fromInteger l * ln2+          in if l == 0+               then expBounded x+               else expBounded y `shiftL` fromInteger l++  -- | 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)+                | l == 0 -> logBounded x+                | l > 0  -> logBounded a + fromIntegral l * ln2++  sqrt (CR x) = CR (\p -> let n = x (2 * p)+                          in isqrt n)++  -- | This will diverge when the base is not positive+  x ** y = exp (log x * y)++  logBase x y = log y / log x++  sin x = cos (x - pi / 2)++  cos x = let CR o = x / piBy4+              s = o 1 /. 2+              octant = fromInteger $ s `mod` 8+              offset = x - (fromIntegral s * piBy4)+              fs = [          cosBounded+                   , negate . sinBounded . subtract piBy4+                   , negate . sinBounded+                   , negate . cosBounded . (piBy4-)+                   , negate . cosBounded+                   ,          sinBounded . subtract piBy4+                   ,          sinBounded+                   ,          cosBounded . (piBy4-)]+          in (fs !! octant) offset++  -- TODO: use multiplyBounded here+  tan x = sin x / cos x++  asin x = 2 * atan (x / (1 + sqrt (1 - x*x)))++  acos x = pi/2 - asin x++  atan x = let -- q is x to the nearest 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)+                 | otherwise -> atanBounded x++  -- TODO: benchmark replacing these with their series expansion+  sinh x = (exp x - exp (-x)) / 2+  cosh x = (exp x + exp (-x)) / 2+  tanh x = let e2x = exp (2 * x)+           in (e2x - 1) / (e2x + 1)++  asinh x = log (x + sqrt (x * x + 1))+  acosh x = log (x + sqrt (x + 1) * sqrt (x - 1))+  atanh x = (log (1 + x) - log (1 - x)) / 2+ -- | Values of type @CReal p@ are compared for equality at precision @p@. This -- may cause values which differ by less than 2^-p to compare as equal. --@@ -128,15 +216,70 @@ --------------------------------------------------------------------------------  --+-- Constants+--++piBy4 :: CReal n+piBy4 = 4 * atanBounded (1/5) - atanBounded (1 / 239) -- Machin Formula++ln2 :: CReal n+ln2 = logBounded 2++--+-- Bounded exponential functions+--++-- | The input to expBounded must be in the range (-1..1)+expBounded :: CReal n -> CReal n+expBounded x = let q = [1 % (n!) | n <- [0..]]+               in powerSeries q (max 5) x++-- | The input must be in [1..2]+logBounded :: CReal n -> CReal n+logBounded x = let q = [1 % n | n <- [1..]]+                   y = (x - 1) / x+               in y * powerSeries q (*2) y++--+-- Bounded trigonometric functions+--++-- | 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)++-- | 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)++-- | 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)++-- -- Multiplication with powers of two -- +-- | @x \`shiftR\` n@ is equal to @x@ divided by 2^@n@+--+-- @n@ can be negative or zero+--+-- 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')) +-- | @x \`shiftL\` n@ is equal to @x@ multiplied by 2^@n@+--+-- @n@ can be negative or zero+--+-- This can be faster than doing the multiplication shiftL :: CReal n -> Int -> CReal n shiftL x = shiftR x . negate @@ -189,6 +332,20 @@ log10 :: Integer -> Int log10 x = I# (integerLogBase# 10 x) +-- | @isqrt x@ returns the square root of @x@ rounded towards zero.+isqrt :: Integer -> Integer+isqrt x | x < 0     = error "Sqrt applied to negative Integer"+        | x == 0    = 0+        | otherwise = until satisfied improve initialGuess+  where improve r    = (r + (x `div` r)) `div` 2+        satisfied r  = sq r <= x && sq (r + 1) > x+        initialGuess = 2 ^ (log2 x `div` 2)+        sq r         = r * r++-- | Factorial function+(!) :: Integer -> Integer+(!) x = product [2..x]+ -- -- Searching --@@ -203,4 +360,33 @@                            in if | l+1 == u  -> l                                  | p m       -> binarySearch l m                                  | otherwise -> binarySearch m u+++--+-- Power series+--++-- | Apply 'negate' to every other element, starting with the second+--+-- >>> alternateSign [1..5]+-- [1, -2, 3, -4, 5]+alternateSign :: Num a => [a] -> [a]+alternateSign = zipWith ($) (cycle [id, negate])++-- | @powerSeries q f x `atPrecision` p@ will evaluate the power series with+-- coefficients @q@ at precision @f p@ at @x@+--+-- @f@ should be a function such that the CReal invariant is maintained+--+-- See any of the trig functions for an example+powerSeries :: [Rational] -> (Int -> Int) -> CReal n -> CReal n+powerSeries q termsAtPrecision (CR x) =+  CR (\p -> let t = termsAtPrecision p+                d = log2 (toInteger t) + 2+                p' = p + d+                p'' = p' + d+                m = x p''+                xs = (%1) <$> iterate (\e -> m * e /. 2^p'') (2^p')+                r = sum . take (t + 1) . fmap (round . (* (2^d))) $ zipWith (*) q xs+            in r /. 4^d) 

@@ -4,10 +4,29 @@   ( module Data.CReal   ) where -import Test.QuickCheck.Checkers (EqProp(..), eq) import Data.CReal+import Data.CReal.Internal (log2)+import Data.Ratio ((%)) import GHC.TypeLits+import System.Random (Random(..))+import Test.QuickCheck (Arbitrary(..), choose)+import Test.QuickCheck.Checkers (EqProp(..), eq)  instance KnownNat n => EqProp (CReal n) where   (=-=) = eq++instance KnownNat n => Arbitrary (CReal n) where+  arbitrary = do+    integralPart <- fromInteger <$> arbitrary+    fractionalPart <- choose (-0.5, 0.5)+    pure (integralPart + fractionalPart)++instance KnownNat n => Random (CReal n) where+  randomR (lo, hi) g = let d = hi - lo+                           l = 1 + log2 (abs d `atPrecision` 0)+                           p = l + crealPrecision lo+                           (n, g') = randomR (0, 2^p) g+                           r = fromRational (n % 2^p)+                       in (r * d + lo, g')+  random = randomR (0, 1) 

@@ -0,0 +1,66 @@+{-# LANGUAGE NoMonomorphismRestriction #-}+{-# LANGUAGE ScopedTypeVariables #-}++module Floating+  ( floating+  ) where++import Fractional (fractional)+import System.Random (Random)+import Test.QuickCheck.Checkers (EqProp, (=-=), inverseL)+import Test.QuickCheck.Extra (UnitInterval(..), Tiny(..), BiunitInterval)+import Test.Tasty (testGroup, TestTree)+import Test.Tasty.QuickCheck (testProperty, NonNegative(..), Positive(..), Arbitrary, (==>))+import Test.Tasty.HUnit (testCase, (@?=))++floating :: forall a. (Arbitrary a, EqProp a, Show a, Floating a, Ord a, Random a) =>+            a -> TestTree+floating _ = testGroup "Test Floating instance" ts+  where e = exp 1+        ts = [ fractional (undefined :: a)+             , testCase "π/4 = atan 1" ((pi::a) @?= 4 * atan 1)+             , testProperty "log == logBase e"+                            (log =-= logBase (e :: Positive a))+             , testProperty "exp == (e **)" (exp =-= ((e::a) **))+             , testProperty "sqrt x * sqrt x = x"+                            (\(NonNegative (x :: a)) -> let r = sqrt x+                                                        in r * r == x)+             , testProperty "law of exponents"+                            (\(Positive (base :: a)) x y ->+                              base ** (x + y) =-= base ** x * base ** y)+             , testProperty "logarithm definition"+                            (\(Positive (b :: a)) (Tiny c) ->+                              let x = b ** c+                              in b /= 1 ==> c =-= logBase b x)+             , testProperty "sine cosine definition"+                            (\x (y :: a) ->+                              cos (x - y) =-= cos x * cos y + sin x * sin y)+               -- TODO: Use open interval+             , testProperty "0 < x cos x"+                            (\(x::UnitInterval a) -> 0 <= x * cos x)+               -- Use <= here because of precision issues :(+             , testProperty "x cos x < sin x"+                            (\(x::UnitInterval a) -> x * cos x <= sin x)+             , testProperty "sin x < x" (\(x::UnitInterval a) -> sin x <= x)+             , testProperty "tangent definition"+                            (\(x::a) -> cos x /= 0 ==> tan x =-= sin x / cos x)+             , testProperty "asin left inverse"+                            (inverseL sin (asin :: BiunitInterval a -> BiunitInterval a))+             , testProperty "acos left inverse"+                            (inverseL cos (acos :: BiunitInterval a -> BiunitInterval a))+             , testProperty "atan left inverse" (inverseL tan (atan :: a -> a))+             , testProperty "sinh definition"+                            (\(x::a) -> sinh x =-= (exp x - exp (-x)) / 2)+             , testProperty "cosh definition"+                            (\(x::a) -> cosh x =-= (exp x + exp (-x)) / 2)+             , testProperty "tanh definition"+                            (\(x::a) -> tanh x =-= sinh x / cosh x)+             , testProperty "sinh left inverse"+                            (inverseL asinh (sinh :: a -> a))+             , testProperty "cosh left inverse"+                            (acosh . cosh =-= (abs :: a -> a))+             , testProperty "tanh left inverse"+                            (inverseL atanh (tanh :: Tiny a -> Tiny a))+             ]++

@@ -14,7 +14,11 @@ import Test.Tasty.QuickCheck (testProperty) import Num (numAuxTests) -fractional :: forall a. (Arbitrary a, EqProp a, Show a, Fractional a, Ord a) => a -> TestTree+-- TODO: Reduce Ord to Eq on the new quickcheck release+-- TODO: Write a program to email me for todo's like that when the conditions+-- are met+fractional :: forall a. (Arbitrary a, EqProp a, Show a, Fractional a, Ord a) =>+              a -> TestTree fractional _ = testGroup "Test Fractional instance" ts   where     ts = [ field "field" (undefined :: a)

@@ -6,29 +6,26 @@  import Test.Tasty (testGroup, TestTree) import Test.Tasty.TH (defaultMainGenerator)-import Test.Tasty.QuickCheck (Arbitrary(..), Positive(..), testProperty, (===), Property, NonNegative(..))+import Test.Tasty.QuickCheck (Positive(..), testProperty, (===), Property)  import Data.CReal.Internal import Data.CReal.Extra () -import Fractional (fractional)+import Floating (floating) import Ord (ord)  -- How many binary digits to use for comparisons TODO: Test with many different -- precisions type Precision = 10 -instance Arbitrary (CReal n) where-  arbitrary = fromInteger <$> arbitrary- infixr 1 ==> (==>) :: Bool -> Bool -> Bool False ==> _ = True True ==> b = b -{-# ANN test_fractional "HLint: ignore Use camelCase" #-}-test_fractional :: [TestTree]-test_fractional = [fractional (undefined :: CReal Precision)]+{-# ANN test_floating "HLint: ignore Use camelCase" #-}+test_floating :: [TestTree]+test_floating = [floating (undefined :: CReal Precision)]  {-# ANN test_ord "HLint: ignore Use camelCase" #-} test_ord :: [TestTree]@@ -42,12 +39,12 @@ prop_showIntegral :: Integer -> Property prop_showIntegral i = show i === show (fromInteger i :: CReal 0) --- TODO: Drop the NonNegative constraint when Floating is implemented and use **-prop_shiftL :: CReal Precision -> NonNegative Int -> Property-prop_shiftL x (NonNegative s) = x `shiftL` s === x * 2^s+prop_shiftL :: CReal Precision -> Int -> Property+prop_shiftL x s = x `shiftL` s === x * 2 ** fromIntegral s -prop_shiftR :: CReal Precision -> NonNegative Int -> Property-prop_shiftR x (NonNegative s) = x `shiftR` s === x / 2^s+prop_shiftR :: CReal Precision -> Int -> Property+prop_shiftR x s = x `shiftR` s === x / 2 ** fromIntegral s  main :: IO () main = $(defaultMainGenerator)+

@@ -64,13 +64,15 @@   where ts = [ring "ring" (undefined :: a),              testProperty "* commutes" (commutes ((*) :: a -> a -> a))] +-- TODO: Reduce the Ord constraint to an Eq constraint on the new quickcheck+-- release field :: forall a. (Arbitrary a, EqProp a, Fractional a, Show a, Ord a) => String -> a -> TestTree field s _ = testGroup s ts   where ts = [abelian "Abelian under Sum" (undefined :: Sum a),               abelian "Abelian under Product NonZero" (undefined :: Product (NonZero a)),               distributes "* distributes over +" (*) ((+) :: a -> a -> a)] -complement :: forall a. (Arbitrary a, EqProp a, Show a, Ord a) =>+complement :: forall a. (Arbitrary a, EqProp a, Show a) =>               String -> (a -> Gen a) -> BinRel a -> BinRel a -> TestTree complement s gen r1 r2 = testGroup s ts   where ts = [testProperty "strictOrd"

@@ -12,15 +12,25 @@   , (<=>)   ) where -import Test.QuickCheck (Arbitrary(..), choose, suchThat)+import Test.QuickCheck import Test.QuickCheck.Checkers (EqProp)-import Test.QuickCheck.Modifiers (NonZero(..))+import Test.QuickCheck.Modifiers (NonZero(..), Positive(..)) import System.Random (Random)  deriving instance Num a => Num (NonZero a) deriving instance Fractional a => Fractional (NonZero a) deriving instance EqProp a => EqProp (NonZero a) +deriving instance Num a => Num (Positive a)+deriving instance Fractional a => Fractional (Positive a)+deriving instance Floating a => Floating (Positive a)+deriving instance EqProp a => EqProp (Positive a)++deriving instance Num a => Num (NonNegative a)+deriving instance Fractional a => Fractional (NonNegative a)+deriving instance Floating a => Floating (NonNegative a)+deriving instance EqProp a => EqProp (NonNegative a)+ newtype UnitInterval a = UnitInterval a   deriving(Eq, Ord, Show, Read, Num, Integral, Fractional, Floating, Real, Enum, Functor, Random, EqProp) @@ -36,7 +46,7 @@   shrink (BiunitInterval a) = BiunitInterval <$> shrink a  newtype Tiny a = Tiny a-  deriving(Eq, Ord, Show, Read, Num, Integral, Real, Enum, Functor)+  deriving(Eq, Ord, Show, Read, Num, Integral, Fractional, Floating, Real, Enum, Functor, Random, EqProp)  -- | Chosen rather arbitrarily just so the tests involving exponentiation don't take too long tinyBound :: Num a => a