diff --git a/exact-real.cabal b/exact-real.cabal --- a/exact-real.cabal +++ b/exact-real.cabal @@ -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, diff --git a/src/Data/CReal.hs b/src/Data/CReal.hs --- a/src/Data/CReal.hs +++ b/src/Data/CReal.hs @@ -1,3 +1,7 @@ +----------------------------------------------------------------------------- +-- | This module exports everything you need to use exact real numbers +---------------------------------------------------------------------------- + module Data.CReal ( CReal , atPrecision 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 @@ -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) diff --git a/test/Data/CReal/Extra.hs b/test/Data/CReal/Extra.hs --- a/test/Data/CReal/Extra.hs +++ b/test/Data/CReal/Extra.hs @@ -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) diff --git a/test/Floating.hs b/test/Floating.hs new file mode 100644 --- /dev/null +++ b/test/Floating.hs @@ -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)) + ] + + diff --git a/test/Fractional.hs b/test/Fractional.hs --- a/test/Fractional.hs +++ b/test/Fractional.hs @@ -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) diff --git a/test/Test.hs b/test/Test.hs --- a/test/Test.hs +++ b/test/Test.hs @@ -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) + diff --git a/test/Test/QuickCheck/Classes/Extra.hs b/test/Test/QuickCheck/Classes/Extra.hs --- a/test/Test/QuickCheck/Classes/Extra.hs +++ b/test/Test/QuickCheck/Classes/Extra.hs @@ -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" diff --git a/test/Test/QuickCheck/Extra.hs b/test/Test/QuickCheck/Extra.hs --- a/test/Test/QuickCheck/Extra.hs +++ b/test/Test/QuickCheck/Extra.hs @@ -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