diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright (c) 2015-2017 Michal Konecny +Copyright (c) 2015-2021 Michal Konecny All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,87 @@ +# aern2-mp + +Variable-precision interval arithmetic + +## Numeric data types + +This package provides the following two data types: + +* `Dyadic`: variable-precision floats with exact ring operations + +* `MPBall`: variable-precision interval (float centre ± error bound) with field & elementary interval operations + +The type `MPBall` has instances of both [mixed-types-num](https://hackage.haskell.org/package/mixed-types-num) type classes such as `CanAdd`, `CanSqrt` as well as with traditional Prelude type classes such as `Ord`, `Num` and `Floating`. +The type `Dyadic` also has an appropriate subset of such instances. + +Package [aern2-real](../aern2-real/README.md) provides an arithmetic of exact real numbers as converging lazy sequences of `MPBall`s of increasing precision. Exact real numbers offer additional convenience and readability to validated numeric programming. + +### Examples + +First, let us test interval arithmetic with **Prelude** operations: + + $ stack ghci aern2-mp:lib --no-load --ghci-options AERN2.MP + *AERN2.MP> import Prelude + *AERN2.MP Prelude> + + ...> pi100 = piBallP (prec 100) + ...> pi10000 = piBallP (prec 10000) + + ...> pi100 ^ 2 + [9.8696044010893586188344909998725639610631902560... ± ~8.1120e-30 ~2^(-96)] + + ...> pi100 ^ pi100 + :18:1: error: + • No instance for (Integral MPBall) arising from a use of ‘^’ + + ...> sin pi100 + [0.0000000000000000000000000000001694818351060767... ± ~7.8925e-31 ~2^(-99)] + + ...> sin pi10000 + [0.0000000000000000000000000000000000000000000000... ± ~0.0000 ~2^(-9999)] + (0.07 secs, 64,466,432 bytes) + + ...> pi100 > 0 + True + + ...> pi100 == pi100 + *** Exception: Failed to decide equality of MPBalls. If you switch to MixedTypesNumPrelude instead of Prelude, comparison of MPBalls returns Kleenean instead of Bool. + +Some things do not work with Prelude. Let us try using **MixedTypesNumPrelude** operations: + + $ stack ghci aern2-mp:lib --no-load --ghci-options AERN2.MP + *AERN2.MP> import MixedTypesNumPrelude + *AERN2.MP MixedTypesNumPrelude> + + ...> pi100 = piBallP (prec 100) + ...> pi10000 = piBallP (prec 10000) + + ...> pi100 ^ pi100 + [36.4621596072079117709908260226198218149834948802... ± ~1.8696e-28 ~2^(-92)] + + ...> pi10000 ^ pi10000 + [36.4621596072079117709908260226921236663655084022... ± ~0.0000 ~2^(-9992)] + (0.27 secs, 204,657,248 bytes) + + ...> pi100 > 0 + CertainTrue + + ...> pi100 == pi100 + TrueOrFalse + +### Internal types and backends + +The type `MPBall` internally uses the type: + +* `MPFloat`: arbitrary-precision floats with both upwards and downwards-rounded arithmetic operations such as `*^` and `*.` + +The package uses [cdar-mBound](https://hackage.haskell.org/package/cdar-mBound), a fork of [cdar](https://github.com/jensblanck/cdar) as its backend for `Dyadic` and `MPFloat`. + +In previous versions, there was an MPFR backend via [rounded](https://hackage.haskell.org/package/rounded). This may be added again in future. + +## Specifications and tests + +This package also provides a fairly complete hspec/QuickCheck specification of algebraic properties for the above types. + +For `MPFloat`, the properties are given mostly as approximate versions of algebraic equalities with a small rounding error tolerance. + +For `MPBall`, the properties are given mostly as (interval) set over-approximations of the usual algebraic equalities. diff --git a/aern2-mp.cabal b/aern2-mp.cabal --- a/aern2-mp.cabal +++ b/aern2-mp.cabal @@ -1,128 +1,106 @@ +cabal-version: 1.12 + +-- This file has been generated from package.yaml by hpack version 0.33.0. +-- +-- see: https://github.com/sol/hpack +-- +-- hash: 210a1027cd95f3ead272fe18229ed64036a8b1408d9ed7558c82f9655f6b8248 + name: aern2-mp -version: 0.1.4 -cabal-version: >= 1.9.2 -build-type: Simple -homepage: https://github.com/michalkonecny/aern2 +version: 0.2.0.0 +synopsis: Multi-precision ball (interval) arithmetic +description: Please see the README on GitHub at +category: Math +homepage: https://github.com/michalkonecny/aern2#readme +bug-reports: https://github.com/michalkonecny/aern2/issues author: Michal Konecny -maintainer: Michal Konecny -copyright: (c) 2015-2019 Michal Konecny +maintainer: mikkonecny@gmail.com +copyright: 2015-2021 Michal Konecny license: BSD3 license-file: LICENSE -extra-source-files: changelog.md -stability: experimental -category: Math -synopsis: Multi-precision ball (interval) arithmetic -Description: - This package provides the following types: - . - * Dyadic: variable-precision floats with exact ring operations - . - * MPBall: float ± error bound with field & elementary interval-like operations - . - The types have instances of both - type classes as well as with traditional Prelude type classes. - . - There is a plan to add an Integer-only backend so that aern2-mp can - be used without MPFR. +build-type: Simple +extra-source-files: + README.md + changelog.md source-repository head - type: git - location: https://github.com/mikkonecny/aern2.git - subdir: aern2-mp - -flag UseCDAR - Description: Use CDAR (mBound branch) as an Integer-only backend instead of MPFR - Default: False + type: git + location: https://github.com/michalkonecny/aern2 library - hs-source-dirs: src + exposed-modules: + AERN2.Limit + AERN2.MP + AERN2.MP.Accuracy + AERN2.MP.Ball + AERN2.MP.Ball.Comparisons + AERN2.MP.Ball.Conversions + AERN2.MP.Ball.Elementary + AERN2.MP.Ball.Field + AERN2.MP.Ball.Limit + AERN2.MP.Ball.PreludeOps + AERN2.MP.Ball.Tests + AERN2.MP.Ball.Type + AERN2.MP.Dyadic + AERN2.MP.Enclosure + AERN2.MP.ErrorBound + AERN2.MP.Float + AERN2.MP.Float.Arithmetic + AERN2.MP.Float.Auxi + AERN2.MP.Float.Conversions + AERN2.MP.Float.Operators + AERN2.MP.Float.Tests + AERN2.MP.Float.Type + AERN2.MP.Precision + AERN2.MP.WithCurrentPrec + AERN2.Norm + AERN2.Normalize + AERN2.Utils.Bench + other-modules: + Paths_aern2_mp + hs-source-dirs: + src + default-extensions: RebindableSyntax, ScopedTypeVariables, DeriveGeneric, TypeFamilies, TypeOperators, ConstraintKinds, DefaultSignatures, MultiParamTypeClasses, FlexibleContexts, FlexibleInstances, UndecidableInstances + other-extensions: TemplateHaskell + ghc-options: -Wall build-depends: - base == 4.* + QuickCheck + , base ==4.* + , cdar-mBound >=0.1.0.0 + , collect-errors >=0.1.3 + , deepseq + , hspec , integer-logarithms - , convertible + , mixed-types-num >=0.5.1 + , reflection , regex-tdfa - , hspec - -- , hspec-smallcheck >= 0.3 - , QuickCheck - , lens , template-haskell - , mixed-types-num >= 0.3.2 - if flag(UseCDAR) - hs-source-dirs: src-cdar - build-depends: - cdar - else - hs-source-dirs: src-rounded - build-depends: - rounded == 0.1.* - ghc-options: -Wall -fno-warn-orphans - extensions: - RebindableSyntax, - PostfixOperators, - ScopedTypeVariables, - DeriveGeneric, - GeneralizedNewtypeDeriving, - TypeFamilies, - TypeOperators, - ConstraintKinds, - DefaultSignatures, - MultiParamTypeClasses, - FlexibleContexts, - FlexibleInstances, - UndecidableInstances - if !flag(UseCDAR) - exposed-modules: - AERN2.MP.Float.RoundedAdaptor - exposed-modules: - -- modules that depend on backend choice: - AERN2.MP.Float.Type - AERN2.MP.Float.Arithmetic - AERN2.MP.Float.Conversions - -- modules common to all backends: - AERN2.Utils.Bench - AERN2.Normalize - AERN2.Norm - AERN2.MP.Precision - AERN2.MP.Accuracy - AERN2.MP.Enclosure - AERN2.MP.ErrorBound - AERN2.MP.Float.Auxi - AERN2.MP.Float.Operators - 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 + default-language: Haskell2010 -test-suite spec - type: - exitcode-stdio-1.0 - ghc-options: - -Wall - extensions: - RebindableSyntax, - PostfixOperators, - ScopedTypeVariables, - FlexibleContexts +test-suite aern2-mp-test + type: exitcode-stdio-1.0 + main-is: Spec.hs + other-modules: + AERN2.MP.BallSpec + AERN2.MP.DyadicSpec + AERN2.MP.FloatSpec + Paths_aern2_mp hs-source-dirs: test - main-is: - Spec.hs - other-modules: - AERN2.MP.BallSpec - AERN2.MP.DyadicSpec - AERN2.MP.FloatSpec + default-extensions: RebindableSyntax, ScopedTypeVariables, DeriveGeneric, TypeFamilies, TypeOperators, ConstraintKinds, DefaultSignatures, MultiParamTypeClasses, FlexibleContexts, FlexibleInstances, UndecidableInstances + other-extensions: TemplateHaskell + ghc-options: -threaded -rtsopts -with-rtsopts=-N -Wall build-depends: - base == 4.* - -- , mixed-types-num >= 0.3.1 && < 0.4 + QuickCheck , aern2-mp - , hspec >= 2.1 - -- , hspec-smallcheck >= 0.3 - , QuickCheck >= 2.7 + , base ==4.* + , cdar-mBound >=0.1.0.0 + , collect-errors >=0.1.3 + , deepseq + , hspec + , integer-logarithms + , mixed-types-num >=0.5.1 + , reflection + , regex-tdfa + , template-haskell + default-language: Haskell2010 diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -1,5 +1,11 @@ # Change log for aern2-mp +* v 0.2.0 2021-05-17 + * switch to new simplified collect-errors, mixed-types-num 0.5.0 + * got rid of EnsureCE etc. + * not introducing CN wrapper unless at least one parameter is already CN + * using CDAR backend only, no MPFR for now + * WithCurrentPrec for specifying default precision via types * v 0.1.4 2019-03-19 * CDAR-based Integer-only backend * needs the mBound branch of CDAR diff --git a/src-cdar/AERN2/MP/Float/Arithmetic.hs b/src-cdar/AERN2/MP/Float/Arithmetic.hs deleted file mode 100644 --- a/src-cdar/AERN2/MP/Float/Arithmetic.hs +++ /dev/null @@ -1,104 +0,0 @@ -{-| - Module : AERN2.MP.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. --} - -module AERN2.MP.Float.Arithmetic - ( - -- * MPFloat basic arithmetic - addCEDU, subCEDU - , mulCEDU, divCEDU, recipCEDU - -- * MPFloat selected constants and operations - , piCEDU - , cosCEDU, sinCEDU - , sqrtCEDU, expCEDU, logCEDU - ) -where - -import MixedTypesNumPrelude -import qualified Prelude as P - -import AERN2.MP.Precision - -import qualified Data.CDAR as MPLow - -import AERN2.MP.Float.Auxi -import AERN2.MP.Float.Type - -{- common functions -} - -instance CanNeg MPFloat where - negate = P.negate - -instance CanAbs MPFloat where - abs = P.abs - -addCEDU :: MPFloat -> MPFloat -> BoundsCEDU MPFloat -addCEDU = binaryCEDU (P.+) - -subCEDU :: MPFloat -> MPFloat -> BoundsCEDU MPFloat -subCEDU = binaryCEDU (P.-) - -mulCEDU :: MPFloat -> MPFloat -> BoundsCEDU MPFloat -mulCEDU = binaryCEDU (P.*) - -divCEDU :: MPFloat -> MPFloat -> BoundsCEDU MPFloat -divCEDU x y - | y P.== (P.fromInteger 0) = getBoundsCEDU MPLow.Bottom - | otherwise = binaryCEDU (P./) x y - -recipCEDU :: MPFloat -> BoundsCEDU MPFloat -recipCEDU = unaryCEDU P.recip - -{- special constants and functions -} - -piCEDU :: Precision -> BoundsCEDU MPFloat -piCEDU pp = - getBoundsCEDU $ MPLow.piA (p2cdarPrec pp) - -cosCEDU :: MPFloat -> BoundsCEDU MPFloat -cosCEDU = unaryPrecCEDU 0 MPLow.cosA - -sinCEDU :: MPFloat -> BoundsCEDU MPFloat -sinCEDU = unaryPrecCEDU 0 MPLow.sinA - -sqrtCEDU :: MPFloat -> BoundsCEDU MPFloat -sqrtCEDU = unaryCEDU MPLow.sqrtA - -expCEDU :: MPFloat -> BoundsCEDU MPFloat -expCEDU = unaryCEDU MPLow.expA - -logCEDU :: MPFloat -> BoundsCEDU MPFloat -logCEDU = unaryCEDU MPLow.logA - -{- auxiliary functions to automatically determine result precision from operand precisions -} - -binaryCEDU :: - (MPFloat -> MPFloat -> MPFloat) -> - (MPFloat -> MPFloat -> BoundsCEDU MPFloat) -binaryCEDU op x y = - getBoundsCEDU $ op x y - -unaryCEDU :: - (MPFloat -> MPFloat) -> - (MPFloat -> BoundsCEDU MPFloat) -unaryCEDU op x = - getBoundsCEDU $ op x - -unaryPrecCEDU :: - Integer -> - (MPLow.Precision -> MPFloat -> MPFloat) -> - (MPFloat -> BoundsCEDU MPFloat) -unaryPrecCEDU addPrec op x@(MPLow.Approx mb _ _ s) = - getBoundsCEDU $ op ((-s P.+ mb) P.+ (int addPrec)) x -unaryPrecCEDU addPrec op MPLow.Bottom = - getBoundsCEDU $ op ((int $ integer defaultPrecision) P.+ (int addPrec)) MPLow.Bottom - diff --git a/src-cdar/AERN2/MP/Float/Conversions.hs b/src-cdar/AERN2/MP/Float/Conversions.hs deleted file mode 100644 --- a/src-cdar/AERN2/MP/Float/Conversions.hs +++ /dev/null @@ -1,152 +0,0 @@ -{-| - Module : AERN2.MP.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 --} - -module AERN2.MP.Float.Conversions - ( - -- * MPFloat to other types (see also instances) - toDouble - -- * MPFloat constructors (see also instances) - , CanBeMPFloat, mpFloat - , fromIntegerCEDU - , fromRationalCEDU - -- * comparisons and constants (see also instances) - , zero, one, two - , nan, infinity - ) -where - -import MixedTypesNumPrelude -import qualified Prelude as P - -import Data.Ratio -import Data.Convertible - --- import AERN2.Norm -import AERN2.MP.Precision - -import qualified Data.CDAR as MPLow - -import AERN2.MP.Float.Auxi -import AERN2.MP.Float.Type -import AERN2.MP.Float.Arithmetic - - -{- conversions to MPFloat -} - -type CanBeMPFloat t = ConvertibleExactly t MPFloat -mpFloat :: (CanBeMPFloat t) => t -> MPFloat -mpFloat = convertExactly - -instance ConvertibleExactly Integer MPFloat where - safeConvertExactly = - Right . P.fromInteger - -instance ConvertibleExactly Int MPFloat where - safeConvertExactly = safeConvertExactly . integer - -fromIntegerCEDU :: Precision -> Integer -> BoundsCEDU MPFloat -fromIntegerCEDU pp = - setPrecisionCEDU pp . P.fromInteger - -fromRationalCEDU :: Precision -> Rational -> BoundsCEDU MPFloat -fromRationalCEDU pp = - setPrecisionCEDU pp . (MPLow.toApprox (p2cdarPrec pp)) - -{- conversions from MPFloat -} - -instance ConvertibleExactly MPFloat Rational where - safeConvertExactly = Right . P.toRational - -toDouble :: MPFloat -> Double -toDouble = P.fromRational . rational - -instance Convertible MPFloat Double where - safeConvert x - | isFinite dbl = Right dbl - | otherwise = convError "conversion to double: out of bounds" x - where - dbl = toDouble x - - -instance CanRound MPFloat where - properFraction x = (n,f) - where - r = rational x - n = (numerator r) `P.quot` (denominator r) - f = ceduCentre $ x `subCEDU` (P.fromInteger n) - -{- 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 - -{- constants -} - -zero, one, two :: MPFloat -zero = mpFloat 0 -one = mpFloat 1 -two = mpFloat 2 - -nan, infinity :: MPFloat -nan = MPLow.Bottom -infinity = nan - -itisNaN :: MPFloat -> Bool -itisNaN MPLow.Bottom = True -itisNaN _ = False - -instance CanTestFinite MPFloat where - isInfinite = itisNaN - isFinite = not . itisNaN - -instance CanTestNaN MPFloat where - isNaN = itisNaN diff --git a/src-cdar/AERN2/MP/Float/Type.hs b/src-cdar/AERN2/MP/Float/Type.hs deleted file mode 100644 --- a/src-cdar/AERN2/MP/Float/Type.hs +++ /dev/null @@ -1,79 +0,0 @@ -{-# LANGUAGE DeriveGeneric, DeriveDataTypeable, StandaloneDeriving #-} -{-| - Module : AERN2.MP.Float.Type - Description : Arbitrary precision floating point numbers (via cdar) - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mikkonecny@gmail.com - Stability : experimental - Portability : portable - - Arbitrary precision floating-point numbers, re-using CDAR Approx type. --} - -module AERN2.MP.Float.Type - ( - -- * MPFloat numbers and their basic operations - MPFloat - , showMPFloat - , setPrecisionCEDU - , p2cdarPrec - , getBoundsCEDU - ) -where - -import MixedTypesNumPrelude -import qualified Prelude as P - --- import Data.Bits (unsafeShiftL) -import Data.Typeable - -import AERN2.Norm -import AERN2.MP.Precision -import AERN2.MP.Float.Auxi - -import qualified Data.CDAR as MPLow - -{-| Multiple-precision floating-point type based on CDAR.Approx with 0 radius. -} -type MPFloat = MPLow.Approx - -showMPFloat :: MPFloat -> String -showMPFloat x = MPLow.showA x - -deriving instance (Typeable MPFloat) - -p2cdarPrec :: Precision -> MPLow.Precision -p2cdarPrec = P.fromInteger . integer - -getBoundsCEDU :: MPFloat -> BoundsCEDU MPFloat -getBoundsCEDU (MPLow.Approx mb m e s) = - BoundsCEDU - (MPLow.Approx mb m 0 s) (MPLow.approxMB eb_mb e 0 s) - (MPLow.Approx mb (m-e) 0 s) (MPLow.Approx mb (m+e) 0 s) -getBoundsCEDU MPLow.Bottom = - BoundsCEDU - MPLow.Bottom MPLow.Bottom MPLow.Bottom MPLow.Bottom - -{-| The bit-size bound for the error bound in CEDU -} -eb_prec :: Precision -eb_prec = prec 63 - -{-| The bit-size bound for the error bound in CEDU -} -eb_mb :: Int -eb_mb = int $ integer eb_prec - -instance HasPrecision MPFloat where - getPrecision (MPLow.Approx mb _ _ _) = prec (P.toInteger $ mb) - getPrecision MPLow.Bottom = error "illegal MPFloat (Bottom)" - - -instance CanSetPrecision MPFloat where - setPrecision p = ceduCentre . setPrecisionCEDU p - -setPrecisionCEDU :: Precision -> MPFloat -> BoundsCEDU MPFloat -setPrecisionCEDU pp = getBoundsCEDU . MPLow.enforceMB . MPLow.setMB (p2cdarPrec pp) - -instance HasNorm MPFloat where - getNormLog (MPLow.Approx _ m _ s) = (getNormLog m) + (integer s) - getNormLog MPLow.Bottom = error "getNormLog undefined for Bottom" diff --git a/src-rounded/AERN2/MP/Float/Arithmetic.hs b/src-rounded/AERN2/MP/Float/Arithmetic.hs deleted file mode 100644 --- a/src-rounded/AERN2/MP/Float/Arithmetic.hs +++ /dev/null @@ -1,117 +0,0 @@ -{-| - Module : AERN2.MP.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. --} - -module AERN2.MP.Float.Arithmetic - ( - -- * MPFloat basic arithmetic - addCEDU, subCEDU - , mulCEDU, divCEDU, recipCEDU - -- * MPFloat selected constants and operations - , piCEDU - , cosCEDU, sinCEDU - , sqrtCEDU, expCEDU, logCEDU - -- * auxiliary functions - , constCEDU, unaryCEDU, binaryCEDU - ) -where - -import MixedTypesNumPrelude -import qualified Prelude as P - -import AERN2.MP.Precision - -import qualified AERN2.MP.Float.RoundedAdaptor as MPLow - -import AERN2.MP.Float.Auxi -import AERN2.MP.Float.Type - -one :: MPFloat -one = MPLow.one - -{- common functions -} - -instance CanNeg MPFloat where - negate = ceduUp . unaryCEDU MPLow.neg - -instance CanAbs MPFloat where - abs x - | x P.< MPLow.zero = negate x - | otherwise = x - - -addCEDU :: MPFloat -> MPFloat -> BoundsCEDU MPFloat -addCEDU = binaryCEDU MPLow.add - -subCEDU :: MPFloat -> MPFloat -> BoundsCEDU MPFloat -subCEDU = binaryCEDU MPLow.sub - -mulCEDU :: MPFloat -> MPFloat -> BoundsCEDU MPFloat -mulCEDU = binaryCEDU MPLow.mul - -divCEDU :: MPFloat -> MPFloat -> BoundsCEDU MPFloat -divCEDU = binaryCEDU MPLow.div - -recipCEDU :: MPFloat -> BoundsCEDU MPFloat -recipCEDU x = divCEDU one x - -{- special constants and functions -} - -piCEDU :: Precision -> BoundsCEDU MPFloat -piCEDU pp = - constCEDU MPLow.pi (p2mpfrPrec pp) - -cosCEDU :: MPFloat -> BoundsCEDU MPFloat -cosCEDU = unaryCEDU MPLow.cos - -sinCEDU :: MPFloat -> BoundsCEDU MPFloat -sinCEDU = unaryCEDU MPLow.sin - -sqrtCEDU :: MPFloat -> BoundsCEDU MPFloat -sqrtCEDU = unaryCEDU MPLow.sqrt - -expCEDU :: MPFloat -> BoundsCEDU MPFloat -expCEDU = unaryCEDU MPLow.exp - -logCEDU :: MPFloat -> BoundsCEDU MPFloat -logCEDU = unaryCEDU MPLow.log - -{- auxiliary functions to automatically determine result precision from operand precisions -} - -binaryCEDU :: - (MPLow.RoundMode -> MPLow.Precision -> MPFloat -> MPFloat -> MPFloat) -> - MPFloat -> MPFloat -> BoundsCEDU MPFloat -binaryCEDU op x y = - getCEDU d u - where - d = op MPLow.Down p x y - u = op MPLow.Up p x y - p = p2mpfrPrec $ (getPrecision x) `max` (getPrecision y) - -unaryCEDU :: - (MPLow.RoundMode -> MPLow.Precision -> MPFloat -> MPFloat) -> - MPFloat -> BoundsCEDU MPFloat -unaryCEDU op x = - getCEDU d u - where - d = op MPLow.Down p x - u = op MPLow.Up p x - p = p2mpfrPrec $ getPrecision x - -constCEDU :: - (MPLow.RoundMode -> MPLow.Precision -> MPFloat) -> - MPLow.Precision -> BoundsCEDU MPFloat -constCEDU op p = - getCEDU d u - where - d = op MPLow.Down p - u = op MPLow.Up p diff --git a/src-rounded/AERN2/MP/Float/Conversions.hs b/src-rounded/AERN2/MP/Float/Conversions.hs deleted file mode 100644 --- a/src-rounded/AERN2/MP/Float/Conversions.hs +++ /dev/null @@ -1,178 +0,0 @@ -{-| - Module : AERN2.MP.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 --} - -module AERN2.MP.Float.Conversions - ( - -- * MPFloat to other types (see also instances) - toDouble - -- * MPFloat constructors (see also instances) - , CanBeMPFloat, mpFloat - , fromIntegerCEDU - , fromRationalCEDU - -- * comparisons and constants (see also instances) - , zero, one, two - , nan, infinity - ) -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.Auxi -import AERN2.MP.Float.Type -import AERN2.MP.Float.Arithmetic - -import qualified AERN2.MP.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 - -{- conversions to MPFloat -} - -type CanBeMPFloat t = ConvertibleExactly t MPFloat -mpFloat :: (CanBeMPFloat t) => t -> MPFloat -mpFloat = convertExactly - -instance ConvertibleExactly Integer MPFloat where - safeConvertExactly n = - findExact $ map (flip fromIntegerCEDU n) $ standardPrecisions initPrec - where - initPrec = - case getNormLog n of - NormBits b -> prec (b + 8) - _ -> prec 8 - findExact [] = - convError "integer too high to represent exactly" n - findExact (cedu : rest) - | ceduErr cedu P.> zero = findExact rest - | otherwise = Right (ceduCentre cedu) - -instance ConvertibleExactly Int MPFloat where - safeConvertExactly = safeConvertExactly . integer - -fromIntegerCEDU :: Precision -> Integer -> BoundsCEDU MPFloat -fromIntegerCEDU pp n = - constCEDU (\r p -> MPLow.fromIntegerA r p n) (p2mpfrPrec pp) - -fromRationalCEDU :: Precision -> Rational -> BoundsCEDU MPFloat -fromRationalCEDU pp q = - constCEDU (\r p -> mpFromRationalA r p q) (p2mpfrPrec pp) - -{- conversions from MPFloat -} - -instance ConvertibleExactly MPFloat Rational where - safeConvertExactly = Right . mpToRational - -toDouble :: MPFloat -> Double -toDouble = mpToDouble MPLow.Up - -instance Convertible MPFloat Double where - safeConvert x - | isFinite dbl = Right dbl - | otherwise = convError "conversion to double: out of bounds" x - where - dbl = toDouble x - -instance CanRound MPFloat where - properFraction x = (n,f) - where - r = rational x - n = (numerator r) `P.quot` (denominator r) - f = ceduCentre $ x `subCEDU` (mpFloat n) - -{- 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 - -{- constants -} - -zero, one, two :: MPFloat -zero = MPLow.zero -one = MPLow.one -two = MPLow.add MPLow.Up (MPLow.getPrec one) one one - -nan, infinity :: MPFloat -nan = ceduCentre $ divCEDU zero zero -infinity = ceduCentre $ divCEDU one zero - -itisNaN :: MPFloat -> Bool -itisNaN x = not $ x P.== x - -itisInfinite :: MPFloat -> Bool -itisInfinite x = - ceduCentre (mulCEDU x two) P.== x - && - x P./= zero - -instance CanTestFinite MPFloat where - isInfinite = itisInfinite - isFinite x = not (itisInfinite x || itisNaN x) - -instance CanTestNaN MPFloat where - isNaN = itisNaN - - diff --git a/src-rounded/AERN2/MP/Float/RoundedAdaptor.hs b/src-rounded/AERN2/MP/Float/RoundedAdaptor.hs deleted file mode 100644 --- a/src-rounded/AERN2/MP/Float/RoundedAdaptor.hs +++ /dev/null @@ -1,84 +0,0 @@ -{-# 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.RoundedAdaptor -( - module AERN2.MP.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 = 53 - -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-rounded/AERN2/MP/Float/Type.hs b/src-rounded/AERN2/MP/Float/Type.hs deleted file mode 100644 --- a/src-rounded/AERN2/MP/Float/Type.hs +++ /dev/null @@ -1,76 +0,0 @@ -{-# LANGUAGE DeriveGeneric, DeriveDataTypeable, StandaloneDeriving #-} -{-| - Module : AERN2.MP.Float.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.Type - ( - -- * MPFloat numbers and their basic operations - MPFloat - , showMPFloat - , setPrecisionCEDU - , getCEDU - , p2mpfrPrec - ) -where - -import MixedTypesNumPrelude -import qualified Prelude as P - -import AERN2.Norm -import AERN2.MP.Precision -import AERN2.MP.Float.Auxi - -import qualified AERN2.MP.Float.RoundedAdaptor as MPLow -import Data.Typeable - -{-| Multiple-precision floating-point type based on MPFR via rounded. -} -type MPFloat = MPLow.Rounded - -showMPFloat :: MPFloat -> String -showMPFloat = show - -deriving instance (Typeable MPFloat) - -p2mpfrPrec :: Precision -> MPLow.Precision -p2mpfrPrec = P.fromInteger . integer - -getCEDU :: MPFloat -> MPFloat -> BoundsCEDU MPFloat -getCEDU d u = BoundsCEDU c e d u - where - c = u - e = MPLow.sub MPLow.Up (MPLow.getPrec c) u d - -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) - -setPrecisionCEDU :: Precision -> MPFloat -> BoundsCEDU MPFloat -setPrecisionCEDU p x = - getCEDU d u - where - d = setPrecisionDown p x - u = setPrecisionUp p x - -instance HasNorm MPFloat where - getNormLog x - | x P.== MPLow.zero = NormZero - | otherwise = NormBits (P.toInteger $ MPLow.getExp x) - diff --git a/src/AERN2/Limit.hs b/src/AERN2/Limit.hs new file mode 100644 --- /dev/null +++ b/src/AERN2/Limit.hs @@ -0,0 +1,26 @@ +{-| + Module : AERN2.Limit + Description : limit operation + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Generic limit operation +-} +module AERN2.Limit where + +-- import MixedTypesNumPrelude + +-- import qualified Numeric.CollectErrors as CN + +--------- +-- limit +--------- + +class HasLimits ix s where + type LimitType ix s + limit :: (ix -> s) -> LimitType ix s + diff --git a/src/AERN2/MP.hs b/src/AERN2/MP.hs --- a/src/AERN2/MP.hs +++ b/src/AERN2/MP.hs @@ -18,6 +18,8 @@ , module AERN2.MP.ErrorBound , module AERN2.MP.Enclosure , MPBall(..), CanBeMPBall, mpBall, CanBeMPBallP, mpBallP +, reducePrecionIfInaccurate +, giveUpIfVeryInaccurate ) where diff --git a/src/AERN2/MP/Accuracy.hs b/src/AERN2/MP/Accuracy.hs --- a/src/AERN2/MP/Accuracy.hs +++ b/src/AERN2/MP/Accuracy.hs @@ -8,20 +8,19 @@ Stability : experimental Portability : portable - A type for roughly measuring the accuracy of an enclosure. + A type for measuring the accuracy of an enclosing set, + roughly corresponding to the maximum absolute error in some distance metric + approximately measured in bits. -} module AERN2.MP.Accuracy (Accuracy(NoInformation, Exact), bits, fromAccuracy, HasAccuracy(..), - HasAccuracyGuide(..), CanSetAccuracyGuide(..), adjustAccuracyGuide, - getFiniteAccuracy, ac2prec, - CanReduceSizeUsingAccuracyGuide(..), - specCanReduceSizeUsingAccuracyGuide, iterateUntilAccurate, - convergentList2CauchySeq, - seqByPrecision2CauchySeq, + convergentList2seqByAccuracy, + seqByPrecision2seqByAccuracy, setPrecisionAtLeastAccuracy, + ShowWithAccuracy(..), HasApproximate(..)) where @@ -32,8 +31,8 @@ import Data.Complex -import Test.Hspec -import Test.QuickCheck +-- import Test.Hspec +import Test.QuickCheck ( Arbitrary(arbitrary), frequency ) import AERN2.Norm import AERN2.MP.Precision @@ -174,52 +173,48 @@ class HasAccuracy a where getAccuracy :: a -> Accuracy + {-| Return accuracy, except when the element is Exact, return its nominal Precision dressed as Accuracy. + This function is useful when we have a convergent sequence where all elements happen to be + actually equal to the limit and we need the property that the sequence elements keep improving. + -} + getFiniteAccuracy :: a -> Accuracy + default getFiniteAccuracy :: (HasPrecision a) => a -> Accuracy + getFiniteAccuracy b = + case getAccuracy b of + Exact -> bits $ getPrecision b + a -> a -instance (HasAccuracy a, SuitableForCE es) => HasAccuracy (CollectErrors es a) where +instance (HasAccuracy a, CanBeErrors es) => HasAccuracy (CollectErrors es a) where getAccuracy (CollectErrors ma es) = case ma of Just a | not (hasCertainError es) -> getAccuracy a _ -> NoInformation + getFiniteAccuracy (CollectErrors ma es) = + case ma of + Just a | not (hasCertainError es) -> getFiniteAccuracy a + _ -> NoInformation -instance HasAccuracy Int where getAccuracy _ = Exact -instance HasAccuracy Integer where getAccuracy _ = Exact -instance HasAccuracy Rational where getAccuracy _ = Exact -instance HasAccuracy Bool where getAccuracy _ = Exact +instance HasAccuracy Int where getAccuracy _ = Exact; getFiniteAccuracy _ = NoInformation +instance HasAccuracy Integer where getAccuracy _ = Exact; getFiniteAccuracy _ = NoInformation +instance HasAccuracy Rational where getAccuracy _ = Exact; getFiniteAccuracy _ = NoInformation +instance HasAccuracy Bool where getAccuracy _ = Exact; getFiniteAccuracy _ = NoInformation +instance HasAccuracy Kleenean where getAccuracy _ = Exact; getFiniteAccuracy _ = NoInformation instance HasAccuracy t => HasAccuracy (Complex t) where getAccuracy (a :+ i) = (getAccuracy a) `min` (getAccuracy i) - + getFiniteAccuracy (a :+ i) = + (getFiniteAccuracy a) `min` (getFiniteAccuracy i) + instance HasAccuracy t => HasAccuracy [t] where getAccuracy xs = foldl min Exact $ map getAccuracy xs + getFiniteAccuracy xs = foldl min Exact $ map getFiniteAccuracy xs instance HasAccuracy t => HasAccuracy (Maybe t) where getAccuracy (Just x) = getAccuracy x getAccuracy _ = NoInformation - -class HasAccuracyGuide a where - getAccuracyGuide :: a -> Accuracy - -class HasAccuracyGuide a => CanSetAccuracyGuide a where - setAccuracyGuide :: Accuracy -> a -> a - -adjustAccuracyGuide :: - (CanSetAccuracyGuide a) => - (Accuracy -> Accuracy) -> a -> a -adjustAccuracyGuide adj_acG a = - setAccuracyGuide (adj_acG (getAccuracyGuide a)) a - -{-| Return accuracy, except when the element is Exact, return its nominal Precision dressed as Accuracy. - This function is useful when we have a convergent sequence where all elements happen to be - actually equal to the limit and we need the property that the sequence elements keep improving. --} -getFiniteAccuracy :: - (HasAccuracy t, HasPrecision t) => - t -> Accuracy -getFiniteAccuracy b = - case getAccuracy b of - Exact -> bits $ getPrecision b - a -> a + getFiniteAccuracy (Just x) = getFiniteAccuracy x + getFiniteAccuracy _ = NoInformation iterateUntilAccurate :: (HasAccuracy t) => @@ -238,22 +233,22 @@ Bits b -> prec (max 2 $ b + 50) _ -> prec 100 -seqByPrecision2CauchySeq :: +seqByPrecision2seqByAccuracy :: (HasAccuracy t) => (Precision -> t) -> (Accuracy -> t) -seqByPrecision2CauchySeq seqByPrecision ac = - convergentList2CauchySeq list ac +seqByPrecision2seqByAccuracy seqByPrecision ac = + convergentList2seqByAccuracy list ac where list = map seqByPrecision $ dropWhile (lowPrec ac) (standardPrecisions (ac2prec ac)) lowPrec Exact _ = False lowPrec _ p = bits p < ac -convergentList2CauchySeq :: (HasAccuracy t) => [t] -> (Accuracy -> t) -convergentList2CauchySeq list ac = findAccurate list +convergentList2seqByAccuracy :: (HasAccuracy t) => [t] -> (Accuracy -> t) +convergentList2seqByAccuracy list ac = findAccurate list where findAccurate [] = - error "convergentList2CauchySeq: the sequence either converges too slowly or it does not converge" + error "convergentList2seqByAccuracy: the sequence either converges too slowly or it does not converge" findAccurate (b : rest) | getAccuracy b >= ac = b | otherwise = findAccurate rest @@ -263,7 +258,7 @@ it is at least as high as the supplied accuracy (assuming the accuracy is finite). -} -setPrecisionAtLeastAccuracy :: (CanSetPrecision t) => Accuracy -> t -> t +setPrecisionAtLeastAccuracy :: (HasPrecision t, CanSetPrecision t) => Accuracy -> t -> t setPrecisionAtLeastAccuracy acc b | p_b < p_acc = setPrecision p_acc b | otherwise = b @@ -275,22 +270,8 @@ _ -> prec $ max 2 (fromAccuracy acc) p_b = getPrecision b - -class CanReduceSizeUsingAccuracyGuide t where - reduceSizeUsingAccuracyGuide :: Accuracy -> t -> t - -specCanReduceSizeUsingAccuracyGuide :: - ( CanReduceSizeUsingAccuracyGuide t - , HasEqCertainly t t - , Arbitrary t, Show t) - => - (T t) -> Spec -specCanReduceSizeUsingAccuracyGuide (T tName :: T t) = - describe ("CanReduceSizeUsingAccuracyGuide " ++ tName) $ do - it "is safe" $ - property $ - \ (t :: t) (ac :: Accuracy) -> - reduceSizeUsingAccuracyGuide ac t ?==? t +class ShowWithAccuracy t where + showWithAccuracy :: Accuracy -> t -> String {-| An unsafe approximation of an enclosure or exact value, useful mainly for showing something brief and readable to humans. 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,3 +1,4 @@ +{-# OPTIONS_GHC -Wno-orphans #-} {-| Module : AERN2.MP.Ball Description : Arbitrary precision ball arithmetic @@ -19,12 +20,17 @@ , module AERN2.MP.ErrorBound , module AERN2.MP.Enclosure -- * The Ball type - , MPBall(..), CanBeMPBall, mpBall, CanBeMPBallP, mpBallP + , MPBall(..), CanBeMPBall, mpBall, cnMPBall, CanBeMPBallP, mpBallP, cnMPBallP , reducePrecionIfInaccurate + , giveUpIfVeryInaccurate + {- -- * Ball construction/extraction functions - -- , endpointsMP, fromEndpointsMP + , endpointsMP, fromEndpointsMP + -} -- * Ball operations (see also instances) , piBallP + , intersectCNMPBall + , hullMPBall -- * Helpers for constructing ball functions , byEndpointsMP , fromApproxWithLipschitz @@ -45,6 +51,7 @@ import AERN2.MP.Ball.Conversions () import AERN2.MP.Ball.Comparisons import AERN2.MP.Ball.Field () +import AERN2.MP.Ball.Limit () import AERN2.MP.Ball.Elementary import AERN2.MP.Ball.PreludeOps () diff --git a/src/AERN2/MP/Ball/Comparisons.hs b/src/AERN2/MP/Ball/Comparisons.hs --- a/src/AERN2/MP/Ball/Comparisons.hs +++ b/src/AERN2/MP/Ball/Comparisons.hs @@ -1,3 +1,4 @@ +{-# OPTIONS_GHC -Wno-orphans #-} {-| Module : AERN2.MP.Ball.Comparisons Description : Comparisons of arbitrary precision dyadic balls @@ -18,13 +19,19 @@ , reducePrecionIfInaccurate -- * Helpers for constructing ball functions , byEndpointsMP + -- * intersection and hull + , intersectCNMPBall + , hullMPBall ) where import MixedTypesNumPrelude -- import qualified Prelude as P +import qualified Control.CollectErrors as CE import Control.CollectErrors + ( CollectErrors(getMaybeValue), CanBeErrors ) +import qualified Numeric.CollectErrors as CN import AERN2.Norm import AERN2.MP.Dyadic (Dyadic) @@ -38,176 +45,168 @@ {- comparisons -} instance HasEqAsymmetric MPBall MPBall where - type EqCompareType MPBall MPBall = Maybe Bool + type EqCompareType MPBall MPBall = Kleenean b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 instance HasEqAsymmetric MPBall Integer where - type EqCompareType MPBall Integer = Maybe Bool + type EqCompareType MPBall Integer = Kleenean b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 instance HasEqAsymmetric Integer MPBall where - type EqCompareType Integer MPBall = Maybe Bool + type EqCompareType Integer MPBall = Kleenean b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 instance HasEqAsymmetric MPBall Int where - type EqCompareType MPBall Int = Maybe Bool + type EqCompareType MPBall Int = Kleenean b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 instance HasEqAsymmetric Int MPBall where - type EqCompareType Int MPBall = Maybe Bool + type EqCompareType Int MPBall = Kleenean b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 instance HasEqAsymmetric MPBall Rational where - type EqCompareType MPBall Rational = Maybe Bool + type EqCompareType MPBall Rational = Kleenean b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 instance HasEqAsymmetric Rational MPBall where - type EqCompareType Rational MPBall = Maybe Bool + type EqCompareType Rational MPBall = Kleenean b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 instance HasEqAsymmetric MPBall Dyadic where - type EqCompareType MPBall Dyadic = Maybe Bool + type EqCompareType MPBall Dyadic = Kleenean b1 `equalTo` b2 = b1 >= b2 && b1 <= b2 instance HasEqAsymmetric Dyadic MPBall where - type EqCompareType Dyadic MPBall = Maybe Bool + type EqCompareType Dyadic MPBall = Kleenean 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) + , IsBool (EqCompareType MPBall b) + , CanBeErrors es) => - HasEqAsymmetric MPBall (CollectErrors es b) + HasEqAsymmetric MPBall (CollectErrors es b) where - type EqCompareType MPBall (CollectErrors es b) = - EnsureCE es (EqCompareType MPBall b) - equalTo = lift2TLCE equalTo + type EqCompareType MPBall (CollectErrors es b) = + CollectErrors es (EqCompareType MPBall b) + equalTo = CE.liftT1 equalTo instance (HasEqAsymmetric a MPBall - , CanEnsureCE es a - , CanEnsureCE es (EqCompareType a MPBall) - , IsBool (EnsureCE es (EqCompareType a MPBall)) - , SuitableForCE es) + , IsBool (EqCompareType a MPBall) + , CanBeErrors es) => HasEqAsymmetric (CollectErrors es a) MPBall where - type EqCompareType (CollectErrors es a) MPBall = - EnsureCE es (EqCompareType a MPBall) - equalTo = lift2TCE equalTo + type EqCompareType (CollectErrors es a) MPBall = + CollectErrors es (EqCompareType a MPBall) + equalTo = CE.lift1T equalTo instance HasOrderAsymmetric MPBall MPBall where - type OrderCompareType MPBall MPBall = Maybe Bool + type OrderCompareType MPBall MPBall = Kleenean lessThan b1 b2 - | r1 < l2 = Just True - | r2 <= l1 = Just False - | otherwise = Nothing + | r1 < l2 = CertainTrue + | r2 <= l1 = CertainFalse + | otherwise = TrueOrFalse where - (l1, r1) = endpointsMP b1 - (l2, r2) = endpointsMP b2 + (l1, r1) = endpoints b1 + (l2, r2) = endpoints b2 leq b1 b2 - | r1 <= l2 = Just True - | r2 < l1 = Just False - | otherwise = Nothing + | r1 <= l2 = CertainTrue + | r2 < l1 = CertainFalse + | otherwise = TrueOrFalse where - (l1, r1) = endpointsMP b1 - (l2, r2) = endpointsMP b2 + (l1, r1) = endpoints b1 + (l2, r2) = endpoints b2 instance HasOrderAsymmetric Integer MPBall where - type OrderCompareType Integer MPBall = Maybe Bool + type OrderCompareType Integer MPBall = Kleenean lessThan = convertFirst lessThan leq = convertFirst leq instance HasOrderAsymmetric MPBall Integer where - type OrderCompareType MPBall Integer = Maybe Bool + type OrderCompareType MPBall Integer = Kleenean lessThan = convertSecond lessThan leq = convertSecond leq instance HasOrderAsymmetric Int MPBall where - type OrderCompareType Int MPBall = Maybe Bool + type OrderCompareType Int MPBall = Kleenean lessThan = convertFirst lessThan leq = convertFirst leq instance HasOrderAsymmetric MPBall Int where - type OrderCompareType MPBall Int = Maybe Bool + type OrderCompareType MPBall Int = Kleenean lessThan = convertSecond lessThan leq = convertSecond leq instance HasOrderAsymmetric Dyadic MPBall where - type OrderCompareType Dyadic MPBall = Maybe Bool + type OrderCompareType Dyadic MPBall = Kleenean lessThan = convertFirst lessThan leq = convertFirst leq instance HasOrderAsymmetric MPBall Dyadic where - type OrderCompareType MPBall Dyadic = Maybe Bool + type OrderCompareType MPBall Dyadic = Kleenean lessThan = convertSecond lessThan leq = convertSecond leq instance HasOrderAsymmetric MPBall Rational where - type OrderCompareType MPBall Rational = Maybe Bool + type OrderCompareType MPBall Rational = Kleenean lessThan b1 q2 - | r1 < l2 = Just True - | r2 <= l1 = Just False - | otherwise = Nothing + | r1 < l2 = CertainTrue + | r2 <= l1 = CertainFalse + | otherwise = TrueOrFalse where - (l1, r1) = endpointsMP b1 + (l1, r1) = endpoints b1 l2 = q2 r2 = q2 leq b1 q2 - | r1 <= l2 = Just True - | r2 < l1 = Just False - | otherwise = Nothing + | r1 <= l2 = CertainTrue + | r2 < l1 = CertainFalse + | otherwise = TrueOrFalse where - (l1, r1) = endpointsMP b1 + (l1, r1) = endpoints b1 l2 = q2 r2 = q2 instance HasOrderAsymmetric Rational MPBall where - type OrderCompareType Rational MPBall = Maybe Bool + type OrderCompareType Rational MPBall = Kleenean lessThan q1 b2 - | r1 < l2 = Just True - | r2 <= l1 = Just False - | otherwise = Nothing + | r1 < l2 = CertainTrue + | r2 <= l1 = CertainFalse + | otherwise = TrueOrFalse where - (l2, r2) = endpointsMP b2 + (l2, r2) = endpoints b2 l1 = q1 r1 = q1 leq q1 b2 - | r1 <= l2 = Just True - | r2 < l1 = Just False - | otherwise = Nothing + | r1 <= l2 = CertainTrue + | r2 < l1 = CertainFalse + | otherwise = TrueOrFalse where - (l2, r2) = endpointsMP b2 + (l2, r2) = endpoints 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) + , IsBool (OrderCompareType MPBall b) + , CanBeErrors 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 + CollectErrors es (OrderCompareType MPBall b) + lessThan = CE.liftT1 lessThan + leq = CE.liftT1 leq + greaterThan = CE.liftT1 greaterThan + geq = CE.liftT1 geq instance (HasOrderAsymmetric a MPBall - , CanEnsureCE es a - , CanEnsureCE es (OrderCompareType a MPBall) - , IsBool (EnsureCE es (OrderCompareType a MPBall)) - , SuitableForCE es) + , IsBool (OrderCompareType a MPBall) + , CanBeErrors 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 + CollectErrors es (OrderCompareType a MPBall) + lessThan = CE.lift1T lessThan + leq = CE.lift1T leq + greaterThan = CE.lift1T greaterThan + geq = CE.lift1T geq instance CanTestZero MPBall instance CanTestPosNeg MPBall @@ -265,117 +264,119 @@ instance (CanMinMaxAsymmetric MPBall b - , CanEnsureCE es b - , CanEnsureCE es (MinMaxType MPBall b) - , SuitableForCE es) + , CanBeErrors es) => CanMinMaxAsymmetric MPBall (CollectErrors es b) where type MinMaxType MPBall (CollectErrors es b) = - EnsureCE es (MinMaxType MPBall b) - min = lift2TLCE min - max = lift2TLCE max + CollectErrors es (MinMaxType MPBall b) + min = CE.liftT1 min + max = CE.liftT1 max instance (CanMinMaxAsymmetric a MPBall - , CanEnsureCE es a - , CanEnsureCE es (MinMaxType a MPBall) - , SuitableForCE es) + , CanBeErrors es) => CanMinMaxAsymmetric (CollectErrors es a) MPBall where type MinMaxType (CollectErrors es a) MPBall = - EnsureCE es (MinMaxType a MPBall) - min = lift2TCE min - max = lift2TCE max + CollectErrors es (MinMaxType a MPBall) + min = CE.lift1T min + max = CE.lift1T 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 + | l > r = + CN.noValueNumErrorCertain $ CN.NumError $ "intersect: empty intersection: " ++ show a ++ "; " ++ show b + | otherwise = cn $ setPrecision p $ fromMPFloatEndpoints l r where - rL = max aL bL - rR = min aR bR - (aL,aR) = endpointsMP a - (bL,bR) = endpointsMP b + p = getPrecision a + l = max aL bL + r = min aR bR + (aL,aR) = endpoints a + (bL,bR) = endpoints b +intersectCNMPBall :: CN MPBall -> CN MPBall -> CN MPBall +intersectCNMPBall = intersect + -- case (fst $ ensureNoCN x, fst $ ensureNoCN y) of + -- (Nothing, Nothing) -> x + -- (Just _ , Nothing) -> x + -- (Nothing, Just _ ) -> y + -- (Just _ , Just _ ) -> lift2CE intersect x y + instance (CanIntersectAsymmetric MPBall b - , CanEnsureCE es b - , CanEnsureCE es (IntersectionType MPBall b) - , SuitableForCE es) + , CanBeErrors es) => CanIntersectAsymmetric MPBall (CollectErrors es b) where type IntersectionType MPBall (CollectErrors es b) = - EnsureCE es (IntersectionType MPBall b) - intersect = lift2TLCE intersect + CollectErrors es (IntersectionType MPBall b) + intersect = CE.liftT1 intersect instance (CanIntersectAsymmetric a MPBall - , CanEnsureCE es a - , CanEnsureCE es (IntersectionType a MPBall) - , SuitableForCE es) + , CanBeErrors es) => CanIntersectAsymmetric (CollectErrors es a) MPBall where type IntersectionType (CollectErrors es a) MPBall = - EnsureCE es (IntersectionType a MPBall) - intersect = lift2TCE intersect + CollectErrors es (IntersectionType a MPBall) + intersect = CE.lift1T intersect +{- hull -} + +hullMPBall :: MPBall -> MPBall -> MPBall +hullMPBall a b = + fromEndpoints rL rR + where + rL = min aL bL + rR = max aR bR + (aL,aR) = endpoints a + (bL,bR) = endpoints b + {- union -} instance CanUnionAsymmetric MPBall MPBall where union a b = - case getMaybeValueCN (a `intersect` b) of - Just _ -> prependErrorsCN [(ErrorPotential, err)] r - _ -> prependErrorsCN [(ErrorCertain, err)] r + case getMaybeValue (a `intersect` b) of + Just _ -> r + _ -> CN.prependErrorCertain 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 - + err = CN.NumError $ "union of enclosures: not enclosing the same value" + r = cn $ hullMPBall a b instance (CanUnionAsymmetric MPBall b - , CanEnsureCE es b - , CanEnsureCE es (UnionType MPBall b) - , SuitableForCE es) + , CanBeErrors es) => CanUnionAsymmetric MPBall (CollectErrors es b) where type UnionType MPBall (CollectErrors es b) = - EnsureCE es (UnionType MPBall b) - union = lift2TLCE union + CollectErrors es (UnionType MPBall b) + union = CE.liftT1 union instance (CanUnionAsymmetric a MPBall - , CanEnsureCE es a - , CanEnsureCE es (UnionType a MPBall) - , SuitableForCE es) + , CanBeErrors es) => CanUnionAsymmetric (CollectErrors es a) MPBall where type UnionType (CollectErrors es a) MPBall = - EnsureCE es (UnionType a MPBall) - union = lift2TCE union + CollectErrors es (UnionType a MPBall) + union = CE.lift1T union {-| - Computes an *increasing* ball fucntion @f@ from *exact* MPFR operations. + Compute an MPBall function from *exact* MPFloat operations on interval endpoints. + This works only for *non-decreasing* operations, eg addition, min, max. -} byEndpointsMP :: (MPFloat -> MPFloat -> MPFloat) -> (MPBall -> MPBall -> MPBall) byEndpointsMP op b1 b2 = - fromEndpointsMP (l1 `op` l2) (r1 `op` r2) + fromEndpoints (l1 `op` l2) (r1 `op` r2) where - (l1,r1) = endpointsMP b1 - (l2,r2) = endpointsMP b2 - -{- random generation -} + (l1,r1) = endpoints b1 + (l2,r2) = endpoints b2 diff --git a/src/AERN2/MP/Ball/Conversions.hs b/src/AERN2/MP/Ball/Conversions.hs --- a/src/AERN2/MP/Ball/Conversions.hs +++ b/src/AERN2/MP/Ball/Conversions.hs @@ -1,3 +1,4 @@ +{-# OPTIONS_GHC -Wno-orphans #-} {-| Module : AERN2.MP.Ball.Conversions Description : Conversions of arbitrary precision dyadic balls @@ -19,16 +20,19 @@ import MixedTypesNumPrelude -- import qualified Prelude as P +-- import qualified Numeric.CollectErrors as CN + import Data.Typeable -import Data.Convertible +-- 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 (MPFloat, mpFloat) -- import AERN2.MP.Float.Operators import AERN2.MP.Precision -- import qualified AERN2.MP.ErrorBound as EB -import AERN2.MP.ErrorBound (ErrorBound, errorBound) +-- import qualified AERN2.MP.ErrorBound as EB +import AERN2.MP.ErrorBound (ErrorBound, errorBound, CanBeErrorBound) import AERN2.MP.Ball.Type @@ -38,13 +42,13 @@ integerBounds b = (floor l, ceiling r) where - (l,r) = endpointsMP b + (l,r) = endpoints b instance Convertible MPBall ErrorBound where safeConvert b = Right (errorBound (max (abs l) (abs r))) where - (l,r) = endpointsMP b + (l,r) = endpoints b {--- constructing an exact ball ---} @@ -115,3 +119,33 @@ b = MPBall xFlt (xe + eUp) -- beware, precision may be too high relative to accuracy (MPBall xFlt xe) = mpBallP p x eUp = errorBound e + +{--- constructing a fat ball ---} + +instance (CanBeErrorBound t) => CanPlusMinus MPBall t where + plusMinus b e = updateRadius (+ (errorBound e)) b + +instance (CanBeErrorBound t) => CanPlusMinus (CN MPBall) t where + plusMinus b e = updateRadius (+ (errorBound e)) b + +instance (CanBeErrorBound t) => CanPlusMinus MPFloat t where + type PlusMinusType MPFloat t = MPBall + plusMinus b e = MPBall b (errorBound e) + +instance (CanBeErrorBound t) => CanPlusMinus Dyadic t where + type PlusMinusType Dyadic t = MPBall + plusMinus b e = MPBall (mpFloat b) (errorBound e) + +instance (CanBeErrorBound t) => CanPlusMinus Integer t where + type PlusMinusType Integer t = MPBall + plusMinus b e = MPBall (mpFloat b) (errorBound e) + +instance (CanBeErrorBound t) => CanPlusMinus Int t where + type PlusMinusType Int t = MPBall + plusMinus b e = MPBall (mpFloat b) (errorBound e) + +instance (CanBeErrorBound t) => CanPlusMinus Rational t where + type PlusMinusType Rational t = MPBall + plusMinus b e = (mpBallP p b) +- e + where + p = prec 100 diff --git a/src/AERN2/MP/Ball/Elementary.hs b/src/AERN2/MP/Ball/Elementary.hs --- a/src/AERN2/MP/Ball/Elementary.hs +++ b/src/AERN2/MP/Ball/Elementary.hs @@ -1,3 +1,5 @@ +{-# LANGUAGE TemplateHaskell #-} +{-# OPTIONS_GHC -Wno-orphans #-} {-| Module : AERN2.MP.Ball.Elementary Description : Elementary operations on arbitrary precision dyadic balls @@ -52,20 +54,22 @@ sinB :: Integer -> MPBall -> MPBall sinB i x = -- increasingPrecisionUntilNotImproving (fromApproxWithLipschitz MPFloat.sinDown MPFloat.sinUp lip) x + max (-1) $ min 1 $ fromApproxWithLipschitz MPFloat.sinCEDU lip x where lip | i == 0 = mpFloat 1 - | otherwise = snd $ endpointsMP $ abs $ cosB (i - 1) x + | otherwise = endpointR $ abs $ cosB (i - 1) x cosB :: Integer -> MPBall -> MPBall cosB i x = -- increasingPrecisionUntilNotImproving (fromApproxWithLipschitz MPFloat.cosDown MPFloat.cosUp lip) x + max (-1) $ min 1 $ fromApproxWithLipschitz MPFloat.cosCEDU lip x where lip | i == 0 = mpFloat 1 - | otherwise = snd $ endpointsMP $ abs $ sinB (i - 1) x + | otherwise = endpointR $ abs $ sinB (i - 1) x -- increasingPrecisionUntilNotImproving :: (MPBall -> MPBall) -> (MPBall -> MPBall) -- increasingPrecisionUntilNotImproving f x = @@ -89,49 +93,45 @@ exp = intervalFunctionByEndpointsUpDown MPFloat.expDown MPFloat.expUp instance CanLog MPBall where - type LogType MPBall = CN MPBall + type LogType MPBall = MPBall log x | x_!>! 1 = - cn $ setPrecision p $ ballFunctionUsingLipschitz log_ logLip x_ + setPrecision p $ ballFunctionUsingLipschitz log_ logLip x_ | x_!>! 0 = - cn $ setPrecision p $ intervalFunctionByEndpoints log_ x_ + setPrecision p $ intervalFunctionByEndpoints log_ x_ | x !>! 0 = - cn $ setPrecision p $ intervalFunctionByEndpoints log_ x - | x !<=! 0 = noValueNumErrorCertainCN err - | otherwise = noValueNumErrorPotentialCN err + intervalFunctionByEndpoints log_ x + | otherwise = err where p = getPrecision x x_ = reducePrecionIfInaccurate x - err = OutOfRange $ "log: argument must be > 0: " ++ show x + err = error $ "log: argument must be > 0: " ++ show x log_ (MPBall c e) = MPBall lc (e + (errorBound le)) where (lc, le) = ceduCentreErr $ MPFloat.logCEDU c - logLip y = errorBound $ (1/!y) + logLip y = errorBound $ (1/y) instance CanPow MPBall MPBall where - powNoCN b e = (~!) $ pow b e - pow = powUsingExpLog (mpBall 0) (mpBall 1) + pow = powUsingExpLog (mpBall 1) instance CanPow MPBall Dyadic where - powNoCN b e = (~!) $ pow b e - pow b e = powUsingExpLog (mpBall 0) (mpBall 1) b (mpBall e) + pow b e = powUsingExpLog (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) + pow b e = powUsingExpLog (mpBall 1) b (mpBallP (getPrecision b) e) instance CanSqrt MPBall where - type SqrtType MPBall = CN MPBall + type SqrtType MPBall = MPBall sqrt x - | x !>=! 0 = cn $ aux x - | x !=! 0 = aux x + | x ?>=? 0 = aux $ max 0 x + | otherwise = err where aux = intervalFunctionByEndpointsUpDown (\ e -> MPFloat.sqrtDown (P.max (mpFloat 0) e)) (\ e -> MPFloat.sqrtUp (P.max (mpFloat 0) e)) - err = OutOfRange $ "sqrt: argument must be >= 0: " ++ show x + err = error $ "sqrt: argument must be >= 0: " ++ show x {- generic methods for computing real functions from MPFR-approximations -} @@ -151,3 +151,15 @@ setPrecision (getPrecision xc) $ -- beware, some MPFloat functions may increase precision, eg sine and cosine (MPBall fxC (errorBound fxErr)) err = (errorBound lip) * xe + fxe + +$(declForTypes + [[t| Integer |], [t| Int |], [t| Rational |]] + (\ b -> [d| + + instance + CanPow $b MPBall + where + type PowType $b MPBall = MPBall + pow x e = pow (mpBallP (getPrecision e) x) e + |])) + diff --git a/src/AERN2/MP/Ball/Field.hs b/src/AERN2/MP/Ball/Field.hs --- a/src/AERN2/MP/Ball/Field.hs +++ b/src/AERN2/MP/Ball/Field.hs @@ -1,4 +1,7 @@ {-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE PartialTypeSignatures #-} +{-# OPTIONS_GHC -Wno-orphans #-} +{-# OPTIONS_GHC -Wno-partial-type-signatures #-} {-| Module : AERN2.MP.Ball.Field Description : Field operations on arbitrary precision dyadic balls @@ -18,7 +21,7 @@ import MixedTypesNumPrelude -- import qualified Prelude as P -import Control.CollectErrors +import qualified Numeric.CollectErrors as CN import AERN2.Normalize @@ -71,28 +74,20 @@ add = convertPFirst add instance - (CanAddAsymmetric MPBall b - , CanEnsureCE es b - , CanEnsureCE es (AddType MPBall b) - , SuitableForCE es) + (CanAddAsymmetric MPBall b) => - CanAddAsymmetric MPBall (CollectErrors es b) + CanAddAsymmetric MPBall (CN b) where - type AddType MPBall (CollectErrors es b) = - EnsureCE es (AddType MPBall b) - add = lift2TLCE add + type AddType MPBall (CN b) = CN (AddType MPBall b) + add = CN.liftT1 add instance - (CanAddAsymmetric a MPBall - , CanEnsureCE es a - , CanEnsureCE es (AddType a MPBall) - , SuitableForCE es) + (CanAddAsymmetric a MPBall) => - CanAddAsymmetric (CollectErrors es a) MPBall + CanAddAsymmetric (CN a) MPBall where - type AddType (CollectErrors es a) MPBall = - EnsureCE es (AddType a MPBall) - add = lift2TCE add + type AddType (CN a) MPBall = CN (AddType a MPBall) + add = CN.lift1T add {- subtraction -} @@ -111,28 +106,20 @@ instance CanSub Dyadic MPBall instance - (CanSub MPBall b - , CanEnsureCE es b - , CanEnsureCE es (SubType MPBall b) - , SuitableForCE es) + (CanSub MPBall b) => - CanSub MPBall (CollectErrors es b) + CanSub MPBall (CN b) where - type SubType MPBall (CollectErrors es b) = - EnsureCE es (SubType MPBall b) - sub = lift2TLCE sub + type SubType MPBall (CN b) = CN (SubType MPBall b) + sub = CN.liftT1 sub instance - (CanSub a MPBall - , CanEnsureCE es a - , CanEnsureCE es (SubType a MPBall) - , SuitableForCE es) + (CanSub a MPBall) => - CanSub (CollectErrors es a) MPBall + CanSub (CN a) MPBall where - type SubType (CollectErrors es a) MPBall = - EnsureCE es (SubType a MPBall) - sub = lift2TCE sub + type SubType (CN a) MPBall = CN (SubType a MPBall) + sub = CN.lift1T sub {- multiplication -} @@ -173,43 +160,27 @@ mul = convertPFirst mul instance - (CanMulAsymmetric MPBall b - , CanEnsureCE es b - , CanEnsureCE es (MulType MPBall b) - , SuitableForCE es) + (CanMulAsymmetric MPBall b) => - CanMulAsymmetric MPBall (CollectErrors es b) + CanMulAsymmetric MPBall (CN b) where - type MulType MPBall (CollectErrors es b) = - EnsureCE es (MulType MPBall b) - mul = lift2TLCE mul + type MulType MPBall (CN b) = CN (MulType MPBall b) + mul = CN.liftT1 mul instance - (CanMulAsymmetric a MPBall - , CanEnsureCE es a - , CanEnsureCE es (MulType a MPBall) - , SuitableForCE es) + (CanMulAsymmetric a MPBall) => - CanMulAsymmetric (CollectErrors es a) MPBall + CanMulAsymmetric (CN a) MPBall where - type MulType (CollectErrors es a) MPBall = - EnsureCE es (MulType a MPBall) - mul = lift2TCE mul + type MulType (CN a) MPBall = CN (MulType a MPBall) + mul = CN.lift1T 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 x12C err - | isCertainlyZero b2 = - noValueNumErrorCertainCN DivByZero - | otherwise = - noValueNumErrorPotentialCN DivByZero + type DivType MPBall MPBall = MPBall + divide (MPBall x1 e1) (MPBall x2 e2) = normalize $ MPBall x12C err where (x12C, e12) = MPFloat.ceduCentreErr $ MPFloat.divCEDU x1 x2 x12AbsUp = (abs x12C) +^ e12 @@ -241,113 +212,77 @@ [[t| Integer |], [t| Int |], [t| Dyadic |]] (\ t -> [d| instance CanDiv MPBall $t where - type DivType MPBall $t = CN MPBall + type DivType MPBall $t = MPBall divide = convertSecond divide - type DivTypeNoCN MPBall $t = MPBall - divideNoCN = convertSecond divideNoCN instance CanDiv $t MPBall where - type DivType $t MPBall = CN MPBall + type DivType $t MPBall = 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) + type DivType Dyadic Dyadic = MPBall divide a b = divide (mpBall a) (mpBall b) instance CanDiv MPBall Rational where - type DivTypeNoCN MPBall Rational = MPBall - divideNoCN = convertPSecond divideNoCN + type DivType MPBall Rational = MPBall divide = convertPSecond divide instance CanDiv Rational MPBall where - type DivTypeNoCN Rational MPBall = MPBall - divideNoCN = convertPFirst divideNoCN + type DivType Rational MPBall = MPBall divide = convertPFirst divide instance - (CanDiv MPBall b - , CanEnsureCE es b - , CanEnsureCE es (DivType MPBall b) - , CanEnsureCE es (DivTypeNoCN MPBall b) - , SuitableForCE es) + (CanDiv MPBall b, CanTestZero b) => - CanDiv MPBall (CollectErrors es b) + CanDiv MPBall (CN 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 + type DivType MPBall (CN b) = CN (DivType MPBall b) + divide a b = divide (cn a) b instance - (CanDiv a MPBall - , CanEnsureCE es a - , CanEnsureCE es (DivType a MPBall) - , CanEnsureCE es (DivTypeNoCN a MPBall) - , SuitableForCE es) + (CanDiv a MPBall) => - CanDiv (CollectErrors es a) MPBall + CanDiv (CN 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 + type DivType (CN a) MPBall = CN (DivType a MPBall) + divide a b = divide a (cn b) {- integer power -} instance CanPow MPBall Integer where - powNoCN b e = (~!) $ powUsingMulRecip (mpBall 1) b e - pow = powUsingMulRecip (mpBall 1) + pow = powUsingMulRecipCutNeg (mpBall 1) instance CanPow MPBall Int where - powNoCN b e = (~!) $ powUsingMulRecip (mpBall 1) b e - pow = powUsingMulRecip (mpBall 1) + pow = powUsingMulRecipCutNeg (mpBall 1) +powUsingMulRecipCutNeg :: _ => t -> t -> e -> DivType Integer t +powUsingMulRecipCutNeg one x e + | even e = + max 0 $ powUsingMulRecip one x e + | otherwise = powUsingMulRecip one x e + instance - (CanPow MPBall b - , CanEnsureCE es b - , CanEnsureCE es (PowType MPBall b) - , CanEnsureCE es (PowTypeNoCN MPBall b) - , SuitableForCE es) + (CanPow MPBall b) => - CanPow MPBall (CollectErrors es b) + CanPow MPBall (CN 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 + type PowType MPBall (CN b) = CN (PowType MPBall b) + pow = CN.liftT1 pow instance - (CanPow a MPBall - , CanEnsureCE es a - , CanEnsureCE es (PowType a MPBall) - , CanEnsureCE es (PowTypeNoCN a MPBall) - , SuitableForCE es) + (CanPow a MPBall) => - CanPow (CollectErrors es a) MPBall + CanPow (CN 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 + type PowType (CN a) MPBall = CN (PowType a MPBall) + pow = CN.lift1T pow instance CanDivIMod MPBall MPBall where + type DivIType MPBall MPBall = Integer divIMod x m - | m !>! 0 = (cn d, cn xm) - | otherwise = (err (0 :: Integer), err xm) + | m !>! 0 = (d, xm) + | otherwise = error $ "modulus not positive: " ++ show m where - d = floor $ centre $ (centreAsBall x) /! (centreAsBall m) + d = floor $ centre $ (centreAsBall x) / (centreAsBall m) xm = x - m*d - err :: (CanEnsureCN t) => t -> EnsureCN t - err s = noValueNumErrorCertainECN (Just s) $ OutOfRange $ "modulus not positive: " ++ show m diff --git a/src/AERN2/MP/Ball/Limit.hs b/src/AERN2/MP/Ball/Limit.hs new file mode 100644 --- /dev/null +++ b/src/AERN2/MP/Ball/Limit.hs @@ -0,0 +1,155 @@ +{-# LANGUAGE CPP #-} +-- #define DEBUG + +{-# OPTIONS_GHC -Wno-orphans #-} +{-| + Module : AERN2.MP.Ball.Limit + Description : limits of MPBall sequences + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Limits of MPBall sequences. +-} +module AERN2.MP.Ball.Limit 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 Numeric.CollectErrors as CN + +-- import Math.NumberTheory.Logarithms (integerLog2) + +import AERN2.Limit +import AERN2.MP.Ball.Type +import AERN2.MP.Ball.Field () + +{- + The following strategies are inspired by + Mueller: The iRRAM: Exact Arithmetic in C++, Section 10.1 + https://link.springer.com/chapter/10.1007/3-540-45335-0_14 +-} + +instance HasLimits Rational (CN MPBall -> CN MPBall) where + type LimitType Rational (CN MPBall -> CN MPBall) = (CN MPBall -> CN MPBall) + limit = limitMPBall (\ac -> 0.5^(fromAccuracy ac)) + +instance HasLimits Integer (CN MPBall -> CN MPBall) where + type LimitType Integer (CN MPBall -> CN MPBall) = (CN MPBall -> CN MPBall) + limit = limitMPBall (\ac -> (fromAccuracy ac)) + -- limit s = limit (s . getN) + -- where + -- -- q > 0 => 0 < 0.5^(getN q) <= q + -- getN :: Rational -> Integer + -- getN r = integer $ integerLog2 (max 1 $ (ceiling $ 1/r) - 1) + 1 + +instance HasLimits Int (CN MPBall -> CN MPBall) where + type LimitType Int (CN MPBall -> CN MPBall) = (CN MPBall -> CN MPBall) + limit s = limit (s . c) + where + c :: Integer -> Int + c = int + + +limitMPBall :: (Accuracy -> ix) -> (ix -> (CN MPBall -> CN MPBall)) -> (CN MPBall -> CN MPBall) +limitMPBall ac2ix fs x = + maybeTrace ("limit (CN MPBall -> CN MPBall): x = " ++ show x) $ + maybeTrace ("limit (CN MPBall -> CN MPBall): xPNext = " ++ show xPNext) $ + maybeTrace ("limit (CN MPBall -> CN MPBall): accuracies = " ++ show accuracies) $ + tryAccuracies accuracies + where + acX = getFiniteAccuracy x + accuracies = aux (fromAccuracy acX) + where + aux a + | a >= 4 = bits a : aux ((100 * a) `divI` 105) + | otherwise = [bits a] + xPNext = setPrecision (increaseP $ getPrecision x) x + increaseP p = + prec $ (integer p) + 10 + -- prec $ ((101 * (integer p)) `div` 100) + 1 + + tryAccuracies [] = + CN.noValueNumErrorPotential $ CN.NumError "limit (CN MPBall -> CN MPBall) failed" + tryAccuracies (ac : rest) = + let result = withAccuracy ac in + case CN.hasError result of + False -> result + _ -> tryAccuracies rest + + withAccuracy ac = + maybeTrace ("limit (CN MPBall -> CN MPBall): withAccuracy: ac = " ++ show ac) $ + (fs ix xPNext) +- e + where + ix = ac2ix ac + e = 0.5^(fromAccuracy ac) + +-- data WithLipschitz f = WithLipschitz f f + +-- instance HasLimits Rational (WithLipschitz (MPBall -> CN MPBall)) where +-- type LimitType Rational (WithLipschitz (MPBall -> CN MPBall)) = (MPBall -> CN MPBall) +-- limit ffs xPre = +-- maybeTrace ("limit (MPBall -> CN MPBall): x = " ++ show x) $ +-- maybeTrace ("limit (MPBall -> CN MPBall): xC = " ++ show xC) $ +-- maybeTrace ("limit (MPBall -> CN MPBall): xE = " ++ show xE) $ +-- maybeTrace ("limit (MPBall -> CN MPBall): accuracies = " ++ show accuracies) $ +-- tryAccuracies accuracies +-- where +-- acX = getFiniteAccuracy x +-- accuracies = aux (fromAccuracy acX) +-- where +-- aux a +-- | a >= 4 = bits a : aux ((100 * a) `divI` 105) +-- | otherwise = [bits a] +-- x = increasePrec xPre +-- where +-- increasePrec z = setPrecision (inc $ getPrecision xPre) z +-- inc p = +-- prec $ (integer p) + 10 +-- -- prec $ ((101 * (integer p)) `div` 100) + 1 +-- xC = centreAsBall x +-- xE = radius x + +-- tryAccuracies [] = +-- noValueNumErrorPotentialECN (Nothing :: Maybe MPBall) $ +-- NumError "limit (MPBall -> CN MPBall) failed" +-- tryAccuracies (ac : rest) = +-- let result = withAccuracy ac in +-- case getErrorsCN result of +-- [] -> result +-- _ -> tryAccuracies rest + +-- withAccuracy ac = +-- maybeTrace ("limit (MPBall -> CN MPBall): withAccuracy: ac = " ++ show ac) $ +-- lift1CE (updateRadius (+ (errorBound e))) fx +-- where +-- e = 0.5^!(fromAccuracy ac) +-- WithLipschitz f f' = ffs e +-- fxC_CN = f xC +-- f'x_CN = f' x +-- fx :: CN MPBall +-- fx = +-- case (ensureNoCN fxC_CN, ensureNoCN f'x_CN) of +-- ((Just fxC, []), (Just f'x, [])) -> +-- cn $ updateRadius (+ (xE * (errorBound f'x))) fxC +-- _ -> +-- f x -- fallback + +-- instance HasLimits Rational (WithLipschitz (CReal -> CReal)) where +-- type LimitType Rational (WithLipschitz (CReal -> CReal)) = (CReal -> CReal) +-- limit ffs x = limit fs x +-- where +-- fs e = f +-- where +-- WithLipschitz f _ = ffs e diff --git a/src/AERN2/MP/Ball/PreludeOps.hs b/src/AERN2/MP/Ball/PreludeOps.hs --- a/src/AERN2/MP/Ball/PreludeOps.hs +++ b/src/AERN2/MP/Ball/PreludeOps.hs @@ -1,3 +1,4 @@ +{-# OPTIONS_GHC -Wno-orphans #-} {-# LANGUAGE CPP #-} {-| Module : AERN2.MP.Ball.PreludeOps @@ -33,18 +34,19 @@ instance P.Eq MPBall where a == b = case a == b of - Just t -> t + CertainTrue -> True + CertainFalse -> False _ -> - error "Failed to decide equality of MPBalls. If you switch to MixedTypesNumPrelude instead of Prelude, comparison of MPBalls returns Maybe Bool instead of Bool." + error "Failed to decide equality of MPBalls. If you switch to MixedTypesNumPrelude instead of Prelude, comparison of MPBalls returns Kleenean 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 + (CertainTrue, _, _) -> P.LT + (_, CertainTrue, _) -> P.EQ + (_, _, CertainTrue) -> 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." + error "Failed to decide order of MPBalls. If you switch to MixedTypesNumPrelude instead of Prelude, comparison of MPBalls returns Kleenean instead of Bool." instance P.Num MPBall where fromInteger = convertExactly @@ -56,16 +58,16 @@ instance P.Fractional MPBall where fromRational = convertExactly . dyadic -- will work only for dyadic rationals - recip = (~!) . recip - (/) = (/!) + recip = recip + (/) = (/) instance P.Floating MPBall where - pi = error "MPBall: no pi :: MPBall, use pi ? (bitsS n) instead" - sqrt = (~!) . sqrt + pi = error "There is no pi :: MPBall, use pi :: Real instead" + sqrt = sqrt exp = exp sin = sin cos = cos - log = (~!) . log + log = log atan = error "MPBall: atan not implemented yet" atanh = error "MPBall: atanh not implemented yet" asin = error "MPBall: asin 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 @@ -1,4 +1,5 @@ {-# LANGUAGE CPP #-} +{-# OPTIONS_GHC -Wno-orphans #-} {-| Module : AERN2.MP.Ball.Tests Description : Tests for operations on arbitrary precision balls @@ -17,7 +18,6 @@ stack test aern2-mp --test-arguments "-a 1000 -m MPBall" @ -} - module AERN2.MP.Ball.Tests ( specMPBall, tMPBall @@ -69,11 +69,13 @@ tMPBall :: T MPBall tMPBall = T "MPBall" +-- tCNMPBall :: T (CN MPBall) +-- tCNMPBall = T "(CN MPBall)" + specMPBall :: Spec specMPBall = describe ("MPBall") $ do specCanSetPrecision tMPBall (printArgsIfFails2 "`contains`" contains) - specCanReduceSizeUsingAccuracyGuide tMPBall specConversion tInteger tMPBall mpBall (fst . integerBounds) describe "order" $ do specHasEqNotMixed tMPBall diff --git a/src/AERN2/MP/Ball/Type.hs b/src/AERN2/MP/Ball/Type.hs --- a/src/AERN2/MP/Ball/Type.hs +++ b/src/AERN2/MP/Ball/Type.hs @@ -1,4 +1,5 @@ {-# LANGUAGE TemplateHaskell #-} +{-# OPTIONS_GHC -Wno-orphans #-} {-| Module : AERN2.MP.Ball.Type Description : Arbitrary precision dyadic balls @@ -18,20 +19,25 @@ , module AERN2.MP.Accuracy , module AERN2.MP.Enclosure -- * The Ball type - , MPBall(..), CanBeMPBall, mpBall, CanBeMPBallP, mpBallP + , MPBall(..), CanBeMPBall, mpBall, cnMPBall + , CanBeMPBallP, mpBallP, cnMPBallP , reducePrecionIfInaccurate -- * Ball construction/extraction functions - , endpointsMP, fromEndpointsMP + , fromMPFloatEndpoints + , mpBallEndpoints, fromMPBallEndpoints ) where import MixedTypesNumPrelude -- import qualified Prelude as P -import Control.CollectErrors +import qualified Numeric.CollectErrors as CN import GHC.Generics (Generic) +import Control.DeepSeq +import qualified Data.List as List + import Text.Printf import AERN2.Normalize @@ -56,20 +62,43 @@ -- } 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]" (showMPFloat 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 NFData MPBall -instance (SuitableForCE es) => CanEnsureCE es MPBall where +instance Show MPBall where + show = showWithAccuracy (bits 50) +instance ShowWithAccuracy MPBall where + showWithAccuracy displayAC b@(MPBall x e) = + -- printf "[%s ± %s](prec=%s)" (show x) (showAC $ getAccuracy b) (show $ integer $ getPrecision b) + printf "[%s ± %s%s]" (dropSomeDigits $ showMPFloat x) eDS (showAC $ getAccuracy b) + -- "[" ++ show x ++ " ± " ++ show e ++ "](prec=" ++ (show $ integer $ getPrecision x) ++ ")" + where + eDS + | e == 0 = "0" + | otherwise = + case safeConvert (dyadic e) of + Right (eD :: Double) -> printf "~%.4g" $ eD + _ -> "" + dropSomeDigits s = + case List.findIndex (== '.') s of + Nothing -> s + Just ix -> withDotIx ix + where + withDotIx ix = + let maxLength = ix + displayAC_n in + let sTrimmed = take maxLength s in + if length sTrimmed < maxLength + then sTrimmed + else (take (maxLength - 3) sTrimmed) <> "..." + displayAC_n = + case displayAC of + Exact -> 1000000000 + NoInformation -> 0 + _ -> integer $ ac2prec displayAC + showAC Exact = "" + showAC NoInformation = "(oo)" + showAC ac = " ~2^(" ++ show (negate $ fromAccuracy ac) ++ ")" + -- instance CanTestValid MPBall where -- isValid = isFinite @@ -108,6 +137,13 @@ p_e_nb = prec $ max 2 (10 + nb + fromAccuracy bAcc) (NormBits nb) = bNorm +instance CanGiveUpIfVeryInaccurate MPBall where + giveUpIfVeryInaccurate = (aux =<<) + where + aux b@(MPBall _ e) + | e > 1000 = CN.noValueNumErrorCertain $ numErrorVeryInaccurate "MPBall" "" + | otherwise = cn b + instance CanTestContains MPBall MPBall where contains (MPBall xLarge eLarge) (MPBall xSmall eSmall) = xLargeDy - eLargeDy <= xSmallDy - eSmallDy @@ -134,42 +170,29 @@ {- ball construction/extraction functions -} -instance IsInterval MPBall MPFloat where +instance IsInterval MPBall where + type IntervalEndpoint MPBall = MPFloat fromEndpoints l u | u < l = fromEndpoints u l | otherwise = - MPBall (mpFloat cDy) (errorBound $ mpFloat eDy) + MPBall c (errorBound e) 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) + c = (l +. u) *. (mpFloat $ dyadic 0.5) + e = (u -^ c) `max` (c -^ l) + endpoints (MPBall x e) = (l, u) where - xDy = dyadic x - eDy = dyadic e - lDy = xDy - eDy - uDy = xDy + eDy + eFl = mpFloat e + l = x -. eFl + u = x +^ eFl -fromEndpointsMP :: MPFloat -> MPFloat -> MPBall -fromEndpointsMP = fromEndpoints +fromMPFloatEndpoints :: MPFloat -> MPFloat -> MPBall +fromMPFloatEndpoints = fromEndpoints -endpointsMP :: MPBall -> (MPFloat, MPFloat) -endpointsMP = endpoints +fromMPBallEndpoints :: MPBall -> MPBall -> MPBall +fromMPBallEndpoints = fromEndpointsAsIntervals -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 +mpBallEndpoints :: MPBall -> (MPBall, MPBall) +mpBallEndpoints = endpointsAsIntervals instance IsBall MPBall where type CentreType MPBall = Dyadic @@ -188,6 +211,8 @@ mpBallP :: (CanBeMPBallP t) => Precision -> t -> MPBall mpBallP = convertP +cnMPBallP :: (CanBeMPBallP a) => Precision -> CN a -> CN MPBall +cnMPBallP p = fmap (mpBallP p) {--- constructing an exact ball ---} @@ -196,41 +221,19 @@ mpBall :: (CanBeMPBall t) => t -> MPBall mpBall = convertExactly +cnMPBall :: (CanBeMPBall a) => CN a -> CN MPBall +cnMPBall = fmap mpBall + + {-- 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 + (_, MPBall boundMP _) = mpBallEndpoints $ absRaw ball instance HasApproximate MPBall where type Approximate MPBall = (MPFloat, Bool) @@ -267,8 +270,8 @@ absRaw :: MPBall -> MPBall absRaw b | l < 0 && 0 < r = - fromEndpointsMP (mpFloat 0) (max (-l) r) + fromEndpoints (mpFloat 0) (max (-l) r) | 0 <= l = b | otherwise = -b where - (l,r) = endpointsMP b + (l,r) = endpoints 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,6 +1,6 @@ {-# LANGUAGE CPP #-} -- #define DEBUG -{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} {-| Module : AERN2.MP.Dyadic Description : Dyadics with exact ring operations @@ -16,7 +16,6 @@ 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 ( -- * Dyadic numbers and their basic operations @@ -40,13 +39,15 @@ import MixedTypesNumPrelude import qualified Prelude as P -import Control.CollectErrors +import Control.CollectErrors (CollectErrors(..), CanBeErrors) +import qualified Control.CollectErrors as CE +-- import qualified Numeric.CollectErrors as CN import Text.Printf import Text.Regex.TDFA import Data.Typeable -import Data.Convertible +-- import Data.Convertible import Test.Hspec import Test.QuickCheck @@ -75,6 +76,7 @@ instance OrderedCertainlyRing (CN Dyadic) instance HasAccuracy Dyadic where getAccuracy _ = Exact +instance CanGiveUpIfVeryInaccurate Dyadic -- ie, never give up instance Show Dyadic where show (Dyadic x) @@ -106,7 +108,7 @@ [nS,eS] -> case (reads nS, reads eS) of ([(n,"")],[(e,"")]) -> - [((n :: Integer)*(dyadic 0.5)^!(e :: Integer), afterS)] + [((n :: Integer)*(dyadic 0.5)^(e :: Integer), afterS)] _ -> tryNext _ -> tryNext where @@ -114,9 +116,7 @@ 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 @@ -151,7 +151,7 @@ | isDyadic = Right $ Dyadic (ceduCentre $ fromRationalCEDU (prec $ max 2 (dp + np + 1)) q) | otherwise = convError "this number is not dyadic" q where - isDyadic = d == 2^!dp + isDyadic = d == 2^dp dp = integerLog2 d d = denominator q np = integerLog2 (max 1 $ abs $ numerator q) @@ -180,29 +180,25 @@ instance (HasEqAsymmetric Dyadic b - , CanEnsureCE es b - , CanEnsureCE es (EqCompareType Dyadic b) - , IsBool (EnsureCE es (EqCompareType Dyadic b)) - , SuitableForCE es) + , IsBool (CollectErrors es (EqCompareType Dyadic b)) + , CanBeErrors es) => - HasEqAsymmetric Dyadic (CollectErrors es b) + HasEqAsymmetric Dyadic (CollectErrors es b) where - type EqCompareType Dyadic (CollectErrors es b) = - EnsureCE es (EqCompareType Dyadic b) - equalTo = lift2TLCE equalTo + type EqCompareType Dyadic (CollectErrors es b) = + CollectErrors es (EqCompareType Dyadic b) + equalTo = CE.liftT1 equalTo instance (HasEqAsymmetric a Dyadic - , CanEnsureCE es a - , CanEnsureCE es (EqCompareType a Dyadic) - , IsBool (EnsureCE es (EqCompareType a Dyadic)) - , SuitableForCE es) + , IsBool (CollectErrors es (EqCompareType a Dyadic)) + , CanBeErrors es) => HasEqAsymmetric (CollectErrors es a) Dyadic where type EqCompareType (CollectErrors es a) Dyadic = - EnsureCE es (EqCompareType a Dyadic) - equalTo = lift2TCE equalTo + CollectErrors es (EqCompareType a Dyadic) + equalTo = CE.lift1T equalTo instance CanTestZero Dyadic @@ -228,35 +224,31 @@ instance (HasOrderAsymmetric Dyadic b - , CanEnsureCE es b - , CanEnsureCE es (OrderCompareType Dyadic b) - , IsBool (EnsureCE es (OrderCompareType Dyadic b)) - , SuitableForCE es) + , IsBool (CollectErrors es (OrderCompareType Dyadic b)) + , CanBeErrors 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 + CollectErrors es (OrderCompareType Dyadic b) + lessThan = CE.liftT1 lessThan + leq = CE.liftT1 leq + greaterThan = CE.liftT1 greaterThan + geq = CE.liftT1 geq instance (HasOrderAsymmetric a Dyadic - , CanEnsureCE es a - , CanEnsureCE es (OrderCompareType a Dyadic) - , IsBool (EnsureCE es (OrderCompareType a Dyadic)) - , SuitableForCE es) + , IsBool (CollectErrors es (OrderCompareType a Dyadic)) + , CanBeErrors 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 + CollectErrors es (OrderCompareType a Dyadic) + lessThan = CE.lift1T lessThan + leq = CE.lift1T leq + greaterThan = CE.lift1T greaterThan + geq = CE.lift1T geq instance CanTestPosNeg Dyadic @@ -306,29 +298,25 @@ instance (CanMinMaxAsymmetric Dyadic b - , CanEnsureCE es b - , CanEnsureCE es (MinMaxType Dyadic b) - , SuitableForCE es) + , CanBeErrors es) => CanMinMaxAsymmetric Dyadic (CollectErrors es b) where type MinMaxType Dyadic (CollectErrors es b) = - EnsureCE es (MinMaxType Dyadic b) - min = lift2TLCE min - max = lift2TLCE max + CollectErrors es (MinMaxType Dyadic b) + min = CE.liftT1 min + max = CE.liftT1 max instance (CanMinMaxAsymmetric a Dyadic - , CanEnsureCE es a - , CanEnsureCE es (MinMaxType a Dyadic) - , SuitableForCE es) + , CanBeErrors es) => CanMinMaxAsymmetric (CollectErrors es a) Dyadic where type MinMaxType (CollectErrors es a) Dyadic = - EnsureCE es (MinMaxType a Dyadic) - min = lift2TCE min - max = lift2TCE max + CollectErrors es (MinMaxType a Dyadic) + min = CE.lift1T min + max = CE.lift1T max {- addition -} @@ -358,27 +346,23 @@ instance (CanAddAsymmetric Dyadic b - , CanEnsureCE es b - , CanEnsureCE es (AddType Dyadic b) - , SuitableForCE es) + , CanBeErrors es) => CanAddAsymmetric Dyadic (CollectErrors es b) where type AddType Dyadic (CollectErrors es b) = - EnsureCE es (AddType Dyadic b) - add = lift2TLCE add + CollectErrors es (AddType Dyadic b) + add = CE.liftT1 add instance (CanAddAsymmetric a Dyadic - , CanEnsureCE es a - , CanEnsureCE es (AddType a Dyadic) - , SuitableForCE es) + , CanBeErrors es) => CanAddAsymmetric (CollectErrors es a) Dyadic where type AddType (CollectErrors es a) Dyadic = - EnsureCE es (AddType a Dyadic) - add = lift2TCE add + CollectErrors es (AddType a Dyadic) + add = CE.lift1T add {- subtraction -} @@ -408,27 +392,23 @@ instance (CanSub Dyadic b - , CanEnsureCE es b - , CanEnsureCE es (SubType Dyadic b) - , SuitableForCE es) + , CanBeErrors es) => CanSub Dyadic (CollectErrors es b) where type SubType Dyadic (CollectErrors es b) = - EnsureCE es (SubType Dyadic b) - sub = lift2TLCE sub + CollectErrors es (SubType Dyadic b) + sub = CE.liftT1 sub instance (CanSub a Dyadic - , CanEnsureCE es a - , CanEnsureCE es (SubType a Dyadic) - , SuitableForCE es) + , CanBeErrors es) => CanSub (CollectErrors es a) Dyadic where type SubType (CollectErrors es a) Dyadic = - EnsureCE es (SubType a Dyadic) - sub = lift2TCE sub + CollectErrors es (SubType a Dyadic) + sub = CE.lift1T sub {- multiplication -} @@ -459,119 +439,89 @@ instance (CanMulAsymmetric Dyadic b - , CanEnsureCE es b - , CanEnsureCE es (MulType Dyadic b) - , SuitableForCE es) + , CanBeErrors es) => CanMulAsymmetric Dyadic (CollectErrors es b) where type MulType Dyadic (CollectErrors es b) = - EnsureCE es (MulType Dyadic b) - mul = lift2TLCE mul + CollectErrors es (MulType Dyadic b) + mul = CE.liftT1 mul instance (CanMulAsymmetric a Dyadic - , CanEnsureCE es a - , CanEnsureCE es (MulType a Dyadic) - , SuitableForCE es) + , CanBeErrors es) => CanMulAsymmetric (CollectErrors es a) Dyadic where type MulType (CollectErrors es a) Dyadic = - EnsureCE es (MulType a Dyadic) - mul = lift2TCE mul + CollectErrors es (MulType a Dyadic) + mul = CE.lift1T mul instance CanPow Dyadic Integer where - powNoCN = powUsingMul (dyadic 1) - pow = integerPowCN (powUsingMul (dyadic 1)) + pow = powUsingMul (dyadic 1) instance CanPow Dyadic Int where - powNoCN = powUsingMul (dyadic 1) - pow = integerPowCN (powUsingMul (dyadic 1)) + pow = powUsingMul (dyadic 1) instance (CanDiv a Dyadic - , CanEnsureCE es a - , CanEnsureCE es (DivType a Dyadic) - , CanEnsureCE es (DivTypeNoCN a Dyadic) - , SuitableForCE es) + , CanBeErrors 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 + CollectErrors es (DivType a Dyadic) + divide = CE.lift1T divide instance CanDiv Integer Dyadic where - type DivTypeNoCN Integer Dyadic = Rational - divideNoCN a b = divideNoCN a (rational b) + type DivType Integer Dyadic = Rational + divide a b = divide a (rational b) instance CanDiv Dyadic Integer where - type DivTypeNoCN Dyadic Integer = Rational - divideNoCN a b = divideNoCN (rational a) b + type DivType Dyadic Integer = Rational + divide a b = divide (rational a) b instance CanDiv Int Dyadic where - type DivTypeNoCN Int Dyadic = Rational - divideNoCN a b = divideNoCN a (rational b) + type DivType Int Dyadic = Rational + divide a b = divide a (rational b) instance CanDiv Dyadic Int where - type DivTypeNoCN Dyadic Int = Rational - divideNoCN a b = divideNoCN (rational a) b + type DivType Dyadic Int = Rational + divide a b = divide (rational a) b instance CanDiv Rational Dyadic where - type DivTypeNoCN Rational Dyadic = Rational - divideNoCN = convertSecond divideNoCN + type DivType Rational Dyadic = Rational + divide = convertSecond divide instance CanDiv Dyadic Rational where - type DivTypeNoCN Dyadic Rational = Rational - divideNoCN = convertFirst divideNoCN + type DivType Dyadic Rational = Rational + divide = convertFirst divide instance (CanDiv Dyadic b - , CanEnsureCE es b - , CanEnsureCE es (DivType Dyadic b) - , CanEnsureCE es (DivTypeNoCN Dyadic b) - , SuitableForCE es) + , CanBeErrors es) => - CanDiv Dyadic (CollectErrors es b) + 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 + type DivType Dyadic (CollectErrors es b) = + CollectErrors es (DivType Dyadic b) + divide = CE.liftT1 divide instance (CanPow Dyadic b - , CanEnsureCE es b - , CanEnsureCE es (PowTypeNoCN Dyadic b) - , CanEnsureCE es (PowType Dyadic b) - , SuitableForCE es) + , CanBeErrors 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 + CollectErrors es (PowType Dyadic b) + pow = CE.liftT1 pow instance (CanPow a Dyadic - , CanEnsureCE es a - , CanEnsureCE es (PowType a Dyadic) - , CanEnsureCE es (PowTypeNoCN a Dyadic) - , SuitableForCE es) + , CanBeErrors 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 + CollectErrors es (PowType a Dyadic) + pow = CE.lift1T pow instance CanTestFinite Dyadic where isFinite = isFinite . dyadicMPFloat 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 @@ -1,4 +1,4 @@ -{-# LANGUAGE Arrows #-} +{-# OPTIONS_GHC -Wno-orphans #-} {-| Module : AERN2.MP.Enclosure Description : Enclosure operations @@ -14,11 +14,16 @@ module AERN2.MP.Enclosure ( IsBall(..), ballFunctionUsingLipschitz - , IsInterval(..), intervalFunctionByEndpoints, intervalFunctionByEndpointsUpDown + , IsInterval(..), endpointL, endpointR + , fromEndpointsAsIntervals, endpointsAsIntervals, endpointLAsInterval, endpointRAsInterval + , intervalFunctionByEndpoints, intervalFunctionByEndpointsUpDown + , CanPlusMinus(..), (+-) , CanTestContains(..), CanMapInside(..), specCanMapInside , CanIntersectAsymmetric(..), CanIntersect + , CanIntersectBy, CanIntersectSameType , CanIntersectCNBy, CanIntersectCNSameType , CanUnionAsymmetric(..), CanUnion + , CanUnionBy, CanUnionSameType , CanUnionCNBy, CanUnionCNSameType ) where @@ -26,15 +31,15 @@ import MixedTypesNumPrelude -- import qualified Prelude as P -import Control.Arrow +-- import Control.Arrow import Test.Hspec import Test.QuickCheck -import Control.CollectErrors +import qualified Numeric.CollectErrors as CN import AERN2.MP.ErrorBound -import AERN2.MP.Accuracy +-- import AERN2.MP.Accuracy {- ball-specific operations -} @@ -52,12 +57,18 @@ to the explicit radius, making the centre exact. This may lose some information, but as a ball is equivalent to the original. For MPBall this function is pointless because it is equivalent to the identity. -} - makeExactCentre :: (IsBall t) => t -> t + makeExactCentre :: t -> t makeExactCentre v = updateRadius (+r) c where (c, r) = centreAsBallAndRadius v +instance (IsBall t => IsBall (CN t)) where + type CentreType (CN t) = CN (CentreType t) + centre = fmap centre + updateRadius f = fmap (updateRadius f) + centreAsBallAndRadius = error $ "centreAsBallAndRadius not defined for CN types" + {-| Computes a ball function @f@ on the centre and updating the error bound using a Lipschitz constant. -} @@ -75,33 +86,84 @@ {- interval-specific operations -} +class IsInterval i where + type IntervalEndpoint i + endpoints :: i -> (IntervalEndpoint i, IntervalEndpoint i) + fromEndpoints :: IntervalEndpoint i -> IntervalEndpoint i -> i -class IsInterval i e where - fromEndpoints :: e -> e -> i - endpoints :: i -> (e,e) +instance (IsInterval t) => (IsInterval (CN t)) where + type (IntervalEndpoint (CN t)) = CN (IntervalEndpoint t) + fromEndpoints l u = CN.lift2 fromEndpoints l u + endpoints = CN.liftPair endpoints +endpointL :: (IsInterval i) => i -> IntervalEndpoint i +endpointL = fst . endpoints + +endpointR :: (IsInterval i) => i -> IntervalEndpoint i +endpointR = snd . endpoints + +endpointsAsIntervals :: + (IsInterval i) => i -> (i,i) +endpointsAsIntervals x = (lI,rI) + where + lI = fromEndpoints l l + rI = fromEndpoints r r + (l,r) = endpoints x + +endpointLAsInterval :: (IsInterval i) => i -> i +endpointLAsInterval = fst . endpointsAsIntervals + +endpointRAsInterval :: (IsInterval i) => i -> i +endpointRAsInterval = snd . endpointsAsIntervals + +fromEndpointsAsIntervals :: + (IsInterval i, CanMinMaxSameType (IntervalEndpoint i)) => + i -> i -> i +fromEndpointsAsIntervals l r = + fromEndpoints lMP uMP + where + lMP = min llMP rlMP + uMP = max luMP ruMP + (llMP, luMP) = endpoints l + (rlMP, ruMP) = endpoints r + +{- plusMinus (+-) operator -} + +class CanPlusMinus t1 t2 where + type PlusMinusType t1 t2 + type PlusMinusType t1 t2 = t1 + {-| Operator for constructing or enlarging enclosures such as balls or intervals -} + plusMinus :: t1 -> t2 -> PlusMinusType t1 t2 + +infixl 6 +- + +{-| Operator for constructing or enlarging enclosures such as balls or intervals -} +(+-) :: (CanPlusMinus t1 t2) => t1 -> t2 -> PlusMinusType t1 t2 +(+-) = plusMinus + + {-| Computes a *monotone* ball function @f@ on intervals using the interval endpoints. -} intervalFunctionByEndpoints :: - (IsInterval t t, HasEqCertainly t t) + (IsInterval t, CanMinMaxSameType (IntervalEndpoint t), HasEqCertainly t t) => (t -> t) {-^ @fThin@: a version of @f@ that works well on thin intervals -} -> (t -> t) {-^ @f@ on *large* intervals -} intervalFunctionByEndpoints fThin x | l !==! u = fThin l - | otherwise = fromEndpoints (fThin l) (fThin u) + | otherwise = fromEndpointsAsIntervals (fThin l) (fThin u) where - (l,u) = endpoints x + (l,u) = endpointsAsIntervals x {-| Computes a *monotone* ball function @f@ on intervals using the interval endpoints. -} intervalFunctionByEndpointsUpDown :: - (IsInterval t e) + (IsInterval t) => - (e -> e) {-^ @fDown@: a version of @f@ working on endpoints, rounded down -} -> - (e -> e) {-^ @fUp@: a version of @f@ working on endpoints, rounded up -} -> + (IntervalEndpoint t -> IntervalEndpoint t) {-^ @fDown@: a version of @f@ working on endpoints, rounded down -} -> + (IntervalEndpoint t -> IntervalEndpoint t) {-^ @fUp@: a version of @f@ working on endpoints, rounded up -} -> (t -> t) {-^ @f@ on intervals rounding *outwards* -} intervalFunctionByEndpointsUpDown fDown fUp x = fromEndpoints (fDown l) (fUp u) @@ -144,170 +206,135 @@ (CanIntersectAsymmetric e1 e2, CanIntersectAsymmetric e1 e2 , IntersectionType e1 e2 ~ IntersectionType e2 e1) +{-| A set intersection (usually partial) -} class CanIntersectAsymmetric e1 e2 where type IntersectionType e1 e2 - type IntersectionType e1 e2 = EnsureCN e1 + type IntersectionType e1 e2 = CN e1 intersect :: e1 -> e2 -> IntersectionType e1 e2 +type CanIntersectBy e1 e2 = + (CanIntersect e1 e2, IntersectionType e1 e2 ~ e1) + +type CanIntersectSameType e1 = + (CanIntersectBy e1 e1) + type CanIntersectCNBy e1 e2 = - (CanIntersect e1 e2, IntersectionType e1 e2 ~ EnsureCN e1 - , CanIntersect (EnsureCN e1) e2, IntersectionType (EnsureCN e1) e2 ~ EnsureCN e1) + (CanIntersect e1 e2, IntersectionType e1 e2 ~ CN e1) + type CanIntersectCNSameType e1 = - (CanIntersectCNBy e1 e1 - , CanIntersect (EnsureCN e1) (EnsureCN e1), IntersectionType (EnsureCN e1) (EnsureCN e1) ~ EnsureCN e1) + (CanIntersectCNBy e1 e1) -instance CanIntersectAsymmetric Bool Bool where +instance + CanIntersectAsymmetric Bool Bool + where intersect b1 b2 | b1 == b2 = cn b1 | otherwise = - noValueNumErrorCertainCN $ NumError "empty Boolean intersection" + CN.noValueNumErrorCertain $ CN.NumError "empty Boolean intersection" instance - (CanIntersectAsymmetric Bool b - , CanEnsureCE es b - , CanEnsureCE es (IntersectionType Bool b) - , SuitableForCE es) - => - CanIntersectAsymmetric Bool (CollectErrors es b) + CanIntersectAsymmetric Kleenean Kleenean where - type IntersectionType Bool (CollectErrors es b) = - EnsureCE es (IntersectionType Bool b) - intersect = lift2TLCE intersect + intersect CertainTrue CertainFalse = + CN.noValueNumErrorCertain $ CN.NumError "empty Kleenean intersection" + intersect CertainFalse CertainTrue = + CN.noValueNumErrorCertain $ CN.NumError "empty Kleenean intersection" + intersect TrueOrFalse k2 = cn k2 + intersect k1 _ = cn k1 -instance - (CanIntersectAsymmetric a Bool - , CanEnsureCE es a - , CanEnsureCE es (IntersectionType a Bool) - , SuitableForCE es) +instance + (CanIntersectAsymmetric a b, IntersectionType a b ~ CN c) => - CanIntersectAsymmetric (CollectErrors es a) Bool + CanIntersectAsymmetric (CN a) (CN b) where - type IntersectionType (CollectErrors es a) Bool = - EnsureCE es (IntersectionType a Bool) - intersect = lift2TCE intersect + type IntersectionType (CN a) (CN b) = IntersectionType a b + intersect = CN.lift2CN intersect instance - (CanIntersectAsymmetric (Maybe a) b - , CanEnsureCE es b - , CanEnsureCE es (IntersectionType (Maybe a) b) - , SuitableForCE es) + (CanIntersectAsymmetric (CN Bool) (CN b)) => - CanIntersectAsymmetric (Maybe a) (CollectErrors es b) + CanIntersectAsymmetric Bool (CN b) where - type IntersectionType (Maybe a) (CollectErrors es b) = - EnsureCE es (IntersectionType (Maybe a) b) - intersect = lift2TLCE intersect + type IntersectionType Bool (CN b) = IntersectionType (CN Bool) (CN b) + intersect b1 = intersect (cn b1) instance - (CanIntersectAsymmetric a (Maybe b) - , CanEnsureCE es a - , CanEnsureCE es (IntersectionType a (Maybe b)) - , SuitableForCE es) + (CanIntersectAsymmetric (CN a) (CN Bool)) => - CanIntersectAsymmetric (CollectErrors es a) (Maybe b) + CanIntersectAsymmetric (CN a) Bool where - type IntersectionType (CollectErrors es a) (Maybe b) = - EnsureCE es (IntersectionType a (Maybe b)) - intersect = lift2TCE intersect - + type IntersectionType (CN a) Bool = IntersectionType (CN a) (CN Bool) + intersect b1 b2 = intersect b1 (cn b2) instance - (CanIntersectAsymmetric a b - , CanEnsureCN a, IntersectionType a b ~ EnsureCN a - , CanEnsureCN (EnsureCN a) - , CanEnsureCN b, EnsureCN b ~ EnsureCN a) + (CanIntersectAsymmetric (CN Kleenean) (CN b)) => - CanIntersectAsymmetric (Maybe a) (Maybe b) + CanIntersectAsymmetric Kleenean (CN b) where - type IntersectionType (Maybe a) (Maybe b) = EnsureCN (Maybe (IntersectionType a b)) - intersect (ma :: Maybe a) (mb :: Maybe b) = - case (ma, mb) of - (Just a, Just b) -> justCN sample_r (intersect a b) - (Just a, Nothing) -> justCN sample_r (ensureCN a) - (Nothing, Just b) -> justCN sample_r (ensureCN b) - _ -> cn (Nothing :: Maybe a) - where - sample_r = Nothing :: EnsureCN (Maybe (IntersectionType a b)) - -justCN :: (CanEnsureCN a) => Maybe a -> EnsureCN a -> EnsureCN (Maybe a) -justCN (_sample_a :: Maybe a) aCN = - case deEnsureCN aCN of - Right a -> cn (Just (a :: a)) - _ -> cn (Nothing :: Maybe a) - + type IntersectionType Kleenean (CN b) = IntersectionType (CN Kleenean) (CN b) + intersect k1 = intersect (cn k1) --- --- Version that removes inner CN: --- instance --- (CanIntersectCNSameType a, CanEnsureCN a) --- => --- CanIntersectAsymmetric (Maybe a) (Maybe a) --- where --- type IntersectionType (Maybe a) (Maybe a) = CN (Maybe (WithoutCN (IntersectionType a a))) --- intersect (Just a) (Just b) = fmap Just (intersect a b) --- intersect (Just a) Nothing = fmap Just (ensureCN a) --- intersect Nothing (Just b) = fmap Just (ensureCN b) --- intersect Nothing Nothing = cn Nothing --- instance - (CanIntersectAsymmetric e1 e2, SuitableForCE es - , CanEnsureCE es e1, CanEnsureCE es e2 - , CanEnsureCE es (IntersectionType e1 e2)) + (CanIntersectAsymmetric (CN a) (CN Kleenean)) => - CanIntersectAsymmetric (CollectErrors es e1) (CollectErrors es e2) + CanIntersectAsymmetric (CN a) Kleenean where - type IntersectionType (CollectErrors es e1) (CollectErrors es e2) = - EnsureCE es (IntersectionType e1 e2) - intersect = lift2CE intersect + type IntersectionType (CN a) Kleenean = IntersectionType (CN a) (CN Kleenean) + intersect k1 k2 = intersect k1 (cn k2) -{- union -} +{- set union -} type CanUnion e1 e2 = (CanUnionAsymmetric e1 e2, CanUnionAsymmetric e1 e2 , UnionType e1 e2 ~ UnionType e2 e1) +{-| A set union (usually partial) -} class CanUnionAsymmetric e1 e2 where type UnionType e1 e2 - type UnionType e1 e2 = EnsureCN e1 + type UnionType e1 e2 = CN e1 union :: e1 -> e2 -> UnionType e1 e2 +type CanUnionBy e1 e2 = + (CanUnion e1 e2, UnionType e1 e2 ~ e1) + +type CanUnionSameType e1 = + (CanUnionBy e1 e1) + type CanUnionCNBy e1 e2 = - (CanUnion e1 e2, UnionType e1 e2 ~ EnsureCN e1 - , CanUnion (EnsureCN e1) e2, UnionType (EnsureCN e1) e2 ~ EnsureCN e1) + (CanUnion e1 e2, UnionType e1 e2 ~ CN e1) type CanUnionCNSameType e1 = - (CanUnionCNBy e1 e1 - , CanUnion (EnsureCN e1) (EnsureCN e1), UnionType (EnsureCN e1) (EnsureCN e1) ~ EnsureCN e1) + (CanUnionCNBy e1 e1) -instance - (CanUnionAsymmetric e1 e2 - , CanEnsureCN e1, CanEnsureCN e2 - , CanEnsureCN (UnionType e1 e2)) +instance + (CanUnionAsymmetric a b, UnionType a b ~ CN c) => - CanUnionAsymmetric (CN e1) (CN e2) - -- a more general `CollectErrors` instance conflicts with the Arrow instance below + CanUnionAsymmetric (CN a) (CN b) where - type UnionType (CN e1) (CN e2) = - EnsureCN (UnionType e1 e2) - union = lift2CE union + type UnionType (CN a) (CN b) = UnionType a b + union = CN.lift2CN union -instance - (Arrow to, CanUnionAsymmetric e1 e2) - => - CanUnionAsymmetric (to Accuracy e1) (to Accuracy e2) - -- this instance is important for "parallel if" - where - type UnionType (to Accuracy e1) (to Accuracy e2) = - to Accuracy (UnionType e1 e2) - union xA yA = - proc ac -> - do - x <- xA -< ac - y <- yA -< ac - returnA -< union x y +-- TODO: move to aern2-real (or aern2-net) +-- instance +-- (Arrow to, CanUnionAsymmetric e1 e2) +-- => +-- CanUnionAsymmetric (to Accuracy e1) (to Accuracy e2) +-- -- this instance is important for "parallel if" +-- where +-- type UnionType (to Accuracy e1) (to Accuracy e2) = +-- to Accuracy (UnionType e1 e2) +-- union xA yA = +-- proc ac -> +-- do +-- x <- xA -< ac +-- y <- yA -< ac +-- returnA -< union x y -instance (CanUnionCNSameType t, CanEnsureCN t) => - HasIfThenElse (Maybe Bool) t +instance (CanUnionSameType t, CN.CanTakeCNErrors t) => + HasIfThenElse Kleenean t where - type IfThenElseType (Maybe Bool) t = EnsureCN t - ifThenElse (Just b) e1 e2 = cn $ if b then e1 else e2 - ifThenElse Nothing e1 e2 = e1 `union` e2 + type IfThenElseType Kleenean t = t + ifThenElse CertainTrue e1 _ = e1 + ifThenElse CertainFalse _ e2 = e2 + ifThenElse TrueOrFalse e1 e2 = e1 `union` e2 + 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 @@ -22,11 +22,15 @@ import qualified Prelude as P import Data.Typeable +import GHC.Generics (Generic) +import Control.DeepSeq import Test.QuickCheck -import Data.Convertible +-- import Data.Convertible +import Data.Ratio (numerator) + import Math.NumberTheory.Logarithms (integerLog2) import AERN2.MP.Precision @@ -45,8 +49,10 @@ {-| 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) + deriving (P.Eq, P.Ord, Typeable, Generic) +instance NFData ErrorBound + instance Show ErrorBound where show (ErrorBound d) = show d @@ -64,9 +70,16 @@ where eN = floor $ rational e eRecipN = ceiling $ rational $ one /. e + getFiniteAccuracy eb@(ErrorBound e) + | e == 0 = bits errorBoundPrecision + | otherwise = getAccuracy eb + {- conversions -} +instance ConvertibleExactly ErrorBound ErrorBound where + safeConvertExactly = Right + instance ConvertibleExactly ErrorBound MPFloat where safeConvertExactly = Right . er2mp @@ -197,12 +210,11 @@ | 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" + divide = 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 = @@ -216,6 +228,6 @@ | otherwise = do (s :: Integer) <- arbitrary - let resultR = ((abs s) `P.mod` (2^!35))/!(2^!32) + let resultR = ((abs s) `P.mod` (numerator $ 2^35))/(2^32) let result = convert resultR return result diff --git a/src/AERN2/MP/Float/Arithmetic.hs b/src/AERN2/MP/Float/Arithmetic.hs new file mode 100644 --- /dev/null +++ b/src/AERN2/MP/Float/Arithmetic.hs @@ -0,0 +1,104 @@ +{-# OPTIONS_GHC -Wno-orphans #-} +{-| + Module : AERN2.MP.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. +-} +module AERN2.MP.Float.Arithmetic + ( + -- * MPFloat basic arithmetic + addCEDU, subCEDU + , mulCEDU, divCEDU, recipCEDU + -- * MPFloat selected constants and operations + , piCEDU + , cosCEDU, sinCEDU + , sqrtCEDU, expCEDU, logCEDU + ) +where + +import MixedTypesNumPrelude +import qualified Prelude as P + +import AERN2.MP.Precision + +import qualified Data.CDAR as MPLow + +import AERN2.MP.Float.Auxi +import AERN2.MP.Float.Type + +{- common functions -} + +instance CanNeg MPFloat where + negate = P.negate + +instance CanAbs MPFloat where + abs = P.abs + +addCEDU :: MPFloat -> MPFloat -> BoundsCEDU MPFloat +addCEDU = binaryCEDU (P.+) + +subCEDU :: MPFloat -> MPFloat -> BoundsCEDU MPFloat +subCEDU = binaryCEDU (P.-) + +mulCEDU :: MPFloat -> MPFloat -> BoundsCEDU MPFloat +mulCEDU = binaryCEDU (P.*) + +divCEDU :: MPFloat -> MPFloat -> BoundsCEDU MPFloat +divCEDU x y + | y P.== (P.fromInteger 0) = getBoundsCEDU MPLow.Bottom + | otherwise = binaryCEDU (P./) x y + +recipCEDU :: MPFloat -> BoundsCEDU MPFloat +recipCEDU = unaryCEDU P.recip + +{- special constants and functions -} + +piCEDU :: Precision -> BoundsCEDU MPFloat +piCEDU pp = + getBoundsCEDU $ MPLow.piA (p2cdarPrec pp) + +cosCEDU :: MPFloat -> BoundsCEDU MPFloat +cosCEDU = unaryCEDU MPLow.cosA + +sinCEDU :: MPFloat -> BoundsCEDU MPFloat +sinCEDU = unaryCEDU MPLow.sinA + +sqrtCEDU :: MPFloat -> BoundsCEDU MPFloat +sqrtCEDU = unaryCEDU MPLow.sqrtA + +expCEDU :: MPFloat -> BoundsCEDU MPFloat +expCEDU = unaryCEDU MPLow.expA + +logCEDU :: MPFloat -> BoundsCEDU MPFloat +logCEDU = unaryCEDU MPLow.logA + +{- auxiliary functions to automatically determine result precision from operand precisions -} + +binaryCEDU :: + (MPFloat -> MPFloat -> MPFloat) -> + (MPFloat -> MPFloat -> BoundsCEDU MPFloat) +binaryCEDU op x y = + getBoundsCEDU $ op x y + +unaryCEDU :: + (MPFloat -> MPFloat) -> + (MPFloat -> BoundsCEDU MPFloat) +unaryCEDU op x = + getBoundsCEDU $ op x + +-- unaryPrecCEDU :: +-- Integer -> +-- (MPLow.Precision -> MPFloat -> MPFloat) -> +-- (MPFloat -> BoundsCEDU MPFloat) +-- unaryPrecCEDU addPrec op x@(MPLow.Approx mb _ _ s) = +-- getBoundsCEDU $ op ((-s P.+ mb) P.+ (int addPrec)) x +-- unaryPrecCEDU addPrec op MPLow.Bottom = +-- getBoundsCEDU $ op ((int $ integer defaultPrecision) P.+ (int addPrec)) MPLow.Bottom + diff --git a/src/AERN2/MP/Float/Conversions.hs b/src/AERN2/MP/Float/Conversions.hs new file mode 100644 --- /dev/null +++ b/src/AERN2/MP/Float/Conversions.hs @@ -0,0 +1,152 @@ +{-# OPTIONS_GHC -Wno-orphans #-} +{-| + Module : AERN2.MP.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 +-} +module AERN2.MP.Float.Conversions + ( + -- * MPFloat to other types (see also instances) + toDouble + -- * MPFloat constructors (see also instances) + , CanBeMPFloat, mpFloat + , fromIntegerCEDU + , fromRationalCEDU + -- * comparisons and constants (see also instances) + , zero, one, two + , nan, infinity + ) +where + +import MixedTypesNumPrelude +import qualified Prelude as P + +import Data.Ratio +-- import Data.Convertible + +-- import AERN2.Norm +import AERN2.MP.Precision + +import qualified Data.CDAR as MPLow + +import AERN2.MP.Float.Auxi +import AERN2.MP.Float.Type +import AERN2.MP.Float.Arithmetic + + +{- conversions to MPFloat -} + +type CanBeMPFloat t = ConvertibleExactly t MPFloat +mpFloat :: (CanBeMPFloat t) => t -> MPFloat +mpFloat = convertExactly + +instance ConvertibleExactly Integer MPFloat where + safeConvertExactly = + Right . P.fromInteger + +instance ConvertibleExactly Int MPFloat where + safeConvertExactly = safeConvertExactly . integer + +fromIntegerCEDU :: Precision -> Integer -> BoundsCEDU MPFloat +fromIntegerCEDU pp = + setPrecisionCEDU pp . P.fromInteger + +fromRationalCEDU :: Precision -> Rational -> BoundsCEDU MPFloat +fromRationalCEDU pp = + setPrecisionCEDU pp . (MPLow.toApproxMB (p2cdarPrec pp)) + +{- conversions from MPFloat -} + +instance ConvertibleExactly MPFloat Rational where + safeConvertExactly = Right . P.toRational + +toDouble :: MPFloat -> Double +toDouble = P.fromRational . rational + +instance Convertible MPFloat Double where + safeConvert x + | isFinite dbl = Right dbl + | otherwise = convError "conversion to double: out of bounds" x + where + dbl = toDouble x + + +instance CanRound MPFloat where + properFraction x = (n,f) + where + r = rational x + n = (numerator r) `P.quot` (denominator r) + f = ceduCentre $ x `subCEDU` (P.fromInteger n) + +{- 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 + +{- constants -} + +zero, one, two :: MPFloat +zero = mpFloat 0 +one = mpFloat 1 +two = mpFloat 2 + +nan, infinity :: MPFloat +nan = MPLow.Bottom +infinity = nan + +itisNaN :: MPFloat -> Bool +itisNaN MPLow.Bottom = True +itisNaN _ = False + +instance CanTestFinite MPFloat where + isInfinite = itisNaN + isFinite = not . itisNaN + +instance CanTestNaN MPFloat where + isNaN = itisNaN diff --git a/src/AERN2/MP/Float/Tests.hs b/src/AERN2/MP/Float/Tests.hs --- a/src/AERN2/MP/Float/Tests.hs +++ b/src/AERN2/MP/Float/Tests.hs @@ -1,3 +1,4 @@ +{-# OPTIONS_GHC -Wno-orphans #-} {-| Module : AERN2.MP.Float.Tests Description : Tests for operations on arbitrary precision floats @@ -16,7 +17,6 @@ stack test aern2-mp --test-arguments "-a 1000 -m MPFloat" @ -} - module AERN2.MP.Float.Tests ( specMPFloat, tMPFloat @@ -36,14 +36,14 @@ import Test.QuickCheck -- import qualified Test.Hspec.SmallCheck as SC -import Control.CollectErrors +-- import qualified Numeric.CollectErrors as CN import AERN2.Norm import AERN2.MP.Precision import AERN2.MP.Float.Auxi import AERN2.MP.Float.Type -import AERN2.MP.Float.Arithmetic +-- import AERN2.MP.Float.Arithmetic import AERN2.MP.Float.Conversions import AERN2.MP.Float.Operators @@ -62,7 +62,7 @@ (p :: Precision) <- arbitrary (s :: Integer) <- arbitrary ex <- choose (-20,10) - let resultR = s * (10.0^!ex) + let resultR = s * (10.0^ex) let result = ceduCentre $ fromRationalCEDU p resultR return result @@ -88,7 +88,7 @@ where l = mpFloat l_ u = mpFloat u_ - b = l +^ ((abs a) `modNoCN` (u-^l)) + b = l +^ ((abs a) `mod` (u-^l)) enforceRangeMP (Just l_, _) a | isInfinite a = abs a | l < a = a @@ -105,21 +105,20 @@ u = mpFloat u_ enforceRangeMP _ a = a -instance CanEnsureCE NumErrors MPFloat - instance CanDivIMod MPFloat MPFloat where - divIMod x m - | (not (isFinite m)) = (errM (d :: Integer), errM xm) - | (not (isFinite x)) = (errX (d :: Integer), errX xm) - | m > zero = (cn d, cn xm) - | otherwise = (errM (d :: Integer), errM xm) + type DivIType MPFloat MPFloat = Integer + divIMod x m + | (not (isFinite m)) = ((d :: Integer), xm) + | (not (isFinite x)) = ((d :: Integer), xm) + | m > zero = (d, xm) + | otherwise = ((d :: Integer), xm) where d = floor (x /^ m) xm = x -^ (mpFloat d)*^m - errM :: (CanEnsureCN t) => t -> EnsureCN t - errM s = noValueNumErrorCertainECN (Just s) $ OutOfRange $ "modulus not finite and positive: " ++ show m - errX :: (CanEnsureCN t) => t -> EnsureCN t - errX s = noValueNumErrorCertainECN (Just s) $ OutOfRange $ "modulus input not finite: " ++ show x + -- errM :: (CanEnsureCN t) => t -> EnsureCN t + -- errM s = CN.noValueNumErrorCertain $ CN.OutOfDomain $ "modulus not finite and positive: " ++ show m + -- errX :: (CanEnsureCN t) => t -> EnsureCN t + -- errX s = CN.noValueNumErrorCertain$ CN.OutOfDomain $ "modulus input not finite: " ++ show x {- approximate comparison -} @@ -145,7 +144,7 @@ | isNaN x || isNaN y = False | isInfinite x || isInfinite y = x == y | otherwise = - abs (x -. y) <= 0.5^!e + abs (x -. y) <= 0.5^e {-| Assert equality of two MPFloat's with tolerance derived from the size and precision @@ -492,14 +491,15 @@ cosDown x <=% 1 && cosUp x >=% -1 - it "cos(pi)=-1" $ do - property $ \ (p :: Precision) -> - let - piA = ceduCentre $ piCEDU p - (=~~=) = approxEqualWithArgs 1 [(piA,"pi")] - infix 4 =~~= - in - cosUp(piA) =~~= (-one) + -- TODO: fix accuracy of CDAR mBounds cosine + -- it "cos(pi)=-1" $ do + -- property $ \ (p :: Precision) -> + -- let + -- piA = ceduCentre $ piCEDU p + -- (=~~=) = approxEqualWithArgs 2 [(piA,"pi")] + -- infix 4 =~~= + -- in + -- cosUp(piA) =~~= (-one) it "cos(x)^2 + sin(x)^2 = 1" $ do property $ \ (x_ :: MPFloat) -> let x = enforceRangeMP (Just (-1000000), Just 1000000) x_ in diff --git a/src/AERN2/MP/Float/Type.hs b/src/AERN2/MP/Float/Type.hs new file mode 100644 --- /dev/null +++ b/src/AERN2/MP/Float/Type.hs @@ -0,0 +1,79 @@ +{-# OPTIONS_GHC -Wno-orphans #-} +{-# LANGUAGE StandaloneDeriving #-} +{-| + Module : AERN2.MP.Float.Type + Description : Arbitrary precision floating point numbers (via cdar) + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Arbitrary precision floating-point numbers, re-using CDAR Approx type. +-} +module AERN2.MP.Float.Type + ( + -- * MPFloat numbers and their basic operations + MPFloat + , showMPFloat + , setPrecisionCEDU + , p2cdarPrec + , getBoundsCEDU + ) +where + +import MixedTypesNumPrelude +import qualified Prelude as P + +-- import Data.Bits (unsafeShiftL) +import Data.Typeable + +import AERN2.Norm +import AERN2.MP.Precision +import AERN2.MP.Float.Auxi + +import qualified Data.CDAR as MPLow + +{-| Multiple-precision floating-point type based on CDAR.Approx with 0 radius. -} +type MPFloat = MPLow.Approx + +showMPFloat :: MPFloat -> String +showMPFloat x = MPLow.showA x + +deriving instance (Typeable MPFloat) + +p2cdarPrec :: Precision -> MPLow.Precision +p2cdarPrec = P.fromInteger . integer + +getBoundsCEDU :: MPFloat -> BoundsCEDU MPFloat +getBoundsCEDU (MPLow.Approx mb m e s) = + BoundsCEDU + (MPLow.Approx mb m 0 s) (MPLow.approxMB eb_mb e 0 s) + (MPLow.Approx mb (m-e) 0 s) (MPLow.Approx mb (m+e) 0 s) +getBoundsCEDU MPLow.Bottom = + BoundsCEDU + MPLow.Bottom MPLow.Bottom MPLow.Bottom MPLow.Bottom + +{-| The bit-size bound for the error bound in CEDU -} +eb_prec :: Precision +eb_prec = prec 63 + +{-| The bit-size bound for the error bound in CEDU -} +eb_mb :: Int +eb_mb = int $ integer eb_prec + +instance HasPrecision MPFloat where + getPrecision (MPLow.Approx mb _ _ _) = prec (P.toInteger $ mb) + getPrecision MPLow.Bottom = error "illegal MPFloat (Bottom)" + + +instance CanSetPrecision MPFloat where + setPrecision p = ceduCentre . setPrecisionCEDU p + +setPrecisionCEDU :: Precision -> MPFloat -> BoundsCEDU MPFloat +setPrecisionCEDU pp = getBoundsCEDU . MPLow.enforceMB . MPLow.setMB (p2cdarPrec pp) + +instance HasNorm MPFloat where + getNormLog (MPLow.Approx _ m _ s) = (getNormLog m) + (integer s) + getNormLog MPLow.Bottom = error "getNormLog undefined for Bottom" diff --git a/src/AERN2/MP/Precision.hs b/src/AERN2/MP/Precision.hs --- a/src/AERN2/MP/Precision.hs +++ b/src/AERN2/MP/Precision.hs @@ -1,8 +1,8 @@ {-# LANGUAGE TemplateHaskell #-} -{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} {-| Module : AERN2.MP.Precision - Description : Floating-point precision + Description : number of significant binary digits Copyright : (c) Michal Konecny License : BSD3 @@ -10,7 +10,8 @@ Stability : experimental Portability : portable - Floating-point precision type and its operations + Precision type and its operations. Precision expresses a limit on the size + of an approximation, roughly corresponding to the number of significant bits for a floating-point number and its generalisations to other types. -} module AERN2.MP.Precision ( @@ -27,7 +28,8 @@ import qualified Prelude as P import Text.Printf -import Control.CollectErrors +import Control.CollectErrors (CollectErrors(..), CanBeErrors) +-- import qualified Control.CollectErrors as CE import Data.Complex @@ -100,7 +102,7 @@ class HasPrecision t where getPrecision :: t -> Precision -class (HasPrecision t) => CanSetPrecision t where +class CanSetPrecision t where setPrecision :: Precision -> t -> t instance HasPrecision t => HasPrecision (Complex t) where @@ -124,24 +126,24 @@ instance HasPrecision t => HasPrecision (CollectErrors es t) where getPrecision vCE = - case getMaybeValueCE vCE of + case getMaybeValue vCE of Just v -> getPrecision v _ -> defaultPrecision instance CanSetPrecision t => CanSetPrecision (CollectErrors es t) where setPrecision p = fmap (setPrecision p) -lowerPrecisionIfAbove :: (CanSetPrecision t) => Precision -> t -> t +lowerPrecisionIfAbove :: (HasPrecision t, CanSetPrecision t) => Precision -> t -> t lowerPrecisionIfAbove p x | getPrecision x > p = setPrecision p x | otherwise = x -raisePrecisionIfBelow :: (CanSetPrecision t) => Precision -> t -> t +raisePrecisionIfBelow :: (HasPrecision t, CanSetPrecision t) => Precision -> t -> t raisePrecisionIfBelow p x | getPrecision x < p = setPrecision p x | otherwise = x specCanSetPrecision :: - (CanSetPrecision t, CanTestFinite t, Arbitrary t, Show t, Testable prop) + (HasPrecision t, CanSetPrecision t, CanTestFinite t, Arbitrary t, Show t, Testable prop) => (T t) -> (t -> t -> prop) -> Spec specCanSetPrecision (T typeName :: T t) check = @@ -213,7 +215,6 @@ (t1 -> t1 -> c) -> (t1 -> t2 -> c) convertPSecond = convertSecondUsing (\ b q -> convertP (getPrecision b) q) - instance Arbitrary Precision where arbitrary = sized $ \size -> choose (4*(size+1),10*(size+1)) >>= return . prec @@ -222,6 +223,6 @@ [[t| Bool |], [t| Integer |], [t| Int |], [t| Rational |], [t| Double |]] (\ t -> [d| - instance (ConvertibleWithPrecision $t t, Monoid es) => ConvertibleWithPrecision $t (CollectErrors es t) where - safeConvertP p = fmap (\v -> CollectErrors (Just v) mempty) . safeConvertP p + instance (ConvertibleWithPrecision $t t, CanBeErrors es) => ConvertibleWithPrecision $t (CollectErrors es t) where + safeConvertP p = fmap pure . safeConvertP p |])) diff --git a/src/AERN2/MP/WithCurrentPrec.hs b/src/AERN2/MP/WithCurrentPrec.hs new file mode 100644 --- /dev/null +++ b/src/AERN2/MP/WithCurrentPrec.hs @@ -0,0 +1,181 @@ +{-# LANGUAGE EmptyDataDecls #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE StandaloneDeriving #-} +{-| + Module : AERN2.MP.WithCurrentPrec + Description : Type wrapper setting default precision + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Type wrapper setting default precision. + + Not suitable for use with MixedTypesNumPrelude since we need binary operators to enforce + the same precision on both operands via the equality of their types. + + Borrowed some tricks from https://github.com/ekmett/rounded/blob/master/src/Numeric/Rounded/Precision.hs +-} +module AERN2.MP.WithCurrentPrec +( + WithCurrentPrec(..), runWithPrec, HasCurrentPrecision(..) + -- , _example1 , _example2 , _example3 +) +where + +import qualified MixedTypesNumPrelude as MxP +import Prelude +-- import Text.Printf + +-- import Text.Printf +import Numeric.CollectErrors (CN, cn, NumErrors, CanTakeErrors(..)) +-- import qualified Numeric.CollectErrors as CN + +import Data.Proxy +import Data.Reflection +import GHC.TypeLits + +-- import Data.Complex + +import AERN2.Limit +import AERN2.MP.Precision +import AERN2.MP.Ball + +class HasCurrentPrecision p where + getCurrentPrecision :: proxy p -> Precision + +instance KnownNat n => HasCurrentPrecision n where + getCurrentPrecision p = max (prec 2) . min maximumPrecision $ prec (natVal p) + +-- data PrecAdd10 (p :: *) + +-- instance (HasCurrentPrecision p) => HasCurrentPrecision (PrecAdd10 p) where +-- isPrecision (_ :: proxy _) = 10 + isPrecision (undefined :: proxy p) + +newtype WithCurrentPrec t p = WithCurrentPrec { unWithCurrentPrec :: t } + deriving (Show) + +deriving instance (CanTakeErrors NumErrors t) => (CanTakeErrors NumErrors (WithCurrentPrec t p)) + +runWithPrec :: Precision -> (forall n. (KnownNat n) => WithCurrentPrec t n) -> t +runWithPrec p (wfp :: (forall n. (KnownNat n) => WithCurrentPrec t n)) = + reifyNat (MxP.integer p) withNat + where + withNat :: KnownNat n => Proxy n -> t + withNat (_ :: Proxy n) = + unWithCurrentPrec (wfp :: WithCurrentPrec t n) + +-- -- The following does not work: +-- instance (CanAddAsymmetric t1 t2) => (CanAddAsymmetric (WithCurrentPrec t1 p) (WithCurrentPrec t2 p)) where +-- type AddType (WithCurrentPrec t1 p) (WithCurrentPrec t2 p) = WithCurrentPrec (AddType t1 t2) p +-- add (WithCurrentPrec a1) (WithCurrentPrec a2) = WithCurrentPrec $ a1 + a2 + +instance + (MxP.HasOrderAsymmetric t1 t2) + => + MxP.HasOrderAsymmetric (WithCurrentPrec t1 p1) (WithCurrentPrec t2 p2) + where + type OrderCompareType (WithCurrentPrec t1 p1) (WithCurrentPrec t2 p2) = MxP.OrderCompareType t1 t2 + greaterThan (WithCurrentPrec v1) (WithCurrentPrec v2) = MxP.greaterThan v1 v2 + lessThan (WithCurrentPrec v1) (WithCurrentPrec v2) = MxP.lessThan v1 v2 + geq (WithCurrentPrec v1) (WithCurrentPrec v2) = MxP.geq v1 v2 + leq (WithCurrentPrec v1) (WithCurrentPrec v2) = MxP.leq v1 v2 + +instance + (MxP.CanMinMaxAsymmetric t1 t2, p1 ~ p2) + => + MxP.CanMinMaxAsymmetric (WithCurrentPrec t1 p1) (WithCurrentPrec t2 p2) + where + type MinMaxType (WithCurrentPrec t1 p1) (WithCurrentPrec t2 p2) = WithCurrentPrec (MxP.MinMaxType t1 t2) p1 + min (WithCurrentPrec v1) (WithCurrentPrec v2) = WithCurrentPrec $ MxP.min v1 v2 + max (WithCurrentPrec v1) (WithCurrentPrec v2) = WithCurrentPrec $ MxP.max v1 v2 + +instance Eq t => Eq (WithCurrentPrec t p) where + (==) = lift2P (==) +instance Ord t => Ord (WithCurrentPrec t p) where + compare = lift2P compare + +instance + (HasCurrentPrecision p, Num t, ConvertibleWithPrecision Integer t) + => + Num (WithCurrentPrec t p) + where + fromInteger n = r + where + r = WithCurrentPrec $ convertP (getCurrentPrecision r) n + negate = lift1 negate + abs = lift1 abs + (+) = lift2 (+) + (*) = lift2 (*) + signum = lift1 signum + +instance + (HasCurrentPrecision p, Fractional t + , ConvertibleWithPrecision Integer t, ConvertibleWithPrecision Rational t) + => + Fractional (WithCurrentPrec t p) + where + fromRational q = r + where + r = WithCurrentPrec $ convertP (getCurrentPrecision r) q + recip = lift1 recip + (/) = lift2 (/) + +instance (HasCurrentPrecision p) => Floating (WithCurrentPrec (CN MPBall) p) where + pi = r + where + r = WithCurrentPrec $ cn $ piBallP (getCurrentPrecision r) + sqrt = lift1 sqrt + exp = lift1 exp + log = lift1 log + sin = lift1 sin + cos = lift1 cos + asin = lift1 asin + acos = lift1 acos + atan = lift1 atan + sinh = lift1 sinh + cosh = lift1 cosh + asinh = lift1 asinh + acosh = lift1 acosh + atanh = lift1 atanh + +instance + (HasLimits ix (CN MPBall -> CN MPBall) + , LimitType ix (CN MPBall -> CN MPBall) ~ (CN MPBall -> CN MPBall) + ,HasCurrentPrecision p) + => + HasLimits ix (WithCurrentPrec (CN MPBall) p) + where + type LimitType ix (WithCurrentPrec (CN MPBall) p) = WithCurrentPrec (CN MPBall) p + limit (s :: ix -> (WithCurrentPrec (CN MPBall) p)) = + WithCurrentPrec $ limit (snop) $ sample + where + sample :: CN MPBall + sample = setPrecision (getCurrentPrecision sampleP) 0 + sampleP :: WithCurrentPrec MPBall p + sampleP = error "sampleP is not defined, it is only a type proxy" + snop :: ix -> (CN MPBall -> CN MPBall) + snop ix _sample = unWithCurrentPrec $ s ix + +lift1 :: (t1 -> t2) -> (WithCurrentPrec t1 p) -> (WithCurrentPrec t2 p) +lift1 f (WithCurrentPrec v1) = WithCurrentPrec (f v1) + +lift2 :: (t1 -> t2 -> t3) -> (WithCurrentPrec t1 p) -> (WithCurrentPrec t2 p) -> (WithCurrentPrec t3 p) +lift2 f (WithCurrentPrec v1) (WithCurrentPrec v2) = WithCurrentPrec (f v1 v2) + +lift2P :: (t1 -> t2 -> t3) -> (WithCurrentPrec t1 p) -> (WithCurrentPrec t2 p) -> t3 +lift2P f (WithCurrentPrec v1) (WithCurrentPrec v2) = f v1 v2 + +_example1 :: CN MPBall +_example1 = runWithPrec (prec 1000) pi + +_example2 :: CN MPBall +_example2 = runWithPrec (prec 1000) $ pi - pi + +_example3 :: CN MPBall +_example3 = runWithPrec (prec 1000) $ sqrt 2 diff --git a/src/AERN2/Norm.hs b/src/AERN2/Norm.hs --- a/src/AERN2/Norm.hs +++ b/src/AERN2/Norm.hs @@ -56,7 +56,7 @@ getNormLog x | x == 0.0 = NormZero | abs x >= 1.0 = getNormLog $ ceiling $ abs x - | otherwise = NormBits $ negate $ integer $ integerLog2 $ floor $ (1 /! (abs x)) + | otherwise = NormBits $ negate $ integer $ integerLog2 $ floor $ (1 / (abs x)) instance (HasNorm t)