diff --git a/aern2-mp.cabal b/aern2-mp.cabal --- a/aern2-mp.cabal +++ b/aern2-mp.cabal @@ -1,36 +1,25 @@ name: aern2-mp -version: 0.1.2.0 +version: 0.1.3.0 cabal-version: >= 1.9.2 build-type: Simple homepage: https://github.com/michalkonecny/aern2 author: Michal Konecny maintainer: Michal Konecny -copyright: (c) 2015-2017 Michal Konecny +copyright: (c) 2015-2018 Michal Konecny license: BSD3 license-file: LICENSE extra-source-files: changelog.md stability: experimental -tested-with: GHC==7.10.3, GHC==8.0.2 category: Math -synopsis: Multi-precision floats via MPFR +synopsis: Multi-precision floats via rounded (MPFR) Description: - Arbitrary-precision dyadic numbers and interval arithmetic, via the types: + Arbitrary-precision dyadic numbers and safely-rounded interval arithmetic: . * Dyadic: variable-precision floats with exact ring operations . * MPBall: float ± error bound with field & elementary (interval-like) operations . - By default, using hmpfr when compiling with ghc >=7.10. Older versions of ghc are no longer tested. - . - Set flag MPFRRounded to bind MPFR not via hmpfr but via a tweaked version of package - rounded: - . - * https://github.com/michalkonecny/rounded - . - This backend is considerably faster than hmpfr. - It is likely to become the default in future. - . - There is a plan to add an Integer-only backend so that this library can + There is a plan to add an Integer-only backend so that aern2-mp can be used without MPFR. source-repository head @@ -38,13 +27,8 @@ location: https://github.com/mikkonecny/aern2.git subdir: aern2-mp -flag MPFR - Description: Use a MPFR backend (default); Integer backend is planned - Default: True - -flag MPFRRounded - Description: Use rounded (Numeric.RoundedSimple) as MPFR backend (not default) - -- Default: True +flag UseCDAR + Description: Use an integer-only backend (work in progress, not default) Default: False library @@ -60,26 +44,14 @@ , lens , template-haskell , mixed-types-num - if flag(MPFR) - cpp-options: -DUSEMPFR - if impl(ghc >= 7.10) - if flag(MPFRRounded) - build-depends: - rounded == 0.1.* - cpp-options: -DMPFRRounded - else - build-depends: - hmpfr == 0.4.*, - integer-gmp >= 1.0 && < 1.1 - cpp-options: -DHMPFR - else - build-depends: - haskell-mpfr == 0.1.*, - integer-gmp < 0.6 - cpp-options: -DHaskellMPFR + if flag(UseCDAR) + cpp-options: -DUseCDAR + cdar else - cpp-options: -DUseIReal - -- TODO + cpp-options: -DMPFRRounded + build-depends: + rounded == 0.1.* +-- TODO ghc-options: -Wall -fno-warn-orphans extensions: RebindableSyntax, @@ -95,25 +67,14 @@ FlexibleContexts, FlexibleInstances, UndecidableInstances - if flag(MPFR) + if flag(UseCDAR) exposed-modules: - AERN2.MP.UseMPFR.Float.Type - AERN2.MP.UseMPFR.Float.RoundedAdaptor - AERN2.MP.UseMPFR.Float.Arithmetic - AERN2.MP.UseMPFR.Float.Conversions - AERN2.MP.UseMPFR.Float.Operators - AERN2.MP.UseMPFR.Float.Constants - AERN2.MP.UseMPFR.Float.Tests - AERN2.MP.UseMPFR.Float - AERN2.MP.UseMPFR.Dyadic - AERN2.MP.UseMPFR.ErrorBound - AERN2.MP.UseMPFR.Ball.Type - AERN2.MP.UseMPFR.Ball.Conversions - AERN2.MP.UseMPFR.Ball.Comparisons - AERN2.MP.UseMPFR.Ball.Field - AERN2.MP.UseMPFR.Ball.Elementary - AERN2.MP.UseMPFR.Ball.PreludeOps - AERN2.MP.UseMPFR.Ball + else + exposed-modules: + AERN2.MP.Float.UseRounded.Type + AERN2.MP.Float.UseRounded.RoundedAdaptor + AERN2.MP.Float.UseRounded.Arithmetic + AERN2.MP.Float.UseRounded.Conversions exposed-modules: AERN2.Utils.Bench AERN2.Normalize @@ -121,10 +82,20 @@ AERN2.MP.Precision AERN2.MP.Accuracy AERN2.MP.Enclosure - AERN2.MP.Dyadic AERN2.MP.ErrorBound - AERN2.MP.Ball.Tests + AERN2.MP.Float.Operators + AERN2.MP.Float.Constants + AERN2.MP.Float.Tests + AERN2.MP.Float + AERN2.MP.Dyadic + AERN2.MP.Ball.Type + AERN2.MP.Ball.Conversions + AERN2.MP.Ball.Comparisons + AERN2.MP.Ball.Field + AERN2.MP.Ball.Elementary + AERN2.MP.Ball.PreludeOps AERN2.MP.Ball + AERN2.MP.Ball.Tests AERN2.MP test-suite spec @@ -144,7 +115,8 @@ other-modules: AERN2.MP.BallSpec AERN2.MP.DyadicSpec - AERN2.MP.UseMPFR.FloatSpec + AERN2.MP.FloatSpec + -- AERN2.MP.UseCDAR.FloatSpec build-depends: base == 4.* -- , mixed-types-num >= 0.3.1 && < 0.4 diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -1,3 +1,8 @@ +# Change log for aern2-mp + +* v 0.1.3.0 2018-11-20 + * only one MPFR backend - rounded + * reduce backend-specific code * v 0.1.2.0 2017-11-14 * fix compilation with haskell-mpfr * v 0.1.1.0 2017-11-14 diff --git a/src/AERN2/MP.hs b/src/AERN2/MP.hs --- a/src/AERN2/MP.hs +++ b/src/AERN2/MP.hs @@ -15,7 +15,7 @@ ( module AERN2.Norm , module AERN2.MP.Precision , module AERN2.MP.Accuracy -, module AERN2.MP.UseMPFR.ErrorBound +, module AERN2.MP.ErrorBound , module AERN2.MP.Enclosure , MPBall(..), CanBeMPBall, mpBall, CanBeMPBallP, mpBallP ) @@ -28,6 +28,6 @@ import AERN2.Norm import AERN2.MP.Precision import AERN2.MP.Accuracy -import AERN2.MP.UseMPFR.ErrorBound +import AERN2.MP.ErrorBound import AERN2.MP.Enclosure -import AERN2.MP.UseMPFR.Ball +import AERN2.MP.Ball diff --git a/src/AERN2/MP/Ball.hs b/src/AERN2/MP/Ball.hs --- a/src/AERN2/MP/Ball.hs +++ b/src/AERN2/MP/Ball.hs @@ -1,4 +1,3 @@ -{-# LANGUAGE CPP #-} {-| Module : AERN2.MP.Ball Description : Arbitrary precision ball arithmetic @@ -11,18 +10,55 @@ Arbitrary precision ball arithmetic -} - module AERN2.MP.Ball - ( -#ifdef USEMPFR - module AERN2.MP.UseMPFR.Ball -#endif - ) +( + -- * Auxiliary types + module AERN2.Norm + , module AERN2.MP.Precision + , module AERN2.MP.Accuracy + , module AERN2.MP.ErrorBound + , module AERN2.MP.Enclosure + -- * The Ball type + , MPBall(..), CanBeMPBall, mpBall, CanBeMPBallP, mpBallP + , reducePrecionIfInaccurate + -- * Ball construction/extraction functions + -- , endpointsMP, fromEndpointsMP + -- * Ball operations (see also instances) + , piBallP + -- * Helpers for constructing ball functions + , byEndpointsMP + , fromApproxWithLipschitz +) where --- import MixedTypesNumPrelude +import MixedTypesNumPrelude -- import qualified Prelude as P -#ifdef USEMPFR -import AERN2.MP.UseMPFR.Ball -#endif +import AERN2.Norm +import AERN2.MP.Precision +import AERN2.MP.Accuracy +import AERN2.MP.Enclosure + +import AERN2.MP.ErrorBound (ErrorBound, CanBeErrorBound, errorBound) + +import AERN2.MP.Ball.Type +import AERN2.MP.Ball.Conversions () +import AERN2.MP.Ball.Comparisons +import AERN2.MP.Ball.Field () +import AERN2.MP.Ball.Elementary +import AERN2.MP.Ball.PreludeOps () + +instance Ring MPBall +instance Ring (CN MPBall) +instance Field MPBall +instance Field (CN MPBall) + +instance OrderedRing MPBall +instance OrderedRing (CN MPBall) +instance OrderedField MPBall +instance OrderedField (CN MPBall) + +instance OrderedCertainlyRing MPBall +instance OrderedCertainlyRing (CN MPBall) +instance OrderedCertainlyField MPBall +instance OrderedCertainlyField (CN MPBall) diff --git a/src/AERN2/MP/Ball/Comparisons.hs b/src/AERN2/MP/Ball/Comparisons.hs new file mode 100644 --- /dev/null +++ b/src/AERN2/MP/Ball/Comparisons.hs @@ -0,0 +1,381 @@ +{-| + Module : AERN2.MP.Ball.Comparisons + Description : Comparisons of arbitrary precision dyadic balls + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Comparisons of arbitrary precision dyadic balls +-} +module AERN2.MP.Ball.Comparisons +( + -- * Auxiliary types + module AERN2.Norm + -- * Ball operations (see also instances) + , reducePrecionIfInaccurate + -- * Helpers for constructing ball functions + , byEndpointsMP +) +where + +import MixedTypesNumPrelude +-- import qualified Prelude as P + +import Control.CollectErrors + +import AERN2.Norm +import AERN2.MP.Dyadic (Dyadic) +import AERN2.MP.Float (MPFloat) +-- import AERN2.MP.Float.Operators +import AERN2.MP.Precision + +import AERN2.MP.Ball.Type +import AERN2.MP.Ball.Conversions () + +{- comparisons -} + +instance HasEqAsymmetric MPBall MPBall where + type EqCompareType MPBall MPBall = Maybe Bool + b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 + +instance HasEqAsymmetric MPBall Integer where + type EqCompareType MPBall Integer = Maybe Bool + b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 +instance HasEqAsymmetric Integer MPBall where + type EqCompareType Integer MPBall = Maybe Bool + b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 + +instance HasEqAsymmetric MPBall Int where + type EqCompareType MPBall Int = Maybe Bool + b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 +instance HasEqAsymmetric Int MPBall where + type EqCompareType Int MPBall = Maybe Bool + b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 + +instance HasEqAsymmetric MPBall Rational where + type EqCompareType MPBall Rational = Maybe Bool + b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 +instance HasEqAsymmetric Rational MPBall where + type EqCompareType Rational MPBall = Maybe Bool + b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 + +instance HasEqAsymmetric MPBall Dyadic where + type EqCompareType MPBall Dyadic = Maybe Bool + b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 +instance HasEqAsymmetric Dyadic MPBall where + type EqCompareType Dyadic MPBall = Maybe Bool + b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 + +instance + (HasEqAsymmetric MPBall b + , CanEnsureCE es b + , CanEnsureCE es (EqCompareType MPBall b) + , IsBool (EnsureCE es (EqCompareType MPBall b)) + , SuitableForCE es) + => + HasEqAsymmetric MPBall (CollectErrors es b) + where + type EqCompareType MPBall (CollectErrors es b) = + EnsureCE es (EqCompareType MPBall b) + equalTo = lift2TLCE equalTo + +instance + (HasEqAsymmetric a MPBall + , CanEnsureCE es a + , CanEnsureCE es (EqCompareType a MPBall) + , IsBool (EnsureCE es (EqCompareType a MPBall)) + , SuitableForCE es) + => + HasEqAsymmetric (CollectErrors es a) MPBall + where + type EqCompareType (CollectErrors es a) MPBall = + EnsureCE es (EqCompareType a MPBall) + equalTo = lift2TCE equalTo + +instance HasOrderAsymmetric MPBall MPBall where + type OrderCompareType MPBall MPBall = Maybe Bool + lessThan b1 b2 + | r1 < l2 = Just True + | r2 <= l1 = Just False + | otherwise = Nothing + where + (l1, r1) = endpointsMP b1 + (l2, r2) = endpointsMP b2 + leq b1 b2 + | r1 <= l2 = Just True + | r2 < l1 = Just False + | otherwise = Nothing + where + (l1, r1) = endpointsMP b1 + (l2, r2) = endpointsMP b2 + +instance HasOrderAsymmetric Integer MPBall where + type OrderCompareType Integer MPBall = Maybe Bool + lessThan = convertFirst lessThan + leq = convertFirst leq +instance HasOrderAsymmetric MPBall Integer where + type OrderCompareType MPBall Integer = Maybe Bool + lessThan = convertSecond lessThan + leq = convertSecond leq + +instance HasOrderAsymmetric Int MPBall where + type OrderCompareType Int MPBall = Maybe Bool + lessThan = convertFirst lessThan + leq = convertFirst leq +instance HasOrderAsymmetric MPBall Int where + type OrderCompareType MPBall Int = Maybe Bool + lessThan = convertSecond lessThan + leq = convertSecond leq + +instance HasOrderAsymmetric Dyadic MPBall where + type OrderCompareType Dyadic MPBall = Maybe Bool + lessThan = convertFirst lessThan + leq = convertFirst leq +instance HasOrderAsymmetric MPBall Dyadic where + type OrderCompareType MPBall Dyadic = Maybe Bool + lessThan = convertSecond lessThan + leq = convertSecond leq + +instance HasOrderAsymmetric MPBall Rational where + type OrderCompareType MPBall Rational = Maybe Bool + lessThan b1 q2 + | r1 < l2 = Just True + | r2 <= l1 = Just False + | otherwise = Nothing + where + (l1, r1) = endpointsMP b1 + l2 = q2 + r2 = q2 + leq b1 q2 + | r1 <= l2 = Just True + | r2 < l1 = Just False + | otherwise = Nothing + where + (l1, r1) = endpointsMP b1 + l2 = q2 + r2 = q2 + +instance HasOrderAsymmetric Rational MPBall where + type OrderCompareType Rational MPBall = Maybe Bool + lessThan q1 b2 + | r1 < l2 = Just True + | r2 <= l1 = Just False + | otherwise = Nothing + where + (l2, r2) = endpointsMP b2 + l1 = q1 + r1 = q1 + leq q1 b2 + | r1 <= l2 = Just True + | r2 < l1 = Just False + | otherwise = Nothing + where + (l2, r2) = endpointsMP b2 + l1 = q1 + r1 = q1 + +instance + (HasOrderAsymmetric MPBall b + , CanEnsureCE es b + , CanEnsureCE es (OrderCompareType MPBall b) + , IsBool (EnsureCE es (OrderCompareType MPBall b)) + , SuitableForCE es) + => + HasOrderAsymmetric MPBall (CollectErrors es b) + where + type OrderCompareType MPBall (CollectErrors es b) = + EnsureCE es (OrderCompareType MPBall b) + lessThan = lift2TLCE lessThan + leq = lift2TLCE leq + greaterThan = lift2TLCE greaterThan + geq = lift2TLCE geq + +instance + (HasOrderAsymmetric a MPBall + , CanEnsureCE es a + , CanEnsureCE es (OrderCompareType a MPBall) + , IsBool (EnsureCE es (OrderCompareType a MPBall)) + , SuitableForCE es) + => + HasOrderAsymmetric (CollectErrors es a) MPBall + where + type OrderCompareType (CollectErrors es a) MPBall = + EnsureCE es (OrderCompareType a MPBall) + lessThan = lift2TCE lessThan + leq = lift2TCE leq + greaterThan = lift2TCE greaterThan + geq = lift2TCE geq + +instance CanTestZero MPBall +instance CanTestPosNeg MPBall + +instance CanTestInteger MPBall where + certainlyNotInteger b = + (rN - lN) == 1 && lN ! + CanMinMaxAsymmetric MPBall (CollectErrors es b) + where + type MinMaxType MPBall (CollectErrors es b) = + EnsureCE es (MinMaxType MPBall b) + min = lift2TLCE min + max = lift2TLCE max + +instance + (CanMinMaxAsymmetric a MPBall + , CanEnsureCE es a + , CanEnsureCE es (MinMaxType a MPBall) + , SuitableForCE es) + => + CanMinMaxAsymmetric (CollectErrors es a) MPBall + where + type MinMaxType (CollectErrors es a) MPBall = + EnsureCE es (MinMaxType a MPBall) + min = lift2TCE min + max = lift2TCE max + +{- intersection -} + +instance CanIntersectAsymmetric MPBall MPBall where + intersect a b + | rL > rR = + noValueNumErrorCertainCN $ NumError $ "intersect: empty intersection: " ++ show a ++ "; " ++ show b + | otherwise = cn $ fromEndpointsMP rL rR + where + rL = max aL bL + rR = min aR bR + (aL,aR) = endpointsMP a + (bL,bR) = endpointsMP b + +instance + (CanIntersectAsymmetric MPBall b + , CanEnsureCE es b + , CanEnsureCE es (IntersectionType MPBall b) + , SuitableForCE es) + => + CanIntersectAsymmetric MPBall (CollectErrors es b) + where + type IntersectionType MPBall (CollectErrors es b) = + EnsureCE es (IntersectionType MPBall b) + intersect = lift2TLCE intersect + +instance + (CanIntersectAsymmetric a MPBall + , CanEnsureCE es a + , CanEnsureCE es (IntersectionType a MPBall) + , SuitableForCE es) + => + CanIntersectAsymmetric (CollectErrors es a) MPBall + where + type IntersectionType (CollectErrors es a) MPBall = + EnsureCE es (IntersectionType a MPBall) + intersect = lift2TCE intersect + +{- union -} + +instance CanUnionAsymmetric MPBall MPBall where + union a b = + case getMaybeValueCN (a `intersect` b) of + Just _ -> prependErrorsCN [(ErrorPotential, err)] r + _ -> prependErrorsCN [(ErrorCertain, err)] r + where + err = NumError $ "union of enclosures: not enclosing the same value" + r = cn $ fromEndpointsMP rL rR + rL = min aL bL + rR = max aR bR + (aL,aR) = endpointsMP a + (bL,bR) = endpointsMP b + + +instance + (CanUnionAsymmetric MPBall b + , CanEnsureCE es b + , CanEnsureCE es (UnionType MPBall b) + , SuitableForCE es) + => + CanUnionAsymmetric MPBall (CollectErrors es b) + where + type UnionType MPBall (CollectErrors es b) = + EnsureCE es (UnionType MPBall b) + union = lift2TLCE union + +instance + (CanUnionAsymmetric a MPBall + , CanEnsureCE es a + , CanEnsureCE es (UnionType a MPBall) + , SuitableForCE es) + => + CanUnionAsymmetric (CollectErrors es a) MPBall + where + type UnionType (CollectErrors es a) MPBall = + EnsureCE es (UnionType a MPBall) + union = lift2TCE union + +{-| + Computes an *increasing* ball fucntion @f@ from *exact* MPFR operations. +-} +byEndpointsMP :: + (MPFloat -> MPFloat -> MPFloat) -> + (MPBall -> MPBall -> MPBall) +byEndpointsMP op b1 b2 = + fromEndpointsMP (l1 `op` l2) (r1 `op` r2) + where + (l1,r1) = endpointsMP b1 + (l2,r2) = endpointsMP b2 + +{- random generation -} diff --git a/src/AERN2/MP/Ball/Conversions.hs b/src/AERN2/MP/Ball/Conversions.hs new file mode 100644 --- /dev/null +++ b/src/AERN2/MP/Ball/Conversions.hs @@ -0,0 +1,119 @@ +{-| + Module : AERN2.MP.Ball.Conversions + Description : Conversions of arbitrary precision dyadic balls + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Conversions of arbitrary precision dyadic balls +-} +module AERN2.MP.Ball.Conversions +( + integerBounds +) +where + +import MixedTypesNumPrelude +-- import qualified Prelude as P + +import Data.Typeable +import Data.Convertible + +import AERN2.MP.Dyadic (Dyadic, dyadic) +import qualified AERN2.MP.Float as MPFloat +import AERN2.MP.Float (mpFloat) +-- import AERN2.MP.Float.Operators +import AERN2.MP.Precision +import qualified AERN2.MP.ErrorBound as EB +import AERN2.MP.ErrorBound (errorBound) + +import AERN2.MP.Ball.Type + +{--- extracting from a ball ---} + +instance HasIntegerBounds MPBall where + integerBounds b = + (floor l, ceiling r) + where + (l,r) = endpointsMP b + +instance Convertible MPBall EB.ErrorBound where + safeConvert b = + Right (errorBound (max (abs l) (abs r))) + where + (l,r) = endpointsMP b + +{--- constructing an exact ball ---} + +instance ConvertibleExactly MPBall MPBall where + safeConvertExactly = Right + +instance ConvertibleExactly Dyadic MPBall where + safeConvertExactly x = Right $ MPBall (mpFloat x) (errorBound 0) + +instance ConvertibleExactly EB.ErrorBound MPBall where + safeConvertExactly eb = Right $ MPBall (mpFloat eb) (errorBound 0) + +instance + (ConvertibleExactly c Dyadic, ConvertibleExactly e Dyadic + , Show c, Show e, Typeable c, Typeable e) + => + ConvertibleExactly (c, e) MPBall + where + safeConvertExactly (c,e) + | isFinite b = Right b + | otherwise = convError "too large to convert to MPBall" (c,e) + where + b = MPBall (mpFloat $ dyadic c) (errorBound $ mpFloat $ dyadic e) + +instance ConvertibleExactly Integer MPBall where + safeConvertExactly x + | isFinite b = Right b + | otherwise = convError "too large to convert to MPBall" x + where + b = MPBall (mpFloat x) (errorBound 0) + +instance ConvertibleExactly Int MPBall where + safeConvertExactly x = Right $ MPBall (mpFloat x) (errorBound 0) + +{--- constructing a ball with a given precision ---} + +instance ConvertibleWithPrecision Integer MPBall where + safeConvertP p x + | isFinite b = Right b + | otherwise = convError ("too large to convert to MPBall with precision " ++ show p) x + where + b = MPBall xUp (xUp `EB.subMP` xDn) + xUp = MPFloat.fromIntegerUp p x + xDn = MPFloat.fromIntegerDown p x + +instance ConvertibleWithPrecision Int MPBall where + safeConvertP p = safeConvertP p . integer + +instance ConvertibleWithPrecision Dyadic MPBall where + safeConvertP p x + | isFinite b = Right b + | otherwise = convError ("too large to convert to MPBall with precision " ++ show p) x + where + b = mpBall x + +instance ConvertibleWithPrecision Rational MPBall where + safeConvertP p x + | isFinite b = Right b + | otherwise = convError ("too large to convert to MPBall with precision " ++ show p) x + where + b = MPBall xUp (xUp `EB.subMP` xDn) + xUp = MPFloat.fromRationalUp p x + xDn = MPFloat.fromRationalDown p x + +instance ConvertibleWithPrecision (Rational, Rational) MPBall where + safeConvertP p (x,e) + | isFinite b = Right b + | otherwise = convError ("too large to convert to MPBall with precision " ++ show p) x + where + b = MPBall xFlt (xe + eUp) -- beware, precision may be too high relative to accuracy + (MPBall xFlt xe) = mpBallP p x + eUp = errorBound e diff --git a/src/AERN2/MP/Ball/Elementary.hs b/src/AERN2/MP/Ball/Elementary.hs new file mode 100644 --- /dev/null +++ b/src/AERN2/MP/Ball/Elementary.hs @@ -0,0 +1,146 @@ +{-| + Module : AERN2.MP.Ball.Elementary + Description : Elementary operations on arbitrary precision dyadic balls + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Elementary operations on arbitrary precision dyadic balls +-} +module AERN2.MP.Ball.Elementary +( + -- * Ball operations (see also instances) + piBallP + -- * Helpers for constructing ball functions + , fromApproxWithLipschitz +) +where + +import MixedTypesNumPrelude +import qualified Prelude as P + +import AERN2.Normalize + +import AERN2.MP.Dyadic (Dyadic) +import qualified AERN2.MP.Float as MPFloat +import AERN2.MP.Float (MPFloat, mpFloat) +-- import AERN2.MP.Float.Operators +import AERN2.MP.Precision +import qualified AERN2.MP.ErrorBound as EB +import AERN2.MP.ErrorBound (errorBound) + +import AERN2.MP.Ball.Type +import AERN2.MP.Ball.Conversions () +import AERN2.MP.Ball.Comparisons () +import AERN2.MP.Ball.Field () + + +{- trigonometrics -} + +piBallP :: Precision -> MPBall +piBallP p = MPBall piUp (piUp `EB.subMP` piDown) + where + piUp = MPFloat.piUp p + piDown = MPFloat.piDown p + +instance CanSinCos MPBall where + sin = sinB 1 + cos = cosB 1 + +sinB :: Integer -> MPBall -> MPBall +sinB i x = + -- increasingPrecisionUntilNotImproving (fromApproxWithLipschitz MPFloat.sinDown MPFloat.sinUp lip) x + fromApproxWithLipschitz MPFloat.sinDown MPFloat.sinUp lip x + where + lip + | i == 0 = mpFloat 1 + | otherwise = snd $ endpointsMP $ abs $ cosB (i - 1) x + +cosB :: Integer -> MPBall -> MPBall +cosB i x = + -- increasingPrecisionUntilNotImproving (fromApproxWithLipschitz MPFloat.cosDown MPFloat.cosUp lip) x + fromApproxWithLipschitz MPFloat.cosDown MPFloat.cosUp lip x + where + lip + | i == 0 = mpFloat 1 + | otherwise = snd $ endpointsMP $ abs $ sinB (i - 1) x + +-- increasingPrecisionUntilNotImproving :: (MPBall -> MPBall) -> (MPBall -> MPBall) +-- increasingPrecisionUntilNotImproving f x = +-- waitUntilNotImproving $ map aux (precisions xPrec (xPrec*2)) +-- where +-- xPrec = getPrecision x +-- precisions p1 p2 = p1 : (precisions p2 (p1 + p2)) +-- aux p = f $ setPrecision p x +-- waitUntilNotImproving xx@(x1:_) = aux2 (getAccuracy x1) xx +-- waitUntilNotImproving _ = error "AERN2.MP.Ball.Elementary: internal error in increasingPrecisionUntilNotImproving" +-- aux2 x1AC (x1:x2:rest) +-- | x1AC < x2AC = aux2 x2AC (x2:rest) +-- | otherwise = x1 +-- where +-- x2AC = getAccuracy x2 +-- aux2 _ _ = error "AERN2.MP.Ball.Elementary: internal error in increasingPrecisionUntilNotImproving" + +{- exp, log, power -} + +instance CanExp MPBall where + exp = intervalFunctionByEndpointsUpDown MPFloat.expDown MPFloat.expUp + +instance CanLog MPBall where + type LogType MPBall = CN MPBall + log x + | x !>! 0 = + cn $ intervalFunctionByEndpointsUpDown MPFloat.logDown MPFloat.logUp x + | x !<=! 0 = noValueNumErrorCertainCN err + | otherwise = noValueNumErrorPotentialCN err + where + err = OutOfRange $ "log: argument must be > 0: " ++ show x + +instance CanPow MPBall MPBall where + powNoCN b e = (~!) $ pow b e + pow = powUsingExpLog (mpBall 0) (mpBall 1) + +instance CanPow MPBall Dyadic where + powNoCN b e = (~!) $ pow b e + pow b e = powUsingExpLog (mpBall 0) (mpBall 1) b (mpBall e) + +instance CanPow MPBall Rational where + powNoCN b e = (~!) $ pow b e + pow b e = powUsingExpLog (mpBall 0) (mpBall 1) b (mpBallP (getPrecision b) e) + +instance CanSqrt MPBall where + type SqrtType MPBall = CN MPBall + sqrt x + | x !>=! 0 = cn $ aux x + | x ! MPFloat.sqrtDown (P.max (mpFloat 0) e)) + (\ e -> MPFloat.sqrtUp (P.max (mpFloat 0) e)) + err = OutOfRange $ "sqrt: argument must be >= 0: " ++ show x + +{- generic methods for computing real functions from MPFR-approximations -} + +{-| + Computes a real function @f@ from correctly rounded MPFR-approximations and a number @lip@ which is a + Lipschitz constant for @f@, i.e. @|f(x) - f(y)| <= lip * |x - y|@ for all @x@,@y@. +-} +fromApproxWithLipschitz :: + (MPFloat -> MPFloat) {-^ @fDown@: a version of @f@ on MPFloat rounding *downwards* -} -> + (MPFloat -> MPFloat) {-^ @fUp@: a version of @f@ on MPFloat rounding *upwards* -} -> + MPFloat {-^ @lip@ a Lipschitz constant for @f@, @lip > 0@ -} -> + (MPBall -> MPBall) {-^ @f@ on MPBall rounding *outwards* -} +fromApproxWithLipschitz fDown fUp lip _x@(MPBall xc xe) = + normalize $ MPBall fxc err + where + fxl = fDown xc + fxu = fUp xc + (MPBall fxc fxe) = + setPrecision (getPrecision xc) $ -- beware, some MPFR functions increase precision, eg sine and cosine + fromEndpointsMP fxl fxu + err = (errorBound lip) * xe + fxe diff --git a/src/AERN2/MP/Ball/Field.hs b/src/AERN2/MP/Ball/Field.hs new file mode 100644 --- /dev/null +++ b/src/AERN2/MP/Ball/Field.hs @@ -0,0 +1,344 @@ +{-# LANGUAGE TemplateHaskell #-} +{-| + Module : AERN2.MP.Ball.Field + Description : Field operations on arbitrary precision dyadic balls + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Field operations on arbitrary precision dyadic balls +-} +module AERN2.MP.Ball.Field +() +where + +import MixedTypesNumPrelude +-- import qualified Prelude as P + +import Control.CollectErrors + +import AERN2.Normalize + +import AERN2.MP.Dyadic (Dyadic) +import AERN2.MP.Float (mpFloat) +import AERN2.MP.Float.Operators +import AERN2.MP.Precision +import qualified AERN2.MP.ErrorBound as EB + +import AERN2.MP.Ball.Type +import AERN2.MP.Ball.Conversions () +import AERN2.MP.Ball.Comparisons () + +{- addition -} + +instance CanAddAsymmetric MPBall MPBall where + type AddType MPBall MPBall = MPBall + add (MPBall x1 e1) (MPBall x2 e2) = + normalize $ MPBall sumUp ((sumUp `EB.subMP` sumDn) + e1 + e2) + where + sumUp = x1 +^ x2 + sumDn = x1 +. x2 + +instance CanAddAsymmetric MPBall Int where + type AddType MPBall Int = MPBall + add = convertSecond add +instance CanAddAsymmetric Int MPBall where + type AddType Int MPBall = MPBall + add = convertFirst add + +instance CanAddAsymmetric MPBall Integer where + type AddType MPBall Integer = MPBall + add = convertSecond add +instance CanAddAsymmetric Integer MPBall where + type AddType Integer MPBall = MPBall + add = convertFirst add + +instance CanAddAsymmetric MPBall Dyadic where + type AddType MPBall Dyadic = MPBall + add = convertSecond add +instance CanAddAsymmetric Dyadic MPBall where + type AddType Dyadic MPBall = MPBall + add = convertFirst add + +instance CanAddAsymmetric MPBall Rational where + type AddType MPBall Rational = MPBall + add = convertPSecond add +instance CanAddAsymmetric Rational MPBall where + type AddType Rational MPBall = MPBall + add = convertPFirst add + +instance + (CanAddAsymmetric MPBall b + , CanEnsureCE es b + , CanEnsureCE es (AddType MPBall b) + , SuitableForCE es) + => + CanAddAsymmetric MPBall (CollectErrors es b) + where + type AddType MPBall (CollectErrors es b) = + EnsureCE es (AddType MPBall b) + add = lift2TLCE add + +instance + (CanAddAsymmetric a MPBall + , CanEnsureCE es a + , CanEnsureCE es (AddType a MPBall) + , SuitableForCE es) + => + CanAddAsymmetric (CollectErrors es a) MPBall + where + type AddType (CollectErrors es a) MPBall = + EnsureCE es (AddType a MPBall) + add = lift2TCE add + +{- subtraction -} + +instance CanSub MPBall MPBall + +instance CanSub MPBall Integer +instance CanSub Integer MPBall + +instance CanSub MPBall Int +instance CanSub Int MPBall + +instance CanSub MPBall Rational +instance CanSub Rational MPBall + +instance CanSub MPBall Dyadic +instance CanSub Dyadic MPBall + +instance + (CanSub MPBall b + , CanEnsureCE es b + , CanEnsureCE es (SubType MPBall b) + , SuitableForCE es) + => + CanSub MPBall (CollectErrors es b) + where + type SubType MPBall (CollectErrors es b) = + EnsureCE es (SubType MPBall b) + sub = lift2TLCE sub + +instance + (CanSub a MPBall + , CanEnsureCE es a + , CanEnsureCE es (SubType a MPBall) + , SuitableForCE es) + => + CanSub (CollectErrors es a) MPBall + where + type SubType (CollectErrors es a) MPBall = + EnsureCE es (SubType a MPBall) + sub = lift2TCE sub + +{- multiplication -} + +instance CanMulAsymmetric MPBall MPBall where + mul (MPBall x1 e1) (MPBall x2 e2) = + normalize $ MPBall x12Up (e12 + e1*(abs x2) + e2*(abs x1) + e1*e2) + -- the mixed operations above automatically convert + -- MPFloat to ErrorBound, checking non-negativity + where + x12Up = x1 *^ x2 + x12Down = x1 *. x2 + e12 = x12Up -^ x12Down + +instance CanMulAsymmetric MPBall Int where + type MulType MPBall Int = MPBall + mul = convertSecond mul +instance CanMulAsymmetric Int MPBall where + type MulType Int MPBall = MPBall + mul = convertFirst mul + +instance CanMulAsymmetric MPBall Integer where + type MulType MPBall Integer = MPBall + mul = convertSecond mul +instance CanMulAsymmetric Integer MPBall where + type MulType Integer MPBall = MPBall + mul = convertFirst mul + +instance CanMulAsymmetric MPBall Dyadic where + type MulType MPBall Dyadic = MPBall + mul = convertSecond mul +instance CanMulAsymmetric Dyadic MPBall where + type MulType Dyadic MPBall = MPBall + mul = convertFirst mul + +instance CanMulAsymmetric MPBall Rational where + type MulType MPBall Rational = MPBall + mul = convertPSecond mul +instance CanMulAsymmetric Rational MPBall where + type MulType Rational MPBall = MPBall + mul = convertPFirst mul + +instance + (CanMulAsymmetric MPBall b + , CanEnsureCE es b + , CanEnsureCE es (MulType MPBall b) + , SuitableForCE es) + => + CanMulAsymmetric MPBall (CollectErrors es b) + where + type MulType MPBall (CollectErrors es b) = + EnsureCE es (MulType MPBall b) + mul = lift2TLCE mul + +instance + (CanMulAsymmetric a MPBall + , CanEnsureCE es a + , CanEnsureCE es (MulType a MPBall) + , SuitableForCE es) + => + CanMulAsymmetric (CollectErrors es a) MPBall + where + type MulType (CollectErrors es a) MPBall = + EnsureCE es (MulType a MPBall) + mul = lift2TCE mul + + +{- division -} + +instance CanDiv MPBall MPBall where + type DivTypeNoCN MPBall MPBall = MPBall + divideNoCN b1 b2 = (~!) (divide b1 b2) + type DivType MPBall MPBall = CN MPBall + divide (MPBall x1 e1) b2@(MPBall x2 e2) + | isCertainlyNonZero b2 = + cn $ normalize $ MPBall x12Up err + | isCertainlyZero b2 = + noValueNumErrorCertainCN DivByZero + | otherwise = + noValueNumErrorPotentialCN DivByZero + where + x12Up = x1 /^ x2 + x12Down = x1 /. x2 + x12AbsUp = (abs x12Up) `max` (abs x12Down) + e12 = x12Up -^ x12Down + err = + ((e12 *^ (abs x2)) -- e12 * |x2| + + + e1 + + + (e2 * x12AbsUp) -- e2 * |x| + ) + * + ((mpFloat 1) /^ ((abs x2) -. (mpFloat e2))) + -- 1/(|x2| - e2) rounded upwards +{- +A derivation of the above formula for an upper bound on the error: + + * e = + * = max ( (x1 ± e1) / (x2 ± e2) - x ) + * = max ( ( x1 ± e1 - (x*(x2 ± e2) ) / (x2 ± e2) ) + * ≤ max ( ( x1 ± e1 - ((x1/x2) ± e12)x2 ± x*e2 ) / (x2 ± e2) ) + * = max ( ( x1 ± e1 - x1 ± e12*x2 ± x*e2 ) / (x2 ± e2) ) + * = max ( ( ± e1 ± e12*x2 ± x*e2 ) / (x2 ± e2) ) + * ≤ (e1 + e12*|x2| + |x|*e2 ) / (|x2| - e2) + * ≤ (e1 +^ e12*^|x2| +^ |x|*^e2 ) /^ (|x2| -. e2) +-} + +$(declForTypes + [[t| Integer |], [t| Int |], [t| Dyadic |]] + (\ t -> [d| + instance CanDiv MPBall $t where + type DivType MPBall $t = CN MPBall + divide = convertSecond divide + type DivTypeNoCN MPBall $t = MPBall + divideNoCN = convertSecond divideNoCN + instance CanDiv $t MPBall where + type DivType $t MPBall = CN MPBall + divide = convertFirst divide + type DivTypeNoCN $t MPBall = MPBall + divideNoCN = convertFirst divideNoCN + |])) + +instance CanDiv Dyadic Dyadic where + type DivTypeNoCN Dyadic Dyadic = MPBall + divideNoCN a b = divideNoCN (mpBall a) (mpBall b) + divide a b = divide (mpBall a) (mpBall b) + +instance CanDiv MPBall Rational where + type DivTypeNoCN MPBall Rational = MPBall + divideNoCN = convertPSecond divideNoCN + divide = convertPSecond divide +instance CanDiv Rational MPBall where + type DivTypeNoCN Rational MPBall = MPBall + divideNoCN = convertPFirst divideNoCN + divide = convertPFirst divide + +instance + (CanDiv MPBall b + , CanEnsureCE es b + , CanEnsureCE es (DivType MPBall b) + , CanEnsureCE es (DivTypeNoCN MPBall b) + , SuitableForCE es) + => + CanDiv MPBall (CollectErrors es b) + where + type DivType MPBall (CollectErrors es b) = + EnsureCE es (DivType MPBall b) + divide = lift2TLCE divide + type DivTypeNoCN MPBall (CollectErrors es b) = + EnsureCE es (DivTypeNoCN MPBall b) + divideNoCN = lift2TLCE divideNoCN + +instance + (CanDiv a MPBall + , CanEnsureCE es a + , CanEnsureCE es (DivType a MPBall) + , CanEnsureCE es (DivTypeNoCN a MPBall) + , SuitableForCE es) + => + CanDiv (CollectErrors es a) MPBall + where + type DivType (CollectErrors es a) MPBall = + EnsureCE es (DivType a MPBall) + divide = lift2TCE divide + type DivTypeNoCN (CollectErrors es a) MPBall = + EnsureCE es (DivTypeNoCN a MPBall) + divideNoCN = lift2TCE divideNoCN + +{- integer power -} + +instance CanPow MPBall Integer where + powNoCN b e = (~!) $ powUsingMulRecip (mpBall 1) b e + pow = powUsingMulRecip (mpBall 1) + +instance CanPow MPBall Int where + powNoCN b e = (~!) $ powUsingMulRecip (mpBall 1) b e + pow = powUsingMulRecip (mpBall 1) + +instance + (CanPow MPBall b + , CanEnsureCE es b + , CanEnsureCE es (PowType MPBall b) + , CanEnsureCE es (PowTypeNoCN MPBall b) + , SuitableForCE es) + => + CanPow MPBall (CollectErrors es b) + where + type PowTypeNoCN MPBall (CollectErrors es b) = + EnsureCE es (PowTypeNoCN MPBall b) + powNoCN = lift2TLCE powNoCN + type PowType MPBall (CollectErrors es b) = + EnsureCE es (PowType MPBall b) + pow = lift2TLCE pow + +instance + (CanPow a MPBall + , CanEnsureCE es a + , CanEnsureCE es (PowType a MPBall) + , CanEnsureCE es (PowTypeNoCN a MPBall) + , SuitableForCE es) + => + CanPow (CollectErrors es a) MPBall + where + type PowTypeNoCN (CollectErrors es a) MPBall = + EnsureCE es (PowTypeNoCN a MPBall) + powNoCN = lift2TCE powNoCN + type PowType (CollectErrors es a) MPBall = + EnsureCE es (PowType a MPBall) + pow = lift2TCE pow diff --git a/src/AERN2/MP/Ball/PreludeOps.hs b/src/AERN2/MP/Ball/PreludeOps.hs new file mode 100644 --- /dev/null +++ b/src/AERN2/MP/Ball/PreludeOps.hs @@ -0,0 +1,76 @@ +{-# LANGUAGE CPP #-} +{-| + Module : AERN2.MP.Ball.PreludeOps + Description : Instances of Prelude.Num etc + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Instances of Prelude classes Eq, Ord, Num etc +-} +module AERN2.MP.Ball.PreludeOps +( +) +where + +import MixedTypesNumPrelude +import qualified Prelude as P + +import AERN2.MP.Dyadic (dyadic) + +import AERN2.MP.Ball.Type +import AERN2.MP.Ball.Conversions () +import AERN2.MP.Ball.Comparisons () +import AERN2.MP.Ball.Field () +import AERN2.MP.Ball.Elementary () + +{- Instances of Prelude numerical classes provided for convenient use outside AERN2 + and also because Template Haskell translates (-x) to (Prelude.negate x) -} + +instance P.Eq MPBall where + a == b = + case a == b of + Just t -> t + _ -> + error "Failed to decide equality of MPBalls. If you switch to MixedTypesNumPrelude instead of Prelude, comparison of MPBalls returns Maybe Bool instead of Bool." + +instance P.Ord MPBall where + compare a b = + case (a < b, a == b, a > b) of + (Just True, _, _) -> P.LT + (_, Just True, _) -> P.EQ + (_, _, Just True) -> P.GT + _ -> + error "Failed to decide order of MPBalls. If you switch to MixedTypesNumPrelude instead of Prelude, comparison of MPBalls returns Maybe Bool instead of Bool." + +instance P.Num MPBall where + fromInteger = convertExactly + negate = negate + (+) = (+) + (*) = (*) + abs = abs + signum = error "Prelude.signum not implemented for MPBall" + +instance P.Fractional MPBall where + fromRational = convertExactly . dyadic -- will work only for dyadic rationals + recip = (~!) . recip + (/) = (/!) + +instance P.Floating MPBall where + pi = error "MPBall: no pi :: MPBall, use pi ? (bitsS n) instead" + sqrt = (~!) . sqrt + exp = exp + sin = sin + cos = cos + log = (~!) . log + atan = error "MPBall: atan not implemented yet" + atanh = error "MPBall: atanh not implemented yet" + asin = error "MPBall: asin not implemented yet" + acos = error "MPBall: acos not implemented yet" + sinh = error "MPBall: sinh not implemented yet" + cosh = error "MPBall: cosh not implemented yet" + asinh = error "MPBall: asinh not implemented yet" + acosh = error "MPBall: acosh not implemented yet" diff --git a/src/AERN2/MP/Ball/Tests.hs b/src/AERN2/MP/Ball/Tests.hs --- a/src/AERN2/MP/Ball/Tests.hs +++ b/src/AERN2/MP/Ball/Tests.hs @@ -36,11 +36,11 @@ -- import AERN2.Norm import AERN2.MP.Precision -import AERN2.MP.UseMPFR.Ball.Type --- import AERN2.MP.UseMPFR.Ball.Conversions () -import AERN2.MP.UseMPFR.Ball.Comparisons () -import AERN2.MP.UseMPFR.Ball.Field () -import AERN2.MP.UseMPFR.Ball.Elementary () +import AERN2.MP.Ball.Type +-- import AERN2.MP.Ball.Conversions () +import AERN2.MP.Ball.Comparisons () +import AERN2.MP.Ball.Field () +import AERN2.MP.Ball.Elementary () instance Arbitrary MPBall where arbitrary = diff --git a/src/AERN2/MP/Ball/Type.hs b/src/AERN2/MP/Ball/Type.hs new file mode 100644 --- /dev/null +++ b/src/AERN2/MP/Ball/Type.hs @@ -0,0 +1,276 @@ +{-# LANGUAGE TemplateHaskell #-} +{-| + Module : AERN2.MP.Ball.Type + Description : Arbitrary precision dyadic balls + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Arbitrary precision dyadic balls +-} +module AERN2.MP.Ball.Type +( + -- * Auxiliary types + module AERN2.MP.Precision + , module AERN2.MP.Accuracy + , module AERN2.MP.Enclosure + -- * The Ball type + , MPBall(..), CanBeMPBall, mpBall, CanBeMPBallP, mpBallP + , reducePrecionIfInaccurate + -- * Ball construction/extraction functions + , endpointsMP, fromEndpointsMP +) +where + +import MixedTypesNumPrelude +-- import qualified Prelude as P + +import Control.CollectErrors + +import GHC.Generics (Generic) + +import Text.Printf + +import AERN2.Normalize + +import AERN2.Norm + +import AERN2.MP.Dyadic +import qualified AERN2.MP.Float as MPFloat +import AERN2.MP.Float (MPFloat, mpFloat) +import AERN2.MP.Float.Operators +import AERN2.MP.Precision +import AERN2.MP.Accuracy +import qualified AERN2.MP.ErrorBound as EB +import AERN2.MP.ErrorBound (ErrorBound, errorBound) +import AERN2.MP.Enclosure + +data MPBall = MPBall + { ball_value :: MPFloat + , ball_error :: ErrorBound + } + -- { ball_value :: {-# UNPACK #-} ! MPFloat + -- , ball_error :: {-# UNPACK #-} ! ErrorBound + -- } + deriving (Generic) + +instance Show MPBall + where + show b@(MPBall x _e) = + -- printf "[%s ± %s](prec=%s)" (show x) (showAC $ getAccuracy b) (show $ integer $ getPrecision b) + printf "[%s ± %s]" (show x) (showAC $ getAccuracy b) + -- "[" ++ show x ++ " ± " ++ show e ++ "](prec=" ++ (show $ integer $ getPrecision x) ++ ")" + where + showAC Exact = "0" + showAC NoInformation = "oo" + showAC ac = "<2^(" ++ show (negate $ fromAccuracy ac) ++ ")" + + +instance (SuitableForCE es) => CanEnsureCE es MPBall where + +-- instance CanTestValid MPBall where +-- isValid = isFinite + +instance CanTestNaN MPBall where + isNaN = not . isFinite +instance CanTestFinite MPBall where + isInfinite = const False + isFinite (MPBall x e) = isFinite x && isFinite (mpFloat e) + +instance CanNormalize MPBall where + normalize b + | isFinite b = + b + -- reducePrecionIfInaccurate b + | otherwise = error $ "invalid MPBall: " ++ show b + +{-| + Reduce the precision of the ball centre if the + accuracy of the ball is poor. + + More precisely, reduce the precision of the centre + so that the ulp is approximately (radius / 1024), + unless the ulp is already lower than this. +-} +reducePrecionIfInaccurate :: MPBall -> MPBall +reducePrecionIfInaccurate b@(MPBall x _) = + case (bAcc, bNorm) of + (Exact, _) -> b + (_, NormZero) -> b + _ | p_e_nb < p_x -> setPrecision p_e_nb b + _ -> b + where + bAcc = getAccuracy b + bNorm = getNormLog b + p_x = getPrecision x + p_e_nb = prec $ max 2 (10 + nb + fromAccuracy bAcc) + (NormBits nb) = bNorm + +instance CanTestContains MPBall MPBall where + contains (MPBall xLarge eLarge) (MPBall xSmall eSmall) = + xLargeDy - eLargeDy <= xSmallDy - eSmallDy + && + xSmallDy + eSmallDy <= xLargeDy + eLargeDy + where + xLargeDy = dyadic xLarge + eLargeDy = dyadic eLarge + xSmallDy = dyadic xSmall + eSmallDy = dyadic eSmall + +$(declForTypes + [[t| Integer |], [t| Int |], [t| Rational |], [t| Dyadic |]] + (\ t -> [d| + instance CanTestContains MPBall $t where + contains (MPBall c e) x = + l <= x && x <= r + where + l = cDy - eDy + r = cDy + eDy + cDy = dyadic c + eDy = dyadic e + |])) + +{- ball construction/extraction functions -} + +instance IsInterval MPBall MPFloat where + fromEndpoints l u + | u < l = fromEndpoints u l + | otherwise = + MPBall (mpFloat cDy) (errorBound $ mpFloat eDy) + where + lDy = dyadic l + uDy = dyadic u + cDy = (lDy + uDy) * (dyadic 0.5) + eDy = (uDy - cDy) `max` (cDy - lDy) + endpoints (MPBall x e) = (mpFloat lDy, mpFloat uDy) + where + xDy = dyadic x + eDy = dyadic e + lDy = xDy - eDy + uDy = xDy + eDy + +fromEndpointsMP :: MPFloat -> MPFloat -> MPBall +fromEndpointsMP = fromEndpoints + +endpointsMP :: MPBall -> (MPFloat, MPFloat) +endpointsMP = endpoints + +instance IsInterval MPBall MPBall where + fromEndpoints l r = -- works as union even when r < l + fromEndpointsMP lMP uMP + where + lMP = min llMP rlMP + uMP = max luMP ruMP + (llMP, luMP) = endpointsMP l + (rlMP, ruMP) = endpointsMP r + endpoints x = (l,u) + where + l = MPBall lMP (errorBound 0) + u = MPBall uMP (errorBound 0) + (lMP, uMP) = endpointsMP x + +instance IsBall MPBall where + type CentreType MPBall = Dyadic + centre (MPBall cMP _e) = dyadic cMP + centreAsBallAndRadius x = (cB,e) + where + (MPBall cMP e) = x + cB = MPBall cMP (errorBound 0) + radius (MPBall _ e) = e + updateRadius updateFn (MPBall c e) = MPBall c (updateFn e) + +{--- constructing a ball with a given precision ---} + +type CanBeMPBallP t = (ConvertibleWithPrecision t MPBall) + +mpBallP :: (CanBeMPBallP t) => Precision -> t -> MPBall +mpBallP = convertP + + +{--- constructing an exact ball ---} + +type CanBeMPBall t = ConvertibleExactly t MPBall + +mpBall :: (CanBeMPBall t) => t -> MPBall +mpBall = convertExactly + +{-- extracting approximate information about a ball --} + +instance HasAccuracy MPBall where + getAccuracy = getAccuracy . ball_error + +instance CanReduceSizeUsingAccuracyGuide MPBall where + reduceSizeUsingAccuracyGuide acGuide b@(MPBall x _e) = + case acGuide of + NoInformation -> lowerPrecisionIfAbove (prec 2) b + _ | getAccuracy b > acGuide -> tryPrec newPrec + _ -> b + where + tryPrec p + | getAccuracy bP >= acGuide = bP + | otherwise = tryPrec (p + 10) + where + bP = lowerPrecisionIfAbove p b + queryBits = 1 + fromAccuracy acGuide + newPrec = + case (getNormLog x) of + NormBits xNormBits -> + prec (max 2 (queryBits + xNormBits + 2)) + NormZero -> + prec $ max 2 queryBits + -- bWithLowAC = + -- case acGuide of + -- Exact -> b + -- NoInformation -> b + -- _ -> normalize $ + -- MPBall x (errorBound ((0.5^(fromAccuracy acGuide))⚡)) + +instance HasNorm MPBall where + getNormLog ball = getNormLog boundMP + where + (_, MPBall boundMP _) = endpoints $ absRaw ball + +instance HasApproximate MPBall where + type Approximate MPBall = (MPFloat, Bool) + getApproximate ac b@(MPBall x e) = + (approx, isAccurate) + where + isAccurate = getAccuracy b < ac + approx + | closeToN = n + | otherwise = MPFloat.setPrecisionUp (prec (fromAccuracy ac)) x + where + n = mpFloat $ round $ rational x + closeToN = ((abs $ x -^ n) <= e) + +instance HasPrecision MPBall where + getPrecision = getPrecision . ball_value + +instance CanSetPrecision MPBall where + setPrecision p (MPBall x e) + | p >= pPrev = MPBall xUp e + | otherwise = MPBall xUp (e + (xUp `EB.subMP` xDown)) + where + pPrev = MPFloat.getPrecision x + xUp = MPFloat.setPrecisionUp p x + xDown = MPFloat.setPrecisionDown p x + +{- negation & abs -} + +instance CanNeg MPBall where + negate (MPBall x e) = MPBall (-x) e + +instance CanAbs MPBall where + abs = normalize . absRaw + +absRaw :: MPBall -> MPBall +absRaw b + | l < 0 && 0 < r = + fromEndpointsMP (mpFloat 0) (max (-l) r) + | 0 <= l = b + | otherwise = -b + where + (l,r) = endpointsMP b diff --git a/src/AERN2/MP/Dyadic.hs b/src/AERN2/MP/Dyadic.hs --- a/src/AERN2/MP/Dyadic.hs +++ b/src/AERN2/MP/Dyadic.hs @@ -1,4 +1,6 @@ {-# LANGUAGE CPP #-} +-- #define DEBUG +{-# LANGUAGE DeriveDataTypeable #-} {-| Module : AERN2.MP.Dyadic Description : Dyadics with exact ring operations @@ -10,19 +12,662 @@ Portability : portable Arbitrary precision floating-point numbers with exact ring operations. + + Currently, we use hmpfr when compiling with ghc 7.10 and higher + and haskell-mpfr when compiling with ghc 7.8. -} module AERN2.MP.Dyadic - ( -#ifdef USEMPFR - module AERN2.MP.UseMPFR.Dyadic -#endif - ) +( + -- * Dyadic numbers and their basic operations + Dyadic, HasDyadics + -- * Dyadic constructors + , CanBeDyadic, dyadic + -- * tests + , specDyadic, tDyadic +) where --- import MixedTypesNumPrelude --- import qualified Prelude as P - -#ifdef USEMPFR -import AERN2.MP.UseMPFR.Dyadic +#ifdef DEBUG +import Debug.Trace (trace) +#define maybeTrace trace +#define maybeTraceIO putStrLn +#else +#define maybeTrace (\ (_ :: String) t -> t) +#define maybeTraceIO (\ (_ :: String) -> return ()) #endif + +import MixedTypesNumPrelude +import qualified Prelude as P + +import Control.CollectErrors + +import Text.Printf +import Text.Regex.TDFA + +import Data.Typeable +import Data.Convertible + +import Test.Hspec +import Test.QuickCheck +-- import qualified Test.Hspec.SmallCheck as SC + +import Data.Ratio (denominator, numerator) + +import Math.NumberTheory.Logarithms (integerLog2) + +import AERN2.Norm +import AERN2.MP.Precision +import AERN2.MP.Accuracy +import AERN2.MP.Float + +{-| Exact dyadic type based on MPFloat. -} +newtype Dyadic = Dyadic { dyadicMPFloat :: MPFloat } + deriving (P.Eq, P.Ord, CanRound, HasPrecision, HasNorm, Typeable) + +instance Ring Dyadic +instance Ring (CN Dyadic) + +instance OrderedRing Dyadic +instance OrderedRing (CN Dyadic) + +instance OrderedCertainlyRing Dyadic +instance OrderedCertainlyRing (CN Dyadic) + +instance HasAccuracy Dyadic where getAccuracy _ = Exact + +instance Show Dyadic where + show (Dyadic x) + | e == 0 = printf "dyadic (%d)" n + | e > 0 = printf "dyadic (%d*0.5^%d)" n e + | otherwise = error "in show Dyadic" + where + xR = rational x + NormBits e = getNormLog (denominator xR) + n = numerator xR + +instance Read Dyadic where + readsPrec _pr dyadicS = + tryInt $ tryWithExp [] + where + tryInt tryNext = + case groups of + [nS] -> + case reads nS of + [(n,"")] -> [(dyadic (n :: Integer), afterS)] + _ -> tryNext + _ -> tryNext + where + (_,_,afterS,groups) = + dyadicS =~ "\\`dyadic \\(([-0-9]*)\\)" + :: (String, String, String, [String]) + tryWithExp tryNext = + case groups of + [nS,eS] -> + case (reads nS, reads eS) of + ([(n,"")],[(e,"")]) -> + [((n :: Integer)*(dyadic 0.5)^!(e :: Integer), afterS)] + _ -> tryNext + _ -> tryNext + where + (_,_,afterS,groups) = + dyadicS =~ "\\`dyadic \\(([-0-9]*)\\*0.5\\^([0-9]*)\\)" + :: (String, String, String, [String]) + +instance (SuitableForCE es) => CanEnsureCE es Dyadic + + +{-- conversions --} + +type HasDyadics t = ConvertibleExactly Dyadic t + +instance ConvertibleExactly Dyadic Dyadic where + safeConvertExactly = Right + +instance ConvertibleExactly Dyadic MPFloat where + safeConvertExactly = Right . dyadicMPFloat + +instance ConvertibleExactly Dyadic Rational where + safeConvertExactly = safeConvertExactly . dyadicMPFloat + +type CanBeDyadic t = ConvertibleExactly t Dyadic +dyadic :: (CanBeDyadic t) => t -> Dyadic +dyadic = convertExactly + +instance ConvertibleExactly MPFloat Dyadic where + safeConvertExactly = Right . Dyadic + +instance HasIntegerBounds Dyadic where + integerBounds d = (floor d, ceiling d) + +instance ConvertibleExactly Integer Dyadic where + safeConvertExactly = fmap Dyadic . safeConvertExactly + +instance ConvertibleExactly Int Dyadic where + safeConvertExactly = fmap Dyadic . safeConvertExactly + +instance ConvertibleExactly Rational Dyadic where + safeConvertExactly q + | isDyadic = Right $ Dyadic (fromRationalUp (prec $ max 2 (dp + np + 1)) q) + | otherwise = convError "this number is not dyadic" q + where + isDyadic = d == 2^!dp + dp = integerLog2 d + d = denominator q + np = integerLog2 (max 1 $ abs $ numerator q) + +instance Convertible Dyadic Double where + safeConvert = safeConvert . dyadicMPFloat + +instance (ConvertibleExactly Dyadic t, Monoid es) => ConvertibleExactly Dyadic (CollectErrors es t) where + safeConvertExactly = fmap (\v -> CollectErrors (Just v) mempty) . safeConvertExactly + +{-- comparisons --} + +instance HasEqAsymmetric Dyadic Dyadic +instance HasEqAsymmetric Dyadic Integer where + equalTo = convertSecond equalTo +instance HasEqAsymmetric Integer Dyadic where + equalTo = convertFirst equalTo +instance HasEqAsymmetric Dyadic Int where + equalTo = convertSecond equalTo +instance HasEqAsymmetric Int Dyadic where + equalTo = convertFirst equalTo +instance HasEqAsymmetric Dyadic Rational where + equalTo = convertFirst equalTo +instance HasEqAsymmetric Rational Dyadic where + equalTo = convertSecond equalTo + +instance + (HasEqAsymmetric Dyadic b + , CanEnsureCE es b + , CanEnsureCE es (EqCompareType Dyadic b) + , IsBool (EnsureCE es (EqCompareType Dyadic b)) + , SuitableForCE es) + => + HasEqAsymmetric Dyadic (CollectErrors es b) + where + type EqCompareType Dyadic (CollectErrors es b) = + EnsureCE es (EqCompareType Dyadic b) + equalTo = lift2TLCE equalTo + +instance + (HasEqAsymmetric a Dyadic + , CanEnsureCE es a + , CanEnsureCE es (EqCompareType a Dyadic) + , IsBool (EnsureCE es (EqCompareType a Dyadic)) + , SuitableForCE es) + => + HasEqAsymmetric (CollectErrors es a) Dyadic + where + type EqCompareType (CollectErrors es a) Dyadic = + EnsureCE es (EqCompareType a Dyadic) + equalTo = lift2TCE equalTo + +instance CanTestZero Dyadic + +instance HasOrderAsymmetric Dyadic Dyadic +instance HasOrderAsymmetric Dyadic Integer where + lessThan = convertSecond lessThan + leq = convertSecond leq +instance HasOrderAsymmetric Integer Dyadic where + lessThan = convertFirst lessThan + leq = convertFirst leq +instance HasOrderAsymmetric Dyadic Int where + lessThan = convertSecond lessThan + leq = convertSecond leq +instance HasOrderAsymmetric Int Dyadic where + lessThan = convertFirst lessThan + leq = convertFirst leq +instance HasOrderAsymmetric Rational Dyadic where + lessThan = convertSecond lessThan + leq = convertSecond leq +instance HasOrderAsymmetric Dyadic Rational where + lessThan = convertFirst lessThan + leq = convertFirst leq + +instance + (HasOrderAsymmetric Dyadic b + , CanEnsureCE es b + , CanEnsureCE es (OrderCompareType Dyadic b) + , IsBool (EnsureCE es (OrderCompareType Dyadic b)) + , SuitableForCE es) + => + HasOrderAsymmetric Dyadic (CollectErrors es b) + where + type OrderCompareType Dyadic (CollectErrors es b) = + EnsureCE es (OrderCompareType Dyadic b) + lessThan = lift2TLCE lessThan + leq = lift2TLCE leq + greaterThan = lift2TLCE greaterThan + geq = lift2TLCE geq + +instance + (HasOrderAsymmetric a Dyadic + , CanEnsureCE es a + , CanEnsureCE es (OrderCompareType a Dyadic) + , IsBool (EnsureCE es (OrderCompareType a Dyadic)) + , SuitableForCE es) + => + HasOrderAsymmetric (CollectErrors es a) Dyadic + where + type OrderCompareType (CollectErrors es a) Dyadic = + EnsureCE es (OrderCompareType a Dyadic) + lessThan = lift2TCE lessThan + leq = lift2TCE leq + greaterThan = lift2TCE greaterThan + geq = lift2TCE geq + + +instance CanTestPosNeg Dyadic + +instance CanTestInteger Dyadic where + certainlyNotInteger = certainlyNotInteger . rational + certainlyIntegerGetIt = certainlyIntegerGetIt . rational + +{- unary functions -} + +instance CanNeg Dyadic where + negate = lift1 negate + +instance CanAbs Dyadic where + abs = lift1 abs + +lift1 :: (MPFloat -> MPFloat) -> (Dyadic -> Dyadic) +lift1 op (Dyadic x) = Dyadic (op x) + +{- min/max -} + +instance CanMinMaxAsymmetric Dyadic Dyadic +instance CanMinMaxAsymmetric Integer Dyadic where + type MinMaxType Integer Dyadic = Dyadic + min = convertFirst min + max = convertFirst max +instance CanMinMaxAsymmetric Dyadic Integer where + type MinMaxType Dyadic Integer = Dyadic + min = convertSecond min + max = convertSecond max +instance CanMinMaxAsymmetric Int Dyadic where + type MinMaxType Int Dyadic = Dyadic + min = convertFirst min + max = convertFirst max +instance CanMinMaxAsymmetric Dyadic Int where + type MinMaxType Dyadic Int = Dyadic + min = convertSecond min + max = convertSecond max +instance CanMinMaxAsymmetric Rational Dyadic where + type MinMaxType Rational Dyadic = Rational + min = convertSecond min + max = convertSecond max +instance CanMinMaxAsymmetric Dyadic Rational where + type MinMaxType Dyadic Rational = Rational + min = convertFirst min + max = convertFirst max + +instance + (CanMinMaxAsymmetric Dyadic b + , CanEnsureCE es b + , CanEnsureCE es (MinMaxType Dyadic b) + , SuitableForCE es) + => + CanMinMaxAsymmetric Dyadic (CollectErrors es b) + where + type MinMaxType Dyadic (CollectErrors es b) = + EnsureCE es (MinMaxType Dyadic b) + min = lift2TLCE min + max = lift2TLCE max + +instance + (CanMinMaxAsymmetric a Dyadic + , CanEnsureCE es a + , CanEnsureCE es (MinMaxType a Dyadic) + , SuitableForCE es) + => + CanMinMaxAsymmetric (CollectErrors es a) Dyadic + where + type MinMaxType (CollectErrors es a) Dyadic = + EnsureCE es (MinMaxType a Dyadic) + min = lift2TCE min + max = lift2TCE max + +{- addition -} + +instance CanAddAsymmetric Dyadic Dyadic where + add = lift2 addDown addUp + +instance CanAddAsymmetric Integer Dyadic where + type AddType Integer Dyadic = Dyadic + add = convertFirst add +instance CanAddAsymmetric Dyadic Integer where + type AddType Dyadic Integer = Dyadic + add = convertSecond add + +instance CanAddAsymmetric Int Dyadic where + type AddType Int Dyadic = Dyadic + add = convertFirst add +instance CanAddAsymmetric Dyadic Int where + type AddType Dyadic Int = Dyadic + add = convertSecond add + +instance CanAddAsymmetric Rational Dyadic where + type AddType Rational Dyadic = Rational + add = convertSecond add +instance CanAddAsymmetric Dyadic Rational where + type AddType Dyadic Rational = Rational + add = convertFirst add + +instance + (CanAddAsymmetric Dyadic b + , CanEnsureCE es b + , CanEnsureCE es (AddType Dyadic b) + , SuitableForCE es) + => + CanAddAsymmetric Dyadic (CollectErrors es b) + where + type AddType Dyadic (CollectErrors es b) = + EnsureCE es (AddType Dyadic b) + add = lift2TLCE add + +instance + (CanAddAsymmetric a Dyadic + , CanEnsureCE es a + , CanEnsureCE es (AddType a Dyadic) + , SuitableForCE es) + => + CanAddAsymmetric (CollectErrors es a) Dyadic + where + type AddType (CollectErrors es a) Dyadic = + EnsureCE es (AddType a Dyadic) + add = lift2TCE add + +{- subtraction -} + +instance CanSub Dyadic Dyadic where + sub = lift2 subDown subUp + +instance CanSub Integer Dyadic where + type SubType Integer Dyadic = Dyadic + sub = convertFirst sub +instance CanSub Dyadic Integer where + type SubType Dyadic Integer = Dyadic + sub = convertSecond sub + +instance CanSub Int Dyadic where + type SubType Int Dyadic = Dyadic + sub = convertFirst sub +instance CanSub Dyadic Int where + type SubType Dyadic Int = Dyadic + sub = convertSecond sub + +instance CanSub Rational Dyadic where + type SubType Rational Dyadic = Rational + sub = convertSecond sub +instance CanSub Dyadic Rational where + type SubType Dyadic Rational = Rational + sub = convertFirst sub + +instance + (CanSub Dyadic b + , CanEnsureCE es b + , CanEnsureCE es (SubType Dyadic b) + , SuitableForCE es) + => + CanSub Dyadic (CollectErrors es b) + where + type SubType Dyadic (CollectErrors es b) = + EnsureCE es (SubType Dyadic b) + sub = lift2TLCE sub + +instance + (CanSub a Dyadic + , CanEnsureCE es a + , CanEnsureCE es (SubType a Dyadic) + , SuitableForCE es) + => + CanSub (CollectErrors es a) Dyadic + where + type SubType (CollectErrors es a) Dyadic = + EnsureCE es (SubType a Dyadic) + sub = lift2TCE sub + + +{- multiplication -} + +instance CanMulAsymmetric Dyadic Dyadic where + mul = lift2 mulDown mulUp + +instance CanMulAsymmetric Integer Dyadic where + type MulType Integer Dyadic = Dyadic + mul = convertFirst mul +instance CanMulAsymmetric Dyadic Integer where + type MulType Dyadic Integer = Dyadic + mul = convertSecond mul + +instance CanMulAsymmetric Int Dyadic where + type MulType Int Dyadic = Dyadic + mul = convertFirst mul +instance CanMulAsymmetric Dyadic Int where + type MulType Dyadic Int = Dyadic + mul = convertSecond mul + +instance CanMulAsymmetric Rational Dyadic where + type MulType Rational Dyadic = Rational + mul = convertSecond mul +instance CanMulAsymmetric Dyadic Rational where + type MulType Dyadic Rational = Rational + mul = convertFirst mul + +instance + (CanMulAsymmetric Dyadic b + , CanEnsureCE es b + , CanEnsureCE es (MulType Dyadic b) + , SuitableForCE es) + => + CanMulAsymmetric Dyadic (CollectErrors es b) + where + type MulType Dyadic (CollectErrors es b) = + EnsureCE es (MulType Dyadic b) + mul = lift2TLCE mul + +instance + (CanMulAsymmetric a Dyadic + , CanEnsureCE es a + , CanEnsureCE es (MulType a Dyadic) + , SuitableForCE es) + => + CanMulAsymmetric (CollectErrors es a) Dyadic + where + type MulType (CollectErrors es a) Dyadic = + EnsureCE es (MulType a Dyadic) + mul = lift2TCE mul + +instance CanPow Dyadic Integer where + powNoCN = powUsingMul (dyadic 1) + pow = integerPowCN (powUsingMul (dyadic 1)) +instance CanPow Dyadic Int where + powNoCN = powUsingMul (dyadic 1) + pow = integerPowCN (powUsingMul (dyadic 1)) + +instance + (CanDiv a Dyadic + , CanEnsureCE es a + , CanEnsureCE es (DivType a Dyadic) + , CanEnsureCE es (DivTypeNoCN a Dyadic) + , SuitableForCE es) + => + CanDiv (CollectErrors es a) Dyadic + where + type DivType (CollectErrors es a) Dyadic = + EnsureCE es (DivType a Dyadic) + divide = lift2TCE divide + type DivTypeNoCN (CollectErrors es a) Dyadic = + EnsureCE es (DivTypeNoCN a Dyadic) + divideNoCN = lift2TCE divideNoCN + +instance CanDiv Integer Dyadic where + type DivTypeNoCN Integer Dyadic = Rational + divideNoCN a b = divideNoCN a (rational b) +instance CanDiv Dyadic Integer where + type DivTypeNoCN Dyadic Integer = Rational + divideNoCN a b = divideNoCN (rational a) b + +instance CanDiv Int Dyadic where + type DivTypeNoCN Int Dyadic = Rational + divideNoCN a b = divideNoCN a (rational b) +instance CanDiv Dyadic Int where + type DivTypeNoCN Dyadic Int = Rational + divideNoCN a b = divideNoCN (rational a) b + +instance CanDiv Rational Dyadic where + type DivTypeNoCN Rational Dyadic = Rational + divideNoCN = convertSecond divideNoCN +instance CanDiv Dyadic Rational where + type DivTypeNoCN Dyadic Rational = Rational + divideNoCN = convertFirst divideNoCN + +instance + (CanDiv Dyadic b + , CanEnsureCE es b + , CanEnsureCE es (DivType Dyadic b) + , CanEnsureCE es (DivTypeNoCN Dyadic b) + , SuitableForCE es) + => + CanDiv Dyadic (CollectErrors es b) + where + type DivType Dyadic (CollectErrors es b) = + EnsureCE es (DivType Dyadic b) + divide = lift2TLCE divide + type DivTypeNoCN Dyadic (CollectErrors es b) = + EnsureCE es (DivTypeNoCN Dyadic b) + divideNoCN = lift2TLCE divideNoCN + +instance + (CanPow Dyadic b + , CanEnsureCE es b + , CanEnsureCE es (PowTypeNoCN Dyadic b) + , CanEnsureCE es (PowType Dyadic b) + , SuitableForCE es) + => + CanPow Dyadic (CollectErrors es b) + where + type PowTypeNoCN Dyadic (CollectErrors es b) = + EnsureCE es (PowTypeNoCN Dyadic b) + powNoCN = lift2TLCE powNoCN + type PowType Dyadic (CollectErrors es b) = + EnsureCE es (PowType Dyadic b) + pow = lift2TLCE pow + +instance + (CanPow a Dyadic + , CanEnsureCE es a + , CanEnsureCE es (PowType a Dyadic) + , CanEnsureCE es (PowTypeNoCN a Dyadic) + , SuitableForCE es) + => + CanPow (CollectErrors es a) Dyadic + where + type PowTypeNoCN (CollectErrors es a) Dyadic = + EnsureCE es (PowTypeNoCN a Dyadic) + powNoCN = lift2TCE powNoCN + type PowType (CollectErrors es a) Dyadic = + EnsureCE es (PowType a Dyadic) + pow = lift2TCE pow + +lift2 :: + (MPFloat -> MPFloat -> MPFloat) -> + (MPFloat -> MPFloat -> MPFloat) -> + (Dyadic -> Dyadic -> Dyadic) +lift2 opDown opUp (Dyadic x0) (Dyadic y0) = Dyadic (opExact x0 y0) + where + opExact x y + | rUp == rDown = rUp + | otherwise = + maybeTrace (printf "Dyadic.lift2: rUp = %s; rDown = %s; p = %s" (show rUp) (show rDown) (show $ integer p)) $ + opExact xH yH + where + rUp = opUp x y + rDown = opDown x y + xH = setPrecision pH x + yH = setPrecision pH y + pH = precisionTimes2 p + p = getPrecision rUp + +instance Arbitrary Dyadic where + arbitrary = + do + c <- finiteMPFloat + return (Dyadic c) + where + finiteMPFloat = + do + x <- arbitrary + if isFinite x + then return x + else finiteMPFloat + +{-| + A runtime representative of type @Dyadic@. + Used for specialising polymorphic tests to concrete types. +-} +tDyadic :: T Dyadic +tDyadic = T "Dyadic" + +specDyadic :: Spec +specDyadic = + describe ("Dyadic") $ do + specConversion tInteger tDyadic dyadic round + specConversion tDyadic tRational rational dyadic + describe "order" $ do + specHasEqNotMixed tDyadic + specHasEq tInt tDyadic tRational + specCanTestZero tDyadic + specHasOrderNotMixed tDyadic + specHasOrder tInt tDyadic tRational + describe "min/max/abs" $ do + specCanNegNum tDyadic + specCanAbs tDyadic + specCanMinMaxNotMixed tDyadic + specCanMinMax tDyadic tInteger tDyadic + it "min Dyadic Rational (dyadic only)" $ do + property $ \ (x :: Dyadic) (y :: Dyadic) -> + x `min` y == x `min` (rational y) + it "max Dyadic Rational (dyadic only)" $ do + property $ \ (x :: Dyadic) (y :: Dyadic) -> + x `max` y == x `max` (rational y) + describe "ring" $ do + specCanAddNotMixed tDyadic + specCanAddSameType tDyadic + specCanAdd tInt tDyadic tInteger + specCanAdd tInteger tDyadic tInt + it "Dyadic + Rational (dyadic only)" $ do + property $ \ (x :: Dyadic) (y :: Dyadic) -> + x + y == x + (rational y) + specCanSubNotMixed tDyadic + specCanSub tDyadic tInteger + specCanSub tInteger tDyadic + specCanSub tDyadic tInt + specCanSub tInt tDyadic + it "Dyadic - Rational (dyadic only)" $ do + property $ \ (x :: Dyadic) (y :: Dyadic) -> + x - y == x - (rational y) + specCanMulNotMixed tDyadic + specCanMulSameType tDyadic + specCanMul tInt tDyadic tInteger + it "Dyadic * Rational (dyadic only)" $ do + property $ \ (x :: Dyadic) (y :: Dyadic) -> + x * y == x * (rational y) + specCanPow tDyadic tInteger + +instance P.Num Dyadic where + fromInteger = convertExactly + negate = negate + (+) = (+) + (*) = (*) + abs = abs + signum d + | d < 0 = dyadic (-1) + | d == 0 = dyadic 0 + | otherwise = dyadic 1 + +instance P.Real Dyadic where + toRational = convertExactly diff --git a/src/AERN2/MP/Enclosure.hs b/src/AERN2/MP/Enclosure.hs --- a/src/AERN2/MP/Enclosure.hs +++ b/src/AERN2/MP/Enclosure.hs @@ -33,7 +33,7 @@ import Control.CollectErrors -import AERN2.MP.UseMPFR.ErrorBound +import AERN2.MP.ErrorBound import AERN2.MP.Accuracy {- ball-specific operations -} diff --git a/src/AERN2/MP/ErrorBound.hs b/src/AERN2/MP/ErrorBound.hs --- a/src/AERN2/MP/ErrorBound.hs +++ b/src/AERN2/MP/ErrorBound.hs @@ -1,4 +1,4 @@ -{-# LANGUAGE CPP #-} +{-# LANGUAGE DeriveDataTypeable #-} {-| Module : AERN2.MP.ErrorBound Description : Fixed precision non-negative up-rounded floating-point numbers @@ -10,19 +10,212 @@ Portability : portable Fixed precision non-negative up-rounded floating-point numbers. --} + Currently using a fixed-precision MPFloat. +-} module AERN2.MP.ErrorBound -#ifdef USEMPFR - ( - module AERN2.MP.UseMPFR.ErrorBound - ) -#endif + (ErrorBound, CanBeErrorBound, errorBound, + absMP, subMP) where --- import MixedTypesNumPrelude --- import qualified Prelude as P +import MixedTypesNumPrelude +import qualified Prelude as P -#ifdef USEMPFR -import AERN2.MP.UseMPFR.ErrorBound -#endif +import Data.Typeable + +import Test.QuickCheck + +import Data.Convertible + +import Math.NumberTheory.Logarithms (integerLog2) + +import AERN2.MP.Precision +import AERN2.MP.Accuracy +import qualified AERN2.MP.Float as MPFloat +import AERN2.MP.Float (MPFloat, mpFloat, frequencyElements) +import AERN2.MP.Float.Operators +import AERN2.MP.Dyadic + +{- example -} + +_example1 :: ErrorBound +_example1 = 2*((errorBound 0.01) + 0.1*(errorBound 0.01)/3) + +{- type -} + +{-| A non-negative Double value to serve as an error bound. Arithmetic is rounded towards +infinity. -} +newtype ErrorBound = ErrorBound { er2mp :: MPFloat } + deriving (P.Eq, P.Ord, Typeable) + +instance Show ErrorBound where + show (ErrorBound d) = show d + +errorBoundPrecision :: Precision +errorBoundPrecision = prec 53 + +instance HasAccuracy ErrorBound where + getAccuracy (ErrorBound e) + | eN > 0 = + bits $ negate $ integerLog2 eN + | e > 0 && eRecipN > 0 = + bits $ integerLog2 eRecipN + | e == 0 = Exact + | otherwise = NoInformation + where + eN = floor $ rational e + eRecipN = ceiling $ rational $ MPFloat.recipDown e + +{- conversions -} + +instance ConvertibleExactly ErrorBound MPFloat where + safeConvertExactly = Right . er2mp + +instance ConvertibleExactly ErrorBound Dyadic where + safeConvertExactly = Right . dyadic . er2mp + +instance ConvertibleExactly ErrorBound Rational where + safeConvertExactly = Right . convertExactly . mpFloat + +type CanBeErrorBound t = Convertible t ErrorBound +errorBound :: (CanBeErrorBound t) => t -> ErrorBound +errorBound = convert + +instance Convertible Rational ErrorBound where + safeConvert x + | x >= 0 = Right $ ErrorBound $ MPFloat.fromRationalUp errorBoundPrecision x + | otherwise = convError "Trying to construct a negative ErrorBound" x + +instance Convertible MPFloat ErrorBound where + safeConvert x + | x >= 0 = Right $ ErrorBound $ MPFloat.setPrecisionUp errorBoundPrecision x + | otherwise = convError "Trying to construct a negative ErrorBound" x + +instance Convertible Integer ErrorBound where + safeConvert x + | x >= 0 = Right $ ErrorBound $ MPFloat.fromIntegerUp errorBoundPrecision x + | otherwise = convError "Trying to construct a negative ErrorBound" x + +instance Convertible Int ErrorBound where + safeConvert = safeConvert . integer + +{- comparisons -} + +instance HasOrderAsymmetric ErrorBound ErrorBound + +instance HasOrderAsymmetric ErrorBound MPFloat where + lessThan = convertFirst lessThan + leq = convertFirst leq +instance HasOrderAsymmetric MPFloat ErrorBound where + lessThan = convertSecond lessThan + leq = convertSecond leq + +instance HasEqAsymmetric ErrorBound Rational where + equalTo = convertFirst equalTo +instance HasEqAsymmetric Rational ErrorBound where + equalTo = convertSecond equalTo +instance HasOrderAsymmetric ErrorBound Rational where + lessThan = convertFirst lessThan + leq = convertFirst leq +instance HasOrderAsymmetric Rational ErrorBound where + lessThan = convertSecond lessThan + leq = convertSecond leq + +instance HasEqAsymmetric ErrorBound Integer where + equalTo a b = equalTo (dyadic a) (dyadic b) +instance HasEqAsymmetric Integer ErrorBound where + equalTo a b = equalTo (dyadic a) (dyadic b) +instance HasOrderAsymmetric ErrorBound Integer where + lessThan a b = lessThan (dyadic a) (dyadic b) + leq a b = leq (dyadic a) (dyadic b) +instance HasOrderAsymmetric Integer ErrorBound where + lessThan a b = lessThan (dyadic a) (dyadic b) + leq a b = leq (dyadic a) (dyadic b) + +instance HasEqAsymmetric ErrorBound Int where + equalTo a b = equalTo (dyadic a) (dyadic b) +instance HasEqAsymmetric Int ErrorBound where + equalTo a b = equalTo (dyadic a) (dyadic b) +instance HasOrderAsymmetric ErrorBound Int where + lessThan a b = lessThan (dyadic a) (dyadic b) + leq a b = leq (dyadic a) (dyadic b) +instance HasOrderAsymmetric Int ErrorBound where + lessThan a b = lessThan (dyadic a) (dyadic b) + leq a b = leq (dyadic a) (dyadic b) + +instance CanMinMaxAsymmetric ErrorBound ErrorBound + +{- converting operations -} + +subMP :: MPFloat -> MPFloat -> ErrorBound +a `subMP` b = errorBound $ a -^ b + +absMP :: MPFloat -> ErrorBound +absMP = errorBound . abs + +{- up-rounded operations -} + +instance CanAddAsymmetric ErrorBound ErrorBound where + add (ErrorBound a) (ErrorBound b) = ErrorBound $ a +^ b + +instance CanAddAsymmetric ErrorBound MPFloat where + type AddType ErrorBound MPFloat = ErrorBound + add = convertSecondUsing (\ _ f -> convert f) add +instance CanAddAsymmetric MPFloat ErrorBound where + type AddType MPFloat ErrorBound = ErrorBound + add = convertFirstUsing (\ f _ -> convert f) add + +instance CanMulAsymmetric ErrorBound ErrorBound where + mul (ErrorBound a) (ErrorBound b) = ErrorBound $ a *^ b + +instance CanMulAsymmetric ErrorBound MPFloat where + type MulType ErrorBound MPFloat = ErrorBound + mul = convertSecondUsing (\ _ f -> convert f) mul +instance CanMulAsymmetric MPFloat ErrorBound where + type MulType MPFloat ErrorBound = ErrorBound + mul = convertFirstUsing (\ f _ -> convert f) mul + +instance CanMulAsymmetric ErrorBound Integer where + type MulType ErrorBound Integer = ErrorBound + mul (ErrorBound a) i + | i >= 0 = ErrorBound $ a *^ (MPFloat.fromIntegerUp errorBoundPrecision i) + | otherwise = error "trying to multiply ErrorBound by a negative integer" +instance CanMulAsymmetric Integer ErrorBound where + type MulType Integer ErrorBound = ErrorBound + mul i (ErrorBound b) + | i >= 0 = ErrorBound $ (MPFloat.fromIntegerUp errorBoundPrecision i) *^ b + | otherwise = error "trying to multiply ErrorBound by a negative integer" + +instance CanMulAsymmetric ErrorBound Rational where + type MulType ErrorBound Rational = ErrorBound + mul (ErrorBound a) r + | r >= 0.0 = ErrorBound $ a *^ (MPFloat.fromRationalUp errorBoundPrecision r) + | otherwise = error "trying to multiply ErrorBound by a negative integer" +instance CanMulAsymmetric Rational ErrorBound where + type MulType Rational ErrorBound = ErrorBound + mul r (ErrorBound b) + | r >= 0.0 = ErrorBound $ (MPFloat.fromRationalUp errorBoundPrecision r) *^ b + | otherwise = error "trying to multiply ErrorBound by a negative integer" + +instance CanDiv ErrorBound Integer where + type DivTypeNoCN ErrorBound Integer = ErrorBound + type DivType ErrorBound Integer = ErrorBound + divideNoCN = divide + divide (ErrorBound a) i + | i > 0 = ErrorBound $ a /^ (MPFloat.fromIntegerUp errorBoundPrecision i) + | otherwise = error "trying to multiply ErrorBound by a non-positive integer" + +instance Arbitrary ErrorBound where + arbitrary = + do + giveSpecialValue <- frequencyElements [(5, False),(1, True)] + aux giveSpecialValue + where + aux giveSpecialValue + | giveSpecialValue = + elements (map convert [0.0,0.0,0.0,10.0,1.0,0.5,0.125]) + | otherwise = + do + (s :: Integer) <- arbitrary + let resultR = ((abs s) `mod` (2^!35))/!(2^!32) + let result = convert resultR + return result diff --git a/src/AERN2/MP/Float.hs b/src/AERN2/MP/Float.hs new file mode 100644 --- /dev/null +++ b/src/AERN2/MP/Float.hs @@ -0,0 +1,78 @@ +{-# LANGUAGE CPP #-} +{-| + Module : AERN2.MP.Float + Description : Arbitrary precision floating point numbers + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Arbitrary precision floating-point numbers with up/down-rounded operations. + + Currently, we use hmpfr when compiling with ghc 7.10 and higher + and haskell-mpfr when compiling with ghc 7.8. +-} + +module AERN2.MP.Float + ( + -- * Precision operations + module AERN2.MP.Precision + -- * The type definition and basic operations +#ifdef UseCDAR + , module AERN2.MP.Float.UseCDAR.Type + -- * Arithmetic operations + , module AERN2.MP.Float.UseCDAR.Arithmetic + , distUp, distDown, avgUp, avgDown + -- * Conversions, comparisons and norm + , module AERN2.MP.Float.UseCDAR.Conversions +#else + , module AERN2.MP.Float.UseRounded.Type + -- * Arithmetic operations + , module AERN2.MP.Float.UseRounded.Arithmetic + , distUp, distDown, avgUp, avgDown + -- * Conversions, comparisons and norm + , module AERN2.MP.Float.UseRounded.Conversions +#endif + -- * Infix operators for up/down-rounded operations + , module AERN2.MP.Float.Operators + -- * Constants such as NaN, infinity + , module AERN2.MP.Float.Constants + -- * Tests + , module AERN2.MP.Float.Tests + ) +where + +import MixedTypesNumPrelude +-- import qualified Prelude as P + +import AERN2.MP.Precision + +#ifdef UseCDAR +import AERN2.MP.Float.UseCDAR.Type +import AERN2.MP.Float.UseCDAR.Arithmetic +import AERN2.MP.Float.UseCDAR.Conversions +#else +import AERN2.MP.Float.UseRounded.Type +import AERN2.MP.Float.UseRounded.Arithmetic +import AERN2.MP.Float.UseRounded.Conversions +#endif + +import AERN2.MP.Float.Operators +import AERN2.MP.Float.Constants +import AERN2.MP.Float.Tests + +-- | Computes an upper bound to the distance @|x - y|@ of @x@ and @y@. +distUp :: MPFloat -> MPFloat -> MPFloat +distUp x y = if x >= y then x -^ y else y -^ x + +-- | Computes a lower bound to the distance @|x - y|@ of @x@ and @y@. +distDown :: MPFloat -> MPFloat -> MPFloat +distDown x y = if x >= y then x -. y else y -. x + +avgUp :: MPFloat -> MPFloat -> MPFloat +avgUp x y = (x +^ y) /^ (mpFloat 2) + +avgDown :: MPFloat -> MPFloat -> MPFloat +avgDown x y = (x +. y) /. (mpFloat 2) diff --git a/src/AERN2/MP/Float/Constants.hs b/src/AERN2/MP/Float/Constants.hs new file mode 100644 --- /dev/null +++ b/src/AERN2/MP/Float/Constants.hs @@ -0,0 +1,58 @@ +{-# LANGUAGE CPP #-} +{-| + Module : AERN2.MP.Float.Constants + Description : Special constants NaN, infinity etc + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Special constants NaN, infinity etc +-} + +module AERN2.MP.Float.Constants + ( + zero, one + , nan, infinity + ) +where + +import MixedTypesNumPrelude +import qualified Prelude as P +-- import Data.Ratio + +#ifdef UseCDAR +import AERN2.MP.Float.UseCDAR.Type +import AERN2.MP.Float.UseCDAR.Conversions +#else +import AERN2.MP.Float.UseRounded.Type +import AERN2.MP.Float.UseRounded.Conversions +#endif + +import AERN2.MP.Float.Operators + +zero, one :: MPFloat +zero = mpFloat 0 +one = mpFloat 1 + +nan, infinity :: MPFloat +nan = zero /. zero +infinity = one /. zero + +itisNaN :: MPFloat -> Bool +itisNaN x = x *^ one /= x + +itisInfinite :: MPFloat -> Bool +itisInfinite x = + x *^ (mpFloat 2) P.== x + && + x P./= (mpFloat 0) + +instance CanTestFinite MPFloat where + isInfinite = itisInfinite + isFinite x = not (itisInfinite x || itisNaN x) + +instance CanTestNaN MPFloat where + isNaN = itisNaN diff --git a/src/AERN2/MP/Float/Operators.hs b/src/AERN2/MP/Float/Operators.hs new file mode 100644 --- /dev/null +++ b/src/AERN2/MP/Float/Operators.hs @@ -0,0 +1,44 @@ +{-# LANGUAGE CPP #-} +{-| + Module : AERN2.MP.Float.Operators + Description : Infix operators for up/down-rounded floating-point numbers + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Infix operators for up/down-rounded floating-point numbers +-} + +module AERN2.MP.Float.Operators where + +#ifdef UseCDAR +import AERN2.MP.Float.UseCDAR.Type +import AERN2.MP.Float.UseCDAR.Arithmetic +#else +import AERN2.MP.Float.UseRounded.Type +import AERN2.MP.Float.UseRounded.Arithmetic +#endif + +infixl 6 +^, -^, +., -. +infixl 7 *^, *., /^, /. + +(+^) :: MPFloat -> MPFloat -> MPFloat +(+^) = addUp +(-^) :: MPFloat -> MPFloat -> MPFloat +(-^) = subUp +(*^) :: MPFloat -> MPFloat -> MPFloat +(*^) = mulUp +(/^) :: MPFloat -> MPFloat -> MPFloat +(/^) = divUp + +(+.) :: MPFloat -> MPFloat -> MPFloat +(+.) = addDown +(-.) :: MPFloat -> MPFloat -> MPFloat +(-.) = subDown +(*.) :: MPFloat -> MPFloat -> MPFloat +(*.) = mulDown +(/.) :: MPFloat -> MPFloat -> MPFloat +(/.) = divDown diff --git a/src/AERN2/MP/Float/Tests.hs b/src/AERN2/MP/Float/Tests.hs new file mode 100644 --- /dev/null +++ b/src/AERN2/MP/Float/Tests.hs @@ -0,0 +1,431 @@ +{-# LANGUAGE CPP #-} +{-| + Module : AERN2.MP.Float.Tests + Description : Tests for operations on arbitrary precision floats + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Tests for operations on arbitrary precision floats. + + To run the tests using stack, execute: + + @ + stack test aern2-mp --test-arguments "-a 1000 -m MPFloat" + @ +-} + +module AERN2.MP.Float.Tests + ( + specMPFloat, tMPFloat + , (=~=), approxEqual, approxEqualWithArgs + , frequencyElements + ) +where + +import MixedTypesNumPrelude +-- import qualified Prelude as P +-- import Data.Ratio +import Text.Printf +import Data.Maybe + +import Test.Hspec +import Test.QuickCheck +-- import qualified Test.Hspec.SmallCheck as SC + + +import AERN2.Norm +import AERN2.MP.Precision + +#ifdef UseCDAR +import AERN2.MP.Float.UseCDAR.Type +import AERN2.MP.Float.UseCDAR.Arithmetic +import AERN2.MP.Float.UseCDAR.Conversions +#else +import AERN2.MP.Float.UseRounded.Type +import AERN2.MP.Float.UseRounded.Arithmetic +import AERN2.MP.Float.UseRounded.Conversions +#endif + +import AERN2.MP.Float.Operators +import AERN2.MP.Float.Constants + +instance Arbitrary MPFloat where + arbitrary = + do + giveSpecialValue <- frequencyElements [(9, False),(1, True)] + aux giveSpecialValue + where + aux giveSpecialValue + | giveSpecialValue = + elements [nan, infinity, -infinity, zero, one, -one] + | otherwise = + do + (p :: Precision) <- arbitrary + (s :: Integer) <- arbitrary + ex <- choose (-20,10) + let resultR = s * (10.0^!ex) + let result = fromRationalUp p resultR + return result + +frequencyElements :: ConvertibleExactly t Int => [(t, a)] -> Gen a +frequencyElements elems = frequency [(int n, return e) | (n,e) <- elems] + +{- approximate comparison -} + +infix 4 =~= + +(=~=) :: MPFloat -> MPFloat -> Property +l =~= r = + approxEqualWithArgs [] l r + +{-| + Assert equality of two MPFloat's with tolerance @1/2^p@. +-} +approxEqual :: + Integer {-^ @p@ precision to guide tolerance -} -> + MPFloat {-^ LHS of equation-} -> + MPFloat {-^ RHS of equation -}-> + Bool +approxEqual e x y + | isNaN x && isNaN y = True + | isNaN x && isInfinite y = True + | isInfinite x && isNaN y = True + | isNaN x || isNaN y = False + | isInfinite x || isInfinite y = x == y + | otherwise = + abs (x -. y) <= 0.5^!e + +{-| + Assert equality of two MPFloat's with tolerance derived from the size and precision + of the given intermediate values. + When the assertion fails, report the given values using the given names. +-} +approxEqualWithArgs :: + [(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 +approxEqualWithArgs argsPre l r = + counterexample description $ approxEqual e l r + where + args = argsPre ++ [(l, "L"), (r, "R"), (abs (l-.r),"|L-R|")] + e = + (foldl min 1000000 $ catMaybes $ map getNminusP args) + - (length argsPre) + getNminusP (x,_) = + case norm of + NormZero -> Nothing -- ideally infinity + NormBits b -> Just (pI-b-1) + where + norm = getNormLog x + pI = integer $ getPrecision x + description = + printf "args:\n%s tolerance: <= %s (e=%d)" argsS (show (double (0.5^!e))) e + argsS = + unlines + [printf " %s = %s (p=%s)" argS (show arg) (show $ getPrecision arg) | (arg, argS) <- args] + +{-| + A runtime representative of type @MPFloat@. + Used for specialising polymorphic tests to concrete types. +-} +tMPFloat :: T MPFloat +tMPFloat = T "MPFloat" + +specMPFloat :: Spec +specMPFloat = + describe ("MPFloat") $ do + specCanSetPrecision tMPFloat (printArgsIfFails2 "=~=" (=~=)) + specCanRound tMPFloat + specCanNegNum tMPFloat + specCanAbs tMPFloat + specCanMinMaxNotMixed tMPFloat + -- specCanMinMax tMPFloat tInteger tMPFloat + describe "special values" $ do + it "0 * infinity = NaN" $ do + isNaN (zero *^ infinity) + && + isNaN (zero *. infinity) + it "infinity / infinity = NaN" $ do + isNaN (infinity /^ infinity) + && + isNaN (infinity /. infinity) + it "infinity - infinity = NaN" $ do + isNaN (infinity -^ infinity) + && + isNaN (infinity -. infinity) + describe "approximate addition" $ do + it "down <= up" $ do + property $ \ (x :: MPFloat) (y :: MPFloat) -> + not (isNaN (x +. y)) + ==> + x +. y <= x +^ y + it "up ~ down" $ do + property $ \ (x :: MPFloat) (y :: MPFloat) -> + x +. y =~= x +^ y + it "absorbs 0" $ do + property $ \ (x :: MPFloat) -> + (not $ isNaN x) ==> + x +. (mpFloat 0) == x + it "approximately commutative" $ do + property $ \ (x :: MPFloat) (y :: MPFloat) -> + (not $ isNaN $ x +. y) ==> + x +. y <= y +^ x + && + x +^ y >= y +. x + it "approximately associative" $ do + property $ \ (x :: MPFloat) (y :: MPFloat) (z :: MPFloat) -> + (not $ isNaN $ x +. y +. z) ==> + (x +. y) +. z <= x +^ (y +^ z) + && + (x +^ y) +^ z >= x +. (y +. z) + describe "approximate subtraction" $ do + it "down <= up" $ do + property $ \ (x :: MPFloat) (y :: MPFloat) -> + not (isNaN (x -. y)) + ==> + x -. y <= x -^ y + it "up ~ down" $ do + property $ \ (x :: MPFloat) (y :: MPFloat) -> + x -. y =~= x -^ y + it "same as negate and add" $ do + property $ \ (x :: MPFloat) (y :: MPFloat) -> + (not $ isNaN $ x -. y) ==> + x -. y <= x +^ (-y) + && + x -^ y >= x +. (-y) + describe "approximate multiplication" $ do + it "down <= up" $ do + property $ \ (x :: MPFloat) (y :: MPFloat) -> + not (isNaN (x *. y)) + ==> + x *. y <= x *^ y + it "up ~ down" $ do + property $ \ (x :: MPFloat) (y :: MPFloat) -> + x *. y =~= x *^ y + it "absorbs 1" $ do + property $ \ (x :: MPFloat) -> + (not $ isNaN x) ==> + x *. (mpFloat 1) == x + it "approximately commutative" $ do + property $ \ (x :: MPFloat) (y :: MPFloat) -> + not (isNaN (x *. y)) ==> + x *. y <= y *^ x + && + x *^ y >= y *. x + it "approximately associative" $ do + property $ \ (x :: MPFloat) (y :: MPFloat) (z :: MPFloat) -> + (x >= 0 && y >= 0 && z >= 0 + && not (isInfinite x) && not (isInfinite y) && not (isInfinite z)) ==> + (x *. y) *. z <= x *^ (y *^ z) + && + (x *^ y) *^ z >= x *. (y *. z) + it "approximately distributes over addition" $ do + property $ \ (x :: MPFloat) (y :: MPFloat) (z :: MPFloat) -> + (x >= 0 && y >= 0 && z >= 0 + && not (isInfinite x) && not (isInfinite y) && not (isInfinite z)) ==> + x *. (y +. z) <= (x *^ y) +^ (x *^ z) + && + x *^ (y +^ z) >= (x *. y) +. (x *. z) + describe "approximate division" $ do + it "down <= up" $ do + property $ \ (x :: MPFloat) (y :: MPFloat) -> + not (isNaN (x /. y)) + ==> + x /. y <= x /^ y + it "up ~ down" $ do + property $ \ (x :: MPFloat) (y :: MPFloat) -> + let + (=~~=) = approxEqualWithArgs [(x /. y,"x/.y")] + infix 4 =~~= + in + not (isNaN (x /. y)) + ==> + x /. y =~~= x /^ y + it "recip(recip x) = x" $ do + property $ \ (x :: MPFloat) -> + (x > 0 || x < 0) ==> + one /. (one /^ x) <= x + && + one /^ (one /. x) >= x + it "x/1 = x" $ do + property $ \ (x :: MPFloat) -> + not (isNaN x) ==> + (x /. one) == x + it "x/x = 1" $ do + property $ \ (x :: MPFloat) -> + (isCertainlyNonZero x && (not $ isNaN $ x /. x)) ==> + (x /. x) <= one + && + (x /^ x) >= one + it "x/y = x*(1/y)" $ do + property $ \ (x :: MPFloat) (y :: MPFloat) -> + (y > 0 && x >= 0 && x/.y >= 0) ==> + (x /. y) <= x *^ (one /^ y) + && + (x /^ y) >= x *. (one /. y) + describe "approximate sqrt" $ do + it "down <= up" $ do + property $ \ (x :: MPFloat) -> + not (isNaN (sqrtDown x)) + ==> + sqrtDown x <= sqrtUp x + it "up ~ down" $ do + property $ \ (x :: MPFloat) -> + (x >= 0) + ==> + sqrtDown x =~= sqrtUp x + it "sqrt(x) >= 0" $ do + property $ \ (x :: MPFloat) -> + (x >= 0) + ==> + sqrtUp x >= 0 + it "sqrt(x)^2 ~ x" $ do + property $ \ (x :: MPFloat) -> + (x >= 0) + ==> + (sqrtDown x) *. (sqrtDown x) <= x + && + (sqrtUp x) *^ (sqrtUp x) >= x + describe "approximate exp" $ do + it "down <= up" $ do + property $ \ (x :: MPFloat) -> + (abs x < 1000000) + ==> + expDown x <= expUp x + it "up ~ down" $ do + property $ \ (x :: MPFloat) -> + (abs x < 1000000) + ==> + let + (=~~=) = approxEqualWithArgs [(x,"x")] + infix 4 =~~= + in + expDown x =~~= expUp x + it "exp(-x) == 1/(exp x)" $ do + property $ \ (x :: MPFloat) -> + (abs x < 1000000) + ==> + one /. (expUp x) <= expUp (-x) + && + one /^ (expDown x) >= expDown (-x) + it "exp(x+y) = exp(x)*exp(y)" $ do + property $ \ (x :: MPFloat) (y :: MPFloat) -> + (abs x < 1000000 && abs y < 1000000) + ==> + expDown (x +. y) <= (expUp x) *^ (expUp y) + && + expUp (x +^ y) >= (expDown x) *. (expDown y) + describe "approximate log" $ do + it "down <= up" $ do + property $ \ (x :: MPFloat) -> + (x > 0) + ==> + logDown x <= logUp x + it "up ~ down" $ do + property $ \ (x :: MPFloat) -> + (x > 0) + ==> + logDown x =~= logUp x + it "log(1/x) == -(log x)" $ do + property $ \ (x :: MPFloat) -> + (x > 0) + ==> + logDown (one /. x) <= -(logDown x) + && + logUp (one /^ x) >= -(logUp x) + it "log(x*y) = log(x)+log(y)" $ do + property $ \ (x :: MPFloat) (y :: MPFloat) -> + (x > 0 && y > 0) + ==> + logDown (x *. y) <= (logUp x) +^ (logUp y) + && + logUp (x *^ y) >= (logDown x) +. (logDown y) + it "log(exp x) == x" $ do + property $ \ (x :: MPFloat) -> + (abs x < 1000000) + ==> + logDown (expDown x) <= x + && + logUp (expUp x) >= x + describe "approximate sine" $ do + it "down <= up" $ do + property $ \ (x :: MPFloat) -> + (abs x < 1000000) + ==> + sinDown x <= sinUp x + it "up ~ down" $ do + property $ \ (x :: MPFloat) -> + (abs x < 1000000) + ==> + let + (=~~=) = approxEqualWithArgs [(x,"x")] + infix 4 =~~= + in + sinDown x =~~= sinUp x + it "sin(pi)=0" $ do + property $ \ (p :: Precision) -> + let + (=~~=) = approxEqualWithArgs [(piDown p,"pi")] + infix 4 =~~= + in + sinUp(piDown p) =~~= (fromIntegerUp p 0) + it "in [-1,1]" $ do + property $ \ (x :: MPFloat) -> + (abs x < 1000000) + ==> + sinDown x <= one + && + sinUp x >= -one + describe "approximate cosine" $ do + it "down <= up" $ do + property $ \ (x :: MPFloat) -> + (abs x < 1000000) + ==> + cosDown x <= cosUp x + it "up ~ down" $ do + property $ \ (x :: MPFloat) -> + (abs x < 1000000) + ==> + let + (=~~=) = approxEqualWithArgs [(x,"x")] + infix 4 =~~= + in + cosDown x =~~= cosUp x + it "in [-1,1]" $ do + property $ \ (x :: MPFloat) -> + (abs x < 1000000) + ==> + cosDown x <= one + && + cosUp x >= -one + it "cos(pi)=-1" $ do + property $ \ (p :: Precision) -> + cosUp(piDown p) =~= (fromIntegerUp p (-1)) + it "cos(x)^2 + sin(x)^2 = 1" $ do + property $ \ (x :: MPFloat) -> + (abs x < 1000000) + ==> + let + cosxU = cosUp x + cosxD = cosDown x + cosx2U = (cosxU *^ cosxU) `max` (cosxD *^ cosxD) + cosx2D + | cosxD > 0 = cosxD *. cosxD + | cosxU < 0 = cosxU *. cosxU + | otherwise = mpFloat 0 + sinxU = sinUp x + sinxD = sinDown x + sinx2U = (sinxU *^ sinxU) `max` (sinxD *^ sinxD) + sinx2D + | sinxD > 0 = sinxD *. sinxD + | sinxU < 0 = sinxU *. sinxU + | otherwise = mpFloat 0 + in + (cosx2D +. sinx2D) <= one + && + (cosx2U +^ sinx2U) >= one diff --git a/src/AERN2/MP/Float/UseRounded/Arithmetic.hs b/src/AERN2/MP/Float/UseRounded/Arithmetic.hs new file mode 100644 --- /dev/null +++ b/src/AERN2/MP/Float/UseRounded/Arithmetic.hs @@ -0,0 +1,151 @@ +{-| + Module : AERN2.MP.Float.UseRounded.Arithmetic + Description : Arbitrary precision floating point numbers + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Arbitrary precision floating-point numbers with up/down-rounded operations. + + Currently, we use hmpfr when compiling with ghc 7.10 and higher + and haskell-mpfr when compiling with ghc 7.8. +-} + +module AERN2.MP.Float.UseRounded.Arithmetic + ( + -- * MPFloat basic arithmetic + addUp, addDown, subUp, subDown + , mulUp, mulDown, divUp, divDown, recipUp, recipDown + -- * MPFloat selected constants and operations + , piUp, piDown + , cosUp, cosDown, sinUp, sinDown + , sqrtUp, sqrtDown, expUp, expDown, logUp, logDown + ) +where + +import MixedTypesNumPrelude +import qualified Prelude as P + +import AERN2.MP.Precision + +import qualified AERN2.MP.Float.UseRounded.RoundedAdaptor as MPLow +import AERN2.MP.Float.UseRounded.Type + +one :: MPFloat +one = MPLow.one + +{- common functions -} + +instance CanNeg MPFloat where + negate = unaryUp MPLow.neg + +instance CanAbs MPFloat where + abs x + | x P.< MPLow.zero = negate x + | otherwise = x + +addUp, addDown :: MPFloat -> MPFloat -> MPFloat +addUp = binaryUp True MPLow.add +addDown = binaryDown True MPLow.add + +subUp, subDown :: MPFloat -> MPFloat -> MPFloat +subUp = binaryUp True MPLow.sub +subDown = binaryDown True MPLow.sub + +mulUp, mulDown :: MPFloat -> MPFloat -> MPFloat +mulUp = binaryUp True MPLow.mul +mulDown = binaryDown True MPLow.mul + +divUp,divDown :: MPFloat -> MPFloat -> MPFloat +divUp = binaryUp False MPLow.div +divDown = binaryDown False MPLow.div + +recipUp :: MPFloat -> MPFloat +recipUp x = divUp one x + +recipDown :: MPFloat -> MPFloat +recipDown x = divDown one x + + +{- special constants and functions -} + +piUp :: Precision -> MPFloat +piUp p = + MPLow.pi MPLow.Up (p2mpfrPrec p) + +piDown :: Precision -> MPFloat +piDown p = + MPLow.pi MPLow.Down (p2mpfrPrec p) + +cosUp :: MPFloat -> MPFloat +cosUp = unaryUp MPLow.cos + +cosDown :: MPFloat -> MPFloat +cosDown = unaryDown MPLow.cos + +sinUp :: MPFloat -> MPFloat +sinUp = unaryUp MPLow.sin + +sinDown :: MPFloat -> MPFloat +sinDown = unaryDown MPLow.sin + +sqrtUp :: MPFloat -> MPFloat +sqrtUp = unaryUp MPLow.sqrt + +sqrtDown :: MPFloat -> MPFloat +sqrtDown = unaryDown MPLow.sqrt + +expUp :: MPFloat -> MPFloat +expUp = unaryUp MPLow.exp + +expDown :: MPFloat -> MPFloat +expDown = unaryDown MPLow.exp + +logUp :: MPFloat -> MPFloat +logUp = unaryUp MPLow.log + +logDown :: MPFloat -> MPFloat +logDown = unaryDown MPLow.log + +{- auxiliary functions to automatically determine result precision from operand precisions -} + +unaryUp :: + (MPLow.RoundMode -> MPLow.Precision -> MPFloat -> MPFloat) -> + (MPFloat -> MPFloat) +unaryUp opRP x = opRP MPLow.Up p x + where + p = MPLow.getPrec x + +unaryDown :: + (MPLow.RoundMode -> MPLow.Precision -> MPFloat -> MPFloat) -> + (MPFloat -> MPFloat) +unaryDown opRP x = opRP MPLow.Down p x + where + p = MPLow.getPrec x + +binaryUp :: + Bool -> + (MPLow.RoundMode -> MPLow.Precision -> MPFloat -> MPFloat -> MPFloat) -> + (MPFloat -> MPFloat -> MPFloat) +binaryUp = binaryApprox True + +binaryDown :: + Bool -> + (MPLow.RoundMode -> MPLow.Precision -> MPFloat -> MPFloat -> MPFloat) -> + (MPFloat -> MPFloat -> MPFloat) +binaryDown = binaryApprox False + +binaryApprox :: + Bool -> Bool -> + (MPLow.RoundMode -> MPLow.Precision -> MPFloat -> MPFloat -> MPFloat) -> + (MPFloat -> MPFloat -> MPFloat) +binaryApprox isUp _canBeExact opRP x y = + withPrec pMax + where + pMax = (getPrecision x) `max` (getPrecision y) + withPrec p + | isUp = opRP MPLow.Up (p2mpfrPrec p) x y + | otherwise = opRP MPLow.Down (p2mpfrPrec p) x y diff --git a/src/AERN2/MP/Float/UseRounded/Conversions.hs b/src/AERN2/MP/Float/UseRounded/Conversions.hs new file mode 100644 --- /dev/null +++ b/src/AERN2/MP/Float/UseRounded/Conversions.hs @@ -0,0 +1,161 @@ +{-| + Module : AERN2.MP.Float.UseRounded.Conversions + Description : Conversions and comparisons of arbitrary precision floats + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Conversions and comparisons of arbitrary precision floating point numbers + + Currently, we use hmpfr when compiling with ghc 7.10 and higher + and haskell-mpfr when compiling with ghc 7.8. +-} + +module AERN2.MP.Float.UseRounded.Conversions + ( + -- * MPFloat to other types (see also instances) + toDoubleUp, toDoubleDown + -- * MPFloat constructors (see also instances) + , CanBeMPFloat, mpFloat + , fromIntegerUp, fromIntegerDown + , fromRationalUp, fromRationalDown + ) +where + +import MixedTypesNumPrelude +import qualified Prelude as P + +import Data.Ratio +import Data.Convertible + +import AERN2.Norm +import AERN2.MP.Precision + +import AERN2.MP.Float.UseRounded.Type +import AERN2.MP.Float.UseRounded.Arithmetic + +import qualified AERN2.MP.Float.UseRounded.RoundedAdaptor as MPLow + +mpToDouble :: MPLow.RoundMode -> MPFloat -> Double +mpToDouble = MPLow.toDoubleA + +mpToRational :: MPFloat -> Rational +mpToRational x + | x == 0 = 0.0 + | otherwise = MPLow.toRationalA x + +mpFromRationalA :: MPLow.RoundMode -> MPLow.Precision -> Rational -> MPFloat +mpFromRationalA = MPLow.fromRationalA + +instance HasNorm MPFloat where + getNormLog x + | x == 0 = NormZero + | otherwise = NormBits (P.toInteger $ MPLow.getExp x) + +{- conversions -} + +instance CanRound MPFloat where + properFraction x = (n,f) + where + r = rational x + n = (numerator r) `quot` (denominator r) + f = x `subUp` (mpFloat n) + +instance ConvertibleExactly MPFloat Rational where + safeConvertExactly = Right . mpToRational + +toDoubleUp :: MPFloat -> Double +toDoubleUp = mpToDouble MPLow.Up + +toDoubleDown :: MPFloat -> Double +toDoubleDown = mpToDouble MPLow.Down + +fromIntegerUp :: Precision -> Integer -> MPFloat +fromIntegerUp p i = MPLow.fromIntegerA MPLow.Up (p2mpfrPrec p) i + +fromIntegerDown :: Precision -> Integer -> MPFloat +fromIntegerDown p i = MPLow.fromIntegerA MPLow.Down (p2mpfrPrec p) i + +type CanBeMPFloat t = ConvertibleExactly t MPFloat +mpFloat :: (CanBeMPFloat t) => t -> MPFloat +mpFloat = convertExactly + +instance ConvertibleExactly Integer MPFloat where + safeConvertExactly n = + findExact $ map upDown $ standardPrecisions initPrec + where + initPrec = + case getNormLog n of + NormBits b -> prec (b + 8) + _ -> prec 8 + upDown p = (fromIntegerDown p n, fromIntegerUp p n) + findExact [] = + convError "integer too high to represent exactly" n + findExact ((nDown, nUp) : rest) + | nDown == nUp = Right nUp + | otherwise = findExact rest + +instance ConvertibleExactly Int MPFloat where + safeConvertExactly = safeConvertExactly . integer + +fromRationalUp :: Precision -> Rational -> MPFloat +fromRationalUp p x = + mpFromRationalA MPLow.Up (p2mpfrPrec p) x + +fromRationalDown :: Precision -> Rational -> MPFloat +fromRationalDown p x = + mpFromRationalA MPLow.Down (p2mpfrPrec p) x + +instance Convertible MPFloat Double where + safeConvert x + | isFinite dbl = Right dbl + | otherwise = convError "conversion to double: out of bounds" x + where + dbl = toDoubleUp x + +{- comparisons -} + +instance HasEqAsymmetric MPFloat MPFloat +instance HasEqAsymmetric MPFloat Integer where + equalTo = convertSecond equalTo +instance HasEqAsymmetric Integer MPFloat where + equalTo = convertFirst equalTo +instance HasEqAsymmetric MPFloat Int where + equalTo = convertSecond equalTo +instance HasEqAsymmetric Int MPFloat where + equalTo = convertFirst equalTo +instance HasEqAsymmetric MPFloat Rational where + equalTo = convertFirst equalTo +instance HasEqAsymmetric Rational MPFloat where + equalTo = convertSecond equalTo + +instance CanTestZero MPFloat + +instance HasOrderAsymmetric MPFloat MPFloat +instance HasOrderAsymmetric MPFloat Integer where + lessThan = convertSecond lessThan + leq = convertSecond leq +instance HasOrderAsymmetric Integer MPFloat where + lessThan = convertFirst lessThan + leq = convertFirst leq +instance HasOrderAsymmetric MPFloat Int where + lessThan = convertSecond lessThan + leq = convertSecond leq +instance HasOrderAsymmetric Int MPFloat where + lessThan = convertFirst lessThan + leq = convertFirst leq +instance HasOrderAsymmetric Rational MPFloat where + lessThan = convertSecond lessThan + leq = convertSecond leq +instance HasOrderAsymmetric MPFloat Rational where + lessThan = convertFirst lessThan + leq = convertFirst leq + +instance CanTestPosNeg MPFloat + +{- min, max -} + +instance CanMinMaxAsymmetric MPFloat MPFloat diff --git a/src/AERN2/MP/Float/UseRounded/RoundedAdaptor.hs b/src/AERN2/MP/Float/UseRounded/RoundedAdaptor.hs new file mode 100644 --- /dev/null +++ b/src/AERN2/MP/Float/UseRounded/RoundedAdaptor.hs @@ -0,0 +1,84 @@ +{-# LANGUAGE DataKinds, ExistentialQuantification, RankNTypes #-} +-- {-# LANGUAGE DeriveGeneric, DeriveDataTypeable, StandaloneDeriving #-} +{-| + Module : AERN2.MP.Float.UseRounded.RoundedAdaptor + Description : Numeric.Rounded + variable precision + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Numeric.Rounded + variable precision +-} +module AERN2.MP.Float.UseRounded.RoundedAdaptor +( + module AERN2.MP.Float.UseRounded.RoundedAdaptor +, module Numeric.Rounded.Simple +) +where + +import Prelude hiding (div, pi) +-- import qualified Prelude as P + +import Numeric.Rounded.Simple +-- import qualified Numeric.RoundedSimple as R + +instance Show Rounded where + show = show' + +getPrec :: Rounded -> Int +getPrec = precision + +getExp :: Rounded -> Int +getExp = exponent' + +data RoundMode = Up | Down + +withRoundMode :: (RoundingMode -> t) -> (RoundMode -> t) +withRoundMode op Up = op TowardInf +withRoundMode op Down = op TowardNegInf +{-# INLINE withRoundMode #-} + +set :: RoundMode -> Precision -> Rounded -> Rounded +set = withRoundMode precRound + +defaultPrecision :: Precision +defaultPrecision = 10 + +pi :: RoundMode -> Precision -> Rounded +pi = withRoundMode kPi + +fromIntegerA :: RoundMode -> Precision -> Integer -> Rounded +fromIntegerA = withRoundMode fromInteger' + +zero, one :: Rounded +zero = fromIntegerA Up defaultPrecision 0 +one = fromIntegerA Up defaultPrecision 1 + +toDoubleA :: RoundMode -> Rounded -> Double +toDoubleA = withRoundMode toDouble + +fromRationalA :: RoundMode -> Precision -> Rational -> Rounded +fromRationalA = withRoundMode fromRational' + +toRationalA :: Rounded -> Rational +toRationalA = toRational' TowardNearest + +add, sub, mul, div, atan2 :: RoundMode -> Precision -> Rounded -> Rounded -> Rounded +add = withRoundMode add_ +sub = withRoundMode sub_ +mul = withRoundMode mul_ +div = withRoundMode div_ +atan2 = withRoundMode atan2_ + +neg, abs, sqrt, exp, log, sin, cos :: RoundMode -> Precision -> Rounded -> Rounded +neg = withRoundMode negate_ +abs = withRoundMode abs_ +sqrt = withRoundMode sqrt_ +exp = withRoundMode exp_ +log = withRoundMode log_ +sin = withRoundMode sin_ +cos = withRoundMode cos_ +-- TODO: add more ops diff --git a/src/AERN2/MP/Float/UseRounded/Type.hs b/src/AERN2/MP/Float/UseRounded/Type.hs new file mode 100644 --- /dev/null +++ b/src/AERN2/MP/Float/UseRounded/Type.hs @@ -0,0 +1,50 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE DeriveGeneric, DeriveDataTypeable, StandaloneDeriving #-} +{-| + Module : AERN2.MP.Float.UseRounded.Type + Description : Arbitrary precision floating point numbers (MPFR) + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Arbitrary precision floating-point numbers using MPFR via package rounded. +-} + +module AERN2.MP.Float.UseRounded.Type + ( + -- * MPFloat numbers and their basic operations + MPFloat, setPrecisionUp, setPrecisionDown + , p2mpfrPrec + ) +where + +import MixedTypesNumPrelude +import qualified Prelude as P + +import AERN2.MP.Precision + +import qualified AERN2.MP.Float.UseRounded.RoundedAdaptor as MPLow +import Data.Typeable + +{-| Multiple-precision floating-point type based on MPFR via rounded. -} +type MPFloat = MPLow.Rounded + +deriving instance (Typeable MPFloat) + +p2mpfrPrec :: Precision -> MPLow.Precision +p2mpfrPrec = P.fromInteger . integer + +instance HasPrecision MPFloat where + getPrecision x = prec (P.toInteger $ MPLow.getPrec x) + +instance CanSetPrecision MPFloat where + setPrecision = setPrecisionUp + +setPrecisionUp :: Precision -> MPFloat -> MPFloat +setPrecisionUp p = MPLow.set MPLow.Up (p2mpfrPrec p) + +setPrecisionDown :: Precision -> MPFloat -> MPFloat +setPrecisionDown p = MPLow.set MPLow.Down (p2mpfrPrec p) diff --git a/src/AERN2/MP/UseMPFR/Ball.hs b/src/AERN2/MP/UseMPFR/Ball.hs deleted file mode 100644 --- a/src/AERN2/MP/UseMPFR/Ball.hs +++ /dev/null @@ -1,64 +0,0 @@ -{-| - Module : AERN2.MP.UseMPFR.Ball - Description : Arbitrary precision ball arithmetic - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mikkonecny@gmail.com - Stability : experimental - Portability : portable - - Arbitrary precision ball arithmetic --} -module AERN2.MP.UseMPFR.Ball -( - -- * Auxiliary types - module AERN2.Norm - , module AERN2.MP.Precision - , module AERN2.MP.Accuracy - , module AERN2.MP.UseMPFR.ErrorBound - , module AERN2.MP.Enclosure - -- * The Ball type - , MPBall(..), CanBeMPBall, mpBall, CanBeMPBallP, mpBallP - , reducePrecionIfInaccurate - -- * Ball construction/extraction functions - -- , endpointsMP, fromEndpointsMP - -- * Ball operations (see also instances) - , piBallP - -- * Helpers for constructing ball functions - , byEndpointsMP - , fromApproxWithLipschitz -) -where - -import MixedTypesNumPrelude --- import qualified Prelude as P - -import AERN2.Norm -import AERN2.MP.Precision -import AERN2.MP.Accuracy -import AERN2.MP.Enclosure - -import AERN2.MP.UseMPFR.ErrorBound (ErrorBound, CanBeErrorBound, errorBound) - -import AERN2.MP.UseMPFR.Ball.Type -import AERN2.MP.UseMPFR.Ball.Conversions () -import AERN2.MP.UseMPFR.Ball.Comparisons -import AERN2.MP.UseMPFR.Ball.Field () -import AERN2.MP.UseMPFR.Ball.Elementary -import AERN2.MP.UseMPFR.Ball.PreludeOps () - -instance Ring MPBall -instance Ring (CN MPBall) -instance Field MPBall -instance Field (CN MPBall) - -instance OrderedRing MPBall -instance OrderedRing (CN MPBall) -instance OrderedField MPBall -instance OrderedField (CN MPBall) - -instance OrderedCertainlyRing MPBall -instance OrderedCertainlyRing (CN MPBall) -instance OrderedCertainlyField MPBall -instance OrderedCertainlyField (CN MPBall) diff --git a/src/AERN2/MP/UseMPFR/Ball/Comparisons.hs b/src/AERN2/MP/UseMPFR/Ball/Comparisons.hs deleted file mode 100644 --- a/src/AERN2/MP/UseMPFR/Ball/Comparisons.hs +++ /dev/null @@ -1,381 +0,0 @@ -{-| - Module : AERN2.MP.UseMPFR.Ball.Comparisons - Description : Comparisons of arbitrary precision dyadic balls - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mikkonecny@gmail.com - Stability : experimental - Portability : portable - - Comparisons of arbitrary precision dyadic balls --} -module AERN2.MP.UseMPFR.Ball.Comparisons -( - -- * Auxiliary types - module AERN2.Norm - -- * Ball operations (see also instances) - , reducePrecionIfInaccurate - -- * Helpers for constructing ball functions - , byEndpointsMP -) -where - -import MixedTypesNumPrelude --- import qualified Prelude as P - -import Control.CollectErrors - -import AERN2.Norm -import AERN2.MP.Dyadic (Dyadic) -import AERN2.MP.UseMPFR.Float (MPFloat) --- import AERN2.MP.UseMPFR.Float.Operators -import AERN2.MP.Precision - -import AERN2.MP.UseMPFR.Ball.Type -import AERN2.MP.UseMPFR.Ball.Conversions () - -{- comparisons -} - -instance HasEqAsymmetric MPBall MPBall where - type EqCompareType MPBall MPBall = Maybe Bool - b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 - -instance HasEqAsymmetric MPBall Integer where - type EqCompareType MPBall Integer = Maybe Bool - b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 -instance HasEqAsymmetric Integer MPBall where - type EqCompareType Integer MPBall = Maybe Bool - b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 - -instance HasEqAsymmetric MPBall Int where - type EqCompareType MPBall Int = Maybe Bool - b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 -instance HasEqAsymmetric Int MPBall where - type EqCompareType Int MPBall = Maybe Bool - b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 - -instance HasEqAsymmetric MPBall Rational where - type EqCompareType MPBall Rational = Maybe Bool - b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 -instance HasEqAsymmetric Rational MPBall where - type EqCompareType Rational MPBall = Maybe Bool - b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 - -instance HasEqAsymmetric MPBall Dyadic where - type EqCompareType MPBall Dyadic = Maybe Bool - b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 -instance HasEqAsymmetric Dyadic MPBall where - type EqCompareType Dyadic MPBall = Maybe Bool - b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 - -instance - (HasEqAsymmetric MPBall b - , CanEnsureCE es b - , CanEnsureCE es (EqCompareType MPBall b) - , IsBool (EnsureCE es (EqCompareType MPBall b)) - , SuitableForCE es) - => - HasEqAsymmetric MPBall (CollectErrors es b) - where - type EqCompareType MPBall (CollectErrors es b) = - EnsureCE es (EqCompareType MPBall b) - equalTo = lift2TLCE equalTo - -instance - (HasEqAsymmetric a MPBall - , CanEnsureCE es a - , CanEnsureCE es (EqCompareType a MPBall) - , IsBool (EnsureCE es (EqCompareType a MPBall)) - , SuitableForCE es) - => - HasEqAsymmetric (CollectErrors es a) MPBall - where - type EqCompareType (CollectErrors es a) MPBall = - EnsureCE es (EqCompareType a MPBall) - equalTo = lift2TCE equalTo - -instance HasOrderAsymmetric MPBall MPBall where - type OrderCompareType MPBall MPBall = Maybe Bool - lessThan b1 b2 - | r1 < l2 = Just True - | r2 <= l1 = Just False - | otherwise = Nothing - where - (l1, r1) = endpointsMP b1 - (l2, r2) = endpointsMP b2 - leq b1 b2 - | r1 <= l2 = Just True - | r2 < l1 = Just False - | otherwise = Nothing - where - (l1, r1) = endpointsMP b1 - (l2, r2) = endpointsMP b2 - -instance HasOrderAsymmetric Integer MPBall where - type OrderCompareType Integer MPBall = Maybe Bool - lessThan = convertFirst lessThan - leq = convertFirst leq -instance HasOrderAsymmetric MPBall Integer where - type OrderCompareType MPBall Integer = Maybe Bool - lessThan = convertSecond lessThan - leq = convertSecond leq - -instance HasOrderAsymmetric Int MPBall where - type OrderCompareType Int MPBall = Maybe Bool - lessThan = convertFirst lessThan - leq = convertFirst leq -instance HasOrderAsymmetric MPBall Int where - type OrderCompareType MPBall Int = Maybe Bool - lessThan = convertSecond lessThan - leq = convertSecond leq - -instance HasOrderAsymmetric Dyadic MPBall where - type OrderCompareType Dyadic MPBall = Maybe Bool - lessThan = convertFirst lessThan - leq = convertFirst leq -instance HasOrderAsymmetric MPBall Dyadic where - type OrderCompareType MPBall Dyadic = Maybe Bool - lessThan = convertSecond lessThan - leq = convertSecond leq - -instance HasOrderAsymmetric MPBall Rational where - type OrderCompareType MPBall Rational = Maybe Bool - lessThan b1 q2 - | r1 < l2 = Just True - | r2 <= l1 = Just False - | otherwise = Nothing - where - (l1, r1) = endpointsMP b1 - l2 = q2 - r2 = q2 - leq b1 q2 - | r1 <= l2 = Just True - | r2 < l1 = Just False - | otherwise = Nothing - where - (l1, r1) = endpointsMP b1 - l2 = q2 - r2 = q2 - -instance HasOrderAsymmetric Rational MPBall where - type OrderCompareType Rational MPBall = Maybe Bool - lessThan q1 b2 - | r1 < l2 = Just True - | r2 <= l1 = Just False - | otherwise = Nothing - where - (l2, r2) = endpointsMP b2 - l1 = q1 - r1 = q1 - leq q1 b2 - | r1 <= l2 = Just True - | r2 < l1 = Just False - | otherwise = Nothing - where - (l2, r2) = endpointsMP b2 - l1 = q1 - r1 = q1 - -instance - (HasOrderAsymmetric MPBall b - , CanEnsureCE es b - , CanEnsureCE es (OrderCompareType MPBall b) - , IsBool (EnsureCE es (OrderCompareType MPBall b)) - , SuitableForCE es) - => - HasOrderAsymmetric MPBall (CollectErrors es b) - where - type OrderCompareType MPBall (CollectErrors es b) = - EnsureCE es (OrderCompareType MPBall b) - lessThan = lift2TLCE lessThan - leq = lift2TLCE leq - greaterThan = lift2TLCE greaterThan - geq = lift2TLCE geq - -instance - (HasOrderAsymmetric a MPBall - , CanEnsureCE es a - , CanEnsureCE es (OrderCompareType a MPBall) - , IsBool (EnsureCE es (OrderCompareType a MPBall)) - , SuitableForCE es) - => - HasOrderAsymmetric (CollectErrors es a) MPBall - where - type OrderCompareType (CollectErrors es a) MPBall = - EnsureCE es (OrderCompareType a MPBall) - lessThan = lift2TCE lessThan - leq = lift2TCE leq - greaterThan = lift2TCE greaterThan - geq = lift2TCE geq - -instance CanTestZero MPBall -instance CanTestPosNeg MPBall - -instance CanTestInteger MPBall where - certainlyNotInteger b = - (rN - lN) == 1 && lN ! - CanMinMaxAsymmetric MPBall (CollectErrors es b) - where - type MinMaxType MPBall (CollectErrors es b) = - EnsureCE es (MinMaxType MPBall b) - min = lift2TLCE min - max = lift2TLCE max - -instance - (CanMinMaxAsymmetric a MPBall - , CanEnsureCE es a - , CanEnsureCE es (MinMaxType a MPBall) - , SuitableForCE es) - => - CanMinMaxAsymmetric (CollectErrors es a) MPBall - where - type MinMaxType (CollectErrors es a) MPBall = - EnsureCE es (MinMaxType a MPBall) - min = lift2TCE min - max = lift2TCE max - -{- intersection -} - -instance CanIntersectAsymmetric MPBall MPBall where - intersect a b - | rL > rR = - noValueNumErrorCertainCN $ NumError $ "intersect: empty intersection: " ++ show a ++ "; " ++ show b - | otherwise = cn $ fromEndpointsMP rL rR - where - rL = max aL bL - rR = min aR bR - (aL,aR) = endpointsMP a - (bL,bR) = endpointsMP b - -instance - (CanIntersectAsymmetric MPBall b - , CanEnsureCE es b - , CanEnsureCE es (IntersectionType MPBall b) - , SuitableForCE es) - => - CanIntersectAsymmetric MPBall (CollectErrors es b) - where - type IntersectionType MPBall (CollectErrors es b) = - EnsureCE es (IntersectionType MPBall b) - intersect = lift2TLCE intersect - -instance - (CanIntersectAsymmetric a MPBall - , CanEnsureCE es a - , CanEnsureCE es (IntersectionType a MPBall) - , SuitableForCE es) - => - CanIntersectAsymmetric (CollectErrors es a) MPBall - where - type IntersectionType (CollectErrors es a) MPBall = - EnsureCE es (IntersectionType a MPBall) - intersect = lift2TCE intersect - -{- union -} - -instance CanUnionAsymmetric MPBall MPBall where - union a b = - case getMaybeValueCN (a `intersect` b) of - Just _ -> prependErrorsCN [(ErrorPotential, err)] r - _ -> prependErrorsCN [(ErrorCertain, err)] r - where - err = NumError $ "union of enclosures: not enclosing the same value" - r = cn $ fromEndpointsMP rL rR - rL = min aL bL - rR = max aR bR - (aL,aR) = endpointsMP a - (bL,bR) = endpointsMP b - - -instance - (CanUnionAsymmetric MPBall b - , CanEnsureCE es b - , CanEnsureCE es (UnionType MPBall b) - , SuitableForCE es) - => - CanUnionAsymmetric MPBall (CollectErrors es b) - where - type UnionType MPBall (CollectErrors es b) = - EnsureCE es (UnionType MPBall b) - union = lift2TLCE union - -instance - (CanUnionAsymmetric a MPBall - , CanEnsureCE es a - , CanEnsureCE es (UnionType a MPBall) - , SuitableForCE es) - => - CanUnionAsymmetric (CollectErrors es a) MPBall - where - type UnionType (CollectErrors es a) MPBall = - EnsureCE es (UnionType a MPBall) - union = lift2TCE union - -{-| - Computes an *increasing* ball fucntion @f@ from *exact* MPFR operations. --} -byEndpointsMP :: - (MPFloat -> MPFloat -> MPFloat) -> - (MPBall -> MPBall -> MPBall) -byEndpointsMP op b1 b2 = - fromEndpointsMP (l1 `op` l2) (r1 `op` r2) - where - (l1,r1) = endpointsMP b1 - (l2,r2) = endpointsMP b2 - -{- random generation -} diff --git a/src/AERN2/MP/UseMPFR/Ball/Conversions.hs b/src/AERN2/MP/UseMPFR/Ball/Conversions.hs deleted file mode 100644 --- a/src/AERN2/MP/UseMPFR/Ball/Conversions.hs +++ /dev/null @@ -1,119 +0,0 @@ -{-| - Module : AERN2.MP.UseMPFR.Ball.Conversions - Description : Conversions of arbitrary precision dyadic balls - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mikkonecny@gmail.com - Stability : experimental - Portability : portable - - Conversions of arbitrary precision dyadic balls --} -module AERN2.MP.UseMPFR.Ball.Conversions -( - integerBounds -) -where - -import MixedTypesNumPrelude --- import qualified Prelude as P - -import Data.Typeable -import Data.Convertible - -import AERN2.MP.Dyadic (Dyadic, dyadic) -import qualified AERN2.MP.UseMPFR.Float as MPFloat -import AERN2.MP.UseMPFR.Float (mpFloat) --- import AERN2.MP.UseMPFR.Float.Operators -import AERN2.MP.Precision -import qualified AERN2.MP.UseMPFR.ErrorBound as EB -import AERN2.MP.UseMPFR.ErrorBound (errorBound) - -import AERN2.MP.UseMPFR.Ball.Type - -{--- extracting from a ball ---} - -instance HasIntegerBounds MPBall where - integerBounds b = - (floor l, ceiling r) - where - (l,r) = endpointsMP b - -instance Convertible MPBall EB.ErrorBound where - safeConvert b = - Right (errorBound (max (abs l) (abs r))) - where - (l,r) = endpointsMP b - -{--- constructing an exact ball ---} - -instance ConvertibleExactly MPBall MPBall where - safeConvertExactly = Right - -instance ConvertibleExactly Dyadic MPBall where - safeConvertExactly x = Right $ MPBall (mpFloat x) (errorBound 0) - -instance ConvertibleExactly EB.ErrorBound MPBall where - safeConvertExactly eb = Right $ MPBall (mpFloat eb) (errorBound 0) - -instance - (ConvertibleExactly c Dyadic, ConvertibleExactly e Dyadic - , Show c, Show e, Typeable c, Typeable e) - => - ConvertibleExactly (c, e) MPBall - where - safeConvertExactly (c,e) - | isFinite b = Right b - | otherwise = convError "too large to convert to MPBall" (c,e) - where - b = MPBall (mpFloat $ dyadic c) (errorBound $ mpFloat $ dyadic e) - -instance ConvertibleExactly Integer MPBall where - safeConvertExactly x - | isFinite b = Right b - | otherwise = convError "too large to convert to MPBall" x - where - b = MPBall (mpFloat x) (errorBound 0) - -instance ConvertibleExactly Int MPBall where - safeConvertExactly x = Right $ MPBall (mpFloat x) (errorBound 0) - -{--- constructing a ball with a given precision ---} - -instance ConvertibleWithPrecision Integer MPBall where - safeConvertP p x - | isFinite b = Right b - | otherwise = convError ("too large to convert to MPBall with precision " ++ show p) x - where - b = MPBall xUp (xUp `EB.subMP` xDn) - xUp = MPFloat.fromIntegerUp p x - xDn = MPFloat.fromIntegerDown p x - -instance ConvertibleWithPrecision Int MPBall where - safeConvertP p = safeConvertP p . integer - -instance ConvertibleWithPrecision Dyadic MPBall where - safeConvertP p x - | isFinite b = Right b - | otherwise = convError ("too large to convert to MPBall with precision " ++ show p) x - where - b = mpBall x - -instance ConvertibleWithPrecision Rational MPBall where - safeConvertP p x - | isFinite b = Right b - | otherwise = convError ("too large to convert to MPBall with precision " ++ show p) x - where - b = MPBall xUp (xUp `EB.subMP` xDn) - xUp = MPFloat.fromRationalUp p x - xDn = MPFloat.fromRationalDown p x - -instance ConvertibleWithPrecision (Rational, Rational) MPBall where - safeConvertP p (x,e) - | isFinite b = Right b - | otherwise = convError ("too large to convert to MPBall with precision " ++ show p) x - where - b = MPBall xFlt (xe + eUp) -- beware, precision may be too high relative to accuracy - (MPBall xFlt xe) = mpBallP p x - eUp = errorBound e diff --git a/src/AERN2/MP/UseMPFR/Ball/Elementary.hs b/src/AERN2/MP/UseMPFR/Ball/Elementary.hs deleted file mode 100644 --- a/src/AERN2/MP/UseMPFR/Ball/Elementary.hs +++ /dev/null @@ -1,147 +0,0 @@ -{-# LANGUAGE CPP #-} -{-| - Module : AERN2.MP.UseMPFR.Ball.Elementary - Description : Elementary operations on arbitrary precision dyadic balls - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mikkonecny@gmail.com - Stability : experimental - Portability : portable - - Elementary operations on arbitrary precision dyadic balls --} -module AERN2.MP.UseMPFR.Ball.Elementary -( - -- * Ball operations (see also instances) - piBallP - -- * Helpers for constructing ball functions - , fromApproxWithLipschitz -) -where - -import MixedTypesNumPrelude -import qualified Prelude as P - -import AERN2.Normalize - -import AERN2.MP.Dyadic (Dyadic) -import qualified AERN2.MP.UseMPFR.Float as MPFloat -import AERN2.MP.UseMPFR.Float (MPFloat, mpFloat) --- import AERN2.MP.UseMPFR.Float.Operators -import AERN2.MP.Precision -import qualified AERN2.MP.UseMPFR.ErrorBound as EB -import AERN2.MP.UseMPFR.ErrorBound (errorBound) - -import AERN2.MP.UseMPFR.Ball.Type -import AERN2.MP.UseMPFR.Ball.Conversions () -import AERN2.MP.UseMPFR.Ball.Comparisons () -import AERN2.MP.UseMPFR.Ball.Field () - - -{- trigonometrics -} - -piBallP :: Precision -> MPBall -piBallP p = MPBall piUp (piUp `EB.subMP` piDown) - where - piUp = MPFloat.piUp p - piDown = MPFloat.piDown p - -instance CanSinCos MPBall where - sin = sinB 1 - cos = cosB 1 - -sinB :: Integer -> MPBall -> MPBall -sinB i x = - -- increasingPrecisionUntilNotImproving (fromApproxWithLipschitz MPFloat.sinDown MPFloat.sinUp lip) x - fromApproxWithLipschitz MPFloat.sinDown MPFloat.sinUp lip x - where - lip - | i == 0 = mpFloat 1 - | otherwise = snd $ endpointsMP $ abs $ cosB (i - 1) x - -cosB :: Integer -> MPBall -> MPBall -cosB i x = - -- increasingPrecisionUntilNotImproving (fromApproxWithLipschitz MPFloat.cosDown MPFloat.cosUp lip) x - fromApproxWithLipschitz MPFloat.cosDown MPFloat.cosUp lip x - where - lip - | i == 0 = mpFloat 1 - | otherwise = snd $ endpointsMP $ abs $ sinB (i - 1) x - --- increasingPrecisionUntilNotImproving :: (MPBall -> MPBall) -> (MPBall -> MPBall) --- increasingPrecisionUntilNotImproving f x = --- waitUntilNotImproving $ map aux (precisions xPrec (xPrec*2)) --- where --- xPrec = getPrecision x --- precisions p1 p2 = p1 : (precisions p2 (p1 + p2)) --- aux p = f $ setPrecision p x --- waitUntilNotImproving xx@(x1:_) = aux2 (getAccuracy x1) xx --- waitUntilNotImproving _ = error "AERN2.MP.UseMPFR.Ball.Elementary: internal error in increasingPrecisionUntilNotImproving" --- aux2 x1AC (x1:x2:rest) --- | x1AC < x2AC = aux2 x2AC (x2:rest) --- | otherwise = x1 --- where --- x2AC = getAccuracy x2 --- aux2 _ _ = error "AERN2.MP.UseMPFR.Ball.Elementary: internal error in increasingPrecisionUntilNotImproving" - -{- exp, log, power -} - -instance CanExp MPBall where - exp = intervalFunctionByEndpointsUpDown MPFloat.expDown MPFloat.expUp - -instance CanLog MPBall where - type LogType MPBall = CN MPBall - log x - | x !>! 0 = - cn $ intervalFunctionByEndpointsUpDown MPFloat.logDown MPFloat.logUp x - | x !<=! 0 = noValueNumErrorCertainCN err - | otherwise = noValueNumErrorPotentialCN err - where - err = OutOfRange $ "log: argument must be > 0: " ++ show x - -instance CanPow MPBall MPBall where - powNoCN b e = (~!) $ pow b e - pow = powUsingExpLog (mpBall 0) (mpBall 1) - -instance CanPow MPBall Dyadic where - powNoCN b e = (~!) $ pow b e - pow b e = powUsingExpLog (mpBall 0) (mpBall 1) b (mpBall e) - -instance CanPow MPBall Rational where - powNoCN b e = (~!) $ pow b e - pow b e = powUsingExpLog (mpBall 0) (mpBall 1) b (mpBallP (getPrecision b) e) - -instance CanSqrt MPBall where - type SqrtType MPBall = CN MPBall - sqrt x - | x !>=! 0 = cn $ aux x - | x ! MPFloat.sqrtDown (P.max (mpFloat 0) e)) - (\ e -> MPFloat.sqrtUp (P.max (mpFloat 0) e)) - err = OutOfRange $ "sqrt: argument must be >= 0: " ++ show x - -{- generic methods for computing real functions from MPFR-approximations -} - -{-| - Computes a real function @f@ from correctly rounded MPFR-approximations and a number @lip@ which is a - Lipschitz constant for @f@, i.e. @|f(x) - f(y)| <= lip * |x - y|@ for all @x@,@y@. --} -fromApproxWithLipschitz :: - (MPFloat -> MPFloat) {-^ @fDown@: a version of @f@ on MPFloat rounding *downwards* -} -> - (MPFloat -> MPFloat) {-^ @fUp@: a version of @f@ on MPFloat rounding *upwards* -} -> - MPFloat {-^ @lip@ a Lipschitz constant for @f@, @lip > 0@ -} -> - (MPBall -> MPBall) {-^ @f@ on MPBall rounding *outwards* -} -fromApproxWithLipschitz fDown fUp lip _x@(MPBall xc xe) = - normalize $ MPBall fxc err - where - fxl = fDown xc - fxu = fUp xc - (MPBall fxc fxe) = - setPrecision (getPrecision xc) $ -- beware, some MPFR functions increase precision, eg sine and cosine - fromEndpointsMP fxl fxu - err = (errorBound lip) * xe + fxe diff --git a/src/AERN2/MP/UseMPFR/Ball/Field.hs b/src/AERN2/MP/UseMPFR/Ball/Field.hs deleted file mode 100644 --- a/src/AERN2/MP/UseMPFR/Ball/Field.hs +++ /dev/null @@ -1,344 +0,0 @@ -{-# LANGUAGE TemplateHaskell #-} -{-| - Module : AERN2.MP.UseMPFR.Ball.Field - Description : Field operations on arbitrary precision dyadic balls - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mikkonecny@gmail.com - Stability : experimental - Portability : portable - - Field operations on arbitrary precision dyadic balls --} -module AERN2.MP.UseMPFR.Ball.Field -() -where - -import MixedTypesNumPrelude --- import qualified Prelude as P - -import Control.CollectErrors - -import AERN2.Normalize - -import AERN2.MP.Dyadic (Dyadic) -import AERN2.MP.UseMPFR.Float (mpFloat) -import AERN2.MP.UseMPFR.Float.Operators -import AERN2.MP.Precision -import qualified AERN2.MP.UseMPFR.ErrorBound as EB - -import AERN2.MP.UseMPFR.Ball.Type -import AERN2.MP.UseMPFR.Ball.Conversions () -import AERN2.MP.UseMPFR.Ball.Comparisons () - -{- addition -} - -instance CanAddAsymmetric MPBall MPBall where - type AddType MPBall MPBall = MPBall - add (MPBall x1 e1) (MPBall x2 e2) = - normalize $ MPBall sumUp ((sumUp `EB.subMP` sumDn) + e1 + e2) - where - sumUp = x1 +^ x2 - sumDn = x1 +. x2 - -instance CanAddAsymmetric MPBall Int where - type AddType MPBall Int = MPBall - add = convertSecond add -instance CanAddAsymmetric Int MPBall where - type AddType Int MPBall = MPBall - add = convertFirst add - -instance CanAddAsymmetric MPBall Integer where - type AddType MPBall Integer = MPBall - add = convertSecond add -instance CanAddAsymmetric Integer MPBall where - type AddType Integer MPBall = MPBall - add = convertFirst add - -instance CanAddAsymmetric MPBall Dyadic where - type AddType MPBall Dyadic = MPBall - add = convertSecond add -instance CanAddAsymmetric Dyadic MPBall where - type AddType Dyadic MPBall = MPBall - add = convertFirst add - -instance CanAddAsymmetric MPBall Rational where - type AddType MPBall Rational = MPBall - add = convertPSecond add -instance CanAddAsymmetric Rational MPBall where - type AddType Rational MPBall = MPBall - add = convertPFirst add - -instance - (CanAddAsymmetric MPBall b - , CanEnsureCE es b - , CanEnsureCE es (AddType MPBall b) - , SuitableForCE es) - => - CanAddAsymmetric MPBall (CollectErrors es b) - where - type AddType MPBall (CollectErrors es b) = - EnsureCE es (AddType MPBall b) - add = lift2TLCE add - -instance - (CanAddAsymmetric a MPBall - , CanEnsureCE es a - , CanEnsureCE es (AddType a MPBall) - , SuitableForCE es) - => - CanAddAsymmetric (CollectErrors es a) MPBall - where - type AddType (CollectErrors es a) MPBall = - EnsureCE es (AddType a MPBall) - add = lift2TCE add - -{- subtraction -} - -instance CanSub MPBall MPBall - -instance CanSub MPBall Integer -instance CanSub Integer MPBall - -instance CanSub MPBall Int -instance CanSub Int MPBall - -instance CanSub MPBall Rational -instance CanSub Rational MPBall - -instance CanSub MPBall Dyadic -instance CanSub Dyadic MPBall - -instance - (CanSub MPBall b - , CanEnsureCE es b - , CanEnsureCE es (SubType MPBall b) - , SuitableForCE es) - => - CanSub MPBall (CollectErrors es b) - where - type SubType MPBall (CollectErrors es b) = - EnsureCE es (SubType MPBall b) - sub = lift2TLCE sub - -instance - (CanSub a MPBall - , CanEnsureCE es a - , CanEnsureCE es (SubType a MPBall) - , SuitableForCE es) - => - CanSub (CollectErrors es a) MPBall - where - type SubType (CollectErrors es a) MPBall = - EnsureCE es (SubType a MPBall) - sub = lift2TCE sub - -{- multiplication -} - -instance CanMulAsymmetric MPBall MPBall where - mul (MPBall x1 e1) (MPBall x2 e2) = - normalize $ MPBall x12Up (e12 + e1*(abs x2) + e2*(abs x1) + e1*e2) - -- the mixed operations above automatically convert - -- MPFloat to ErrorBound, checking non-negativity - where - x12Up = x1 *^ x2 - x12Down = x1 *. x2 - e12 = x12Up -^ x12Down - -instance CanMulAsymmetric MPBall Int where - type MulType MPBall Int = MPBall - mul = convertSecond mul -instance CanMulAsymmetric Int MPBall where - type MulType Int MPBall = MPBall - mul = convertFirst mul - -instance CanMulAsymmetric MPBall Integer where - type MulType MPBall Integer = MPBall - mul = convertSecond mul -instance CanMulAsymmetric Integer MPBall where - type MulType Integer MPBall = MPBall - mul = convertFirst mul - -instance CanMulAsymmetric MPBall Dyadic where - type MulType MPBall Dyadic = MPBall - mul = convertSecond mul -instance CanMulAsymmetric Dyadic MPBall where - type MulType Dyadic MPBall = MPBall - mul = convertFirst mul - -instance CanMulAsymmetric MPBall Rational where - type MulType MPBall Rational = MPBall - mul = convertPSecond mul -instance CanMulAsymmetric Rational MPBall where - type MulType Rational MPBall = MPBall - mul = convertPFirst mul - -instance - (CanMulAsymmetric MPBall b - , CanEnsureCE es b - , CanEnsureCE es (MulType MPBall b) - , SuitableForCE es) - => - CanMulAsymmetric MPBall (CollectErrors es b) - where - type MulType MPBall (CollectErrors es b) = - EnsureCE es (MulType MPBall b) - mul = lift2TLCE mul - -instance - (CanMulAsymmetric a MPBall - , CanEnsureCE es a - , CanEnsureCE es (MulType a MPBall) - , SuitableForCE es) - => - CanMulAsymmetric (CollectErrors es a) MPBall - where - type MulType (CollectErrors es a) MPBall = - EnsureCE es (MulType a MPBall) - mul = lift2TCE mul - - -{- division -} - -instance CanDiv MPBall MPBall where - type DivTypeNoCN MPBall MPBall = MPBall - divideNoCN b1 b2 = (~!) (divide b1 b2) - type DivType MPBall MPBall = CN MPBall - divide (MPBall x1 e1) b2@(MPBall x2 e2) - | isCertainlyNonZero b2 = - cn $ normalize $ MPBall x12Up err - | isCertainlyZero b2 = - noValueNumErrorCertainCN DivByZero - | otherwise = - noValueNumErrorPotentialCN DivByZero - where - x12Up = x1 /^ x2 - x12Down = x1 /. x2 - x12AbsUp = (abs x12Up) `max` (abs x12Down) - e12 = x12Up -^ x12Down - err = - ((e12 *^ (abs x2)) -- e12 * |x2| - + - e1 - + - (e2 * x12AbsUp) -- e2 * |x| - ) - * - ((mpFloat 1) /^ ((abs x2) -. (mpFloat e2))) - -- 1/(|x2| - e2) rounded upwards -{- -A derivation of the above formula for an upper bound on the error: - - * e = - * = max ( (x1 ± e1) / (x2 ± e2) - x ) - * = max ( ( x1 ± e1 - (x*(x2 ± e2) ) / (x2 ± e2) ) - * ≤ max ( ( x1 ± e1 - ((x1/x2) ± e12)x2 ± x*e2 ) / (x2 ± e2) ) - * = max ( ( x1 ± e1 - x1 ± e12*x2 ± x*e2 ) / (x2 ± e2) ) - * = max ( ( ± e1 ± e12*x2 ± x*e2 ) / (x2 ± e2) ) - * ≤ (e1 + e12*|x2| + |x|*e2 ) / (|x2| - e2) - * ≤ (e1 +^ e12*^|x2| +^ |x|*^e2 ) /^ (|x2| -. e2) --} - -$(declForTypes - [[t| Integer |], [t| Int |], [t| Dyadic |]] - (\ t -> [d| - instance CanDiv MPBall $t where - type DivType MPBall $t = CN MPBall - divide = convertSecond divide - type DivTypeNoCN MPBall $t = MPBall - divideNoCN = convertSecond divideNoCN - instance CanDiv $t MPBall where - type DivType $t MPBall = CN MPBall - divide = convertFirst divide - type DivTypeNoCN $t MPBall = MPBall - divideNoCN = convertFirst divideNoCN - |])) - -instance CanDiv Dyadic Dyadic where - type DivTypeNoCN Dyadic Dyadic = MPBall - divideNoCN a b = divideNoCN (mpBall a) (mpBall b) - divide a b = divide (mpBall a) (mpBall b) - -instance CanDiv MPBall Rational where - type DivTypeNoCN MPBall Rational = MPBall - divideNoCN = convertPSecond divideNoCN - divide = convertPSecond divide -instance CanDiv Rational MPBall where - type DivTypeNoCN Rational MPBall = MPBall - divideNoCN = convertPFirst divideNoCN - divide = convertPFirst divide - -instance - (CanDiv MPBall b - , CanEnsureCE es b - , CanEnsureCE es (DivType MPBall b) - , CanEnsureCE es (DivTypeNoCN MPBall b) - , SuitableForCE es) - => - CanDiv MPBall (CollectErrors es b) - where - type DivType MPBall (CollectErrors es b) = - EnsureCE es (DivType MPBall b) - divide = lift2TLCE divide - type DivTypeNoCN MPBall (CollectErrors es b) = - EnsureCE es (DivTypeNoCN MPBall b) - divideNoCN = lift2TLCE divideNoCN - -instance - (CanDiv a MPBall - , CanEnsureCE es a - , CanEnsureCE es (DivType a MPBall) - , CanEnsureCE es (DivTypeNoCN a MPBall) - , SuitableForCE es) - => - CanDiv (CollectErrors es a) MPBall - where - type DivType (CollectErrors es a) MPBall = - EnsureCE es (DivType a MPBall) - divide = lift2TCE divide - type DivTypeNoCN (CollectErrors es a) MPBall = - EnsureCE es (DivTypeNoCN a MPBall) - divideNoCN = lift2TCE divideNoCN - -{- integer power -} - -instance CanPow MPBall Integer where - powNoCN b e = (~!) $ powUsingMulRecip (mpBall 1) b e - pow = powUsingMulRecip (mpBall 1) - -instance CanPow MPBall Int where - powNoCN b e = (~!) $ powUsingMulRecip (mpBall 1) b e - pow = powUsingMulRecip (mpBall 1) - -instance - (CanPow MPBall b - , CanEnsureCE es b - , CanEnsureCE es (PowType MPBall b) - , CanEnsureCE es (PowTypeNoCN MPBall b) - , SuitableForCE es) - => - CanPow MPBall (CollectErrors es b) - where - type PowTypeNoCN MPBall (CollectErrors es b) = - EnsureCE es (PowTypeNoCN MPBall b) - powNoCN = lift2TLCE powNoCN - type PowType MPBall (CollectErrors es b) = - EnsureCE es (PowType MPBall b) - pow = lift2TLCE pow - -instance - (CanPow a MPBall - , CanEnsureCE es a - , CanEnsureCE es (PowType a MPBall) - , CanEnsureCE es (PowTypeNoCN a MPBall) - , SuitableForCE es) - => - CanPow (CollectErrors es a) MPBall - where - type PowTypeNoCN (CollectErrors es a) MPBall = - EnsureCE es (PowTypeNoCN a MPBall) - powNoCN = lift2TCE powNoCN - type PowType (CollectErrors es a) MPBall = - EnsureCE es (PowType a MPBall) - pow = lift2TCE pow diff --git a/src/AERN2/MP/UseMPFR/Ball/PreludeOps.hs b/src/AERN2/MP/UseMPFR/Ball/PreludeOps.hs deleted file mode 100644 --- a/src/AERN2/MP/UseMPFR/Ball/PreludeOps.hs +++ /dev/null @@ -1,76 +0,0 @@ -{-# LANGUAGE CPP #-} -{-| - Module : AERN2.MP.UseMPFR.Ball.PreludeOps - Description : Instances of Prelude.Num etc - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mikkonecny@gmail.com - Stability : experimental - Portability : portable - - Instances of Prelude classes Eq, Ord, Num etc --} -module AERN2.MP.UseMPFR.Ball.PreludeOps -( -) -where - -import MixedTypesNumPrelude -import qualified Prelude as P - -import AERN2.MP.Dyadic (dyadic) - -import AERN2.MP.UseMPFR.Ball.Type -import AERN2.MP.UseMPFR.Ball.Conversions () -import AERN2.MP.UseMPFR.Ball.Comparisons () -import AERN2.MP.UseMPFR.Ball.Field () -import AERN2.MP.UseMPFR.Ball.Elementary - -{- Instances of Prelude numerical classes provided for convenient use outside AERN2 - and also because Template Haskell translates (-x) to (Prelude.negate x) -} - -instance P.Eq MPBall where - a == b = - case a == b of - Just t -> t - _ -> - error "Failed to decide equality of MPBalls. If you switch to MixedTypesNumPrelude instead of Prelude, comparison of MPBalls returns Maybe Bool instead of Bool." - -instance P.Ord MPBall where - compare a b = - case (a < b, a == b, a > b) of - (Just True, _, _) -> P.LT - (_, Just True, _) -> P.EQ - (_, _, Just True) -> P.GT - _ -> - error "Failed to decide order of MPBalls. If you switch to MixedTypesNumPrelude instead of Prelude, comparison of MPBalls returns Maybe Bool instead of Bool." - -instance P.Num MPBall where - fromInteger = convertExactly - negate = negate - (+) = (+) - (*) = (*) - abs = abs - signum = error "Prelude.signum not implemented for MPBall" - -instance P.Fractional MPBall where - fromRational = convertExactly . dyadic -- will work only for dyadic rationals - recip = (~!) . recip - (/) = (/!) - -instance P.Floating MPBall where - pi = error "MPBall: no pi :: MPBall, use pi ? (bitsS n) instead" - sqrt = (~!) . sqrt - exp = exp - sin = sin - cos = cos - log = (~!) . log - atan = error "MPBall: atan not implemented yet" - atanh = error "MPBall: atanh not implemented yet" - asin = error "MPBall: asin not implemented yet" - acos = error "MPBall: acos not implemented yet" - sinh = error "MPBall: sinh not implemented yet" - cosh = error "MPBall: cosh not implemented yet" - asinh = error "MPBall: asinh not implemented yet" - acosh = error "MPBall: acosh not implemented yet" diff --git a/src/AERN2/MP/UseMPFR/Ball/Type.hs b/src/AERN2/MP/UseMPFR/Ball/Type.hs deleted file mode 100644 --- a/src/AERN2/MP/UseMPFR/Ball/Type.hs +++ /dev/null @@ -1,276 +0,0 @@ -{-# LANGUAGE TemplateHaskell #-} -{-| - Module : AERN2.MP.UseMPFR.Ball.Type - Description : Arbitrary precision dyadic balls - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mikkonecny@gmail.com - Stability : experimental - Portability : portable - - Arbitrary precision dyadic balls --} -module AERN2.MP.UseMPFR.Ball.Type -( - -- * Auxiliary types - module AERN2.MP.Precision - , module AERN2.MP.Accuracy - , module AERN2.MP.Enclosure - -- * The Ball type - , MPBall(..), CanBeMPBall, mpBall, CanBeMPBallP, mpBallP - , reducePrecionIfInaccurate - -- * Ball construction/extraction functions - , endpointsMP, fromEndpointsMP -) -where - -import MixedTypesNumPrelude --- import qualified Prelude as P - -import Control.CollectErrors - -import GHC.Generics (Generic) - -import Text.Printf - -import AERN2.Normalize - -import AERN2.Norm - -import AERN2.MP.Dyadic -import qualified AERN2.MP.UseMPFR.Float as MPFloat -import AERN2.MP.UseMPFR.Float (MPFloat, mpFloat) -import AERN2.MP.UseMPFR.Float.Operators -import AERN2.MP.Precision -import AERN2.MP.Accuracy -import qualified AERN2.MP.UseMPFR.ErrorBound as EB -import AERN2.MP.UseMPFR.ErrorBound (ErrorBound, errorBound) -import AERN2.MP.Enclosure - -data MPBall = MPBall - { ball_value :: MPFloat - , ball_error :: ErrorBound - } - -- { ball_value :: {-# UNPACK #-} ! MPFloat - -- , ball_error :: {-# UNPACK #-} ! ErrorBound - -- } - deriving (Generic) - -instance Show MPBall - where - show b@(MPBall x _e) = - -- printf "[%s ± %s](prec=%s)" (show x) (showAC $ getAccuracy b) (show $ integer $ getPrecision b) - printf "[%s ± %s]" (show x) (showAC $ getAccuracy b) - -- "[" ++ show x ++ " ± " ++ show e ++ "](prec=" ++ (show $ integer $ getPrecision x) ++ ")" - where - showAC Exact = "0" - showAC NoInformation = "oo" - showAC ac = "<2^(" ++ show (negate $ fromAccuracy ac) ++ ")" - - -instance (SuitableForCE es) => CanEnsureCE es MPBall where - --- instance CanTestValid MPBall where --- isValid = isFinite - -instance CanTestNaN MPBall where - isNaN = not . isFinite -instance CanTestFinite MPBall where - isInfinite = const False - isFinite (MPBall x e) = isFinite x && isFinite (mpFloat e) - -instance CanNormalize MPBall where - normalize b - | isFinite b = - b - -- reducePrecionIfInaccurate b - | otherwise = error $ "invalid MPBall: " ++ show b - -{-| - Reduce the precision of the ball centre if the - accuracy of the ball is poor. - - More precisely, reduce the precision of the centre - so that the ulp is approximately (radius / 1024), - unless the ulp is already lower than this. --} -reducePrecionIfInaccurate :: MPBall -> MPBall -reducePrecionIfInaccurate b@(MPBall x _) = - case (bAcc, bNorm) of - (Exact, _) -> b - (_, NormZero) -> b - _ | p_e_nb < p_x -> setPrecision p_e_nb b - _ -> b - where - bAcc = getAccuracy b - bNorm = getNormLog b - p_x = getPrecision x - p_e_nb = prec $ max 2 (10 + nb + fromAccuracy bAcc) - (NormBits nb) = bNorm - -instance CanTestContains MPBall MPBall where - contains (MPBall xLarge eLarge) (MPBall xSmall eSmall) = - xLargeDy - eLargeDy <= xSmallDy - eSmallDy - && - xSmallDy + eSmallDy <= xLargeDy + eLargeDy - where - xLargeDy = dyadic xLarge - eLargeDy = dyadic eLarge - xSmallDy = dyadic xSmall - eSmallDy = dyadic eSmall - -$(declForTypes - [[t| Integer |], [t| Int |], [t| Rational |], [t| Dyadic |]] - (\ t -> [d| - instance CanTestContains MPBall $t where - contains (MPBall c e) x = - l <= x && x <= r - where - l = cDy - eDy - r = cDy + eDy - cDy = dyadic c - eDy = dyadic e - |])) - -{- ball construction/extraction functions -} - -instance IsInterval MPBall MPFloat where - fromEndpoints l u - | u < l = fromEndpoints u l - | otherwise = - MPBall (mpFloat cDy) (errorBound $ mpFloat eDy) - where - lDy = dyadic l - uDy = dyadic u - cDy = (lDy + uDy) * (dyadic 0.5) - eDy = (uDy - cDy) `max` (cDy - lDy) - endpoints (MPBall x e) = (mpFloat lDy, mpFloat uDy) - where - xDy = dyadic x - eDy = dyadic e - lDy = xDy - eDy - uDy = xDy + eDy - -fromEndpointsMP :: MPFloat -> MPFloat -> MPBall -fromEndpointsMP = fromEndpoints - -endpointsMP :: MPBall -> (MPFloat, MPFloat) -endpointsMP = endpoints - -instance IsInterval MPBall MPBall where - fromEndpoints l r = -- works as union even when r < l - fromEndpointsMP lMP uMP - where - lMP = min llMP rlMP - uMP = max luMP ruMP - (llMP, luMP) = endpointsMP l - (rlMP, ruMP) = endpointsMP r - endpoints x = (l,u) - where - l = MPBall lMP (errorBound 0) - u = MPBall uMP (errorBound 0) - (lMP, uMP) = endpointsMP x - -instance IsBall MPBall where - type CentreType MPBall = Dyadic - centre (MPBall cMP _e) = dyadic cMP - centreAsBallAndRadius x = (cB,e) - where - (MPBall cMP e) = x - cB = MPBall cMP (errorBound 0) - radius (MPBall _ e) = e - updateRadius updateFn (MPBall c e) = MPBall c (updateFn e) - -{--- constructing a ball with a given precision ---} - -type CanBeMPBallP t = (ConvertibleWithPrecision t MPBall) - -mpBallP :: (CanBeMPBallP t) => Precision -> t -> MPBall -mpBallP = convertP - - -{--- constructing an exact ball ---} - -type CanBeMPBall t = ConvertibleExactly t MPBall - -mpBall :: (CanBeMPBall t) => t -> MPBall -mpBall = convertExactly - -{-- extracting approximate information about a ball --} - -instance HasAccuracy MPBall where - getAccuracy = getAccuracy . ball_error - -instance CanReduceSizeUsingAccuracyGuide MPBall where - reduceSizeUsingAccuracyGuide acGuide b@(MPBall x _e) = - case acGuide of - NoInformation -> lowerPrecisionIfAbove (prec 2) b - _ | getAccuracy b > acGuide -> tryPrec newPrec - _ -> b - where - tryPrec p - | getAccuracy bP >= acGuide = bP - | otherwise = tryPrec (p + 10) - where - bP = lowerPrecisionIfAbove p b - queryBits = 1 + fromAccuracy acGuide - newPrec = - case (getNormLog x) of - NormBits xNormBits -> - prec (max 2 (queryBits + xNormBits + 2)) - NormZero -> - prec $ max 2 queryBits - -- bWithLowAC = - -- case acGuide of - -- Exact -> b - -- NoInformation -> b - -- _ -> normalize $ - -- MPBall x (errorBound ((0.5^(fromAccuracy acGuide))⚡)) - -instance HasNorm MPBall where - getNormLog ball = getNormLog boundMP - where - (_, MPBall boundMP _) = endpoints $ absRaw ball - -instance HasApproximate MPBall where - type Approximate MPBall = (MPFloat, Bool) - getApproximate ac b@(MPBall x e) = - (approx, isAccurate) - where - isAccurate = getAccuracy b < ac - approx - | closeToN = n - | otherwise = MPFloat.setPrecisionUp (prec (fromAccuracy ac)) x - where - n = mpFloat $ round $ rational x - closeToN = ((abs $ x -^ n) <= e) - -instance HasPrecision MPBall where - getPrecision = getPrecision . ball_value - -instance CanSetPrecision MPBall where - setPrecision p (MPBall x e) - | p >= pPrev = MPBall xUp e - | otherwise = MPBall xUp (e + (xUp `EB.subMP` xDown)) - where - pPrev = MPFloat.getPrecision x - xUp = MPFloat.setPrecisionUp p x - xDown = MPFloat.setPrecisionDown p x - -{- negation & abs -} - -instance CanNeg MPBall where - negate (MPBall x e) = MPBall (-x) e - -instance CanAbs MPBall where - abs = normalize . absRaw - -absRaw :: MPBall -> MPBall -absRaw b - | l < 0 && 0 < r = - fromEndpointsMP (mpFloat 0) (max (-l) r) - | 0 <= l = b - | otherwise = -b - where - (l,r) = endpointsMP b diff --git a/src/AERN2/MP/UseMPFR/Dyadic.hs b/src/AERN2/MP/UseMPFR/Dyadic.hs deleted file mode 100644 --- a/src/AERN2/MP/UseMPFR/Dyadic.hs +++ /dev/null @@ -1,673 +0,0 @@ -{-# LANGUAGE CPP #-} --- #define DEBUG -{-# LANGUAGE DeriveDataTypeable #-} -{-| - Module : AERN2.MP.UseMPFR.Dyadic - Description : Dyadics with exact ring operations - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mikkonecny@gmail.com - Stability : experimental - Portability : portable - - Arbitrary precision floating-point numbers with exact ring operations. - - Currently, we use hmpfr when compiling with ghc 7.10 and higher - and haskell-mpfr when compiling with ghc 7.8. --} - -module AERN2.MP.UseMPFR.Dyadic -( - -- * Dyadic numbers and their basic operations - Dyadic, HasDyadics - -- * Dyadic constructors - , CanBeDyadic, dyadic - -- * tests - , specDyadic, tDyadic -) -where - -#ifdef DEBUG -import Debug.Trace (trace) -#define maybeTrace trace -#define maybeTraceIO putStrLn -#else -#define maybeTrace (\ (_ :: String) t -> t) -#define maybeTraceIO (\ (_ :: String) -> return ()) -#endif - -import MixedTypesNumPrelude -import qualified Prelude as P - -import Control.CollectErrors - -import Text.Printf -import Text.Regex.TDFA - -import Data.Typeable -import Data.Convertible - -import Test.Hspec -import Test.QuickCheck --- import qualified Test.Hspec.SmallCheck as SC - -import Data.Ratio (denominator, numerator) - -import Math.NumberTheory.Logarithms (integerLog2) - -import AERN2.Norm -import AERN2.MP.Precision -import AERN2.MP.Accuracy -import AERN2.MP.UseMPFR.Float - -{-| Exact dyadic type based on MPFloat. -} -newtype Dyadic = Dyadic { dyadicMPFloat :: MPFloat } - deriving (P.Eq, P.Ord, CanRound, HasPrecision, HasNorm, Typeable) - -instance Ring Dyadic -instance Ring (CN Dyadic) - -instance OrderedRing Dyadic -instance OrderedRing (CN Dyadic) - -instance OrderedCertainlyRing Dyadic -instance OrderedCertainlyRing (CN Dyadic) - -instance HasAccuracy Dyadic where getAccuracy _ = Exact - -instance Show Dyadic where - show (Dyadic x) - | e == 0 = printf "dyadic (%d)" n - | e > 0 = printf "dyadic (%d*0.5^%d)" n e - | otherwise = error "in show Dyadic" - where - xR = rational x - NormBits e = getNormLog (denominator xR) - n = numerator xR - -instance Read Dyadic where - readsPrec _pr dyadicS = - tryInt $ tryWithExp [] - where - tryInt tryNext = - case groups of - [nS] -> - case reads nS of - [(n,"")] -> [(dyadic (n :: Integer), afterS)] - _ -> tryNext - _ -> tryNext - where - (_,_,afterS,groups) = - dyadicS =~ "\\`dyadic \\(([-0-9]*)\\)" - :: (String, String, String, [String]) - tryWithExp tryNext = - case groups of - [nS,eS] -> - case (reads nS, reads eS) of - ([(n,"")],[(e,"")]) -> - [((n :: Integer)*(dyadic 0.5)^!(e :: Integer), afterS)] - _ -> tryNext - _ -> tryNext - where - (_,_,afterS,groups) = - dyadicS =~ "\\`dyadic \\(([-0-9]*)\\*0.5\\^([0-9]*)\\)" - :: (String, String, String, [String]) - -instance (SuitableForCE es) => CanEnsureCE es Dyadic - - -{-- conversions --} - -type HasDyadics t = ConvertibleExactly Dyadic t - -instance ConvertibleExactly Dyadic Dyadic where - safeConvertExactly = Right - -instance ConvertibleExactly Dyadic MPFloat where - safeConvertExactly = Right . dyadicMPFloat - -instance ConvertibleExactly Dyadic Rational where - safeConvertExactly = safeConvertExactly . dyadicMPFloat - -type CanBeDyadic t = ConvertibleExactly t Dyadic -dyadic :: (CanBeDyadic t) => t -> Dyadic -dyadic = convertExactly - -instance ConvertibleExactly MPFloat Dyadic where - safeConvertExactly = Right . Dyadic - -instance HasIntegerBounds Dyadic where - integerBounds d = (floor d, ceiling d) - -instance ConvertibleExactly Integer Dyadic where - safeConvertExactly = fmap Dyadic . safeConvertExactly - -instance ConvertibleExactly Int Dyadic where - safeConvertExactly = fmap Dyadic . safeConvertExactly - -instance ConvertibleExactly Rational Dyadic where - safeConvertExactly q - | isDyadic = Right $ Dyadic (fromRationalUp (prec $ max 2 (dp + np + 1)) q) - | otherwise = convError "this number is not dyadic" q - where - isDyadic = d == 2^!dp - dp = integerLog2 d - d = denominator q - np = integerLog2 (max 1 $ abs $ numerator q) - -instance Convertible Dyadic Double where - safeConvert = safeConvert . dyadicMPFloat - -instance (ConvertibleExactly Dyadic t, Monoid es) => ConvertibleExactly Dyadic (CollectErrors es t) where - safeConvertExactly = fmap (\v -> CollectErrors (Just v) mempty) . safeConvertExactly - -{-- comparisons --} - -instance HasEqAsymmetric Dyadic Dyadic -instance HasEqAsymmetric Dyadic Integer where - equalTo = convertSecond equalTo -instance HasEqAsymmetric Integer Dyadic where - equalTo = convertFirst equalTo -instance HasEqAsymmetric Dyadic Int where - equalTo = convertSecond equalTo -instance HasEqAsymmetric Int Dyadic where - equalTo = convertFirst equalTo -instance HasEqAsymmetric Dyadic Rational where - equalTo = convertFirst equalTo -instance HasEqAsymmetric Rational Dyadic where - equalTo = convertSecond equalTo - -instance - (HasEqAsymmetric Dyadic b - , CanEnsureCE es b - , CanEnsureCE es (EqCompareType Dyadic b) - , IsBool (EnsureCE es (EqCompareType Dyadic b)) - , SuitableForCE es) - => - HasEqAsymmetric Dyadic (CollectErrors es b) - where - type EqCompareType Dyadic (CollectErrors es b) = - EnsureCE es (EqCompareType Dyadic b) - equalTo = lift2TLCE equalTo - -instance - (HasEqAsymmetric a Dyadic - , CanEnsureCE es a - , CanEnsureCE es (EqCompareType a Dyadic) - , IsBool (EnsureCE es (EqCompareType a Dyadic)) - , SuitableForCE es) - => - HasEqAsymmetric (CollectErrors es a) Dyadic - where - type EqCompareType (CollectErrors es a) Dyadic = - EnsureCE es (EqCompareType a Dyadic) - equalTo = lift2TCE equalTo - -instance CanTestZero Dyadic - -instance HasOrderAsymmetric Dyadic Dyadic -instance HasOrderAsymmetric Dyadic Integer where - lessThan = convertSecond lessThan - leq = convertSecond leq -instance HasOrderAsymmetric Integer Dyadic where - lessThan = convertFirst lessThan - leq = convertFirst leq -instance HasOrderAsymmetric Dyadic Int where - lessThan = convertSecond lessThan - leq = convertSecond leq -instance HasOrderAsymmetric Int Dyadic where - lessThan = convertFirst lessThan - leq = convertFirst leq -instance HasOrderAsymmetric Rational Dyadic where - lessThan = convertSecond lessThan - leq = convertSecond leq -instance HasOrderAsymmetric Dyadic Rational where - lessThan = convertFirst lessThan - leq = convertFirst leq - -instance - (HasOrderAsymmetric Dyadic b - , CanEnsureCE es b - , CanEnsureCE es (OrderCompareType Dyadic b) - , IsBool (EnsureCE es (OrderCompareType Dyadic b)) - , SuitableForCE es) - => - HasOrderAsymmetric Dyadic (CollectErrors es b) - where - type OrderCompareType Dyadic (CollectErrors es b) = - EnsureCE es (OrderCompareType Dyadic b) - lessThan = lift2TLCE lessThan - leq = lift2TLCE leq - greaterThan = lift2TLCE greaterThan - geq = lift2TLCE geq - -instance - (HasOrderAsymmetric a Dyadic - , CanEnsureCE es a - , CanEnsureCE es (OrderCompareType a Dyadic) - , IsBool (EnsureCE es (OrderCompareType a Dyadic)) - , SuitableForCE es) - => - HasOrderAsymmetric (CollectErrors es a) Dyadic - where - type OrderCompareType (CollectErrors es a) Dyadic = - EnsureCE es (OrderCompareType a Dyadic) - lessThan = lift2TCE lessThan - leq = lift2TCE leq - greaterThan = lift2TCE greaterThan - geq = lift2TCE geq - - -instance CanTestPosNeg Dyadic - -instance CanTestInteger Dyadic where - certainlyNotInteger = certainlyNotInteger . rational - certainlyIntegerGetIt = certainlyIntegerGetIt . rational - -{- unary functions -} - -instance CanNeg Dyadic where - negate = lift1 negate - -instance CanAbs Dyadic where - abs = lift1 abs - -lift1 :: (MPFloat -> MPFloat) -> (Dyadic -> Dyadic) -lift1 op (Dyadic x) = Dyadic (op x) - -{- min/max -} - -instance CanMinMaxAsymmetric Dyadic Dyadic -instance CanMinMaxAsymmetric Integer Dyadic where - type MinMaxType Integer Dyadic = Dyadic - min = convertFirst min - max = convertFirst max -instance CanMinMaxAsymmetric Dyadic Integer where - type MinMaxType Dyadic Integer = Dyadic - min = convertSecond min - max = convertSecond max -instance CanMinMaxAsymmetric Int Dyadic where - type MinMaxType Int Dyadic = Dyadic - min = convertFirst min - max = convertFirst max -instance CanMinMaxAsymmetric Dyadic Int where - type MinMaxType Dyadic Int = Dyadic - min = convertSecond min - max = convertSecond max -instance CanMinMaxAsymmetric Rational Dyadic where - type MinMaxType Rational Dyadic = Rational - min = convertSecond min - max = convertSecond max -instance CanMinMaxAsymmetric Dyadic Rational where - type MinMaxType Dyadic Rational = Rational - min = convertFirst min - max = convertFirst max - -instance - (CanMinMaxAsymmetric Dyadic b - , CanEnsureCE es b - , CanEnsureCE es (MinMaxType Dyadic b) - , SuitableForCE es) - => - CanMinMaxAsymmetric Dyadic (CollectErrors es b) - where - type MinMaxType Dyadic (CollectErrors es b) = - EnsureCE es (MinMaxType Dyadic b) - min = lift2TLCE min - max = lift2TLCE max - -instance - (CanMinMaxAsymmetric a Dyadic - , CanEnsureCE es a - , CanEnsureCE es (MinMaxType a Dyadic) - , SuitableForCE es) - => - CanMinMaxAsymmetric (CollectErrors es a) Dyadic - where - type MinMaxType (CollectErrors es a) Dyadic = - EnsureCE es (MinMaxType a Dyadic) - min = lift2TCE min - max = lift2TCE max - -{- addition -} - -instance CanAddAsymmetric Dyadic Dyadic where - add = lift2 addDown addUp - -instance CanAddAsymmetric Integer Dyadic where - type AddType Integer Dyadic = Dyadic - add = convertFirst add -instance CanAddAsymmetric Dyadic Integer where - type AddType Dyadic Integer = Dyadic - add = convertSecond add - -instance CanAddAsymmetric Int Dyadic where - type AddType Int Dyadic = Dyadic - add = convertFirst add -instance CanAddAsymmetric Dyadic Int where - type AddType Dyadic Int = Dyadic - add = convertSecond add - -instance CanAddAsymmetric Rational Dyadic where - type AddType Rational Dyadic = Rational - add = convertSecond add -instance CanAddAsymmetric Dyadic Rational where - type AddType Dyadic Rational = Rational - add = convertFirst add - -instance - (CanAddAsymmetric Dyadic b - , CanEnsureCE es b - , CanEnsureCE es (AddType Dyadic b) - , SuitableForCE es) - => - CanAddAsymmetric Dyadic (CollectErrors es b) - where - type AddType Dyadic (CollectErrors es b) = - EnsureCE es (AddType Dyadic b) - add = lift2TLCE add - -instance - (CanAddAsymmetric a Dyadic - , CanEnsureCE es a - , CanEnsureCE es (AddType a Dyadic) - , SuitableForCE es) - => - CanAddAsymmetric (CollectErrors es a) Dyadic - where - type AddType (CollectErrors es a) Dyadic = - EnsureCE es (AddType a Dyadic) - add = lift2TCE add - -{- subtraction -} - -instance CanSub Dyadic Dyadic where - sub = lift2 subDown subUp - -instance CanSub Integer Dyadic where - type SubType Integer Dyadic = Dyadic - sub = convertFirst sub -instance CanSub Dyadic Integer where - type SubType Dyadic Integer = Dyadic - sub = convertSecond sub - -instance CanSub Int Dyadic where - type SubType Int Dyadic = Dyadic - sub = convertFirst sub -instance CanSub Dyadic Int where - type SubType Dyadic Int = Dyadic - sub = convertSecond sub - -instance CanSub Rational Dyadic where - type SubType Rational Dyadic = Rational - sub = convertSecond sub -instance CanSub Dyadic Rational where - type SubType Dyadic Rational = Rational - sub = convertFirst sub - -instance - (CanSub Dyadic b - , CanEnsureCE es b - , CanEnsureCE es (SubType Dyadic b) - , SuitableForCE es) - => - CanSub Dyadic (CollectErrors es b) - where - type SubType Dyadic (CollectErrors es b) = - EnsureCE es (SubType Dyadic b) - sub = lift2TLCE sub - -instance - (CanSub a Dyadic - , CanEnsureCE es a - , CanEnsureCE es (SubType a Dyadic) - , SuitableForCE es) - => - CanSub (CollectErrors es a) Dyadic - where - type SubType (CollectErrors es a) Dyadic = - EnsureCE es (SubType a Dyadic) - sub = lift2TCE sub - - -{- multiplication -} - -instance CanMulAsymmetric Dyadic Dyadic where - mul = lift2 mulDown mulUp - -instance CanMulAsymmetric Integer Dyadic where - type MulType Integer Dyadic = Dyadic - mul = convertFirst mul -instance CanMulAsymmetric Dyadic Integer where - type MulType Dyadic Integer = Dyadic - mul = convertSecond mul - -instance CanMulAsymmetric Int Dyadic where - type MulType Int Dyadic = Dyadic - mul = convertFirst mul -instance CanMulAsymmetric Dyadic Int where - type MulType Dyadic Int = Dyadic - mul = convertSecond mul - -instance CanMulAsymmetric Rational Dyadic where - type MulType Rational Dyadic = Rational - mul = convertSecond mul -instance CanMulAsymmetric Dyadic Rational where - type MulType Dyadic Rational = Rational - mul = convertFirst mul - -instance - (CanMulAsymmetric Dyadic b - , CanEnsureCE es b - , CanEnsureCE es (MulType Dyadic b) - , SuitableForCE es) - => - CanMulAsymmetric Dyadic (CollectErrors es b) - where - type MulType Dyadic (CollectErrors es b) = - EnsureCE es (MulType Dyadic b) - mul = lift2TLCE mul - -instance - (CanMulAsymmetric a Dyadic - , CanEnsureCE es a - , CanEnsureCE es (MulType a Dyadic) - , SuitableForCE es) - => - CanMulAsymmetric (CollectErrors es a) Dyadic - where - type MulType (CollectErrors es a) Dyadic = - EnsureCE es (MulType a Dyadic) - mul = lift2TCE mul - -instance CanPow Dyadic Integer where - powNoCN = powUsingMul (dyadic 1) - pow = integerPowCN (powUsingMul (dyadic 1)) -instance CanPow Dyadic Int where - powNoCN = powUsingMul (dyadic 1) - pow = integerPowCN (powUsingMul (dyadic 1)) - -instance - (CanDiv a Dyadic - , CanEnsureCE es a - , CanEnsureCE es (DivType a Dyadic) - , CanEnsureCE es (DivTypeNoCN a Dyadic) - , SuitableForCE es) - => - CanDiv (CollectErrors es a) Dyadic - where - type DivType (CollectErrors es a) Dyadic = - EnsureCE es (DivType a Dyadic) - divide = lift2TCE divide - type DivTypeNoCN (CollectErrors es a) Dyadic = - EnsureCE es (DivTypeNoCN a Dyadic) - divideNoCN = lift2TCE divideNoCN - -instance CanDiv Integer Dyadic where - type DivTypeNoCN Integer Dyadic = Rational - divideNoCN a b = divideNoCN a (rational b) -instance CanDiv Dyadic Integer where - type DivTypeNoCN Dyadic Integer = Rational - divideNoCN a b = divideNoCN (rational a) b - -instance CanDiv Int Dyadic where - type DivTypeNoCN Int Dyadic = Rational - divideNoCN a b = divideNoCN a (rational b) -instance CanDiv Dyadic Int where - type DivTypeNoCN Dyadic Int = Rational - divideNoCN a b = divideNoCN (rational a) b - -instance CanDiv Rational Dyadic where - type DivTypeNoCN Rational Dyadic = Rational - divideNoCN = convertSecond divideNoCN -instance CanDiv Dyadic Rational where - type DivTypeNoCN Dyadic Rational = Rational - divideNoCN = convertFirst divideNoCN - -instance - (CanDiv Dyadic b - , CanEnsureCE es b - , CanEnsureCE es (DivType Dyadic b) - , CanEnsureCE es (DivTypeNoCN Dyadic b) - , SuitableForCE es) - => - CanDiv Dyadic (CollectErrors es b) - where - type DivType Dyadic (CollectErrors es b) = - EnsureCE es (DivType Dyadic b) - divide = lift2TLCE divide - type DivTypeNoCN Dyadic (CollectErrors es b) = - EnsureCE es (DivTypeNoCN Dyadic b) - divideNoCN = lift2TLCE divideNoCN - -instance - (CanPow Dyadic b - , CanEnsureCE es b - , CanEnsureCE es (PowTypeNoCN Dyadic b) - , CanEnsureCE es (PowType Dyadic b) - , SuitableForCE es) - => - CanPow Dyadic (CollectErrors es b) - where - type PowTypeNoCN Dyadic (CollectErrors es b) = - EnsureCE es (PowTypeNoCN Dyadic b) - powNoCN = lift2TLCE powNoCN - type PowType Dyadic (CollectErrors es b) = - EnsureCE es (PowType Dyadic b) - pow = lift2TLCE pow - -instance - (CanPow a Dyadic - , CanEnsureCE es a - , CanEnsureCE es (PowType a Dyadic) - , CanEnsureCE es (PowTypeNoCN a Dyadic) - , SuitableForCE es) - => - CanPow (CollectErrors es a) Dyadic - where - type PowTypeNoCN (CollectErrors es a) Dyadic = - EnsureCE es (PowTypeNoCN a Dyadic) - powNoCN = lift2TCE powNoCN - type PowType (CollectErrors es a) Dyadic = - EnsureCE es (PowType a Dyadic) - pow = lift2TCE pow - -lift2 :: - (MPFloat -> MPFloat -> MPFloat) -> - (MPFloat -> MPFloat -> MPFloat) -> - (Dyadic -> Dyadic -> Dyadic) -lift2 opDown opUp (Dyadic x0) (Dyadic y0) = Dyadic (opExact x0 y0) - where - opExact x y - | rUp == rDown = rUp - | otherwise = - maybeTrace (printf "Dyadic.lift2: rUp = %s; rDown = %s; p = %s" (show rUp) (show rDown) (show $ integer p)) $ - opExact xH yH - where - rUp = opUp x y - rDown = opDown x y - xH = setPrecision pH x - yH = setPrecision pH y - pH = precisionTimes2 p - p = getPrecision rUp - -instance Arbitrary Dyadic where - arbitrary = - do - c <- finiteMPFloat - return (Dyadic c) - where - finiteMPFloat = - do - x <- arbitrary - if isFinite x - then return x - else finiteMPFloat - -{-| - A runtime representative of type @Dyadic@. - Used for specialising polymorphic tests to concrete types. --} -tDyadic :: T Dyadic -tDyadic = T "Dyadic" - -specDyadic :: Spec -specDyadic = - describe ("Dyadic") $ do - specConversion tInteger tDyadic dyadic round - specConversion tDyadic tRational rational dyadic - describe "order" $ do - specHasEqNotMixed tDyadic - specHasEq tInt tDyadic tRational - specCanTestZero tDyadic - specHasOrderNotMixed tDyadic - specHasOrder tInt tDyadic tRational - describe "min/max/abs" $ do - specCanNegNum tDyadic - specCanAbs tDyadic - specCanMinMaxNotMixed tDyadic - specCanMinMax tDyadic tInteger tDyadic - it "min Dyadic Rational (dyadic only)" $ do - property $ \ (x :: Dyadic) (y :: Dyadic) -> - x `min` y == x `min` (rational y) - it "max Dyadic Rational (dyadic only)" $ do - property $ \ (x :: Dyadic) (y :: Dyadic) -> - x `max` y == x `max` (rational y) - describe "ring" $ do - specCanAddNotMixed tDyadic - specCanAddSameType tDyadic - specCanAdd tInt tDyadic tInteger - specCanAdd tInteger tDyadic tInt - it "Dyadic + Rational (dyadic only)" $ do - property $ \ (x :: Dyadic) (y :: Dyadic) -> - x + y == x + (rational y) - specCanSubNotMixed tDyadic - specCanSub tDyadic tInteger - specCanSub tInteger tDyadic - specCanSub tDyadic tInt - specCanSub tInt tDyadic - it "Dyadic - Rational (dyadic only)" $ do - property $ \ (x :: Dyadic) (y :: Dyadic) -> - x - y == x - (rational y) - specCanMulNotMixed tDyadic - specCanMulSameType tDyadic - specCanMul tInt tDyadic tInteger - it "Dyadic * Rational (dyadic only)" $ do - property $ \ (x :: Dyadic) (y :: Dyadic) -> - x * y == x * (rational y) - specCanPow tDyadic tInteger - -instance P.Num Dyadic where - fromInteger = convertExactly - negate = negate - (+) = (+) - (*) = (*) - abs = abs - signum d - | d < 0 = dyadic (-1) - | d == 0 = dyadic 0 - | otherwise = dyadic 1 - -instance P.Real Dyadic where - toRational = convertExactly diff --git a/src/AERN2/MP/UseMPFR/ErrorBound.hs b/src/AERN2/MP/UseMPFR/ErrorBound.hs deleted file mode 100644 --- a/src/AERN2/MP/UseMPFR/ErrorBound.hs +++ /dev/null @@ -1,221 +0,0 @@ -{-# LANGUAGE DeriveDataTypeable #-} -{-| - Module : AERN2.MP.UseMPFR.ErrorBound - Description : Fixed precision non-negative up-rounded floating-point numbers - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mikkonecny@gmail.com - Stability : experimental - Portability : portable - - Fixed precision non-negative up-rounded floating-point numbers. - - Currently using a fixed-precision MPFloat. --} -module AERN2.MP.UseMPFR.ErrorBound - (ErrorBound, CanBeErrorBound, errorBound, - absMP, subMP) -where - -import MixedTypesNumPrelude -import qualified Prelude as P - -import Data.Typeable - -import Test.QuickCheck - -import Data.Convertible - -import Math.NumberTheory.Logarithms (integerLog2) - -import AERN2.MP.Precision -import AERN2.MP.Accuracy -import qualified AERN2.MP.UseMPFR.Float as MPFloat -import AERN2.MP.UseMPFR.Float (MPFloat, mpFloat, frequencyElements) -import AERN2.MP.UseMPFR.Float.Operators -import AERN2.MP.Dyadic - -{- example -} - -_example1 :: ErrorBound -_example1 = 2*((errorBound 0.01) + 0.1*(errorBound 0.01)/3) - -{- type -} - -{-| A non-negative Double value to serve as an error bound. Arithmetic is rounded towards +infinity. -} -newtype ErrorBound = ErrorBound { er2mp :: MPFloat } - deriving (P.Eq, P.Ord, Typeable) - -instance Show ErrorBound where - show (ErrorBound d) = show d - -errorBoundPrecision :: Precision -errorBoundPrecision = prec 53 - -instance HasAccuracy ErrorBound where - getAccuracy (ErrorBound e) - | eN > 0 = - bits $ negate $ integerLog2 eN - | e > 0 && eRecipN > 0 = - bits $ integerLog2 eRecipN - | e == 0 = Exact - | otherwise = NoInformation - where - eN = floor $ rational e - eRecipN = ceiling $ rational $ MPFloat.recipDown e - -{- conversions -} - -instance ConvertibleExactly ErrorBound MPFloat where - safeConvertExactly = Right . er2mp - -instance ConvertibleExactly ErrorBound Dyadic where - safeConvertExactly = Right . dyadic . er2mp - -instance ConvertibleExactly ErrorBound Rational where - safeConvertExactly = Right . convertExactly . mpFloat - -type CanBeErrorBound t = Convertible t ErrorBound -errorBound :: (CanBeErrorBound t) => t -> ErrorBound -errorBound = convert - -instance Convertible Rational ErrorBound where - safeConvert x - | x >= 0 = Right $ ErrorBound $ MPFloat.fromRationalUp errorBoundPrecision x - | otherwise = convError "Trying to construct a negative ErrorBound" x - -instance Convertible MPFloat ErrorBound where - safeConvert x - | x >= 0 = Right $ ErrorBound $ MPFloat.setPrecisionUp errorBoundPrecision x - | otherwise = convError "Trying to construct a negative ErrorBound" x - -instance Convertible Integer ErrorBound where - safeConvert x - | x >= 0 = Right $ ErrorBound $ MPFloat.fromIntegerUp errorBoundPrecision x - | otherwise = convError "Trying to construct a negative ErrorBound" x - -instance Convertible Int ErrorBound where - safeConvert = safeConvert . integer - -{- comparisons -} - -instance HasOrderAsymmetric ErrorBound ErrorBound - -instance HasOrderAsymmetric ErrorBound MPFloat where - lessThan = convertFirst lessThan - leq = convertFirst leq -instance HasOrderAsymmetric MPFloat ErrorBound where - lessThan = convertSecond lessThan - leq = convertSecond leq - -instance HasEqAsymmetric ErrorBound Rational where - equalTo = convertFirst equalTo -instance HasEqAsymmetric Rational ErrorBound where - equalTo = convertSecond equalTo -instance HasOrderAsymmetric ErrorBound Rational where - lessThan = convertFirst lessThan - leq = convertFirst leq -instance HasOrderAsymmetric Rational ErrorBound where - lessThan = convertSecond lessThan - leq = convertSecond leq - -instance HasEqAsymmetric ErrorBound Integer where - equalTo a b = equalTo (dyadic a) (dyadic b) -instance HasEqAsymmetric Integer ErrorBound where - equalTo a b = equalTo (dyadic a) (dyadic b) -instance HasOrderAsymmetric ErrorBound Integer where - lessThan a b = lessThan (dyadic a) (dyadic b) - leq a b = leq (dyadic a) (dyadic b) -instance HasOrderAsymmetric Integer ErrorBound where - lessThan a b = lessThan (dyadic a) (dyadic b) - leq a b = leq (dyadic a) (dyadic b) - -instance HasEqAsymmetric ErrorBound Int where - equalTo a b = equalTo (dyadic a) (dyadic b) -instance HasEqAsymmetric Int ErrorBound where - equalTo a b = equalTo (dyadic a) (dyadic b) -instance HasOrderAsymmetric ErrorBound Int where - lessThan a b = lessThan (dyadic a) (dyadic b) - leq a b = leq (dyadic a) (dyadic b) -instance HasOrderAsymmetric Int ErrorBound where - lessThan a b = lessThan (dyadic a) (dyadic b) - leq a b = leq (dyadic a) (dyadic b) - -instance CanMinMaxAsymmetric ErrorBound ErrorBound - -{- converting operations -} - -subMP :: MPFloat -> MPFloat -> ErrorBound -a `subMP` b = errorBound $ a -^ b - -absMP :: MPFloat -> ErrorBound -absMP = errorBound . abs - -{- up-rounded operations -} - -instance CanAddAsymmetric ErrorBound ErrorBound where - add (ErrorBound a) (ErrorBound b) = ErrorBound $ a +^ b - -instance CanAddAsymmetric ErrorBound MPFloat where - type AddType ErrorBound MPFloat = ErrorBound - add = convertSecondUsing (\ _ f -> convert f) add -instance CanAddAsymmetric MPFloat ErrorBound where - type AddType MPFloat ErrorBound = ErrorBound - add = convertFirstUsing (\ f _ -> convert f) add - -instance CanMulAsymmetric ErrorBound ErrorBound where - mul (ErrorBound a) (ErrorBound b) = ErrorBound $ a *^ b - -instance CanMulAsymmetric ErrorBound MPFloat where - type MulType ErrorBound MPFloat = ErrorBound - mul = convertSecondUsing (\ _ f -> convert f) mul -instance CanMulAsymmetric MPFloat ErrorBound where - type MulType MPFloat ErrorBound = ErrorBound - mul = convertFirstUsing (\ f _ -> convert f) mul - -instance CanMulAsymmetric ErrorBound Integer where - type MulType ErrorBound Integer = ErrorBound - mul (ErrorBound a) i - | i >= 0 = ErrorBound $ a *^ (MPFloat.fromIntegerUp errorBoundPrecision i) - | otherwise = error "trying to multiply ErrorBound by a negative integer" -instance CanMulAsymmetric Integer ErrorBound where - type MulType Integer ErrorBound = ErrorBound - mul i (ErrorBound b) - | i >= 0 = ErrorBound $ (MPFloat.fromIntegerUp errorBoundPrecision i) *^ b - | otherwise = error "trying to multiply ErrorBound by a negative integer" - -instance CanMulAsymmetric ErrorBound Rational where - type MulType ErrorBound Rational = ErrorBound - mul (ErrorBound a) r - | r >= 0.0 = ErrorBound $ a *^ (MPFloat.fromRationalUp errorBoundPrecision r) - | otherwise = error "trying to multiply ErrorBound by a negative integer" -instance CanMulAsymmetric Rational ErrorBound where - type MulType Rational ErrorBound = ErrorBound - mul r (ErrorBound b) - | r >= 0.0 = ErrorBound $ (MPFloat.fromRationalUp errorBoundPrecision r) *^ b - | otherwise = error "trying to multiply ErrorBound by a negative integer" - -instance CanDiv ErrorBound Integer where - type DivTypeNoCN ErrorBound Integer = ErrorBound - type DivType ErrorBound Integer = ErrorBound - divideNoCN = divide - divide (ErrorBound a) i - | i > 0 = ErrorBound $ a /^ (MPFloat.fromIntegerUp errorBoundPrecision i) - | otherwise = error "trying to multiply ErrorBound by a non-positive integer" - -instance Arbitrary ErrorBound where - arbitrary = - do - giveSpecialValue <- frequencyElements [(5, False),(1, True)] - aux giveSpecialValue - where - aux giveSpecialValue - | giveSpecialValue = - elements (map convert [0.0,0.0,0.0,10.0,1.0,0.5,0.125]) - | otherwise = - do - (s :: Integer) <- arbitrary - let resultR = ((abs s) `mod` (2^!35))/!(2^!32) - let result = convert resultR - return result diff --git a/src/AERN2/MP/UseMPFR/Float.hs b/src/AERN2/MP/UseMPFR/Float.hs deleted file mode 100644 --- a/src/AERN2/MP/UseMPFR/Float.hs +++ /dev/null @@ -1,60 +0,0 @@ -{-| - Module : AERN2.MP.UseMPFR.Float - Description : Arbitrary precision floating point numbers - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mikkonecny@gmail.com - Stability : experimental - Portability : portable - - Arbitrary precision floating-point numbers with up/down-rounded operations. - - Currently, we use hmpfr when compiling with ghc 7.10 and higher - and haskell-mpfr when compiling with ghc 7.8. --} - -module AERN2.MP.UseMPFR.Float - ( - -- * Precision operations - module AERN2.MP.Precision - -- * The type definition and basic operations - , module AERN2.MP.UseMPFR.Float.Type - -- * Arithmetic operations - , module AERN2.MP.UseMPFR.Float.Arithmetic - , distUp, distDown, avgUp, avgDown - -- * Conversions, comparisons and norm - , module AERN2.MP.UseMPFR.Float.Conversions - -- * Infix operators for up/down-rounded operations - , module AERN2.MP.UseMPFR.Float.Operators - -- * Constants such as NaN, infinity - , module AERN2.MP.UseMPFR.Float.Constants - -- * Tests - , module AERN2.MP.UseMPFR.Float.Tests - ) -where - -import MixedTypesNumPrelude --- import qualified Prelude as P - -import AERN2.MP.Precision -import AERN2.MP.UseMPFR.Float.Type -import AERN2.MP.UseMPFR.Float.Arithmetic -import AERN2.MP.UseMPFR.Float.Conversions -import AERN2.MP.UseMPFR.Float.Operators -import AERN2.MP.UseMPFR.Float.Constants -import AERN2.MP.UseMPFR.Float.Tests - --- | Computes an upper bound to the distance @|x - y|@ of @x@ and @y@. -distUp :: MPFloat -> MPFloat -> MPFloat -distUp x y = if x >= y then x -^ y else y -^ x - --- | Computes a lower bound to the distance @|x - y|@ of @x@ and @y@. -distDown :: MPFloat -> MPFloat -> MPFloat -distDown x y = if x >= y then x -. y else y -. x - -avgUp :: MPFloat -> MPFloat -> MPFloat -avgUp x y = (x +^ y) /^ (mpFloat 2) - -avgDown :: MPFloat -> MPFloat -> MPFloat -avgDown x y = (x +. y) /. (mpFloat 2) diff --git a/src/AERN2/MP/UseMPFR/Float/Arithmetic.hs b/src/AERN2/MP/UseMPFR/Float/Arithmetic.hs deleted file mode 100644 --- a/src/AERN2/MP/UseMPFR/Float/Arithmetic.hs +++ /dev/null @@ -1,170 +0,0 @@ -{-# LANGUAGE CPP #-} -{-| - Module : AERN2.MP.UseMPFR.Float.Arithmetic - Description : Arbitrary precision floating point numbers - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mikkonecny@gmail.com - Stability : experimental - Portability : portable - - Arbitrary precision floating-point numbers with up/down-rounded operations. - - Currently, we use hmpfr when compiling with ghc 7.10 and higher - and haskell-mpfr when compiling with ghc 7.8. --} - -module AERN2.MP.UseMPFR.Float.Arithmetic - ( - -- * MPFloat basic arithmetic - addUp, addDown, subUp, subDown - , mulUp, mulDown, divUp, divDown, recipUp, recipDown - -- * MPFloat selected constants and operations - , piUp, piDown - , cosUp, cosDown, sinUp, sinDown - , sqrtUp, sqrtDown, expUp, expDown, logUp, logDown - ) -where - -import MixedTypesNumPrelude -import qualified Prelude as P - -import AERN2.MP.Precision -import AERN2.MP.UseMPFR.Float.Type - -#ifdef HaskellMPFR -import qualified Data.Approximate.MPFRLowLevel as MPLow - -one :: MPFloat -one = MPLow.fromInt MPLow.Up (P.fromInteger 10) (int 1) - -#endif - -#ifdef MPFRRounded -import qualified AERN2.MP.UseMPFR.Float.RoundedAdaptor as MPLow - -one :: MPFloat -one = MPLow.one - -#endif - -#ifdef HMPFR -import qualified Data.Number.MPFR as MPLow - -one :: MPFloat -one = MPLow.one -#endif - -{- common functions -} - -instance CanNeg MPFloat where - negate = unaryUp MPLow.neg - -instance CanAbs MPFloat where - abs x - | x P.< MPLow.zero = negate x - | otherwise = x - -addUp, addDown :: MPFloat -> MPFloat -> MPFloat -addUp = binaryUp True MPLow.add -addDown = binaryDown True MPLow.add - -subUp, subDown :: MPFloat -> MPFloat -> MPFloat -subUp = binaryUp True MPLow.sub -subDown = binaryDown True MPLow.sub - -mulUp, mulDown :: MPFloat -> MPFloat -> MPFloat -mulUp = binaryUp True MPLow.mul -mulDown = binaryDown True MPLow.mul - -divUp,divDown :: MPFloat -> MPFloat -> MPFloat -divUp = binaryUp False MPLow.div -divDown = binaryDown False MPLow.div - -recipUp :: MPFloat -> MPFloat -recipUp x = divUp one x - -recipDown :: MPFloat -> MPFloat -recipDown x = divDown one x - - -{- special constants and functions -} - -piUp :: Precision -> MPFloat -piUp p = - MPLow.pi MPLow.Up (p2mpfrPrec p) - -piDown :: Precision -> MPFloat -piDown p = - MPLow.pi MPLow.Down (p2mpfrPrec p) - -cosUp :: MPFloat -> MPFloat -cosUp = unaryUp MPLow.cos - -cosDown :: MPFloat -> MPFloat -cosDown = unaryDown MPLow.cos - -sinUp :: MPFloat -> MPFloat -sinUp = unaryUp MPLow.sin - -sinDown :: MPFloat -> MPFloat -sinDown = unaryDown MPLow.sin - -sqrtUp :: MPFloat -> MPFloat -sqrtUp = unaryUp MPLow.sqrt - -sqrtDown :: MPFloat -> MPFloat -sqrtDown = unaryDown MPLow.sqrt - -expUp :: MPFloat -> MPFloat -expUp = unaryUp MPLow.exp - -expDown :: MPFloat -> MPFloat -expDown = unaryDown MPLow.exp - -logUp :: MPFloat -> MPFloat -logUp = unaryUp MPLow.log - -logDown :: MPFloat -> MPFloat -logDown = unaryDown MPLow.log - -{- auxiliary functions to automatically determine result precision from operand precisions -} - -unaryUp :: - (MPLow.RoundMode -> MPLow.Precision -> MPFloat -> MPFloat) -> - (MPFloat -> MPFloat) -unaryUp opRP x = opRP MPLow.Up p x - where - p = MPLow.getPrec x - -unaryDown :: - (MPLow.RoundMode -> MPLow.Precision -> MPFloat -> MPFloat) -> - (MPFloat -> MPFloat) -unaryDown opRP x = opRP MPLow.Down p x - where - p = MPLow.getPrec x - -binaryUp :: - Bool -> - (MPLow.RoundMode -> MPLow.Precision -> MPFloat -> MPFloat -> MPFloat) -> - (MPFloat -> MPFloat -> MPFloat) -binaryUp = binaryApprox True - -binaryDown :: - Bool -> - (MPLow.RoundMode -> MPLow.Precision -> MPFloat -> MPFloat -> MPFloat) -> - (MPFloat -> MPFloat -> MPFloat) -binaryDown = binaryApprox False - -binaryApprox :: - Bool -> Bool -> - (MPLow.RoundMode -> MPLow.Precision -> MPFloat -> MPFloat -> MPFloat) -> - (MPFloat -> MPFloat -> MPFloat) -binaryApprox isUp _canBeExact opRP x y = - withPrec pMax - where - pMax = (getPrecision x) `max` (getPrecision y) - withPrec p - | isUp = opRP MPLow.Up (p2mpfrPrec p) x y - | otherwise = opRP MPLow.Down (p2mpfrPrec p) x y diff --git a/src/AERN2/MP/UseMPFR/Float/Constants.hs b/src/AERN2/MP/UseMPFR/Float/Constants.hs deleted file mode 100644 --- a/src/AERN2/MP/UseMPFR/Float/Constants.hs +++ /dev/null @@ -1,52 +0,0 @@ -{-# LANGUAGE CPP #-} -{-| - Module : AERN2.MP.UseMPFR.Float.Constants - Description : Special constants NaN, infinity etc - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mikkonecny@gmail.com - Stability : experimental - Portability : portable - - Special constants NaN, infinity etc --} - -module AERN2.MP.UseMPFR.Float.Constants - ( - zero, one - , nan, infinity - ) -where - -import MixedTypesNumPrelude -import qualified Prelude as P --- import Data.Ratio - -import AERN2.MP.UseMPFR.Float.Type -import AERN2.MP.UseMPFR.Float.Conversions -import AERN2.MP.UseMPFR.Float.Operators - -zero, one :: MPFloat -zero = mpFloat 0 -one = mpFloat 1 - -nan, infinity :: MPFloat -nan = zero /. zero -infinity = one /. zero - -itisNaN :: MPFloat -> Bool -itisNaN x = x *^ one /= x - -itisInfinite :: MPFloat -> Bool -itisInfinite x = - x *^ (mpFloat 2) P.== x - && - x P./= (mpFloat 0) - -instance CanTestFinite MPFloat where - isInfinite = itisInfinite - isFinite x = not (itisInfinite x || itisNaN x) - -instance CanTestNaN MPFloat where - isNaN = itisNaN diff --git a/src/AERN2/MP/UseMPFR/Float/Conversions.hs b/src/AERN2/MP/UseMPFR/Float/Conversions.hs deleted file mode 100644 --- a/src/AERN2/MP/UseMPFR/Float/Conversions.hs +++ /dev/null @@ -1,208 +0,0 @@ -{-# LANGUAGE CPP #-} -{-| - Module : AERN2.MP.UseMPFR.Float.Conversions - Description : Conversions and comparisons of arbitrary precision floats - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mikkonecny@gmail.com - Stability : experimental - Portability : portable - - Conversions and comparisons of arbitrary precision floating point numbers - - Currently, we use hmpfr when compiling with ghc 7.10 and higher - and haskell-mpfr when compiling with ghc 7.8. --} - -module AERN2.MP.UseMPFR.Float.Conversions - ( - -- * MPFloat to other types (see also instances) - toDoubleUp, toDoubleDown - -- * MPFloat constructors (see also instances) - , CanBeMPFloat, mpFloat - , fromIntegerUp, fromIntegerDown - , fromRationalUp, fromRationalDown - ) -where - -import MixedTypesNumPrelude -import qualified Prelude as P - -import Data.Ratio -import Data.Convertible - -import AERN2.Norm -import AERN2.MP.Precision - -import AERN2.MP.UseMPFR.Float.Type -import AERN2.MP.UseMPFR.Float.Arithmetic - -#ifdef HaskellMPFR -import qualified Data.Approximate.MPFRLowLevel as MPLow - -mpToDouble :: MPLow.RoundMode -> MPFloat -> Double -mpToDouble = MPLow.toDoubleA - -mpToRational :: MPFloat -> Rational -mpToRational x - | x == 0 = 0.0 - | otherwise = MPLow.toRationalA x - -mpFromRationalA :: MPLow.RoundMode -> MPLow.Precision -> Rational -> MPFloat -mpFromRationalA = MPLow.fromRationalA - -#endif -#ifdef MPFRRounded -import qualified AERN2.MP.UseMPFR.Float.RoundedAdaptor as MPLow - -mpToDouble :: MPLow.RoundMode -> MPFloat -> Double -mpToDouble = MPLow.toDoubleA - -mpToRational :: MPFloat -> Rational -mpToRational x - | x == 0 = 0.0 - | otherwise = MPLow.toRationalA x - -mpFromRationalA :: MPLow.RoundMode -> MPLow.Precision -> Rational -> MPFloat -mpFromRationalA = MPLow.fromRationalA - -#endif -#ifdef HMPFR -import qualified Data.Number.MPFR as MPLow - -mpToDouble :: MPLow.RoundMode -> MPLow.MPFR -> Double -mpToDouble = MPLow.toDouble - -mpToRational :: MPFloat -> Rational -mpToRational x - | x == 0 = 0.0 - | otherwise = mantissa * 2.0^!e - where - (mantissa, ePre) = MPLow.decompose x - e = P.toInteger ePre - -mpFromRationalA :: MPLow.RoundMode -> MPLow.Precision -> Rational -> MPFloat -mpFromRationalA dir p q - | q < 0 = - MPLow.fromIntegerA dir p (numerator q) `divDir` MPLow.fromIntegerA dir p (denominator q) - | otherwise = - MPLow.fromIntegerA dir p (numerator q) `divDir` MPLow.fromIntegerA dirOpp p (denominator q) - where - (divDir, dirOpp) = - case dir of - MPLow.Down -> (divDown, MPLow.Up) - MPLow.Up -> (divUp, MPLow.Down) - _ -> error "in mpFromRationalA" - -#endif - -instance HasNorm MPFloat where - getNormLog x - | x == 0 = NormZero - | otherwise = NormBits (P.toInteger $ MPLow.getExp x) - -{- conversions -} - -instance CanRound MPFloat where - properFraction x = (n,f) - where - r = rational x - n = (numerator r) `quot` (denominator r) - f = x `subUp` (mpFloat n) - -instance ConvertibleExactly MPFloat Rational where - safeConvertExactly = Right . mpToRational - -toDoubleUp :: MPFloat -> Double -toDoubleUp = mpToDouble MPLow.Up - -toDoubleDown :: MPFloat -> Double -toDoubleDown = mpToDouble MPLow.Down - -fromIntegerUp :: Precision -> Integer -> MPFloat -fromIntegerUp p i = MPLow.fromIntegerA MPLow.Up (p2mpfrPrec p) i - -fromIntegerDown :: Precision -> Integer -> MPFloat -fromIntegerDown p i = MPLow.fromIntegerA MPLow.Down (p2mpfrPrec p) i - -type CanBeMPFloat t = ConvertibleExactly t MPFloat -mpFloat :: (CanBeMPFloat t) => t -> MPFloat -mpFloat = convertExactly - -instance ConvertibleExactly Integer MPFloat where - safeConvertExactly n = - findExact $ map upDown $ standardPrecisions initPrec - where - initPrec = - case getNormLog n of - NormBits b -> prec (b + 8) - _ -> prec 8 - upDown p = (fromIntegerDown p n, fromIntegerUp p n) - findExact [] = - convError "integer too high to represent exactly" n - findExact ((nDown, nUp) : rest) - | nDown == nUp = Right nUp - | otherwise = findExact rest - -instance ConvertibleExactly Int MPFloat where - safeConvertExactly = safeConvertExactly . integer - -fromRationalUp :: Precision -> Rational -> MPFloat -fromRationalUp p x = - mpFromRationalA MPLow.Up (p2mpfrPrec p) x - -fromRationalDown :: Precision -> Rational -> MPFloat -fromRationalDown p x = - mpFromRationalA MPLow.Down (p2mpfrPrec p) x - -instance Convertible MPFloat Double where - safeConvert x - | isFinite dbl = Right dbl - | otherwise = convError "conversion to double: out of bounds" x - where - dbl = toDoubleUp x - -{- comparisons -} - -instance HasEqAsymmetric MPFloat MPFloat -instance HasEqAsymmetric MPFloat Integer where - equalTo = convertSecond equalTo -instance HasEqAsymmetric Integer MPFloat where - equalTo = convertFirst equalTo -instance HasEqAsymmetric MPFloat Int where - equalTo = convertSecond equalTo -instance HasEqAsymmetric Int MPFloat where - equalTo = convertFirst equalTo -instance HasEqAsymmetric MPFloat Rational where - equalTo = convertFirst equalTo -instance HasEqAsymmetric Rational MPFloat where - equalTo = convertSecond equalTo - -instance CanTestZero MPFloat - -instance HasOrderAsymmetric MPFloat MPFloat -instance HasOrderAsymmetric MPFloat Integer where - lessThan = convertSecond lessThan - leq = convertSecond leq -instance HasOrderAsymmetric Integer MPFloat where - lessThan = convertFirst lessThan - leq = convertFirst leq -instance HasOrderAsymmetric MPFloat Int where - lessThan = convertSecond lessThan - leq = convertSecond leq -instance HasOrderAsymmetric Int MPFloat where - lessThan = convertFirst lessThan - leq = convertFirst leq -instance HasOrderAsymmetric Rational MPFloat where - lessThan = convertSecond lessThan - leq = convertSecond leq -instance HasOrderAsymmetric MPFloat Rational where - lessThan = convertFirst lessThan - leq = convertFirst leq - -instance CanTestPosNeg MPFloat - -{- min, max -} - -instance CanMinMaxAsymmetric MPFloat MPFloat diff --git a/src/AERN2/MP/UseMPFR/Float/Operators.hs b/src/AERN2/MP/UseMPFR/Float/Operators.hs deleted file mode 100644 --- a/src/AERN2/MP/UseMPFR/Float/Operators.hs +++ /dev/null @@ -1,38 +0,0 @@ -{-| - Module : AERN2.MP.UseMPFR.Float.Operators - Description : Infix operators for up/down-rounded floating-point numbers - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mikkonecny@gmail.com - Stability : experimental - Portability : portable - - Infix operators for up/down-rounded floating-point numbers --} - -module AERN2.MP.UseMPFR.Float.Operators where - -import AERN2.MP.UseMPFR.Float.Type -import AERN2.MP.UseMPFR.Float.Arithmetic - -infixl 6 +^, -^, +., -. -infixl 7 *^, *., /^, /. - -(+^) :: MPFloat -> MPFloat -> MPFloat -(+^) = addUp -(-^) :: MPFloat -> MPFloat -> MPFloat -(-^) = subUp -(*^) :: MPFloat -> MPFloat -> MPFloat -(*^) = mulUp -(/^) :: MPFloat -> MPFloat -> MPFloat -(/^) = divUp - -(+.) :: MPFloat -> MPFloat -> MPFloat -(+.) = addDown -(-.) :: MPFloat -> MPFloat -> MPFloat -(-.) = subDown -(*.) :: MPFloat -> MPFloat -> MPFloat -(*.) = mulDown -(/.) :: MPFloat -> MPFloat -> MPFloat -(/.) = divDown diff --git a/src/AERN2/MP/UseMPFR/Float/RoundedAdaptor.hs b/src/AERN2/MP/UseMPFR/Float/RoundedAdaptor.hs deleted file mode 100644 --- a/src/AERN2/MP/UseMPFR/Float/RoundedAdaptor.hs +++ /dev/null @@ -1,90 +0,0 @@ -{-# LANGUAGE CPP #-} -{-# LANGUAGE DataKinds, ExistentialQuantification, RankNTypes #-} --- {-# LANGUAGE DeriveGeneric, DeriveDataTypeable, StandaloneDeriving #-} -{-| - Module : AERN2.MP.UseMPFR.Float.RoundedAdaptor - Description : Numeric.Rounded + variable precision - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mikkonecny@gmail.com - Stability : experimental - Portability : portable - - Numeric.Rounded + variable precision --} -module AERN2.MP.UseMPFR.Float.RoundedAdaptor -#ifndef MPFRRounded -() where -#else -( - module AERN2.MP.UseMPFR.Float.RoundedAdaptor -, module Numeric.Rounded.Simple -) -where - -import Prelude hiding (div, pi) --- import qualified Prelude as P - -import Numeric.Rounded.Simple --- import qualified Numeric.RoundedSimple as R - -instance Show Rounded where - show = show' - -getPrec :: Rounded -> Int -getPrec = precision - -getExp :: Rounded -> Int -getExp = exponent' - -data RoundMode = Up | Down - -withRoundMode :: (RoundingMode -> t) -> (RoundMode -> t) -withRoundMode op Up = op TowardInf -withRoundMode op Down = op TowardNegInf -{-# INLINE withRoundMode #-} - -set :: RoundMode -> Precision -> Rounded -> Rounded -set = withRoundMode precRound - -defaultPrecision :: Precision -defaultPrecision = 10 - -pi :: RoundMode -> Precision -> Rounded -pi = withRoundMode kPi - -fromIntegerA :: RoundMode -> Precision -> Integer -> Rounded -fromIntegerA = withRoundMode fromInteger' - -zero, one :: Rounded -zero = fromIntegerA Up defaultPrecision 0 -one = fromIntegerA Up defaultPrecision 1 - -toDoubleA :: RoundMode -> Rounded -> Double -toDoubleA = withRoundMode toDouble - -fromRationalA :: RoundMode -> Precision -> Rational -> Rounded -fromRationalA = withRoundMode fromRational' - -toRationalA :: Rounded -> Rational -toRationalA = toRational' TowardNearest - -add, sub, mul, div, atan2 :: RoundMode -> Precision -> Rounded -> Rounded -> Rounded -add = withRoundMode add_ -sub = withRoundMode sub_ -mul = withRoundMode mul_ -div = withRoundMode div_ -atan2 = withRoundMode atan2_ - -neg, abs, sqrt, exp, log, sin, cos :: RoundMode -> Precision -> Rounded -> Rounded -neg = withRoundMode negate_ -abs = withRoundMode abs_ -sqrt = withRoundMode sqrt_ -exp = withRoundMode exp_ -log = withRoundMode log_ -sin = withRoundMode sin_ -cos = withRoundMode cos_ --- TODO: add more ops - -#endif diff --git a/src/AERN2/MP/UseMPFR/Float/Tests.hs b/src/AERN2/MP/UseMPFR/Float/Tests.hs deleted file mode 100644 --- a/src/AERN2/MP/UseMPFR/Float/Tests.hs +++ /dev/null @@ -1,424 +0,0 @@ -{-# LANGUAGE CPP #-} -{-| - Module : AERN2.MP.UseMPFR.Float.Tests - Description : Tests for operations on arbitrary precision floats - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mikkonecny@gmail.com - Stability : experimental - Portability : portable - - Tests for operations on arbitrary precision floats. - - To run the tests using stack, execute: - - @ - stack test aern2-mp --test-arguments "-a 1000 -m MPFloat" - @ --} - -module AERN2.MP.UseMPFR.Float.Tests - ( - specMPFloat, tMPFloat - , (=~=), approxEqual, approxEqualWithArgs - , frequencyElements - ) -where - -import MixedTypesNumPrelude --- import qualified Prelude as P --- import Data.Ratio -import Text.Printf -import Data.Maybe - -import Test.Hspec -import Test.QuickCheck --- import qualified Test.Hspec.SmallCheck as SC - - -import AERN2.Norm -import AERN2.MP.Precision - -import AERN2.MP.UseMPFR.Float.Type -import AERN2.MP.UseMPFR.Float.Arithmetic -import AERN2.MP.UseMPFR.Float.Conversions -import AERN2.MP.UseMPFR.Float.Operators -import AERN2.MP.UseMPFR.Float.Constants - -instance Arbitrary MPFloat where - arbitrary = - do - giveSpecialValue <- frequencyElements [(9, False),(1, True)] - aux giveSpecialValue - where - aux giveSpecialValue - | giveSpecialValue = - elements [nan, infinity, -infinity, zero, one, -one] - | otherwise = - do - (p :: Precision) <- arbitrary - (s :: Integer) <- arbitrary - ex <- choose (-20,10) - let resultR = s * (10.0^!ex) - let result = fromRationalUp p resultR - return result - -frequencyElements :: ConvertibleExactly t Int => [(t, a)] -> Gen a -frequencyElements elems = frequency [(int n, return e) | (n,e) <- elems] - -{- approximate comparison -} - -infix 4 =~= - -(=~=) :: MPFloat -> MPFloat -> Property -l =~= r = - approxEqualWithArgs [] l r - -{-| - Assert equality of two MPFloat's with tolerance @1/2^p@. --} -approxEqual :: - Integer {-^ @p@ precision to guide tolerance -} -> - MPFloat {-^ LHS of equation-} -> - MPFloat {-^ RHS of equation -}-> - Bool -approxEqual e x y - | isNaN x && isNaN y = True - | isNaN x && isInfinite y = True - | isInfinite x && isNaN y = True - | isNaN x || isNaN y = False - | isInfinite x || isInfinite y = x == y - | otherwise = - abs (x -. y) <= 0.5^!e - -{-| - Assert equality of two MPFloat's with tolerance derived from the size and precision - of the given intermediate values. - When the assertion fails, report the given values using the given names. --} -approxEqualWithArgs :: - [(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 -approxEqualWithArgs argsPre l r = - counterexample description $ approxEqual e l r - where - args = argsPre ++ [(l, "L"), (r, "R"), (abs (l-.r),"|L-R|")] - e = - (foldl min 1000000 $ catMaybes $ map getNminusP args) - - (length argsPre) - getNminusP (x,_) = - case norm of - NormZero -> Nothing -- ideally infinity - NormBits b -> Just (pI-b-1) - where - norm = getNormLog x - pI = integer $ getPrecision x - description = - printf "args:\n%s tolerance: <= %s (e=%d)" argsS (show (double (0.5^!e))) e - argsS = - unlines - [printf " %s = %s (p=%s)" argS (show arg) (show $ getPrecision arg) | (arg, argS) <- args] - -{-| - A runtime representative of type @MPFloat@. - Used for specialising polymorphic tests to concrete types. --} -tMPFloat :: T MPFloat -tMPFloat = T "MPFloat" - -specMPFloat :: Spec -specMPFloat = - describe ("MPFloat") $ do - specCanSetPrecision tMPFloat (printArgsIfFails2 "=~=" (=~=)) - specCanRound tMPFloat - specCanNegNum tMPFloat - specCanAbs tMPFloat - specCanMinMaxNotMixed tMPFloat - -- specCanMinMax tMPFloat tInteger tMPFloat - describe "special values" $ do - it "0 * infinity = NaN" $ do - isNaN (zero *^ infinity) - && - isNaN (zero *. infinity) - it "infinity / infinity = NaN" $ do - isNaN (infinity /^ infinity) - && - isNaN (infinity /. infinity) - it "infinity - infinity = NaN" $ do - isNaN (infinity -^ infinity) - && - isNaN (infinity -. infinity) - describe "approximate addition" $ do - it "down <= up" $ do - property $ \ (x :: MPFloat) (y :: MPFloat) -> - not (isNaN (x +. y)) - ==> - x +. y <= x +^ y - it "up ~ down" $ do - property $ \ (x :: MPFloat) (y :: MPFloat) -> - x +. y =~= x +^ y - it "absorbs 0" $ do - property $ \ (x :: MPFloat) -> - (not $ isNaN x) ==> - x +. (mpFloat 0) == x - it "approximately commutative" $ do - property $ \ (x :: MPFloat) (y :: MPFloat) -> - (not $ isNaN $ x +. y) ==> - x +. y <= y +^ x - && - x +^ y >= y +. x - it "approximately associative" $ do - property $ \ (x :: MPFloat) (y :: MPFloat) (z :: MPFloat) -> - (not $ isNaN $ x +. y +. z) ==> - (x +. y) +. z <= x +^ (y +^ z) - && - (x +^ y) +^ z >= x +. (y +. z) - describe "approximate subtraction" $ do - it "down <= up" $ do - property $ \ (x :: MPFloat) (y :: MPFloat) -> - not (isNaN (x -. y)) - ==> - x -. y <= x -^ y - it "up ~ down" $ do - property $ \ (x :: MPFloat) (y :: MPFloat) -> - x -. y =~= x -^ y - it "same as negate and add" $ do - property $ \ (x :: MPFloat) (y :: MPFloat) -> - (not $ isNaN $ x -. y) ==> - x -. y <= x +^ (-y) - && - x -^ y >= x +. (-y) - describe "approximate multiplication" $ do - it "down <= up" $ do - property $ \ (x :: MPFloat) (y :: MPFloat) -> - not (isNaN (x *. y)) - ==> - x *. y <= x *^ y - it "up ~ down" $ do - property $ \ (x :: MPFloat) (y :: MPFloat) -> - x *. y =~= x *^ y - it "absorbs 1" $ do - property $ \ (x :: MPFloat) -> - (not $ isNaN x) ==> - x *. (mpFloat 1) == x - it "approximately commutative" $ do - property $ \ (x :: MPFloat) (y :: MPFloat) -> - not (isNaN (x *. y)) ==> - x *. y <= y *^ x - && - x *^ y >= y *. x - it "approximately associative" $ do - property $ \ (x :: MPFloat) (y :: MPFloat) (z :: MPFloat) -> - (x >= 0 && y >= 0 && z >= 0 - && not (isInfinite x) && not (isInfinite y) && not (isInfinite z)) ==> - (x *. y) *. z <= x *^ (y *^ z) - && - (x *^ y) *^ z >= x *. (y *. z) - it "approximately distributes over addition" $ do - property $ \ (x :: MPFloat) (y :: MPFloat) (z :: MPFloat) -> - (x >= 0 && y >= 0 && z >= 0 - && not (isInfinite x) && not (isInfinite y) && not (isInfinite z)) ==> - x *. (y +. z) <= (x *^ y) +^ (x *^ z) - && - x *^ (y +^ z) >= (x *. y) +. (x *. z) - describe "approximate division" $ do - it "down <= up" $ do - property $ \ (x :: MPFloat) (y :: MPFloat) -> - not (isNaN (x /. y)) - ==> - x /. y <= x /^ y - it "up ~ down" $ do - property $ \ (x :: MPFloat) (y :: MPFloat) -> - let - (=~~=) = approxEqualWithArgs [(x /. y,"x/.y")] - infix 4 =~~= - in - not (isNaN (x /. y)) - ==> - x /. y =~~= x /^ y - it "recip(recip x) = x" $ do - property $ \ (x :: MPFloat) -> - (x > 0 || x < 0) ==> - one /. (one /^ x) <= x - && - one /^ (one /. x) >= x - it "x/1 = x" $ do - property $ \ (x :: MPFloat) -> - not (isNaN x) ==> - (x /. one) == x - it "x/x = 1" $ do - property $ \ (x :: MPFloat) -> - (isCertainlyNonZero x && (not $ isNaN $ x /. x)) ==> - (x /. x) <= one - && - (x /^ x) >= one - it "x/y = x*(1/y)" $ do - property $ \ (x :: MPFloat) (y :: MPFloat) -> - (y > 0 && x >= 0 && x/.y >= 0) ==> - (x /. y) <= x *^ (one /^ y) - && - (x /^ y) >= x *. (one /. y) - describe "approximate sqrt" $ do - it "down <= up" $ do - property $ \ (x :: MPFloat) -> - not (isNaN (sqrtDown x)) - ==> - sqrtDown x <= sqrtUp x - it "up ~ down" $ do - property $ \ (x :: MPFloat) -> - (x >= 0) - ==> - sqrtDown x =~= sqrtUp x - it "sqrt(x) >= 0" $ do - property $ \ (x :: MPFloat) -> - (x >= 0) - ==> - sqrtUp x >= 0 - it "sqrt(x)^2 ~ x" $ do - property $ \ (x :: MPFloat) -> - (x >= 0) - ==> - (sqrtDown x) *. (sqrtDown x) <= x - && - (sqrtUp x) *^ (sqrtUp x) >= x - describe "approximate exp" $ do - it "down <= up" $ do - property $ \ (x :: MPFloat) -> - (abs x < 1000000) - ==> - expDown x <= expUp x - it "up ~ down" $ do - property $ \ (x :: MPFloat) -> - (abs x < 1000000) - ==> - let - (=~~=) = approxEqualWithArgs [(x,"x")] - infix 4 =~~= - in - expDown x =~~= expUp x - it "exp(-x) == 1/(exp x)" $ do - property $ \ (x :: MPFloat) -> - (abs x < 1000000) - ==> - one /. (expUp x) <= expUp (-x) - && - one /^ (expDown x) >= expDown (-x) - it "exp(x+y) = exp(x)*exp(y)" $ do - property $ \ (x :: MPFloat) (y :: MPFloat) -> - (abs x < 1000000 && abs y < 1000000) - ==> - expDown (x +. y) <= (expUp x) *^ (expUp y) - && - expUp (x +^ y) >= (expDown x) *. (expDown y) - describe "approximate log" $ do - it "down <= up" $ do - property $ \ (x :: MPFloat) -> - (x > 0) - ==> - logDown x <= logUp x - it "up ~ down" $ do - property $ \ (x :: MPFloat) -> - (x > 0) - ==> - logDown x =~= logUp x - it "log(1/x) == -(log x)" $ do - property $ \ (x :: MPFloat) -> - (x > 0) - ==> - logDown (one /. x) <= -(logDown x) - && - logUp (one /^ x) >= -(logUp x) - it "log(x*y) = log(x)+log(y)" $ do - property $ \ (x :: MPFloat) (y :: MPFloat) -> - (x > 0 && y > 0) - ==> - logDown (x *. y) <= (logUp x) +^ (logUp y) - && - logUp (x *^ y) >= (logDown x) +. (logDown y) - it "log(exp x) == x" $ do - property $ \ (x :: MPFloat) -> - (abs x < 1000000) - ==> - logDown (expDown x) <= x - && - logUp (expUp x) >= x - describe "approximate sine" $ do - it "down <= up" $ do - property $ \ (x :: MPFloat) -> - (abs x < 1000000) - ==> - sinDown x <= sinUp x - it "up ~ down" $ do - property $ \ (x :: MPFloat) -> - (abs x < 1000000) - ==> - let - (=~~=) = approxEqualWithArgs [(x,"x")] - infix 4 =~~= - in - sinDown x =~~= sinUp x - it "sin(pi)=0" $ do - property $ \ (p :: Precision) -> - let - (=~~=) = approxEqualWithArgs [(piDown p,"pi")] - infix 4 =~~= - in - sinUp(piDown p) =~~= (fromIntegerUp p 0) - it "in [-1,1]" $ do - property $ \ (x :: MPFloat) -> - (abs x < 1000000) - ==> - sinDown x <= one - && - sinUp x >= -one - describe "approximate cosine" $ do - it "down <= up" $ do - property $ \ (x :: MPFloat) -> - (abs x < 1000000) - ==> - cosDown x <= cosUp x - it "up ~ down" $ do - property $ \ (x :: MPFloat) -> - (abs x < 1000000) - ==> - let - (=~~=) = approxEqualWithArgs [(x,"x")] - infix 4 =~~= - in - cosDown x =~~= cosUp x - it "in [-1,1]" $ do - property $ \ (x :: MPFloat) -> - (abs x < 1000000) - ==> - cosDown x <= one - && - cosUp x >= -one - it "cos(pi)=-1" $ do - property $ \ (p :: Precision) -> - cosUp(piDown p) =~= (fromIntegerUp p (-1)) - it "cos(x)^2 + sin(x)^2 = 1" $ do - property $ \ (x :: MPFloat) -> - (abs x < 1000000) - ==> - let - cosxU = cosUp x - cosxD = cosDown x - cosx2U = (cosxU *^ cosxU) `max` (cosxD *^ cosxD) - cosx2D - | cosxD > 0 = cosxD *. cosxD - | cosxU < 0 = cosxU *. cosxU - | otherwise = mpFloat 0 - sinxU = sinUp x - sinxD = sinDown x - sinx2U = (sinxU *^ sinxU) `max` (sinxD *^ sinxD) - sinx2D - | sinxD > 0 = sinxD *. sinxD - | sinxU < 0 = sinxU *. sinxU - | otherwise = mpFloat 0 - in - (cosx2D +. sinx2D) <= one - && - (cosx2U +^ sinx2U) >= one diff --git a/src/AERN2/MP/UseMPFR/Float/Type.hs b/src/AERN2/MP/UseMPFR/Float/Type.hs deleted file mode 100644 --- a/src/AERN2/MP/UseMPFR/Float/Type.hs +++ /dev/null @@ -1,91 +0,0 @@ -{-# LANGUAGE CPP #-} -{-# LANGUAGE DeriveGeneric, DeriveDataTypeable, StandaloneDeriving #-} -{-| - Module : AERN2.MP.UseMPFR.Float.Type - Description : Arbitrary precision floating point numbers - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mikkonecny@gmail.com - Stability : experimental - Portability : portable - - Arbitrary precision floating-point numbers - - Currently, we use hmpfr when compiling with ghc 7.10 and higher - and haskell-mpfr when compiling with ghc 7.8. --} - -module AERN2.MP.UseMPFR.Float.Type - ( - -- * MPFloat numbers and their basic operations - MPFloat, setPrecisionUp, setPrecisionDown - , p2mpfrPrec - ) -where - -import MixedTypesNumPrelude -import qualified Prelude as P - -import AERN2.MP.Precision - -#ifdef IntegerBackend -import qualified AERN2.MP.UseMPFR.Float.Native as MPLow - -type MPFloat = MPLow.MPFloat - -p2mpfrPrec :: Precision -> Precision -p2mpfrPrec = id - -#endif - -#ifdef HaskellMPFR -import qualified Data.Approximate.MPFRLowLevel as MPLow -import Data.Typeable - -{-| Multiple-precision floating-point type based on MPFR via haskell-mpfr. -} -type MPFloat = MPLow.Rounded - -deriving instance (Typeable MPFloat) - -p2mpfrPrec :: Precision -> MPLow.Precision -p2mpfrPrec = P.fromInteger . integer - -#endif - -#ifdef MPFRRounded -import qualified AERN2.MP.UseMPFR.Float.RoundedAdaptor as MPLow -import Data.Typeable - -{-| Multiple-precision floating-point type based on MPFR via rounded. -} -type MPFloat = MPLow.Rounded - -deriving instance (Typeable MPFloat) - -p2mpfrPrec :: Precision -> MPLow.Precision -p2mpfrPrec = P.fromInteger . integer - -#endif - -#ifdef HMPFR -import qualified Data.Number.MPFR as MPLow - -{-| Multiple-precision floating-point type based on MPFR via hmpfr. -} -type MPFloat = MPLow.MPFR - -p2mpfrPrec :: Precision -> MPLow.Precision -p2mpfrPrec = P.fromInteger . integer - -#endif - -instance HasPrecision MPFloat where - getPrecision x = prec (P.toInteger $ MPLow.getPrec x) - -instance CanSetPrecision MPFloat where - setPrecision = setPrecisionUp - -setPrecisionUp :: Precision -> MPFloat -> MPFloat -setPrecisionUp p = MPLow.set MPLow.Up (p2mpfrPrec p) - -setPrecisionDown :: Precision -> MPFloat -> MPFloat -setPrecisionDown p = MPLow.set MPLow.Down (p2mpfrPrec p) diff --git a/test/AERN2/MP/FloatSpec.hs b/test/AERN2/MP/FloatSpec.hs new file mode 100644 --- /dev/null +++ b/test/AERN2/MP/FloatSpec.hs @@ -0,0 +1,20 @@ +{-| + Module : AERN2.MP.FloatSpec + Description : hspec tests for MPFloat + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable +-} + +module AERN2.MP.FloatSpec (spec) where + +-- import MixedTypesNumPrelude +import AERN2.MP.Float.Tests + +import Test.Hspec + +spec :: Spec +spec = specMPFloat diff --git a/test/AERN2/MP/UseMPFR/FloatSpec.hs b/test/AERN2/MP/UseMPFR/FloatSpec.hs deleted file mode 100644 --- a/test/AERN2/MP/UseMPFR/FloatSpec.hs +++ /dev/null @@ -1,20 +0,0 @@ -{-| - Module : AERN2.MP.FloatSpec - Description : hspec tests for MPFloat - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mikkonecny@gmail.com - Stability : experimental - Portability : portable --} - -module AERN2.MP.UseMPFR.FloatSpec (spec) where - --- import MixedTypesNumPrelude -import AERN2.MP.UseMPFR.Float.Tests - -import Test.Hspec - -spec :: Spec -spec = specMPFloat