diff --git a/aern2-mp.cabal b/aern2-mp.cabal --- a/aern2-mp.cabal +++ b/aern2-mp.cabal @@ -4,10 +4,10 @@ -- -- see: https://github.com/sol/hpack -- --- hash: 1b1c633ecdc7ad7935778819abf62906d7d1002d0942466df287afb79ee1e593 +-- hash: b173a95426982332154ff9bd4b378842b3f7e36987e6c3cd6105fa8b5381f8c2 name: aern2-mp -version: 0.2.1.0 +version: 0.2.1.1 synopsis: Multi-precision ball (interval) arithmetic description: Please see the README on GitHub at category: Math diff --git a/src/AERN2/MP/Float/Conversions.hs b/src/AERN2/MP/Float/Conversions.hs --- a/src/AERN2/MP/Float/Conversions.hs +++ b/src/AERN2/MP/Float/Conversions.hs @@ -127,7 +127,16 @@ {- min, max -} -instance CanMinMaxAsymmetric MPFloat MPFloat +instance CanMinMaxAsymmetric MPFloat MPFloat where + type MinMaxType MPFloat MPFloat = MPFloat + max x y + | isNaN x = x + | isNaN y = y + | otherwise = P.max x y + min x y + | isNaN x = x + | isNaN y = y + | otherwise = P.min x y {- constants -} diff --git a/src/AERN2/MP/Float/Tests.hs b/src/AERN2/MP/Float/Tests.hs --- a/src/AERN2/MP/Float/Tests.hs +++ b/src/AERN2/MP/Float/Tests.hs @@ -14,6 +14,12 @@ To run the tests using stack, execute: @ + stack test aern2-mp --test-arguments "-a 1000" + @ + + Or to run only MPFloat tests, execute: + + @ stack test aern2-mp --test-arguments "-a 1000 -m MPFloat" @ -} @@ -30,7 +36,7 @@ -- import qualified Prelude as P -- import Data.Ratio import Text.Printf --- import Data.Maybe +import Data.Maybe ( catMaybes ) import Test.Hspec import Test.QuickCheck @@ -43,7 +49,7 @@ import AERN2.MP.Float.Auxi import AERN2.MP.Float.Type --- import AERN2.MP.Float.Arithmetic +import AERN2.MP.Float.Arithmetic import AERN2.MP.Float.Conversions import AERN2.MP.Float.Operators @@ -147,17 +153,17 @@ abs (x -. y) <= 0.5^e {-| - Assert equality of two MPFloat's with tolerance derived from the size and precision - of the given list of input and intermediate values. - The result is expected to have at least as many significant digits - as the (highest) nominal precision of the input and intermediate numbers + Assert equality of two MPFloat's with tolerance derived from the given list of input and intermediate values + and their differentials (ie the derivative of the result wrt this value). + The difference in the left and right results is expected to share at least as many significant digits + as the worst error step among the input and intermediate numbers shifted by the slope of the differential minus the given precision loss parameter. When the assertion fails, report the given values using the given names. -} approxEqualWithArgs :: Integer {-^ bits of extra precision loss allowed -} -> - [(MPFloat, String)] {-^ intermediate values from which to determine tolerance, their names to report when the equality fails -} -> + [(MPFloat, MPFloat, String)] {-^ intermediate values from which to determine tolerance, their names to report when the equality fails -} -> MPFloat {-^ LHS of equation-} -> MPFloat {-^ RHS of equation -}-> Property @@ -168,31 +174,25 @@ printf "args:\n%s tolerance: <= 2^(%d)" argsS (-e) argsS = unlines - [printf " %s = %s (p=%s)" argS (show arg) (show $ getPrecision arg) - | (arg, argS) <- args ++ [(l, "L"), (r, "R"), (abs(r-.l), "|R-L|")] + [printf " %s = %s ~ %s bits (dR/d%s = %s ~ %s bits)" + argS (show arg) (show $ getErrorStepSizeLog arg) + argS (showMPFloat argD) (show $ getNormLog argD) + | (arg, argD, argS) <- argsLR ] - - e = p - resNorm - precLoss - resNorm = - case (getNormLog l, getNormLog r) of - (NormBits nl, NormBits nr) -> nl `max` nr; - (NormBits nl, _) -> nl - (_, NormBits nr) -> nr - _ -> 0 - p = foldl max 2 $ map (integer . getPrecision . fst) args - - {- - args = argsPre ++ [(l, "L"), (r, "R"), (abs (l-.r),"|L-R|")] - e = - (foldl min 1000000 $ catMaybes $ map getAbsPrecBits args) - - (length argsPre) - getAbsPrecBits (x,_) = - case getNormLog x of - NormZero -> Nothing -- ideally infinity - NormBits b -> Just (pI-b-precLoss) - where - pI = integer $ getPrecision x - -} + argsLR = args ++ [(l, one, "L"), (r, one, "R"), (abs(r-.l), zero, "|R-L|")] + e = 0 - maxStepLoss - precLoss + maxStepLoss = sum $ map (max 1) $ (catMaybes $ map stepLoss argsLR) + stepLoss (arg, argD, _argS) = + case (getNormLog argD, getErrorStepSizeLog arg) of + (NormBits dbits, Just stepBits) -> Just $ dbits + stepBits + _ -> Nothing + -- resNorm = + -- case (getNormLog l, getNormLog r) of + -- (NormBits nl, NormBits nr) -> nl `max` nr; + -- (NormBits nl, _) -> nl + -- (_, NormBits nr) -> nr + -- _ -> 0 + -- p = foldl max 2 $ map (integer . getPrecision . fst) args {-| A runtime representative of type @MPFloat@. @@ -227,7 +227,7 @@ in describe ("MPFloat") $ do specCanSetPrecision tMPFloat - (printArgsIfFails2 "=~=" (\xPrec x -> approxEqualWithArgs 1 [(xPrec, "xPrec")] x xPrec)) + (printArgsIfFails2 "=~=" (\xPrec x -> approxEqualWithArgs 1 [(xPrec, one, "xPrec")] x xPrec)) specCanRound tMPFloat specCanNegNum tMPFloat specCanAbs tMPFloat @@ -253,7 +253,7 @@ it "up ~ down" $ do property $ \ (x :: MPFloat) (y :: MPFloat) -> let - (=~~=) = approxEqualWithArgs 1 [(x,"x"), (y,"y")] + (=~~=) = approxEqualWithArgs 1 [(x,one,"x"), (y,one,"y")] infix 4 =~~= in x +. y =~~= x +^ y @@ -278,7 +278,7 @@ it "up ~ down" $ do property $ \ (x :: MPFloat) (y :: MPFloat) -> let - (=~~=) = approxEqualWithArgs 1 [(x,"x"), (y,"y")] + (=~~=) = approxEqualWithArgs 1 [(x,one,"x"), (y,one,"y")] infix 4 =~~= in x -. y =~~= x -^ y @@ -294,7 +294,7 @@ it "up ~ down" $ do property $ \ (x :: MPFloat) (y :: MPFloat) -> let - (=~~=) = approxEqualWithArgs 1 [(x,"x"), (y,"y")] + (=~~=) = approxEqualWithArgs 1 [(x,y,"x"), (y,x,"y")] infix 4 =~~= in x *. y =~~= x *^ y @@ -331,7 +331,7 @@ it "up ~ down" $ do property $ \ (x :: MPFloat) (y :: MPFloat) -> let - (=~~=) = approxEqualWithArgs 10 [(x,"x"), (y,"y"), (x /. y,"x/.y"), (x /^ y,"x/^y")] + (=~~=) = approxEqualWithArgs 10 [(x,one/^y,"x"), (y,-(x/^y)/^y,"y")] infix 4 =~~= in isFinite y && y /= 0 @@ -370,7 +370,7 @@ property $ \ (x_ :: MPFloat) -> let x = enforceRangeMP (Just 0, Nothing) x_ - (=~~=) = approxEqualWithArgs 2 [(x,"x")] + (=~~=) = approxEqualWithArgs 7 [(x,(one/^(sqrtUp x))/^two,"x")] infix 4 =~~= in sqrtDown x =~~= sqrtUp x @@ -381,7 +381,7 @@ it "sqrt(x)^2 ~ x" $ do property $ \ (x_ :: MPFloat) -> let x = enforceRangeMP (Just 0, Nothing) x_ in - (sqrtDown x) *. (sqrtDown x) <=% x + (max zero $ sqrtDown x) *. (max zero $ sqrtDown x) <=% x && (sqrtUp x) *^ (sqrtUp x) >=% x describe "approximate exp" $ do @@ -393,7 +393,7 @@ property $ \ (x_ :: MPFloat) -> let x = enforceRangeMP (Just (-1000000), Just 1000000) x_ in let - (=~~=) = approxEqualWithArgs 3 [(x,"x")] + (=~~=) = approxEqualWithArgs 4 [(x,expUp x,"x")] infix 4 =~~= in expDown x =~~= expUp x @@ -415,15 +415,14 @@ property $ \ (x_ :: MPFloat) -> let x = enforceRangeMP (Just 0, Nothing) x_ in logDown x <=% logUp x - -- TODO: fix accuracy of CDAR mBounds logA x for x near 1 - -- it "up ~ down" $ do - -- property $ \ (x_ :: MPFloat) -> - -- let x = enforceRangeMP (Just 0, Nothing) x_ in - -- let - -- (=~~=) = approxEqualWithArgs 10 [(x,"x")] - -- infix 4 =~~= - -- in - -- logDown x =~~= logUp x + it "up ~ down" $ do + property $ \ (x_ :: MPFloat) -> + let x = enforceRangeMP (Just 0, Nothing) x_ in + let + (=~~=) = approxEqualWithArgs 4 [(x,one/^x,"x")] + infix 4 =~~= + in + logDown x =~~= logUp x it "log(1/x) == -(log x)" $ do property $ \ (x_ :: MPFloat) -> let x = enforceRangeMP (Just 0, Nothing) x_ in @@ -449,22 +448,22 @@ let x = enforceRangeMP (Just (-1000000), Just 1000000) x_ in sinDown x <=% sinUp x -- TODO: fix accuracy of CDAR mBounds sine - -- it "up ~ down" $ do - -- property $ \ (x_ :: MPFloat) -> - -- let x = enforceRangeMP (Just (-1000000), Just 1000000) x_ in - -- let - -- (=~~=) = approxEqualWithArgs 1 [(x,"x")] - -- infix 4 =~~= - -- in - -- sinDown x =~~= sinUp x - -- it "sin(pi/2) ~ 1" $ do - -- property $ \ (p :: Precision) -> - -- let - -- piA = ceduCentre $ piCEDU p - -- (=~~=) = approxEqualWithArgs 1 [(piA,"pi")] - -- infix 4 =~~= - -- in - -- sinUp(piA/.(setPrecision (p+10) $ mpFloat 2)) =~~= one + it "up ~ down" $ do + property $ \ (x_ :: MPFloat) -> + let x = enforceRangeMP (Just (-1000000), Just 1000000) x_ in + let + (=~~=) = approxEqualWithArgs 2 [(x,cosUp x,"x")] + infix 4 =~~= + in + sinDown x =~~= sinUp x + it "sin(pi/2) ~ 1" $ do + property $ \ (p :: Precision) -> + let + piA = ceduCentre $ piCEDU p + (=~~=) = approxEqualWithArgs 2 [] + infix 4 =~~= + in + sinUp(piA/.(setPrecision (p+10) $ mpFloat 2)) =~~= one it "in [-1,1]" $ do property $ \ (x_ :: MPFloat) -> let x = enforceRangeMP (Just (-1000000), Just 1000000) x_ in @@ -476,30 +475,28 @@ property $ \ (x_ :: MPFloat) -> let x = enforceRangeMP (Just (-1000000), Just 1000000) x_ in cosDown x <=% cosUp x - -- TODO: fix accuracy of CDAR mBounds cosine - -- it "up ~ down" $ do - -- property $ \ (x_ :: MPFloat) -> - -- let x = enforceRangeMP (Just (-1000000), Just 1000000) x_ in - -- let - -- (=~~=) = approxEqualWithArgs 1 [(x,"x")] - -- infix 4 =~~= - -- in - -- cosDown x =~~= cosUp x + it "up ~ down" $ do + property $ \ (x_ :: MPFloat) -> + let x = enforceRangeMP (Just (-1000000), Just 1000000) x_ in + let + (=~~=) = approxEqualWithArgs 2 [(x,-(sinUp x),"x")] + infix 4 =~~= + in + cosDown x =~~= cosUp x it "in [-1,1]" $ do property $ \ (x_ :: MPFloat) -> let x = enforceRangeMP (Just (-1000000), Just 1000000) x_ in cosDown x <=% 1 && cosUp x >=% -1 - -- TODO: fix accuracy of CDAR mBounds cosine - -- it "cos(pi)=-1" $ do - -- property $ \ (p :: Precision) -> - -- let - -- piA = ceduCentre $ piCEDU p - -- (=~~=) = approxEqualWithArgs 2 [(piA,"pi")] - -- infix 4 =~~= - -- in - -- cosUp(piA) =~~= (-one) + it "cos(pi)=-1" $ do + property $ \ (p :: Precision) -> + let + piA = ceduCentre $ piCEDU p + (=~~=) = approxEqualWithArgs 2 [] + infix 4 =~~= + in + cosUp(piA) =~~= (-one) it "cos(x)^2 + sin(x)^2 = 1" $ do property $ \ (x_ :: MPFloat) -> let x = enforceRangeMP (Just (-1000000), Just 1000000) x_ in diff --git a/src/AERN2/MP/Float/Type.hs b/src/AERN2/MP/Float/Type.hs --- a/src/AERN2/MP/Float/Type.hs +++ b/src/AERN2/MP/Float/Type.hs @@ -17,6 +17,7 @@ -- * MPFloat numbers and their basic operations MPFloat , showMPFloat + , getErrorStepSizeLog , setPrecisionCEDU , p2cdarPrec , getBoundsCEDU @@ -67,7 +68,6 @@ getPrecision (MPLow.Approx mb _ _ _) = prec (P.toInteger $ mb) getPrecision MPLow.Bottom = error "illegal MPFloat (Bottom)" - instance CanSetPrecision MPFloat where setPrecision p = ceduCentre . setPrecisionCEDU p @@ -77,3 +77,11 @@ instance HasNorm MPFloat where getNormLog (MPLow.Approx _ m _ s) = (getNormLog m) + (integer s) getNormLog MPLow.Bottom = error "getNormLog undefined for Bottom" + +{-| + Returns @s@ such that @2^s@ is the distance to the nearest other number with the same precision. + Returns Nothing for Bottom. +-} +getErrorStepSizeLog :: MPLow.Approx -> Maybe Int +getErrorStepSizeLog (MPLow.Approx _ _ _ s) = Just $ s +getErrorStepSizeLog _ = Nothing -- represents +Infinity