packages feed

aern2-mp (empty) → 0.1.0.0

raw patch · 34 files changed

+4816/−0 lines, 34 filesdep +QuickCheckdep +aern2-mpdep +base

Dependencies added: QuickCheck, aern2-mp, base, convertible, haskell-mpfr, hmpfr, hspec, integer-gmp, integer-logarithms, lens, mixed-types-num, regex-tdfa, rounded, template-haskell

Files

+ LICENSE view
@@ -0,0 +1,12 @@+Copyright (c) 2015-2017 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:++1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.++2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.++3. Neither the name of the software nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ aern2-mp.cabal view
@@ -0,0 +1,153 @@+name:           aern2-mp+version:        0.1.0.0+cabal-version:  >= 1.9.2+build-type:     Simple+homepage:       https://github.com/michalkonecny/aern2/aern2-mp+author:         Michal Konecny+maintainer:     Michal Konecny <mikkonecny@gmail.com>+copyright:      (c) 2015-2017 Michal Konecny+license:        BSD3+license-file:   LICENSE+stability:      experimental+tested-with:    GHC==7.10.3, GHC==8.0.2+category:       Math+synopsis:       Multi-precision floats via MPFR+Description:+  Arbitrary-precision dyadic numbers and interval arithmetic, via the types:+  .+  * Dyadic:  variable-precision floats with exact ring operations+  .+  * MPBall: float ± error bound with field & elementary (interval-like) operations+  .+  By default, using hmpfr when compiling with ghc >=7.10.  Older versions of ghc are no longer tested.+  .+  Set flag MPFRRounded to bind MPFR not via hmpfr but via a tweaked version of package+  rounded:+  .+  * https://github.com/michalkonecny/rounded+  .+  This backend is considerably faster than hmpfr.+  It is likely to become the default in future.+  .+  There is a plan to add an Integer-only backend so that this library can+  be used without MPFR.++source-repository head+  type:     git+  location: https://github.com/mikkonecny/aern2.git+  subdir: aern2-mp++flag MPFR+  Description: Use a MPFR backend (default); Integer backend is planned+  Default:     True++flag MPFRRounded+  Description: Use rounded (Numeric.RoundedSimple) as MPFR backend (not default)+  -- Default:     True+  Default:     False++library+  hs-source-dirs:  src+  build-depends:+    base == 4.*+    , integer-logarithms >= 1.0 && < 1.1+    , convertible >= 1.1.1.0 && < 1.2+    , regex-tdfa >= 1.2.1 && < 1.3+    , hspec >= 2.1 && < 2.5+    -- , hspec-smallcheck >= 0.3 && < 0.5+    , QuickCheck >= 2.7 && < 2.10+    , lens >= 4.13 && < 4.16+    , template-haskell+    , mixed-types-num >= 0.3.1 && < 0.4+  if flag(MPFR)+    cpp-options: -DUSEMPFR+    if impl(ghc >= 7.10)+      if flag(MPFRRounded)+        build-depends:+          rounded == 0.1.*+        cpp-options: -DMPFRRounded+      else+        build-depends:+          hmpfr == 0.4.*,+          integer-gmp >= 1.0 && < 1.1+        cpp-options: -DHMPFR+    else+      build-depends:+        haskell-mpfr == 0.1.*,+        integer-gmp < 0.6+      cpp-options: -DHaskellMPFR+  else+    cpp-options: -DUseIReal+    -- TODO+  ghc-options:     -Wall -fno-warn-orphans+  extensions:+    RebindableSyntax,+    PostfixOperators,+    ScopedTypeVariables,+    DeriveGeneric,+    GeneralizedNewtypeDeriving,+    TypeFamilies,+    TypeOperators,+    ConstraintKinds,+    DefaultSignatures,+    MultiParamTypeClasses,+    FlexibleContexts,+    FlexibleInstances,+    UndecidableInstances+  if flag(MPFR)+    exposed-modules:+      AERN2.MP.UseMPFR.Float.Type+      AERN2.MP.UseMPFR.Float.RoundedAdaptor+      AERN2.MP.UseMPFR.Float.Arithmetic+      AERN2.MP.UseMPFR.Float.Conversions+      AERN2.MP.UseMPFR.Float.Operators+      AERN2.MP.UseMPFR.Float.Constants+      AERN2.MP.UseMPFR.Float.Tests+      AERN2.MP.UseMPFR.Float+      AERN2.MP.UseMPFR.Dyadic+      AERN2.MP.UseMPFR.ErrorBound+      AERN2.MP.UseMPFR.Ball.Type+      AERN2.MP.UseMPFR.Ball.Conversions+      AERN2.MP.UseMPFR.Ball.Comparisons+      AERN2.MP.UseMPFR.Ball.Field+      AERN2.MP.UseMPFR.Ball.Elementary+      AERN2.MP.UseMPFR.Ball.PreludeOps+      AERN2.MP.UseMPFR.Ball+  exposed-modules:+    AERN2.Utils.Bench+    AERN2.Normalize+    AERN2.Norm+    AERN2.MP.Precision+    AERN2.MP.Accuracy+    AERN2.MP.Enclosure+    AERN2.MP.Dyadic+    AERN2.MP.ErrorBound+    AERN2.MP.Ball.Tests+    AERN2.MP.Ball+    AERN2.MP++test-suite spec+  type:+      exitcode-stdio-1.0+  ghc-options:+      -Wall+  extensions:+    RebindableSyntax,+    PostfixOperators,+    ScopedTypeVariables,+    FlexibleContexts+  hs-source-dirs:+      test+  main-is:+      Spec.hs+  other-modules:+    AERN2.MP.BallSpec+    AERN2.MP.DyadicSpec+    AERN2.MP.UseMPFR.FloatSpec+  build-depends:+    base == 4.*+    -- , mixed-types-num >= 0.3.1 && < 0.4+    , aern2-mp+    , hspec >= 2.1 && < 2.5+    -- , hspec-smallcheck >= 0.3 && < 0.5+    , QuickCheck >= 2.7 && < 2.10
+ src/AERN2/MP.hs view
@@ -0,0 +1,33 @@+{-|+    Module      :  AERN2.MP+    Description :  Multiple-precision ball arithmetic+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Multiple-precision ball arithmetic+-}++module AERN2.MP+( module AERN2.Norm+, module AERN2.MP.Precision+, module AERN2.MP.Accuracy+, module AERN2.MP.UseMPFR.ErrorBound+, module AERN2.MP.Enclosure+, MPBall(..), CanBeMPBall, mpBall, CanBeMPBallP, mpBallP+)+where++-- import MixedTypesNumPrelude+-- import qualified Prelude as P+-- import Text.Printf++import AERN2.Norm+import AERN2.MP.Precision+import AERN2.MP.Accuracy+import AERN2.MP.UseMPFR.ErrorBound+import AERN2.MP.Enclosure+import AERN2.MP.UseMPFR.Ball
+ src/AERN2/MP/Accuracy.hs view
@@ -0,0 +1,300 @@+{-|+    Module      :  AERN2.MP.Accuracy+    Description :  Rough accuracy of an enclosure+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    A type for roughly measuring the accuracy of an enclosure.+-}+module AERN2.MP.Accuracy+    (Accuracy(NoInformation, Exact), bits, fromAccuracy,+     HasAccuracy(..),+     HasAccuracyGuide(..), CanSetAccuracyGuide(..), adjustAccuracyGuide,+     getFiniteAccuracy,+     ac2prec,+     CanReduceSizeUsingAccuracyGuide(..),+      specCanReduceSizeUsingAccuracyGuide,+     iterateUntilAccurate,+     convergentList2CauchySeq,+     seqByPrecision2CauchySeq,+     setPrecisionAtLeastAccuracy,+     HasApproximate(..))+where++import MixedTypesNumPrelude+import qualified Prelude as P++import Control.CollectErrors++import Data.Complex++import Test.Hspec+import Test.QuickCheck++import AERN2.Norm+import AERN2.MP.Precision++{- example -}++_example1 :: Accuracy+_example1 = 1 + 2*(bits 100)++{-| A non-negative Double value to serve as an error bound. Arithmetic is rounded towards +infinity. -}+data Accuracy = NoInformation | Bits { fromAccuracy :: Integer } | Exact+  deriving (P.Eq, P.Ord)++instance Arbitrary Accuracy where+  arbitrary =+    frequency+      [(int 1, pure Exact),+       (int 1, pure NoInformation),+       (int 8, Bits <$> arbitrary)]++instance Enum Accuracy where+    fromEnum NoInformation = minBound+    fromEnum (Bits i) = int i+    fromEnum Exact = maxBound+    toEnum i = Bits (integer i)++instance Bounded Accuracy where+    minBound = NoInformation+    maxBound = Exact++instance ConvertibleExactly Integer Accuracy where+  safeConvertExactly = Right . Bits+instance ConvertibleExactly Int Accuracy where+  safeConvertExactly = Right . Bits . integer+instance ConvertibleExactly Precision Accuracy where+  safeConvertExactly = Right . Bits . integer+instance ConvertibleExactly NormLog Accuracy where+  safeConvertExactly (NormBits b) = Right $ bits (-b)+  safeConvertExactly NormZero = Right Exact++bits :: (ConvertibleExactly t Accuracy) => t -> Accuracy+bits = convertExactly++instance Show Accuracy where+    show (NoInformation) = "NoInformation"+    show (Bits a) = "bits " ++ show a+    show (Exact) = "Exact"++instance HasEqAsymmetric Accuracy Accuracy+instance HasOrderAsymmetric Accuracy Accuracy+instance CanMinMaxAsymmetric Accuracy Accuracy++instance HasEqAsymmetric Accuracy Integer where+  equalTo = convertSecond equalTo+instance HasEqAsymmetric Integer Accuracy where+  equalTo = convertFirst equalTo+instance HasEqAsymmetric Accuracy Int where+  equalTo = convertSecond equalTo+instance HasEqAsymmetric Int Accuracy where+  equalTo = convertFirst equalTo++instance HasOrderAsymmetric Accuracy Integer where+  lessThan = convertSecond lessThan+  leq = convertSecond leq+instance HasOrderAsymmetric Integer Accuracy where+  lessThan = convertFirst lessThan+  leq = convertFirst leq+instance HasOrderAsymmetric Accuracy Int where+  lessThan = convertSecond lessThan+  leq = convertSecond leq+instance HasOrderAsymmetric Int Accuracy where+  lessThan = convertFirst lessThan+  leq = convertFirst leq++instance CanMinMaxAsymmetric Accuracy Integer where+    type MinMaxType Accuracy Integer = Accuracy+    min = convertSecond min+    max = convertSecond max+instance CanMinMaxAsymmetric Integer Accuracy where+    type MinMaxType Integer Accuracy = Accuracy+    min = convertFirst min+    max = convertFirst max+instance CanMinMaxAsymmetric Accuracy Int where+    type MinMaxType Accuracy Int = Accuracy+    min = convertSecond min+    max = convertSecond max+instance CanMinMaxAsymmetric Int Accuracy where+    type MinMaxType Int Accuracy = Accuracy+    min = convertFirst min+    max = convertFirst max++instance CanNeg Accuracy where+  negate NoInformation = Exact+  negate Exact = NoInformation+  negate (Bits a) = Bits (-a)++instance CanAddAsymmetric Accuracy Accuracy where+   add NoInformation _ = NoInformation+   add _ NoInformation = NoInformation+   add (Bits a) (Bits b) = Bits $ a + b+   add Exact _ = Exact+   add _ Exact = Exact++instance CanSub Accuracy Accuracy++--instance CanMulAsymmetric Accuracy Accuracy where+--    mulA NoInformation _ = NoInformation+--    mulA _ NoInformation = NoInformation+--    mulA (Bits a) (Bits b) = Bits $ a * b+--    mulA Exact _ = Exact+--    mulA _ Exact = Exact++instance CanMulAsymmetric Accuracy Integer where+    type MulType Accuracy Integer = Accuracy+    mul NoInformation _ = NoInformation+    mul (Bits a) i = Bits $ a * i+    mul Exact _ = Exact++instance CanMulAsymmetric Integer Accuracy where+    type MulType Integer Accuracy = Accuracy+    mul i a = mul a i++instance CanAddAsymmetric Accuracy Integer where+    type AddType Accuracy Integer = Accuracy+    add NoInformation _ = NoInformation+    add (Bits a) i = Bits $ a + i+    add Exact _ = Exact++instance CanAddAsymmetric Integer Accuracy where+    type AddType Integer Accuracy = Accuracy+    add i a = add a i++instance CanSub Accuracy Integer where+    type SubType Accuracy Integer = Accuracy+    sub NoInformation _ = NoInformation+    sub (Bits a) i = Bits $ a - i+    sub Exact _ = Exact++class HasAccuracy a where+  getAccuracy :: a -> Accuracy++instance (HasAccuracy a, SuitableForCE es) => HasAccuracy (CollectErrors es a) where+  getAccuracy (CollectErrors ma es) =+    case ma of+      Just a | not (hasCertainError es) -> getAccuracy 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 t => HasAccuracy (Complex t) where+  getAccuracy (a :+ i) =+    (getAccuracy a) `min` (getAccuracy i)++instance HasAccuracy t => HasAccuracy [t] where+  getAccuracy xs = foldl min Exact $ map getAccuracy 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++iterateUntilAccurate ::+  (HasAccuracy t) =>+  Accuracy ->+  (Precision -> Maybe t) ->+  [(Precision, Maybe t)]+iterateUntilAccurate ac =+  iterateUntilOK (ac2prec ac) $ \maybeResult ->+    case maybeResult of+      Just result -> getAccuracy result >= ac+      _ -> False++ac2prec :: Accuracy -> Precision+ac2prec ac =+  case ac of+    Bits b -> prec (max 2 $ b + 50)+    _ -> prec 100++seqByPrecision2CauchySeq ::+    (HasAccuracy t) =>+    (Precision -> t) -> (Accuracy -> t)+seqByPrecision2CauchySeq seqByPrecision ac =+    convergentList2CauchySeq 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+  where+  findAccurate [] =+    error "convergentList2CauchySeq: the sequence either converges too slowly or it does not converge"+  findAccurate (b : rest)+    | getAccuracy b >= ac = b+    | otherwise = findAccurate rest++{-|+    Change the precision so that+    it is at least as high as the supplied accuracy+    (assuming the accuracy is finite).+-}+setPrecisionAtLeastAccuracy :: (CanSetPrecision t) => Accuracy -> t -> t+setPrecisionAtLeastAccuracy acc b+    | p_b < p_acc = setPrecision p_acc b+    | otherwise = b+    where+    p_acc =+        case acc of+          Exact -> error $ "setPrecisionAtLeastAccuracy: cannot match Exact accuracy"+          NoInformation -> p_b+          _ -> 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++{-| An unsafe approximation of an enclosure or exact value,+    useful mainly for showing something brief and readable to humans.+-}+class HasApproximate t where+    type Approximate t+    getApproximate :: Accuracy -> t -> (Approximate t)
+ src/AERN2/MP/Ball.hs view
@@ -0,0 +1,28 @@+{-# LANGUAGE CPP #-}+{-|+    Module      :  AERN2.MP.Ball+    Description :  Arbitrary precision ball arithmetic+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Arbitrary precision ball arithmetic+-}++module AERN2.MP.Ball+  (+#ifdef USEMPFR+   module AERN2.MP.UseMPFR.Ball+#endif+  )+where++-- import MixedTypesNumPrelude+-- import qualified Prelude as P++#ifdef USEMPFR+import AERN2.MP.UseMPFR.Ball+#endif
+ src/AERN2/MP/Ball/Tests.hs view
@@ -0,0 +1,112 @@+{-# LANGUAGE CPP #-}+{-|+    Module      :  AERN2.MP.Ball.Tests+    Description :  Tests for operations on arbitrary precision balls+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Tests for operations on arbitrary precision balls.++    To run the tests using stack, execute:++    @+    stack test aern2-mp --test-arguments "-a 1000 -m MPBall"+    @+-}++module AERN2.MP.Ball.Tests+  (+    specMPBall, tMPBall+  )+where++import MixedTypesNumPrelude+-- import qualified Prelude as P+-- import Data.Ratio+-- import Text.Printf++import Test.Hspec+import Test.QuickCheck+-- import qualified Test.Hspec.SmallCheck as SC++-- import AERN2.Norm+import AERN2.MP.Precision++import AERN2.MP.UseMPFR.Ball.Type+-- import AERN2.MP.UseMPFR.Ball.Conversions ()+import AERN2.MP.UseMPFR.Ball.Comparisons ()+import AERN2.MP.UseMPFR.Ball.Field ()+import AERN2.MP.UseMPFR.Ball.Elementary ()++instance Arbitrary MPBall where+  arbitrary =+    do+      c <- finiteMPFloat+      e <- smallEB+      return (reducePrecionIfInaccurate $ MPBall c e)+    where+      smallEB =+        do+          e <- arbitrary+          if (mpBall e) !<! 10+            then return e+            else smallEB+      finiteMPFloat =+        do+          x <- arbitrary+          if isFinite x+            then return x+            else finiteMPFloat++{-|+  A runtime representative of type @MPBall@.+  Used for specialising polymorphic tests to concrete types.+-}+tMPBall :: T MPBall+tMPBall = T "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+      specHasEq tInt tMPBall tRational+      specCanTestZero tMPBall+      specHasOrderNotMixed tMPBall+      specHasOrder tInt tMPBall tRational+    describe "min/max/abs" $ do+      specCanNegNum tMPBall+      specCanAbs tMPBall+      specCanMinMaxNotMixed tMPBall+      specCanMinMax tMPBall tInteger tMPBall+    describe "ring" $ do+      specCanAddNotMixed tMPBall+      specCanAddSameType tMPBall+      specCanAdd tInt tMPBall tRational+      specCanAdd tInteger tMPBall tInt+      specCanSubNotMixed tMPBall+      specCanSub tMPBall tInteger+      specCanSub tInteger tMPBall+      specCanSub tMPBall tInt+      specCanSub tInt tMPBall+      specCanMulNotMixed tMPBall+      specCanMulSameType tMPBall+      specCanMul tInt tMPBall tRational+      -- specCanPow tMPBall tInteger+    describe "field" $ do+      specCanDivNotMixed tMPBall+      specCanDiv tInteger tMPBall+      specCanDiv tMPBall tInt+      specCanDiv tMPBall tRational+    describe "elementary" $ do+      specCanExpReal tMPBall+      specCanLogReal tMPBall+      specCanSqrtReal tMPBall+      specCanSinCosReal tMPBall
+ src/AERN2/MP/Dyadic.hs view
@@ -0,0 +1,28 @@+{-# LANGUAGE CPP #-}+{-|+    Module      :  AERN2.MP.Dyadic+    Description :  Dyadics with exact ring operations+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Arbitrary precision floating-point numbers with exact ring operations.+-}++module AERN2.MP.Dyadic+  (+#ifdef USEMPFR+   module AERN2.MP.UseMPFR.Dyadic+#endif+  )+where++-- import MixedTypesNumPrelude+-- import qualified Prelude as P++#ifdef USEMPFR+import AERN2.MP.UseMPFR.Dyadic+#endif
+ src/AERN2/MP/Enclosure.hs view
@@ -0,0 +1,297 @@+{-# LANGUAGE Arrows #-}+{-|+    Module      :  AERN2.MP.Enclosure+    Description :  Enclosure operations+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Enclosure classes and operations.+-}+module AERN2.MP.Enclosure+(+  IsBall(..)+  , IsInterval(..), intervalFunctionByEndpoints, intervalFunctionByEndpointsUpDown+  , CanTestContains(..), CanMapInside(..), specCanMapInside+  , CanIntersectAsymmetric(..), CanIntersect+  , CanIntersectCNBy, CanIntersectCNSameType+  , CanUnionAsymmetric(..), CanUnion+  , CanUnionCNBy, CanUnionCNSameType+  )+where++import MixedTypesNumPrelude+-- import qualified Prelude as P++import Control.Arrow++import Test.Hspec+import Test.QuickCheck++import Control.CollectErrors++import AERN2.MP.UseMPFR.ErrorBound+import AERN2.MP.Accuracy++{- ball-specific operations -}++class IsBall t where+  type CentreType t+  centre :: t -> CentreType t+  centreAsBallAndRadius :: t-> (t,ErrorBound)+  centreAsBall :: t -> t+  centreAsBall = fst . centreAsBallAndRadius+  radius :: t -> ErrorBound+  radius = snd . centreAsBallAndRadius+  updateRadius :: (ErrorBound -> ErrorBound) -> (t -> t)+  {-|  When the radius of the ball is implicitly contributed to by imprecision in the centre+     (eg if the centre is a polynomial with inexact coefficients), move all that imprecision+     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 v =+    updateRadius (+r) c+    where+    (c, r) = centreAsBallAndRadius v++{- interval-specific operations -}++class IsInterval i e where+  fromEndpoints :: e -> e -> i+  endpoints :: i -> (e,e)++{-|+    Computes a *monotone* ball function @f@ on intervals using the interval endpoints.+-}+intervalFunctionByEndpoints ::+  (IsInterval t 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)+  where+  (l,u) = endpoints x++{-|+    Computes a *monotone* ball function @f@ on intervals using the interval endpoints.+-}+intervalFunctionByEndpointsUpDown ::+  (IsInterval t e)+  =>+  (e -> e) {-^ @fDown@: a version of @f@ working on endpoints, rounded down -} ->+  (e -> e) {-^ @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)+  where+  (l,u) = endpoints x+++{- containment -}++class CanTestContains dom e where+  {-| Test if @e@ is inside @dom@. -}+  contains :: dom {-^ @dom@ -} -> e  {-^ @e@ -} -> Bool++class CanMapInside dom e where+  {-| Return some value contained in @dom@.+      The returned value does not have to equal the given @e@+      even if @e@ is already inside @dom@.+      All elements of @dom@ should be covered with roughly the same probability+      when calling this function for evenly distributed @e@'s.++      This function is intended mainly for generating values inside @dom@+      for randomised tests.+  -}+  mapInside :: dom {-^ @dom@ -} -> e  {-^ @e@ -} -> e++specCanMapInside ::+  (CanMapInside d e, CanTestContains d e+  , Arbitrary d, Arbitrary e, Show d, Show e)+  =>+  T d -> T e -> Spec+specCanMapInside (T dName :: T d) (T eName :: T e) =+  it ("CanMapInside " ++ dName ++ " " ++ eName) $ do+    property $+      \ (d :: d) (e :: e) ->+        contains d $ mapInside d e++{- intersection -}++type CanIntersect e1 e2 =+  (CanIntersectAsymmetric e1 e2, CanIntersectAsymmetric e1 e2+  , IntersectionType e1 e2 ~ IntersectionType e2 e1)++class CanIntersectAsymmetric e1 e2 where+  type IntersectionType e1 e2+  type IntersectionType e1 e2 = EnsureCN e1+  intersect :: e1 -> e2 -> IntersectionType e1 e2++type CanIntersectCNBy e1 e2 =+  (CanIntersect e1 e2, IntersectionType e1 e2 ~ EnsureCN e1+  , CanIntersect (EnsureCN e1) e2, IntersectionType (EnsureCN e1) e2 ~ EnsureCN e1)+type CanIntersectCNSameType e1 =+  (CanIntersectCNBy e1 e1+  , CanIntersect (EnsureCN e1) (EnsureCN e1), IntersectionType (EnsureCN e1) (EnsureCN e1) ~ EnsureCN e1)++instance CanIntersectAsymmetric Bool Bool where+  intersect b1 b2+    | b1 == b2 = cn b1+    | otherwise =+      noValueNumErrorCertainCN $ NumError "empty Boolean intersection"++instance+  (CanIntersectAsymmetric Bool b+  , CanEnsureCE es b+  , CanEnsureCE es (IntersectionType Bool b)+  , SuitableForCE es)+  =>+  CanIntersectAsymmetric Bool (CollectErrors es b)+  where+  type IntersectionType Bool (CollectErrors es b) =+    EnsureCE es (IntersectionType Bool b)+  intersect = lift2TLCE intersect++instance+  (CanIntersectAsymmetric a Bool+  , CanEnsureCE es a+  , CanEnsureCE es (IntersectionType a Bool)+  , SuitableForCE es)+  =>+  CanIntersectAsymmetric (CollectErrors es a) Bool+  where+  type IntersectionType (CollectErrors es  a) Bool =+    EnsureCE es (IntersectionType a Bool)+  intersect = lift2TCE intersect++instance+  (CanIntersectAsymmetric (Maybe a) b+  , CanEnsureCE es b+  , CanEnsureCE es (IntersectionType (Maybe a) b)+  , SuitableForCE es)+  =>+  CanIntersectAsymmetric (Maybe a) (CollectErrors es b)+  where+  type IntersectionType (Maybe a) (CollectErrors es b) =+    EnsureCE es (IntersectionType (Maybe a) b)+  intersect = lift2TLCE intersect++instance+  (CanIntersectAsymmetric a (Maybe b)+  , CanEnsureCE es a+  , CanEnsureCE es (IntersectionType a (Maybe b))+  , SuitableForCE es)+  =>+  CanIntersectAsymmetric (CollectErrors es a) (Maybe b)+  where+  type IntersectionType (CollectErrors es  a) (Maybe b) =+    EnsureCE es (IntersectionType a (Maybe b))+  intersect = lift2TCE intersect+++instance+  (CanIntersectAsymmetric a b+  , CanEnsureCN a, IntersectionType a b ~ EnsureCN a+  , CanEnsureCN (EnsureCN a)+  , CanEnsureCN b, EnsureCN b ~ EnsureCN a)+  =>+  CanIntersectAsymmetric (Maybe a) (Maybe 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)+++-- --- 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 (CollectErrors es e1) (CollectErrors es e2)+  where+  type IntersectionType (CollectErrors es e1) (CollectErrors es e2) =+    EnsureCE es (IntersectionType e1 e2)+  intersect = lift2CE intersect++{- union -}++type CanUnion e1 e2 =+  (CanUnionAsymmetric e1 e2, CanUnionAsymmetric e1 e2+  , UnionType e1 e2 ~ UnionType e2 e1)++class CanUnionAsymmetric e1 e2 where+  type UnionType e1 e2+  type UnionType e1 e2 = EnsureCN e1+  union :: e1 -> e2 -> UnionType e1 e2++type CanUnionCNBy e1 e2 =+  (CanUnion e1 e2, UnionType e1 e2 ~ EnsureCN e1+  , CanUnion (EnsureCN e1) e2, UnionType (EnsureCN e1) e2 ~ EnsureCN e1)++type CanUnionCNSameType e1 =+  (CanUnionCNBy e1 e1+  , CanUnion (EnsureCN e1) (EnsureCN e1), UnionType (EnsureCN e1) (EnsureCN e1) ~ EnsureCN e1)++instance+  (CanUnionAsymmetric e1 e2+  , CanEnsureCN e1, CanEnsureCN e2+  , CanEnsureCN (UnionType e1 e2))+  =>+  CanUnionAsymmetric (CN e1) (CN e2)+  -- a more general `CollectErrors` instance conflicts with the Arrow instance below+  where+  type UnionType (CN e1) (CN e2) =+    EnsureCN (UnionType e1 e2)+  union = lift2CE 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++instance (CanUnionCNSameType t, CanEnsureCN t) =>+  HasIfThenElse (Maybe Bool) 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
+ src/AERN2/MP/ErrorBound.hs view
@@ -0,0 +1,28 @@+{-# LANGUAGE CPP #-}+{-|+    Module      :  AERN2.MP.ErrorBound+    Description :  Fixed precision non-negative up-rounded floating-point numbers+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Fixed precision non-negative up-rounded floating-point numbers.+-}++module AERN2.MP.ErrorBound+#ifdef USEMPFR+  (+   module AERN2.MP.UseMPFR.ErrorBound+  )+#endif+where++-- import MixedTypesNumPrelude+-- import qualified Prelude as P++#ifdef USEMPFR+import AERN2.MP.UseMPFR.ErrorBound+#endif
+ src/AERN2/MP/Precision.hs view
@@ -0,0 +1,226 @@+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-|+    Module      :  AERN2.MP.Precision+    Description :  Floating-point precision+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Floating-point precision type and its operations+-}+module AERN2.MP.Precision+(+     Precision, prec+     , HasPrecision(..), CanSetPrecision(..), lowerPrecisionIfAbove, raisePrecisionIfBelow, specCanSetPrecision+     , defaultPrecision, maximumPrecision, standardPrecisions, precisionTimes2+     , iterateUntilOK+     , ConvertibleWithPrecision(..), convertP+     , convertPFirst, convertPSecond+)+where++import MixedTypesNumPrelude+import qualified Prelude as P+import Text.Printf++import Control.CollectErrors++import Data.Complex++import Data.Typeable++import Test.Hspec+import Test.QuickCheck++newtype Precision = Precision Integer+    deriving (P.Eq, P.Ord, P.Show, P.Enum, P.Num, P.Real, P.Integral, Typeable)++instance HasEqAsymmetric Precision Precision+instance HasOrderAsymmetric Precision Precision+instance CanMinMaxAsymmetric Precision Precision++instance ConvertibleExactly Precision Integer where+  safeConvertExactly (Precision p) = Right p++instance ConvertibleExactly Integer Precision where+  safeConvertExactly p+    | p < 2 = convError errmsg p+    | Precision p > maximumPrecision = convError errmsg p+      -- beware: if one removes "Precision" in the line above, it will type-check but loop+    | otherwise = Right $ Precision p+    where+    errmsg =+        "Precision must be between 2 and " ++ show maximumPrecision ++ " (given: p=" ++ show p ++ ")."++prec :: Integer -> Precision+prec = convertExactly++instance HasEqAsymmetric Precision Integer where+  equalTo p i = equalTo p (prec i)+instance HasEqAsymmetric Integer Precision where+  equalTo i p = equalTo (prec i) p+instance HasOrderAsymmetric Precision Integer where+  lessThan p i = lessThan p (prec i)+  leq p i = leq p (prec i)+instance HasOrderAsymmetric Integer Precision where+  lessThan i p = lessThan (prec i) p+  leq i p = leq (prec i) p++instance HasEqAsymmetric Precision Int where+  equalTo p i = equalTo p (prec (integer i))+instance HasEqAsymmetric Int Precision where+  equalTo i p = equalTo (prec (integer i)) p+instance HasOrderAsymmetric Precision Int where+  lessThan p i = lessThan p (prec (integer i))+  leq p i = leq p (prec (integer i))+instance HasOrderAsymmetric Int Precision where+  lessThan i p = lessThan (prec (integer i)) p+  leq i p = leq (prec (integer i)) p++instance CanAddAsymmetric Precision Precision+instance CanAddAsymmetric Integer Precision where+  type AddType Integer Precision = Precision+  add n (Precision p) = prec (n + p)+instance CanAddAsymmetric Precision Integer where+  type AddType Precision Integer = Precision+  add (Precision p) n = prec (n + p)++instance CanMulAsymmetric Precision Precision+instance CanMulAsymmetric Integer Precision where+  type MulType Integer Precision = Precision+  mul n (Precision p) = prec (n * p)+instance CanMulAsymmetric Precision Integer where+  type MulType Precision Integer = Precision+  mul (Precision p) n = prec (n * p)++class HasPrecision t where+    getPrecision :: t -> Precision++class (HasPrecision t) => CanSetPrecision t where+    setPrecision :: Precision -> t -> t++instance HasPrecision t => HasPrecision (Complex t) where+  getPrecision (a :+ i) =+    (getPrecision a) `min` (getPrecision i)++instance CanSetPrecision t => CanSetPrecision (Complex t) where+  setPrecision p (a :+ i) =+    (setPrecision p a) :+ (setPrecision p i)++instance HasPrecision t => HasPrecision (Maybe t) where+  getPrecision (Just v) = getPrecision v+  getPrecision Nothing = defaultPrecision+instance CanSetPrecision t => CanSetPrecision (Maybe t) where+  setPrecision p = fmap (setPrecision p)++instance HasPrecision Bool where+  getPrecision _ = defaultPrecision+instance CanSetPrecision Bool where+  setPrecision _ = id++instance HasPrecision t => HasPrecision (CollectErrors es t) where+  getPrecision vCE =+    case getMaybeValueCE 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 p x+  | getPrecision x > p = setPrecision p x+  | otherwise = x++raisePrecisionIfBelow :: (CanSetPrecision t) => Precision -> t -> t+raisePrecisionIfBelow p x+  | getPrecision x < p = setPrecision p x+  | otherwise = x++specCanSetPrecision ::+  (CanSetPrecision t, Arbitrary t, Show t, Testable prop)+  =>+  (T t) -> (t -> t -> prop) -> Spec+specCanSetPrecision (T typeName :: T t) check =+  describe (printf "CanSetPrecision %s" typeName) $ do+    it "set then get" $ do+      property $ \ (x :: t) (p :: Precision) ->+        let xP = setPrecision p x in+          p == getPrecision xP+    it "setPrecision x ~ x" $ do+      property $ \ (x :: t) (p :: Precision) ->+        let xP = setPrecision p x in+          check xP x++maximumPrecision :: Precision+maximumPrecision = Precision 5000000++defaultPrecision :: Precision+defaultPrecision = Precision 100++standardPrecisions :: Precision -> [Precision]+standardPrecisions (Precision initPrec0) =+    map (Precision . (+ initPrec)) $ aux 0 (max 2 (initPrec `P.div` 16))+    where+    initPrec = max 2 initPrec0+    aux j j'+        | Precision j <= maximumPrecision = j : (aux j' (j+j'))+        | otherwise = []++precisionTimes2 :: Precision -> Precision+precisionTimes2 (Precision p) = Precision (2*p)++iterateUntilOK ::+    Precision ->+    (a -> Bool) ->+    (Precision -> a) ->+    [(Precision, a)]+iterateUntilOK initPrec isOK fn =+    stopWhenAccurate ps+    where+--    fnWrap p =+--        unsafePerformIO $+--            catch (return $! Just $! fn p)+--                (\e -> let _ = e :: SomeException in return Nothing)+    ps = standardPrecisions initPrec+    stopWhenAccurate [] = []+    stopWhenAccurate (p : rest)+      | isOK result = [(p, result)]+      | otherwise = (p, result) : stopWhenAccurate rest+      where+      result = fn p++class ConvertibleWithPrecision t1 t2 where+  safeConvertP :: Precision -> t1 -> ConvertResult t2++convertP :: (ConvertibleWithPrecision t1 t2) => Precision -> t1 -> t2+convertP p a =+  case safeConvertP p a of+    Right v -> v+    Left err -> error (show err)++convertPFirst ::+  (ConvertibleWithPrecision t1 t2, HasPrecision t2) =>+  (t2 -> t2 -> c) -> (t1 -> t2 -> c)+convertPFirst = convertFirstUsing (\ q b -> convertP (getPrecision b) q)++convertPSecond ::+  (ConvertibleWithPrecision t2 t1, HasPrecision t1) =>+  (t1 -> t1 -> c) -> (t1 -> t2 -> c)+convertPSecond = convertSecondUsing (\ b q -> convertP (getPrecision b) q)+++instance Arbitrary Precision where+  arbitrary =+    sized $ \size -> choose (4,10+size) >>= return . prec++$(declForTypes+  [[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+  |]))
+ src/AERN2/MP/UseMPFR/Ball.hs view
@@ -0,0 +1,64 @@+{-|+    Module      :  AERN2.MP.UseMPFR.Ball+    Description :  Arbitrary precision ball arithmetic+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Arbitrary precision ball arithmetic+-}+module AERN2.MP.UseMPFR.Ball+(+  -- * Auxiliary types+  module AERN2.Norm+  , module AERN2.MP.Precision+  , module AERN2.MP.Accuracy+  , module AERN2.MP.UseMPFR.ErrorBound+  , module AERN2.MP.Enclosure+  -- * The Ball type+  , MPBall(..), CanBeMPBall, mpBall, CanBeMPBallP, mpBallP+  , reducePrecionIfInaccurate+  -- * Ball construction/extraction functions+  -- , endpointsMP, fromEndpointsMP+  -- * Ball operations (see also instances)+  , piBallP+  -- * Helpers for constructing ball functions+  , byEndpointsMP+  , fromApproxWithLipschitz+)+where++import MixedTypesNumPrelude+-- import qualified Prelude as P++import AERN2.Norm+import AERN2.MP.Precision+import AERN2.MP.Accuracy+import AERN2.MP.Enclosure++import AERN2.MP.UseMPFR.ErrorBound (ErrorBound, CanBeErrorBound, errorBound)++import AERN2.MP.UseMPFR.Ball.Type+import AERN2.MP.UseMPFR.Ball.Conversions ()+import AERN2.MP.UseMPFR.Ball.Comparisons+import AERN2.MP.UseMPFR.Ball.Field ()+import AERN2.MP.UseMPFR.Ball.Elementary+import AERN2.MP.UseMPFR.Ball.PreludeOps ()++instance Ring MPBall+instance Ring (CN MPBall)+instance Field MPBall+instance Field (CN MPBall)++instance OrderedRing MPBall+instance OrderedRing (CN MPBall)+instance OrderedField MPBall+instance OrderedField (CN MPBall)++instance OrderedCertainlyRing MPBall+instance OrderedCertainlyRing (CN MPBall)+instance OrderedCertainlyField MPBall+instance OrderedCertainlyField (CN MPBall)
+ src/AERN2/MP/UseMPFR/Ball/Comparisons.hs view
@@ -0,0 +1,381 @@+{-|+    Module      :  AERN2.MP.UseMPFR.Ball.Comparisons+    Description :  Comparisons of arbitrary precision dyadic balls+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Comparisons of arbitrary precision dyadic balls+-}+module AERN2.MP.UseMPFR.Ball.Comparisons+(+  -- * Auxiliary types+  module AERN2.Norm+  -- * Ball operations (see also instances)+  , reducePrecionIfInaccurate+  -- * Helpers for constructing ball functions+  , byEndpointsMP+)+where++import MixedTypesNumPrelude+-- import qualified Prelude as P++import Control.CollectErrors++import AERN2.Norm+import AERN2.MP.Dyadic (Dyadic)+import AERN2.MP.UseMPFR.Float (MPFloat)+-- import AERN2.MP.UseMPFR.Float.Operators+import AERN2.MP.Precision++import AERN2.MP.UseMPFR.Ball.Type+import AERN2.MP.UseMPFR.Ball.Conversions ()++{- comparisons -}++instance HasEqAsymmetric MPBall MPBall where+  type EqCompareType MPBall MPBall = Maybe Bool+  b1 `equalTo` b2 =   b1 >= b2 && b1 <= b2++instance HasEqAsymmetric MPBall Integer where+  type EqCompareType MPBall Integer = Maybe Bool+  b1 `equalTo` b2 =   b1 >= b2 && b1 <= b2+instance HasEqAsymmetric Integer MPBall where+  type EqCompareType Integer MPBall = Maybe Bool+  b1 `equalTo` b2 =   b1 >= b2 && b1 <= b2++instance HasEqAsymmetric MPBall Int where+  type EqCompareType MPBall Int = Maybe Bool+  b1 `equalTo` b2 =   b1 >= b2 && b1 <= b2+instance HasEqAsymmetric Int MPBall where+  type EqCompareType Int MPBall = Maybe Bool+  b1 `equalTo` b2 =   b1 >= b2 && b1 <= b2++instance HasEqAsymmetric MPBall Rational where+  type EqCompareType MPBall Rational = Maybe Bool+  b1 `equalTo` b2 =   b1 >= b2 && b1 <= b2+instance HasEqAsymmetric Rational MPBall where+  type EqCompareType Rational MPBall = Maybe Bool+  b1 `equalTo` b2 =   b1 >= b2 && b1 <= b2++instance HasEqAsymmetric MPBall Dyadic where+  type EqCompareType MPBall Dyadic = Maybe Bool+  b1 `equalTo` b2 =   b1 >= b2 && b1 <= b2+instance HasEqAsymmetric Dyadic MPBall where+  type EqCompareType Dyadic MPBall = Maybe Bool+  b1 `equalTo` b2 =   b1 >= b2 && b1 <= b2++instance+  (HasEqAsymmetric MPBall b+  , CanEnsureCE es b+  , CanEnsureCE es (EqCompareType MPBall b)+  , IsBool (EnsureCE es (EqCompareType MPBall b))+  , SuitableForCE es)+  =>+  HasEqAsymmetric MPBall (CollectErrors es  b)+  where+  type EqCompareType MPBall (CollectErrors es  b) =+    EnsureCE es (EqCompareType MPBall b)+  equalTo = lift2TLCE equalTo++instance+  (HasEqAsymmetric a MPBall+  , CanEnsureCE es a+  , CanEnsureCE es (EqCompareType a MPBall)+  , IsBool (EnsureCE es (EqCompareType a MPBall))+  , SuitableForCE es)+  =>+  HasEqAsymmetric (CollectErrors es a) MPBall+  where+  type EqCompareType (CollectErrors es  a) MPBall =+    EnsureCE es (EqCompareType a MPBall)+  equalTo = lift2TCE equalTo++instance HasOrderAsymmetric MPBall MPBall where+  type OrderCompareType MPBall MPBall = Maybe Bool+  lessThan b1 b2+    | r1 < l2 = Just True+    | r2 <= l1 = Just False+    | otherwise = Nothing+    where+    (l1, r1) = endpointsMP b1+    (l2, r2) = endpointsMP b2+  leq b1 b2+    | r1 <= l2 = Just True+    | r2 < l1 = Just False+    | otherwise = Nothing+    where+    (l1, r1) = endpointsMP b1+    (l2, r2) = endpointsMP b2++instance HasOrderAsymmetric Integer MPBall where+  type OrderCompareType Integer MPBall = Maybe Bool+  lessThan = convertFirst lessThan+  leq = convertFirst leq+instance HasOrderAsymmetric MPBall Integer where+  type OrderCompareType MPBall Integer = Maybe Bool+  lessThan = convertSecond lessThan+  leq = convertSecond leq++instance HasOrderAsymmetric Int MPBall where+  type OrderCompareType Int MPBall = Maybe Bool+  lessThan = convertFirst lessThan+  leq = convertFirst leq+instance HasOrderAsymmetric MPBall Int where+  type OrderCompareType MPBall Int = Maybe Bool+  lessThan = convertSecond lessThan+  leq = convertSecond leq++instance HasOrderAsymmetric Dyadic MPBall where+  type OrderCompareType Dyadic MPBall = Maybe Bool+  lessThan = convertFirst lessThan+  leq = convertFirst leq+instance HasOrderAsymmetric MPBall Dyadic where+  type OrderCompareType MPBall Dyadic = Maybe Bool+  lessThan = convertSecond lessThan+  leq = convertSecond leq++instance HasOrderAsymmetric MPBall Rational where+  type OrderCompareType MPBall Rational = Maybe Bool+  lessThan b1 q2+    | r1 < l2 = Just True+    | r2 <= l1 = Just False+    | otherwise = Nothing+    where+    (l1, r1) = endpointsMP b1+    l2 = q2+    r2 = q2+  leq b1 q2+    | r1 <= l2 = Just True+    | r2 < l1 = Just False+    | otherwise = Nothing+    where+    (l1, r1) = endpointsMP b1+    l2 = q2+    r2 = q2++instance HasOrderAsymmetric Rational MPBall where+  type OrderCompareType Rational MPBall = Maybe Bool+  lessThan q1 b2+    | r1 < l2 = Just True+    | r2 <= l1 = Just False+    | otherwise = Nothing+    where+    (l2, r2) = endpointsMP b2+    l1 = q1+    r1 = q1+  leq q1 b2+    | r1 <= l2 = Just True+    | r2 < l1 = Just False+    | otherwise = Nothing+    where+    (l2, r2) = endpointsMP b2+    l1 = q1+    r1 = q1++instance+  (HasOrderAsymmetric MPBall b+  , CanEnsureCE es b+  , CanEnsureCE es (OrderCompareType MPBall b)+  , IsBool (EnsureCE es (OrderCompareType MPBall b))+  , SuitableForCE es)+  =>+  HasOrderAsymmetric MPBall (CollectErrors es  b)+  where+  type OrderCompareType MPBall (CollectErrors es  b) =+    EnsureCE es (OrderCompareType MPBall b)+  lessThan = lift2TLCE lessThan+  leq = lift2TLCE leq+  greaterThan = lift2TLCE greaterThan+  geq = lift2TLCE geq++instance+  (HasOrderAsymmetric a MPBall+  , CanEnsureCE es a+  , CanEnsureCE es (OrderCompareType a MPBall)+  , IsBool (EnsureCE es (OrderCompareType a MPBall))+  , SuitableForCE es)+  =>+  HasOrderAsymmetric (CollectErrors es a) MPBall+  where+  type OrderCompareType (CollectErrors es  a) MPBall =+    EnsureCE es (OrderCompareType a MPBall)+  lessThan = lift2TCE lessThan+  leq = lift2TCE leq+  greaterThan = lift2TCE greaterThan+  geq = lift2TCE geq++instance CanTestZero MPBall+instance CanTestPosNeg MPBall++instance CanTestInteger MPBall where+  certainlyNotInteger b =+    (rN - lN) == 1 && lN !<! b && b !<! rN+    where+      (lN, rN) = integerBounds b+  certainlyIntegerGetIt b+    | rN == lN = Just lN+    | otherwise = Nothing+    where+      (lN, rN) = integerBounds b++instance CanMinMaxAsymmetric MPBall MPBall where+  min = byEndpointsMP min+  max = byEndpointsMP max++instance CanMinMaxAsymmetric MPBall Integer where+  type MinMaxType MPBall Integer = MPBall+  min = convertSecond min+  max = convertSecond max+instance CanMinMaxAsymmetric Integer MPBall where+  type MinMaxType Integer MPBall = MPBall+  min = convertFirst min+  max = convertFirst max++instance CanMinMaxAsymmetric MPBall Int where+  type MinMaxType MPBall Int = MPBall+  min = convertSecond min+  max = convertSecond max+instance CanMinMaxAsymmetric Int MPBall where+  type MinMaxType Int MPBall = MPBall+  min = convertFirst min+  max = convertFirst max++instance CanMinMaxAsymmetric MPBall Dyadic where+  type MinMaxType MPBall Dyadic = MPBall+  min = convertSecond min+  max = convertSecond max+instance CanMinMaxAsymmetric Dyadic MPBall where+  type MinMaxType Dyadic MPBall = MPBall+  min = convertFirst min+  max = convertFirst max++instance CanMinMaxAsymmetric MPBall Rational where+  type MinMaxType MPBall Rational = MPBall+  min = convertPSecond min+  max = convertPSecond max+instance CanMinMaxAsymmetric Rational MPBall where+  type MinMaxType Rational MPBall = MPBall+  min = convertPFirst min+  max = convertPFirst max++instance+  (CanMinMaxAsymmetric MPBall b+  , CanEnsureCE es b+  , CanEnsureCE es (MinMaxType MPBall b)+  , SuitableForCE es)+  =>+  CanMinMaxAsymmetric MPBall (CollectErrors es  b)+  where+  type MinMaxType MPBall (CollectErrors es  b) =+    EnsureCE es (MinMaxType MPBall b)+  min = lift2TLCE min+  max = lift2TLCE max++instance+  (CanMinMaxAsymmetric a MPBall+  , CanEnsureCE es a+  , CanEnsureCE es (MinMaxType a MPBall)+  , SuitableForCE es)+  =>+  CanMinMaxAsymmetric (CollectErrors es a) MPBall+  where+  type MinMaxType (CollectErrors es  a) MPBall =+    EnsureCE es (MinMaxType a MPBall)+  min = lift2TCE min+  max = lift2TCE max++{- intersection -}++instance CanIntersectAsymmetric MPBall MPBall where+  intersect a b+    | rL > rR =+        noValueNumErrorCertainCN $ NumError $ "intersect: empty intersection: " ++ show a ++ "; " ++ show b+    | otherwise = cn $ fromEndpointsMP rL rR+    where+    rL = max aL bL+    rR = min aR bR+    (aL,aR) = endpointsMP a+    (bL,bR) = endpointsMP b++instance+  (CanIntersectAsymmetric MPBall b+  , CanEnsureCE es b+  , CanEnsureCE es (IntersectionType MPBall b)+  , SuitableForCE es)+  =>+  CanIntersectAsymmetric MPBall (CollectErrors es b)+  where+  type IntersectionType MPBall (CollectErrors es b) =+    EnsureCE es (IntersectionType MPBall b)+  intersect = lift2TLCE intersect++instance+  (CanIntersectAsymmetric a MPBall+  , CanEnsureCE es a+  , CanEnsureCE es (IntersectionType a MPBall)+  , SuitableForCE es)+  =>+  CanIntersectAsymmetric (CollectErrors es a) MPBall+  where+  type IntersectionType (CollectErrors es  a) MPBall =+    EnsureCE es (IntersectionType a MPBall)+  intersect = lift2TCE intersect++{- union -}++instance CanUnionAsymmetric MPBall MPBall where+  union a b =+    case getMaybeValueCN (a `intersect` b) of+      Just _ -> prependErrorsCN [(ErrorPotential, err)] r+      _ -> prependErrorsCN [(ErrorCertain, err)] r+    where+    err = NumError $ "union of enclosures: not enclosing the same value"+    r = cn $ fromEndpointsMP rL rR+    rL = min aL bL+    rR = max aR bR+    (aL,aR) = endpointsMP a+    (bL,bR) = endpointsMP b+++instance+  (CanUnionAsymmetric MPBall b+  , CanEnsureCE es b+  , CanEnsureCE es (UnionType MPBall b)+  , SuitableForCE es)+  =>+  CanUnionAsymmetric MPBall (CollectErrors es b)+  where+  type UnionType MPBall (CollectErrors es b) =+    EnsureCE es (UnionType MPBall b)+  union = lift2TLCE union++instance+  (CanUnionAsymmetric a MPBall+  , CanEnsureCE es a+  , CanEnsureCE es (UnionType a MPBall)+  , SuitableForCE es)+  =>+  CanUnionAsymmetric (CollectErrors es a) MPBall+  where+  type UnionType (CollectErrors es  a) MPBall =+    EnsureCE es (UnionType a MPBall)+  union = lift2TCE union++{-|+  Computes an *increasing* ball fucntion @f@ from *exact* MPFR operations.+-}+byEndpointsMP ::+    (MPFloat -> MPFloat -> MPFloat) ->+    (MPBall -> MPBall -> MPBall)+byEndpointsMP op b1 b2 =+    fromEndpointsMP (l1 `op` l2) (r1 `op` r2)+    where+    (l1,r1) = endpointsMP b1+    (l2,r2) = endpointsMP b2++{-  random generation -}
+ src/AERN2/MP/UseMPFR/Ball/Conversions.hs view
@@ -0,0 +1,119 @@+{-|+    Module      :  AERN2.MP.UseMPFR.Ball.Conversions+    Description :  Conversions of arbitrary precision dyadic balls+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Conversions of arbitrary precision dyadic balls+-}+module AERN2.MP.UseMPFR.Ball.Conversions+(+  integerBounds+)+where++import MixedTypesNumPrelude+-- import qualified Prelude as P++import Data.Typeable+import Data.Convertible++import AERN2.MP.Dyadic (Dyadic, dyadic)+import qualified AERN2.MP.UseMPFR.Float as MPFloat+import AERN2.MP.UseMPFR.Float (mpFloat)+-- import AERN2.MP.UseMPFR.Float.Operators+import AERN2.MP.Precision+import qualified AERN2.MP.UseMPFR.ErrorBound as EB+import AERN2.MP.UseMPFR.ErrorBound (errorBound)++import AERN2.MP.UseMPFR.Ball.Type++{--- extracting from a ball ---}++instance HasIntegerBounds MPBall where+  integerBounds b =+    (floor l, ceiling r)+    where+      (l,r) = endpointsMP b++instance Convertible MPBall EB.ErrorBound where+  safeConvert b =+    Right (errorBound (max (abs l) (abs r)))+    where+    (l,r) = endpointsMP b++{--- constructing an exact ball ---}++instance ConvertibleExactly MPBall MPBall where+  safeConvertExactly = Right++instance ConvertibleExactly Dyadic MPBall where+  safeConvertExactly x = Right $ MPBall (mpFloat x) (errorBound 0)++instance ConvertibleExactly EB.ErrorBound MPBall where+  safeConvertExactly eb = Right $ MPBall (mpFloat eb) (errorBound 0)++instance+  (ConvertibleExactly c Dyadic, ConvertibleExactly e Dyadic+  , Show c, Show e, Typeable c, Typeable e)+  =>+  ConvertibleExactly (c, e) MPBall+  where+  safeConvertExactly (c,e)+    | isFinite b = Right b+    | otherwise = convError "too large to convert to MPBall" (c,e)+    where+    b = MPBall (mpFloat $ dyadic c) (errorBound $ mpFloat $ dyadic e)++instance ConvertibleExactly Integer MPBall where+  safeConvertExactly x+    | isFinite b = Right b+    | otherwise = convError "too large to convert to MPBall" x+    where+      b = MPBall (mpFloat x) (errorBound 0)++instance ConvertibleExactly Int MPBall where+  safeConvertExactly x = Right $ MPBall (mpFloat x) (errorBound 0)++{--- constructing a ball with a given precision ---}++instance ConvertibleWithPrecision Integer MPBall where+  safeConvertP p x+    | isFinite b = Right b+    | otherwise = convError ("too large to convert to MPBall with precision " ++ show p) x+    where+    b = MPBall xUp (xUp `EB.subMP` xDn)+    xUp = MPFloat.fromIntegerUp p x+    xDn = MPFloat.fromIntegerDown p x++instance ConvertibleWithPrecision Int MPBall where+  safeConvertP p = safeConvertP p . integer++instance ConvertibleWithPrecision Dyadic MPBall where+  safeConvertP p x+    | isFinite b = Right b+    | otherwise = convError ("too large to convert to MPBall with precision " ++ show p) x+    where+    b = mpBall x++instance ConvertibleWithPrecision Rational MPBall where+  safeConvertP p x+    | isFinite b = Right b+    | otherwise = convError ("too large to convert to MPBall with precision " ++ show p) x+    where+    b = MPBall xUp (xUp `EB.subMP` xDn)+    xUp = MPFloat.fromRationalUp p x+    xDn = MPFloat.fromRationalDown p x++instance ConvertibleWithPrecision (Rational, Rational) MPBall where+  safeConvertP p (x,e)+    | isFinite b = Right b+    | otherwise = convError ("too large to convert to MPBall with precision " ++ show p) x+    where+    b = MPBall xFlt (xe + eUp) -- beware, precision may be too high relative to accuracy+    (MPBall xFlt xe) = mpBallP p x+    eUp = errorBound e
+ src/AERN2/MP/UseMPFR/Ball/Elementary.hs view
@@ -0,0 +1,147 @@+{-# LANGUAGE CPP #-}+{-|+    Module      :  AERN2.MP.UseMPFR.Ball.Elementary+    Description :  Elementary operations on arbitrary precision dyadic balls+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Elementary operations on arbitrary precision dyadic balls+-}+module AERN2.MP.UseMPFR.Ball.Elementary+(+  -- * Ball operations (see also instances)+  piBallP+  -- * Helpers for constructing ball functions+  , fromApproxWithLipschitz+)+where++import MixedTypesNumPrelude+import qualified Prelude as P++import AERN2.Normalize++import AERN2.MP.Dyadic (Dyadic)+import qualified AERN2.MP.UseMPFR.Float as MPFloat+import AERN2.MP.UseMPFR.Float (MPFloat, mpFloat)+-- import AERN2.MP.UseMPFR.Float.Operators+import AERN2.MP.Precision+import qualified AERN2.MP.UseMPFR.ErrorBound as EB+import AERN2.MP.UseMPFR.ErrorBound (errorBound)++import AERN2.MP.UseMPFR.Ball.Type+import AERN2.MP.UseMPFR.Ball.Conversions ()+import AERN2.MP.UseMPFR.Ball.Comparisons ()+import AERN2.MP.UseMPFR.Ball.Field ()+++{- trigonometrics -}++piBallP :: Precision -> MPBall+piBallP p = MPBall piUp (piUp `EB.subMP` piDown)+  where+  piUp = MPFloat.piUp p+  piDown = MPFloat.piDown p++instance CanSinCos MPBall where+  sin = sinB 1+  cos = cosB 1++sinB :: Integer -> MPBall -> MPBall+sinB i x =+    -- increasingPrecisionUntilNotImproving (fromApproxWithLipschitz MPFloat.sinDown MPFloat.sinUp lip) x+    fromApproxWithLipschitz MPFloat.sinDown MPFloat.sinUp lip x+    where+    lip+        | i == 0 = mpFloat 1+        | otherwise = snd $ endpointsMP $ abs $ cosB (i - 1) x++cosB :: Integer -> MPBall -> MPBall+cosB i x =+    -- increasingPrecisionUntilNotImproving (fromApproxWithLipschitz MPFloat.cosDown MPFloat.cosUp lip) x+    fromApproxWithLipschitz MPFloat.cosDown MPFloat.cosUp lip x+    where+    lip+        | i == 0 = mpFloat 1+        | otherwise = snd $ endpointsMP $ abs $ sinB (i - 1) x++-- increasingPrecisionUntilNotImproving :: (MPBall -> MPBall) -> (MPBall -> MPBall)+-- increasingPrecisionUntilNotImproving f x =+--   waitUntilNotImproving $ map aux (precisions xPrec (xPrec*2))+--   where+--   xPrec = getPrecision x+--   precisions p1 p2 = p1 : (precisions p2 (p1 + p2))+--   aux p = f $ setPrecision p x+--   waitUntilNotImproving xx@(x1:_) = aux2 (getAccuracy x1) xx+--   waitUntilNotImproving _ = error "AERN2.MP.UseMPFR.Ball.Elementary: internal error in increasingPrecisionUntilNotImproving"+--   aux2 x1AC (x1:x2:rest)+--     | x1AC < x2AC = aux2 x2AC (x2:rest)+--     | otherwise = x1+--     where+--     x2AC = getAccuracy x2+--   aux2 _ _ = error "AERN2.MP.UseMPFR.Ball.Elementary: internal error in increasingPrecisionUntilNotImproving"++{- exp, log, power -}++instance CanExp MPBall where+  exp = intervalFunctionByEndpointsUpDown MPFloat.expDown MPFloat.expUp++instance CanLog MPBall where+  type LogType MPBall = CN MPBall+  log x+    | x !>! 0 =+        cn $ intervalFunctionByEndpointsUpDown MPFloat.logDown MPFloat.logUp x+    | x !<=! 0 = noValueNumErrorCertainCN err+    | otherwise = noValueNumErrorPotentialCN err+    where+    err = OutOfRange $ "log: argument must be > 0: " ++ show x++instance CanPow MPBall MPBall where+  powNoCN b e = (~!) $ pow b e+  pow = powUsingExpLog (mpBall 0) (mpBall 1)++instance CanPow MPBall Dyadic where+  powNoCN b e = (~!) $ pow b e+  pow b e = powUsingExpLog (mpBall 0) (mpBall 1) b (mpBall e)++instance CanPow MPBall Rational where+  powNoCN b e = (~!) $ pow b e+  pow b e = powUsingExpLog (mpBall 0) (mpBall 1) b (mpBallP (getPrecision b) e)++instance CanSqrt MPBall where+  type SqrtType MPBall = CN MPBall+  sqrt x+    | x !>=! 0 = cn $ aux x+    | x !<! 0 = noValueNumErrorCertainCN err+    | otherwise = prependErrorsCN [(ErrorPotential, err)] $ cn $ aux (max 0 x)+    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++{- generic methods for computing real functions from MPFR-approximations -}++{-|+    Computes a real function @f@ from correctly rounded MPFR-approximations and a number @lip@ which is a+    Lipschitz constant for @f@, i.e. @|f(x) - f(y)| <= lip * |x - y|@ for all @x@,@y@.+-}+fromApproxWithLipschitz ::+    (MPFloat -> MPFloat) {-^ @fDown@: a version of @f@ on MPFloat rounding *downwards* -} ->+    (MPFloat -> MPFloat) {-^ @fUp@: a version of @f@ on MPFloat rounding *upwards* -} ->+    MPFloat {-^ @lip@ a Lipschitz constant for @f@, @lip > 0@ -} ->+    (MPBall -> MPBall) {-^ @f@ on MPBall rounding *outwards* -}+fromApproxWithLipschitz fDown fUp lip _x@(MPBall xc xe) =+    normalize $ MPBall fxc err+    where+    fxl = fDown xc+    fxu = fUp xc+    (MPBall fxc fxe) =+      setPrecision (getPrecision xc) $ -- beware, some MPFR functions increase precision, eg sine and cosine+        fromEndpointsMP fxl fxu+    err = (errorBound lip) * xe  +  fxe
+ src/AERN2/MP/UseMPFR/Ball/Field.hs view
@@ -0,0 +1,344 @@+{-# LANGUAGE TemplateHaskell #-}+{-|+    Module      :  AERN2.MP.UseMPFR.Ball.Field+    Description :  Field operations on arbitrary precision dyadic balls+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Field operations on arbitrary precision dyadic balls+-}+module AERN2.MP.UseMPFR.Ball.Field+()+where++import MixedTypesNumPrelude+-- import qualified Prelude as P++import Control.CollectErrors++import AERN2.Normalize++import AERN2.MP.Dyadic (Dyadic)+import AERN2.MP.UseMPFR.Float (mpFloat)+import AERN2.MP.UseMPFR.Float.Operators+import AERN2.MP.Precision+import qualified AERN2.MP.UseMPFR.ErrorBound as EB++import AERN2.MP.UseMPFR.Ball.Type+import AERN2.MP.UseMPFR.Ball.Conversions ()+import AERN2.MP.UseMPFR.Ball.Comparisons ()++{- addition -}++instance CanAddAsymmetric MPBall MPBall where+  type AddType MPBall MPBall = MPBall+  add (MPBall x1 e1) (MPBall x2 e2) =+    normalize $ MPBall sumUp ((sumUp `EB.subMP` sumDn) + e1 + e2)+    where+    sumUp = x1 +^ x2+    sumDn = x1 +. x2++instance CanAddAsymmetric MPBall Int where+  type AddType MPBall Int = MPBall+  add = convertSecond add+instance CanAddAsymmetric Int MPBall where+  type AddType Int MPBall = MPBall+  add = convertFirst add++instance CanAddAsymmetric MPBall Integer where+  type AddType MPBall Integer = MPBall+  add = convertSecond add+instance CanAddAsymmetric Integer MPBall where+  type AddType Integer MPBall = MPBall+  add = convertFirst add++instance CanAddAsymmetric MPBall Dyadic where+  type AddType MPBall Dyadic = MPBall+  add = convertSecond add+instance CanAddAsymmetric Dyadic MPBall where+  type AddType Dyadic MPBall = MPBall+  add = convertFirst add++instance CanAddAsymmetric MPBall Rational where+  type AddType MPBall Rational = MPBall+  add = convertPSecond add+instance CanAddAsymmetric Rational MPBall where+  type AddType Rational MPBall = MPBall+  add = convertPFirst add++instance+  (CanAddAsymmetric MPBall b+  , CanEnsureCE es b+  , CanEnsureCE es (AddType MPBall b)+  , SuitableForCE es)+  =>+  CanAddAsymmetric MPBall (CollectErrors es  b)+  where+  type AddType MPBall (CollectErrors es  b) =+    EnsureCE es (AddType MPBall b)+  add = lift2TLCE add++instance+  (CanAddAsymmetric a MPBall+  , CanEnsureCE es a+  , CanEnsureCE es (AddType a MPBall)+  , SuitableForCE es)+  =>+  CanAddAsymmetric (CollectErrors es a) MPBall+  where+  type AddType (CollectErrors es  a) MPBall =+    EnsureCE es (AddType a MPBall)+  add = lift2TCE add++{- subtraction -}++instance CanSub MPBall MPBall++instance CanSub MPBall Integer+instance CanSub Integer MPBall++instance CanSub MPBall Int+instance CanSub Int MPBall++instance CanSub MPBall Rational+instance CanSub Rational MPBall++instance CanSub MPBall Dyadic+instance CanSub Dyadic MPBall++instance+  (CanSub MPBall b+  , CanEnsureCE es b+  , CanEnsureCE es (SubType MPBall b)+  , SuitableForCE es)+  =>+  CanSub MPBall (CollectErrors es  b)+  where+  type SubType MPBall (CollectErrors es  b) =+    EnsureCE es (SubType MPBall b)+  sub = lift2TLCE sub++instance+  (CanSub a MPBall+  , CanEnsureCE es a+  , CanEnsureCE es (SubType a MPBall)+  , SuitableForCE es)+  =>+  CanSub (CollectErrors es a) MPBall+  where+  type SubType (CollectErrors es  a) MPBall =+    EnsureCE es (SubType a MPBall)+  sub = lift2TCE sub++{- multiplication -}++instance CanMulAsymmetric MPBall MPBall where+  mul (MPBall x1 e1) (MPBall x2 e2) =+    normalize $ MPBall x12Up (e12 + e1*(abs x2) + e2*(abs x1) + e1*e2)+      -- the mixed operations above automatically convert+      -- MPFloat to ErrorBound, checking non-negativity+    where+    x12Up = x1 *^ x2+    x12Down = x1 *. x2+    e12 = x12Up -^ x12Down++instance CanMulAsymmetric MPBall Int where+  type MulType MPBall Int = MPBall+  mul = convertSecond mul+instance CanMulAsymmetric Int MPBall where+  type MulType Int MPBall = MPBall+  mul = convertFirst mul++instance CanMulAsymmetric MPBall Integer where+  type MulType MPBall Integer = MPBall+  mul = convertSecond mul+instance CanMulAsymmetric Integer MPBall where+  type MulType Integer MPBall = MPBall+  mul = convertFirst mul++instance CanMulAsymmetric MPBall Dyadic where+  type MulType MPBall Dyadic = MPBall+  mul = convertSecond mul+instance CanMulAsymmetric Dyadic MPBall where+  type MulType Dyadic MPBall = MPBall+  mul = convertFirst mul++instance CanMulAsymmetric MPBall Rational where+  type MulType MPBall Rational = MPBall+  mul = convertPSecond mul+instance CanMulAsymmetric Rational MPBall where+  type MulType Rational MPBall = MPBall+  mul = convertPFirst mul++instance+  (CanMulAsymmetric MPBall b+  , CanEnsureCE es b+  , CanEnsureCE es (MulType MPBall b)+  , SuitableForCE es)+  =>+  CanMulAsymmetric MPBall (CollectErrors es  b)+  where+  type MulType MPBall (CollectErrors es  b) =+    EnsureCE es (MulType MPBall b)+  mul = lift2TLCE mul++instance+  (CanMulAsymmetric a MPBall+  , CanEnsureCE es a+  , CanEnsureCE es (MulType a MPBall)+  , SuitableForCE es)+  =>+  CanMulAsymmetric (CollectErrors es a) MPBall+  where+  type MulType (CollectErrors es  a) MPBall =+    EnsureCE es (MulType a MPBall)+  mul = lift2TCE mul+++{- division -}++instance CanDiv MPBall MPBall where+  type DivTypeNoCN MPBall MPBall = MPBall+  divideNoCN b1 b2 = (~!) (divide b1 b2)+  type DivType MPBall MPBall = CN MPBall+  divide (MPBall x1 e1) b2@(MPBall x2 e2)+    | isCertainlyNonZero b2 =+        cn $ normalize $ MPBall x12Up err+    | isCertainlyZero b2 =+        noValueNumErrorCertainCN DivByZero+    | otherwise =+        noValueNumErrorPotentialCN DivByZero+    where+    x12Up = x1 /^ x2+    x12Down = x1 /. x2+    x12AbsUp = (abs x12Up) `max` (abs x12Down)+    e12 = x12Up -^ x12Down+    err =+        ((e12 *^ (abs x2)) -- e12 * |x2|+         ++         e1+         ++         (e2 * x12AbsUp) -- e2 * |x|+        )+        *+        ((mpFloat 1) /^ ((abs x2) -. (mpFloat e2)))+            -- 1/(|x2| - e2) rounded upwards+{-+A derivation of the above formula for an upper bound on the error:++    * e =+        * = max ( (x1 ± e1) / (x2 ± e2) - x )+        * = max ( ( x1 ± e1 - (x*(x2 ± e2) ) / (x2 ± e2) )+        * ≤ max ( ( x1 ± e1 - ((x1/x2) ± e12)x2 ± x*e2 ) / (x2 ± e2) )+        * = max ( ( x1 ± e1 - x1 ± e12*x2 ± x*e2 ) / (x2 ± e2) )+        * = max ( ( ± e1 ± e12*x2 ± x*e2 ) / (x2 ± e2) )+        * ≤ (e1 + e12*|x2| + |x|*e2 ) / (|x2| - e2)+        * ≤ (e1 +^ e12*^|x2| +^ |x|*^e2 ) /^ (|x2| -. e2)+-}++$(declForTypes+  [[t| Integer |], [t| Int |], [t| Dyadic |]]+  (\ t -> [d|+    instance CanDiv MPBall $t where+      type DivType MPBall $t = CN MPBall+      divide = convertSecond divide+      type DivTypeNoCN MPBall $t = MPBall+      divideNoCN = convertSecond divideNoCN+    instance CanDiv $t MPBall where+      type DivType $t MPBall = CN MPBall+      divide = convertFirst divide+      type DivTypeNoCN $t MPBall = MPBall+      divideNoCN = convertFirst divideNoCN+  |]))++instance CanDiv Dyadic Dyadic where+  type DivTypeNoCN Dyadic Dyadic = MPBall+  divideNoCN a b = divideNoCN (mpBall a) (mpBall b)+  divide a b = divide (mpBall a) (mpBall b)++instance CanDiv MPBall Rational where+  type DivTypeNoCN MPBall Rational = MPBall+  divideNoCN = convertPSecond divideNoCN+  divide = convertPSecond divide+instance CanDiv Rational MPBall where+  type DivTypeNoCN Rational MPBall = MPBall+  divideNoCN = convertPFirst divideNoCN+  divide = convertPFirst divide++instance+  (CanDiv MPBall b+  , CanEnsureCE es b+  , CanEnsureCE es (DivType MPBall b)+  , CanEnsureCE es (DivTypeNoCN MPBall b)+  , SuitableForCE es)+  =>+  CanDiv MPBall (CollectErrors es  b)+  where+  type DivType MPBall (CollectErrors es  b) =+    EnsureCE es (DivType MPBall b)+  divide = lift2TLCE divide+  type DivTypeNoCN MPBall (CollectErrors es  b) =+    EnsureCE es (DivTypeNoCN MPBall b)+  divideNoCN = lift2TLCE divideNoCN++instance+  (CanDiv a MPBall+  , CanEnsureCE es a+  , CanEnsureCE es (DivType a MPBall)+  , CanEnsureCE es (DivTypeNoCN a MPBall)+  , SuitableForCE es)+  =>+  CanDiv (CollectErrors es a) MPBall+  where+  type DivType (CollectErrors es  a) MPBall =+    EnsureCE es (DivType a MPBall)+  divide = lift2TCE divide+  type DivTypeNoCN (CollectErrors es  a) MPBall =+    EnsureCE es (DivTypeNoCN a MPBall)+  divideNoCN = lift2TCE divideNoCN++{- integer power -}++instance CanPow MPBall Integer where+  powNoCN b e = (~!) $ powUsingMulRecip (mpBall 1) b e+  pow = powUsingMulRecip (mpBall 1)++instance CanPow MPBall Int where+  powNoCN b e = (~!) $ powUsingMulRecip (mpBall 1) b e+  pow = powUsingMulRecip (mpBall 1)++instance+  (CanPow MPBall b+  , CanEnsureCE es b+  , CanEnsureCE es (PowType MPBall b)+  , CanEnsureCE es (PowTypeNoCN MPBall b)+  , SuitableForCE es)+  =>+  CanPow MPBall (CollectErrors es  b)+  where+  type PowTypeNoCN MPBall (CollectErrors es  b) =+    EnsureCE es (PowTypeNoCN MPBall b)+  powNoCN = lift2TLCE powNoCN+  type PowType MPBall (CollectErrors es  b) =+    EnsureCE es (PowType MPBall b)+  pow = lift2TLCE pow++instance+  (CanPow a MPBall+  , CanEnsureCE es a+  , CanEnsureCE es (PowType a MPBall)+  , CanEnsureCE es (PowTypeNoCN a MPBall)+  , SuitableForCE es)+  =>+  CanPow (CollectErrors es a) MPBall+  where+  type PowTypeNoCN (CollectErrors es  a) MPBall =+    EnsureCE es (PowTypeNoCN a MPBall)+  powNoCN = lift2TCE powNoCN+  type PowType (CollectErrors es  a) MPBall =+    EnsureCE es (PowType a MPBall)+  pow = lift2TCE pow
+ src/AERN2/MP/UseMPFR/Ball/PreludeOps.hs view
@@ -0,0 +1,76 @@+{-# LANGUAGE CPP #-}+{-|+    Module      :  AERN2.MP.UseMPFR.Ball.PreludeOps+    Description :  Instances of Prelude.Num etc+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Instances of Prelude classes Eq, Ord, Num etc+-}+module AERN2.MP.UseMPFR.Ball.PreludeOps+(+)+where++import MixedTypesNumPrelude+import qualified Prelude as P++import AERN2.MP.Dyadic (dyadic)++import AERN2.MP.UseMPFR.Ball.Type+import AERN2.MP.UseMPFR.Ball.Conversions ()+import AERN2.MP.UseMPFR.Ball.Comparisons ()+import AERN2.MP.UseMPFR.Ball.Field ()+import AERN2.MP.UseMPFR.Ball.Elementary++{- Instances of Prelude numerical classes provided for convenient use outside AERN2+   and also because Template Haskell translates (-x) to (Prelude.negate x) -}++instance P.Eq MPBall where+  a == b =+    case a == b of+      Just t -> t+      _ ->+        error "Failed to decide equality of MPBalls.  If you switch to MixedTypesNumPrelude instead of Prelude, comparison of MPBalls returns Maybe Bool instead of Bool."++instance P.Ord MPBall where+  compare a b =+    case (a < b, a == b, a > b) of+      (Just True, _, _) -> P.LT+      (_, Just True, _) -> P.EQ+      (_, _, Just True) -> P.GT+      _ ->+        error "Failed to decide order of MPBalls.  If you switch to MixedTypesNumPrelude instead of Prelude, comparison of MPBalls returns Maybe Bool instead of Bool."++instance P.Num MPBall where+    fromInteger = convertExactly+    negate = negate+    (+) = (+)+    (*) = (*)+    abs = abs+    signum = error "Prelude.signum not implemented for MPBall"++instance P.Fractional MPBall where+    fromRational = convertExactly . dyadic -- will work only for dyadic rationals+    recip = (~!) . recip+    (/) = (/!)++instance P.Floating MPBall where+    pi = error "MPBall: no pi :: MPBall, use pi ? (bitsS n) instead"+    sqrt = (~!) . sqrt+    exp = exp+    sin = sin+    cos = cos+    log = (~!) . log+    atan = error "MPBall: atan not implemented yet"+    atanh = error "MPBall: atanh not implemented yet"+    asin = error "MPBall: asin not implemented yet"+    acos = error "MPBall: acos not implemented yet"+    sinh = error "MPBall: sinh not implemented yet"+    cosh = error "MPBall: cosh not implemented yet"+    asinh = error "MPBall: asinh not implemented yet"+    acosh = error "MPBall: acosh not implemented yet"
+ src/AERN2/MP/UseMPFR/Ball/Type.hs view
@@ -0,0 +1,276 @@+{-# LANGUAGE TemplateHaskell #-}+{-|+    Module      :  AERN2.MP.UseMPFR.Ball.Type+    Description :  Arbitrary precision dyadic balls+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Arbitrary precision dyadic balls+-}+module AERN2.MP.UseMPFR.Ball.Type+(+  -- * Auxiliary types+  module AERN2.MP.Precision+  , module AERN2.MP.Accuracy+  , module AERN2.MP.Enclosure+  -- * The Ball type+  , MPBall(..), CanBeMPBall, mpBall, CanBeMPBallP, mpBallP+  , reducePrecionIfInaccurate+  -- * Ball construction/extraction functions+  , endpointsMP, fromEndpointsMP+)+where++import MixedTypesNumPrelude+-- import qualified Prelude as P++import Control.CollectErrors++import GHC.Generics (Generic)++import Text.Printf++import AERN2.Normalize++import AERN2.Norm++import AERN2.MP.Dyadic+import qualified AERN2.MP.UseMPFR.Float as MPFloat+import AERN2.MP.UseMPFR.Float (MPFloat, mpFloat)+import AERN2.MP.UseMPFR.Float.Operators+import AERN2.MP.Precision+import AERN2.MP.Accuracy+import qualified AERN2.MP.UseMPFR.ErrorBound as EB+import AERN2.MP.UseMPFR.ErrorBound (ErrorBound, errorBound)+import AERN2.MP.Enclosure++data MPBall = MPBall+  { ball_value :: MPFloat+  , ball_error :: ErrorBound+  }+  -- { ball_value :: {-# UNPACK #-} ! MPFloat+  -- , ball_error :: {-# UNPACK #-} ! ErrorBound+  -- }+  deriving (Generic)++instance Show MPBall+    where+    show b@(MPBall x _e) =+      -- printf "[%s ± %s](prec=%s)" (show x) (showAC $ getAccuracy b) (show $ integer $ getPrecision b)+      printf "[%s ± %s]" (show x) (showAC $ getAccuracy b)+      -- "[" ++ show x ++ " ± " ++ show e ++ "](prec=" ++ (show $ integer $ getPrecision x) ++ ")"+      where+      showAC Exact = "0"+      showAC NoInformation = "oo"+      showAC ac = "<2^(" ++ show (negate $ fromAccuracy ac) ++ ")"+++instance (SuitableForCE es) => CanEnsureCE es MPBall where++-- instance CanTestValid MPBall where+--   isValid = isFinite++instance CanTestNaN MPBall where+  isNaN = not . isFinite+instance CanTestFinite MPBall where+  isInfinite = const False+  isFinite (MPBall x e) = isFinite x && isFinite (mpFloat e)++instance CanNormalize MPBall where+  normalize b+    | isFinite b =+        b+        -- reducePrecionIfInaccurate b+    | otherwise = error $ "invalid MPBall: " ++ show b++{-|+    Reduce the precision of the ball centre if the+    accuracy of the ball is poor.++    More precisely, reduce the precision of the centre+    so that the ulp is approximately (radius / 1024),+    unless the ulp is already lower than this.+-}+reducePrecionIfInaccurate :: MPBall -> MPBall+reducePrecionIfInaccurate b@(MPBall x _) =+    case (bAcc, bNorm) of+        (Exact, _) -> b+        (_, NormZero) -> b+        _ | p_e_nb < p_x -> setPrecision p_e_nb b+        _ -> b+    where+    bAcc = getAccuracy b+    bNorm = getNormLog b+    p_x = getPrecision x+    p_e_nb = prec $ max 2 (10 + nb + fromAccuracy bAcc)+    (NormBits nb) = bNorm++instance CanTestContains MPBall MPBall where+  contains (MPBall xLarge eLarge) (MPBall xSmall eSmall) =+    xLargeDy - eLargeDy <= xSmallDy - eSmallDy+    &&+    xSmallDy + eSmallDy <= xLargeDy + eLargeDy+    where+    xLargeDy = dyadic xLarge+    eLargeDy = dyadic eLarge+    xSmallDy = dyadic xSmall+    eSmallDy = dyadic eSmall++$(declForTypes+  [[t| Integer |], [t| Int |], [t| Rational |], [t| Dyadic |]]+  (\ t -> [d|+    instance CanTestContains MPBall $t where+      contains (MPBall c e) x =+        l <= x && x <= r+        where+        l = cDy - eDy+        r = cDy + eDy+        cDy = dyadic c+        eDy = dyadic e+  |]))++{- ball construction/extraction functions -}++instance IsInterval MPBall MPFloat where+  fromEndpoints l u+    | u < l = fromEndpoints u l+    | otherwise =+      MPBall (mpFloat cDy) (errorBound $ mpFloat eDy)+      where+      lDy = dyadic l+      uDy = dyadic u+      cDy = (lDy + uDy) * (dyadic 0.5)+      eDy = (uDy - cDy) `max` (cDy - lDy)+  endpoints (MPBall x e) = (mpFloat lDy, mpFloat uDy)+      where+      xDy = dyadic x+      eDy = dyadic e+      lDy   = xDy - eDy+      uDy   = xDy + eDy++fromEndpointsMP :: MPFloat -> MPFloat -> MPBall+fromEndpointsMP = fromEndpoints++endpointsMP :: MPBall -> (MPFloat, MPFloat)+endpointsMP = endpoints++instance IsInterval MPBall MPBall where+  fromEndpoints l r = -- works as union even when r < l+      fromEndpointsMP lMP uMP+      where+      lMP = min llMP rlMP+      uMP = max luMP ruMP+      (llMP, luMP) = endpointsMP l+      (rlMP, ruMP) = endpointsMP r+  endpoints x = (l,u)+      where+      l = MPBall lMP (errorBound 0)+      u = MPBall uMP (errorBound 0)+      (lMP, uMP) = endpointsMP x++instance IsBall MPBall where+  type CentreType MPBall = Dyadic+  centre (MPBall cMP _e) = dyadic cMP+  centreAsBallAndRadius x = (cB,e)+    where+    (MPBall cMP e) = x+    cB = MPBall cMP (errorBound 0)+  radius (MPBall _ e) = e+  updateRadius updateFn (MPBall c e) = MPBall c (updateFn e)++{--- constructing a ball with a given precision ---}++type CanBeMPBallP t = (ConvertibleWithPrecision t MPBall)++mpBallP :: (CanBeMPBallP t) => Precision -> t -> MPBall+mpBallP = convertP+++{--- constructing an exact ball ---}++type CanBeMPBall t = ConvertibleExactly t MPBall++mpBall :: (CanBeMPBall t) => t -> MPBall+mpBall = convertExactly++{-- extracting approximate information about a ball --}++instance HasAccuracy MPBall where+    getAccuracy = getAccuracy . ball_error++instance CanReduceSizeUsingAccuracyGuide MPBall where+  reduceSizeUsingAccuracyGuide acGuide b@(MPBall x _e) =+    case acGuide of+      NoInformation -> lowerPrecisionIfAbove (prec 2) b+      _ | getAccuracy b > acGuide -> tryPrec newPrec+      _ -> b+    where+    tryPrec p+      | getAccuracy bP >= acGuide = bP+      | otherwise = tryPrec (p + 10)+      where+      bP = lowerPrecisionIfAbove p b+    queryBits = 1 + fromAccuracy acGuide+    newPrec =+      case (getNormLog x) of+        NormBits xNormBits ->+          prec (max 2 (queryBits + xNormBits + 2))+        NormZero ->+          prec $ max 2 queryBits+    -- bWithLowAC =+    --   case acGuide of+    --     Exact -> b+    --     NoInformation -> b+    --     _ -> normalize $+    --           MPBall x (errorBound ((0.5^(fromAccuracy acGuide))⚡))++instance HasNorm MPBall where+    getNormLog ball = getNormLog boundMP+        where+        (_, MPBall boundMP _) = endpoints $ absRaw ball++instance HasApproximate MPBall where+    type Approximate MPBall = (MPFloat, Bool)+    getApproximate ac b@(MPBall x e) =+        (approx, isAccurate)+        where+        isAccurate = getAccuracy b < ac+        approx+            | closeToN = n+            | otherwise = MPFloat.setPrecisionUp (prec (fromAccuracy ac)) x+            where+            n = mpFloat $ round $ rational x+            closeToN = ((abs $ x -^ n) <= e)++instance HasPrecision MPBall where+    getPrecision  = getPrecision . ball_value++instance CanSetPrecision MPBall where+    setPrecision p (MPBall x e)+        | p >= pPrev = MPBall xUp e+        | otherwise  = MPBall xUp (e + (xUp `EB.subMP` xDown))+        where+        pPrev = MPFloat.getPrecision x+        xUp = MPFloat.setPrecisionUp p x+        xDown = MPFloat.setPrecisionDown p x++{- negation & abs -}++instance CanNeg MPBall where+  negate (MPBall x e) = MPBall (-x) e++instance CanAbs MPBall where+  abs = normalize . absRaw++absRaw :: MPBall -> MPBall+absRaw b+  | l < 0 && 0 < r =+    fromEndpointsMP (mpFloat 0) (max (-l) r)+  | 0 <= l = b+  | otherwise = -b+  where+  (l,r) = endpointsMP b
+ src/AERN2/MP/UseMPFR/Dyadic.hs view
@@ -0,0 +1,673 @@+{-# LANGUAGE CPP #-}+-- #define DEBUG+{-# LANGUAGE DeriveDataTypeable #-}+{-|+    Module      :  AERN2.MP.UseMPFR.Dyadic+    Description :  Dyadics with exact ring operations+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Arbitrary precision floating-point numbers with exact ring operations.++    Currently, we use hmpfr when compiling with ghc 7.10 and higher+    and haskell-mpfr when compiling with ghc 7.8.+-}++module AERN2.MP.UseMPFR.Dyadic+(+   -- * Dyadic numbers and their basic operations+   Dyadic, HasDyadics+   -- * Dyadic constructors+   , CanBeDyadic, dyadic+   -- * tests+   , specDyadic, tDyadic+)+where++#ifdef DEBUG+import Debug.Trace (trace)+#define maybeTrace trace+#define maybeTraceIO putStrLn+#else+#define maybeTrace (\ (_ :: String) t -> t)+#define maybeTraceIO (\ (_ :: String) -> return ())+#endif++import MixedTypesNumPrelude+import qualified Prelude as P++import Control.CollectErrors++import Text.Printf+import Text.Regex.TDFA++import Data.Typeable+import Data.Convertible++import Test.Hspec+import Test.QuickCheck+-- import qualified Test.Hspec.SmallCheck as SC++import Data.Ratio (denominator, numerator)++import Math.NumberTheory.Logarithms (integerLog2)++import AERN2.Norm+import AERN2.MP.Precision+import AERN2.MP.Accuracy+import AERN2.MP.UseMPFR.Float++{-| Exact dyadic type based on MPFloat. -}+newtype Dyadic = Dyadic { dyadicMPFloat :: MPFloat }+  deriving (P.Eq, P.Ord, CanRound, HasPrecision, HasNorm, Typeable)++instance Ring Dyadic+instance Ring (CN Dyadic)++instance OrderedRing Dyadic+instance OrderedRing (CN Dyadic)++instance OrderedCertainlyRing Dyadic+instance OrderedCertainlyRing (CN Dyadic)++instance HasAccuracy Dyadic where getAccuracy _ = Exact++instance Show Dyadic where+  show (Dyadic x)+    | e == 0 = printf "dyadic (%d)" n+    | e > 0 = printf "dyadic (%d*0.5^%d)" n e+    | otherwise = error "in show Dyadic"+    where+    xR = rational x+    NormBits e = getNormLog (denominator xR)+    n = numerator xR++instance Read Dyadic where+  readsPrec _pr dyadicS =+    tryInt $ tryWithExp []+    where+    tryInt tryNext =+      case groups of+        [nS] ->+          case reads nS of+            [(n,"")] -> [(dyadic (n :: Integer), afterS)]+            _ -> tryNext+        _ -> tryNext+      where+      (_,_,afterS,groups) =+        dyadicS =~ "\\`dyadic \\(([-0-9]*)\\)"+          :: (String, String, String, [String])+    tryWithExp tryNext =+      case groups of+        [nS,eS] ->+          case (reads nS, reads eS) of+            ([(n,"")],[(e,"")]) ->+              [((n :: Integer)*(dyadic 0.5)^!(e :: Integer), afterS)]+            _ -> tryNext+        _ -> tryNext+      where+      (_,_,afterS,groups) =+        dyadicS =~ "\\`dyadic \\(([-0-9]*)\\*0.5\\^([0-9]*)\\)"+          :: (String, String, String, [String])++instance (SuitableForCE es) => CanEnsureCE es Dyadic+++{-- conversions --}++type HasDyadics t = ConvertibleExactly Dyadic t++instance ConvertibleExactly Dyadic Dyadic where+  safeConvertExactly = Right++instance ConvertibleExactly Dyadic MPFloat where+  safeConvertExactly = Right . dyadicMPFloat++instance ConvertibleExactly Dyadic Rational where+  safeConvertExactly = safeConvertExactly . dyadicMPFloat++type CanBeDyadic t = ConvertibleExactly t Dyadic+dyadic :: (CanBeDyadic t) => t -> Dyadic+dyadic = convertExactly++instance ConvertibleExactly MPFloat Dyadic where+  safeConvertExactly = Right . Dyadic++instance HasIntegerBounds Dyadic where+  integerBounds d = (floor d, ceiling d)++instance ConvertibleExactly Integer Dyadic where+  safeConvertExactly = fmap Dyadic . safeConvertExactly++instance ConvertibleExactly Int Dyadic where+  safeConvertExactly = fmap Dyadic . safeConvertExactly++instance ConvertibleExactly Rational Dyadic where+  safeConvertExactly q+    | isDyadic = Right $ Dyadic (fromRationalUp (prec $ max 2 (dp + np + 1)) q)+    | otherwise = convError "this number is not dyadic" q+    where+    isDyadic = d == 2^!dp+    dp = integerLog2 d+    d = denominator q+    np = integerLog2 (max 1 $ abs $ numerator q)++instance Convertible Dyadic Double where+  safeConvert = safeConvert . dyadicMPFloat++instance (ConvertibleExactly Dyadic t, Monoid es) => ConvertibleExactly Dyadic (CollectErrors es t) where+  safeConvertExactly = fmap (\v -> CollectErrors (Just v) mempty) . safeConvertExactly++{-- comparisons --}++instance HasEqAsymmetric Dyadic Dyadic+instance HasEqAsymmetric Dyadic Integer where+  equalTo = convertSecond equalTo+instance HasEqAsymmetric Integer Dyadic where+  equalTo = convertFirst equalTo+instance HasEqAsymmetric Dyadic Int where+  equalTo = convertSecond equalTo+instance HasEqAsymmetric Int Dyadic where+  equalTo = convertFirst equalTo+instance HasEqAsymmetric Dyadic Rational where+  equalTo = convertFirst equalTo+instance HasEqAsymmetric Rational Dyadic where+  equalTo = convertSecond equalTo++instance+  (HasEqAsymmetric Dyadic b+  , CanEnsureCE es b+  , CanEnsureCE es (EqCompareType Dyadic b)+  , IsBool (EnsureCE es (EqCompareType Dyadic b))+  , SuitableForCE es)+  =>+  HasEqAsymmetric Dyadic (CollectErrors es  b)+  where+  type EqCompareType Dyadic (CollectErrors es  b) =+    EnsureCE es (EqCompareType Dyadic b)+  equalTo = lift2TLCE equalTo++instance+  (HasEqAsymmetric a Dyadic+  , CanEnsureCE es a+  , CanEnsureCE es (EqCompareType a Dyadic)+  , IsBool (EnsureCE es (EqCompareType a Dyadic))+  , SuitableForCE es)+  =>+  HasEqAsymmetric (CollectErrors es a) Dyadic+  where+  type EqCompareType (CollectErrors es  a) Dyadic =+    EnsureCE es (EqCompareType a Dyadic)+  equalTo = lift2TCE equalTo++instance CanTestZero Dyadic++instance HasOrderAsymmetric Dyadic Dyadic+instance HasOrderAsymmetric Dyadic Integer where+  lessThan = convertSecond lessThan+  leq = convertSecond leq+instance HasOrderAsymmetric Integer Dyadic where+  lessThan = convertFirst lessThan+  leq = convertFirst leq+instance HasOrderAsymmetric Dyadic Int where+  lessThan = convertSecond lessThan+  leq = convertSecond leq+instance HasOrderAsymmetric Int Dyadic where+  lessThan = convertFirst lessThan+  leq = convertFirst leq+instance HasOrderAsymmetric Rational Dyadic where+  lessThan = convertSecond lessThan+  leq = convertSecond leq+instance HasOrderAsymmetric Dyadic Rational where+  lessThan = convertFirst lessThan+  leq = convertFirst leq++instance+  (HasOrderAsymmetric Dyadic b+  , CanEnsureCE es b+  , CanEnsureCE es (OrderCompareType Dyadic b)+  , IsBool (EnsureCE es (OrderCompareType Dyadic b))+  , SuitableForCE es)+  =>+  HasOrderAsymmetric Dyadic (CollectErrors es  b)+  where+  type OrderCompareType Dyadic (CollectErrors es  b) =+    EnsureCE es (OrderCompareType Dyadic b)+  lessThan = lift2TLCE lessThan+  leq = lift2TLCE leq+  greaterThan = lift2TLCE greaterThan+  geq = lift2TLCE geq++instance+  (HasOrderAsymmetric a Dyadic+  , CanEnsureCE es a+  , CanEnsureCE es (OrderCompareType a Dyadic)+  , IsBool (EnsureCE es (OrderCompareType a Dyadic))+  , SuitableForCE es)+  =>+  HasOrderAsymmetric (CollectErrors es a) Dyadic+  where+  type OrderCompareType (CollectErrors es  a) Dyadic =+    EnsureCE es (OrderCompareType a Dyadic)+  lessThan = lift2TCE lessThan+  leq = lift2TCE leq+  greaterThan = lift2TCE greaterThan+  geq = lift2TCE geq+++instance CanTestPosNeg Dyadic++instance CanTestInteger Dyadic where+  certainlyNotInteger = certainlyNotInteger . rational+  certainlyIntegerGetIt = certainlyIntegerGetIt . rational++{- unary functions -}++instance CanNeg Dyadic where+  negate = lift1 negate++instance CanAbs Dyadic where+  abs = lift1 abs++lift1 :: (MPFloat -> MPFloat) -> (Dyadic -> Dyadic)+lift1 op (Dyadic x) = Dyadic (op x)++{- min/max -}++instance CanMinMaxAsymmetric Dyadic Dyadic+instance CanMinMaxAsymmetric Integer Dyadic where+  type MinMaxType Integer Dyadic = Dyadic+  min = convertFirst min+  max = convertFirst max+instance CanMinMaxAsymmetric Dyadic Integer where+  type MinMaxType Dyadic Integer = Dyadic+  min = convertSecond min+  max = convertSecond max+instance CanMinMaxAsymmetric Int Dyadic where+  type MinMaxType Int Dyadic = Dyadic+  min = convertFirst min+  max = convertFirst max+instance CanMinMaxAsymmetric Dyadic Int where+  type MinMaxType Dyadic Int = Dyadic+  min = convertSecond min+  max = convertSecond max+instance CanMinMaxAsymmetric Rational Dyadic where+  type MinMaxType Rational Dyadic = Rational+  min = convertSecond min+  max = convertSecond max+instance CanMinMaxAsymmetric Dyadic Rational where+  type MinMaxType Dyadic Rational = Rational+  min = convertFirst min+  max = convertFirst max++instance+  (CanMinMaxAsymmetric Dyadic b+  , CanEnsureCE es b+  , CanEnsureCE es (MinMaxType Dyadic b)+  , SuitableForCE es)+  =>+  CanMinMaxAsymmetric Dyadic (CollectErrors es  b)+  where+  type MinMaxType Dyadic (CollectErrors es  b) =+    EnsureCE es (MinMaxType Dyadic b)+  min = lift2TLCE min+  max = lift2TLCE max++instance+  (CanMinMaxAsymmetric a Dyadic+  , CanEnsureCE es a+  , CanEnsureCE es (MinMaxType a Dyadic)+  , SuitableForCE es)+  =>+  CanMinMaxAsymmetric (CollectErrors es a) Dyadic+  where+  type MinMaxType (CollectErrors es  a) Dyadic =+    EnsureCE es (MinMaxType a Dyadic)+  min = lift2TCE min+  max = lift2TCE max++{- addition -}++instance CanAddAsymmetric Dyadic Dyadic where+  add = lift2 addDown addUp++instance CanAddAsymmetric Integer Dyadic where+  type AddType Integer Dyadic = Dyadic+  add = convertFirst add+instance CanAddAsymmetric Dyadic Integer where+  type AddType Dyadic Integer = Dyadic+  add = convertSecond add++instance CanAddAsymmetric Int Dyadic where+  type AddType Int Dyadic = Dyadic+  add = convertFirst add+instance CanAddAsymmetric Dyadic Int where+  type AddType Dyadic Int = Dyadic+  add = convertSecond add++instance CanAddAsymmetric Rational Dyadic where+  type AddType Rational Dyadic = Rational+  add = convertSecond add+instance CanAddAsymmetric Dyadic Rational where+  type AddType Dyadic Rational = Rational+  add = convertFirst add++instance+  (CanAddAsymmetric Dyadic b+  , CanEnsureCE es b+  , CanEnsureCE es (AddType Dyadic b)+  , SuitableForCE es)+  =>+  CanAddAsymmetric Dyadic (CollectErrors es  b)+  where+  type AddType Dyadic (CollectErrors es  b) =+    EnsureCE es (AddType Dyadic b)+  add = lift2TLCE add++instance+  (CanAddAsymmetric a Dyadic+  , CanEnsureCE es a+  , CanEnsureCE es (AddType a Dyadic)+  , SuitableForCE es)+  =>+  CanAddAsymmetric (CollectErrors es a) Dyadic+  where+  type AddType (CollectErrors es  a) Dyadic =+    EnsureCE es (AddType a Dyadic)+  add = lift2TCE add++{- subtraction -}++instance CanSub Dyadic Dyadic where+  sub = lift2 subDown subUp++instance CanSub Integer Dyadic where+  type SubType Integer Dyadic = Dyadic+  sub = convertFirst sub+instance CanSub Dyadic Integer where+  type SubType Dyadic Integer = Dyadic+  sub = convertSecond sub++instance CanSub Int Dyadic where+  type SubType Int Dyadic = Dyadic+  sub = convertFirst sub+instance CanSub Dyadic Int where+  type SubType Dyadic Int = Dyadic+  sub = convertSecond sub++instance CanSub Rational Dyadic where+  type SubType Rational Dyadic = Rational+  sub = convertSecond sub+instance CanSub Dyadic Rational where+  type SubType Dyadic Rational = Rational+  sub = convertFirst sub++instance+  (CanSub Dyadic b+  , CanEnsureCE es b+  , CanEnsureCE es (SubType Dyadic b)+  , SuitableForCE es)+  =>+  CanSub Dyadic (CollectErrors es  b)+  where+  type SubType Dyadic (CollectErrors es  b) =+    EnsureCE es (SubType Dyadic b)+  sub = lift2TLCE sub++instance+  (CanSub a Dyadic+  , CanEnsureCE es a+  , CanEnsureCE es (SubType a Dyadic)+  , SuitableForCE es)+  =>+  CanSub (CollectErrors es a) Dyadic+  where+  type SubType (CollectErrors es  a) Dyadic =+    EnsureCE es (SubType a Dyadic)+  sub = lift2TCE sub+++{- multiplication -}++instance CanMulAsymmetric Dyadic Dyadic where+  mul = lift2 mulDown mulUp++instance CanMulAsymmetric Integer Dyadic where+  type MulType Integer Dyadic = Dyadic+  mul = convertFirst mul+instance CanMulAsymmetric Dyadic Integer where+  type MulType Dyadic Integer = Dyadic+  mul = convertSecond mul++instance CanMulAsymmetric Int Dyadic where+  type MulType Int Dyadic = Dyadic+  mul = convertFirst mul+instance CanMulAsymmetric Dyadic Int where+  type MulType Dyadic Int = Dyadic+  mul = convertSecond mul++instance CanMulAsymmetric Rational Dyadic where+  type MulType Rational Dyadic = Rational+  mul = convertSecond mul+instance CanMulAsymmetric Dyadic Rational where+  type MulType Dyadic Rational = Rational+  mul = convertFirst mul++instance+  (CanMulAsymmetric Dyadic b+  , CanEnsureCE es b+  , CanEnsureCE es (MulType Dyadic b)+  , SuitableForCE es)+  =>+  CanMulAsymmetric Dyadic (CollectErrors es  b)+  where+  type MulType Dyadic (CollectErrors es  b) =+    EnsureCE es (MulType Dyadic b)+  mul = lift2TLCE mul++instance+  (CanMulAsymmetric a Dyadic+  , CanEnsureCE es a+  , CanEnsureCE es (MulType a Dyadic)+  , SuitableForCE es)+  =>+  CanMulAsymmetric (CollectErrors es a) Dyadic+  where+  type MulType (CollectErrors es  a) Dyadic =+    EnsureCE es (MulType a Dyadic)+  mul = lift2TCE mul++instance CanPow Dyadic Integer where+  powNoCN = powUsingMul (dyadic 1)+  pow = integerPowCN (powUsingMul (dyadic 1))+instance CanPow Dyadic Int where+  powNoCN = powUsingMul (dyadic 1)+  pow = integerPowCN (powUsingMul (dyadic 1))++instance+  (CanDiv a Dyadic+  , CanEnsureCE es a+  , CanEnsureCE es (DivType a Dyadic)+  , CanEnsureCE es (DivTypeNoCN a Dyadic)+  , SuitableForCE es)+  =>+  CanDiv (CollectErrors es a) Dyadic+  where+  type DivType (CollectErrors es  a) Dyadic =+    EnsureCE es (DivType a Dyadic)+  divide = lift2TCE divide+  type DivTypeNoCN (CollectErrors es a) Dyadic =+    EnsureCE es (DivTypeNoCN a Dyadic)+  divideNoCN = lift2TCE divideNoCN++instance CanDiv Integer Dyadic where+  type DivTypeNoCN Integer Dyadic = Rational+  divideNoCN a b = divideNoCN a (rational b)+instance CanDiv Dyadic Integer where+  type DivTypeNoCN Dyadic Integer = Rational+  divideNoCN a b = divideNoCN (rational a) b++instance CanDiv Int Dyadic where+  type DivTypeNoCN Int Dyadic = Rational+  divideNoCN a b = divideNoCN a (rational b)+instance CanDiv Dyadic Int where+  type DivTypeNoCN Dyadic Int = Rational+  divideNoCN a b = divideNoCN (rational a) b++instance CanDiv Rational Dyadic where+  type DivTypeNoCN Rational Dyadic = Rational+  divideNoCN = convertSecond divideNoCN+instance CanDiv Dyadic Rational where+  type DivTypeNoCN Dyadic Rational = Rational+  divideNoCN = convertFirst divideNoCN++instance+  (CanDiv Dyadic b+  , CanEnsureCE es b+  , CanEnsureCE es (DivType Dyadic b)+  , CanEnsureCE es (DivTypeNoCN Dyadic b)+  , SuitableForCE es)+  =>+  CanDiv Dyadic (CollectErrors es  b)+  where+  type DivType Dyadic (CollectErrors es  b) =+    EnsureCE es (DivType Dyadic b)+  divide = lift2TLCE divide+  type DivTypeNoCN Dyadic (CollectErrors es  b) =+    EnsureCE es (DivTypeNoCN Dyadic b)+  divideNoCN = lift2TLCE divideNoCN++instance+  (CanPow Dyadic b+  , CanEnsureCE es b+  , CanEnsureCE es (PowTypeNoCN Dyadic b)+  , CanEnsureCE es (PowType Dyadic b)+  , SuitableForCE es)+  =>+  CanPow Dyadic (CollectErrors es  b)+  where+  type PowTypeNoCN Dyadic (CollectErrors es b) =+    EnsureCE es (PowTypeNoCN Dyadic b)+  powNoCN = lift2TLCE powNoCN+  type PowType Dyadic (CollectErrors es b) =+    EnsureCE es (PowType Dyadic b)+  pow = lift2TLCE pow++instance+  (CanPow a Dyadic+  , CanEnsureCE es a+  , CanEnsureCE es (PowType a Dyadic)+  , CanEnsureCE es (PowTypeNoCN a Dyadic)+  , SuitableForCE es)+  =>+  CanPow (CollectErrors es a) Dyadic+  where+  type PowTypeNoCN (CollectErrors es  a) Dyadic =+    EnsureCE es (PowTypeNoCN a Dyadic)+  powNoCN = lift2TCE powNoCN+  type PowType (CollectErrors es  a) Dyadic =+    EnsureCE es (PowType a Dyadic)+  pow = lift2TCE pow++lift2 ::+  (MPFloat -> MPFloat -> MPFloat) ->+  (MPFloat -> MPFloat -> MPFloat) ->+  (Dyadic -> Dyadic -> Dyadic)+lift2 opDown opUp (Dyadic x0) (Dyadic y0) = Dyadic (opExact x0 y0)+  where+    opExact x y+      | rUp == rDown = rUp+      | otherwise =+          maybeTrace (printf "Dyadic.lift2: rUp = %s; rDown = %s; p = %s" (show rUp) (show rDown) (show $ integer p)) $+          opExact xH yH+      where+      rUp = opUp x y+      rDown = opDown x y+      xH = setPrecision pH x+      yH = setPrecision pH y+      pH = precisionTimes2 p+      p = getPrecision rUp++instance Arbitrary Dyadic where+  arbitrary =+    do+      c <- finiteMPFloat+      return (Dyadic c)+    where+      finiteMPFloat =+        do+          x <- arbitrary+          if isFinite x+            then return x+            else finiteMPFloat++{-|+  A runtime representative of type @Dyadic@.+  Used for specialising polymorphic tests to concrete types.+-}+tDyadic :: T Dyadic+tDyadic = T "Dyadic"++specDyadic :: Spec+specDyadic =+  describe ("Dyadic") $ do+    specConversion tInteger tDyadic dyadic round+    specConversion tDyadic tRational rational dyadic+    describe "order" $ do+      specHasEqNotMixed tDyadic+      specHasEq tInt tDyadic tRational+      specCanTestZero tDyadic+      specHasOrderNotMixed tDyadic+      specHasOrder tInt tDyadic tRational+    describe "min/max/abs" $ do+      specCanNegNum tDyadic+      specCanAbs tDyadic+      specCanMinMaxNotMixed tDyadic+      specCanMinMax tDyadic tInteger tDyadic+      it "min Dyadic Rational (dyadic only)" $ do+        property $ \ (x :: Dyadic) (y :: Dyadic) ->+          x `min` y == x `min` (rational y)+      it "max Dyadic Rational (dyadic only)" $ do+        property $ \ (x :: Dyadic) (y :: Dyadic) ->+          x `max` y == x `max` (rational y)+    describe "ring" $ do+      specCanAddNotMixed tDyadic+      specCanAddSameType tDyadic+      specCanAdd tInt tDyadic tInteger+      specCanAdd tInteger tDyadic tInt+      it "Dyadic + Rational (dyadic only)" $ do+        property $ \ (x :: Dyadic) (y :: Dyadic) ->+          x + y == x + (rational y)+      specCanSubNotMixed tDyadic+      specCanSub tDyadic tInteger+      specCanSub tInteger tDyadic+      specCanSub tDyadic tInt+      specCanSub tInt tDyadic+      it "Dyadic - Rational (dyadic only)" $ do+        property $ \ (x :: Dyadic) (y :: Dyadic) ->+          x - y == x - (rational y)+      specCanMulNotMixed tDyadic+      specCanMulSameType tDyadic+      specCanMul tInt tDyadic tInteger+      it "Dyadic * Rational (dyadic only)" $ do+        property $ \ (x :: Dyadic) (y :: Dyadic) ->+          x * y == x * (rational y)+      specCanPow tDyadic tInteger++instance P.Num Dyadic where+    fromInteger = convertExactly+    negate = negate+    (+) = (+)+    (*) = (*)+    abs = abs+    signum d+      | d < 0 = dyadic (-1)+      | d == 0 = dyadic 0+      | otherwise = dyadic 1++instance P.Real Dyadic where+    toRational = convertExactly
+ src/AERN2/MP/UseMPFR/ErrorBound.hs view
@@ -0,0 +1,221 @@+{-# LANGUAGE DeriveDataTypeable #-}+{-|+    Module      :  AERN2.MP.UseMPFR.ErrorBound+    Description :  Fixed precision non-negative up-rounded floating-point numbers+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Fixed precision non-negative up-rounded floating-point numbers.++    Currently using a fixed-precision MPFloat.+-}+module AERN2.MP.UseMPFR.ErrorBound+    (ErrorBound, CanBeErrorBound, errorBound,+     absMP, subMP)+where++import MixedTypesNumPrelude+import qualified Prelude as P++import Data.Typeable++import Test.QuickCheck++import Data.Convertible++import Math.NumberTheory.Logarithms (integerLog2)++import AERN2.MP.Precision+import AERN2.MP.Accuracy+import qualified AERN2.MP.UseMPFR.Float as MPFloat+import AERN2.MP.UseMPFR.Float (MPFloat, mpFloat, frequencyElements)+import AERN2.MP.UseMPFR.Float.Operators+import AERN2.MP.Dyadic++{- example -}++_example1 :: ErrorBound+_example1 = 2*((errorBound 0.01) + 0.1*(errorBound 0.01)/3)++{- type -}++{-| A non-negative Double value to serve as an error bound. Arithmetic is rounded towards +infinity. -}+newtype ErrorBound = ErrorBound { er2mp :: MPFloat }+  deriving (P.Eq, P.Ord, Typeable)++instance Show ErrorBound where+    show (ErrorBound d) = show d++errorBoundPrecision :: Precision+errorBoundPrecision = prec 53++instance HasAccuracy ErrorBound where+  getAccuracy (ErrorBound e)+      | eN > 0 =+          bits $ negate $ integerLog2 eN+      | e > 0 && eRecipN > 0 =+          bits $ integerLog2 eRecipN+      | e == 0 = Exact+      | otherwise = NoInformation+      where+      eN = floor $ rational e+      eRecipN = ceiling $ rational $ MPFloat.recipDown e++{- conversions -}++instance ConvertibleExactly ErrorBound MPFloat where+  safeConvertExactly = Right . er2mp++instance ConvertibleExactly ErrorBound Dyadic where+  safeConvertExactly = Right . dyadic . er2mp++instance ConvertibleExactly ErrorBound Rational where+  safeConvertExactly = Right . convertExactly . mpFloat++type CanBeErrorBound t = Convertible t ErrorBound+errorBound :: (CanBeErrorBound t) => t -> ErrorBound+errorBound = convert++instance Convertible Rational ErrorBound where+  safeConvert x+    | x >= 0 = Right $ ErrorBound $ MPFloat.fromRationalUp errorBoundPrecision x+    | otherwise = convError "Trying to construct a negative ErrorBound" x++instance Convertible MPFloat ErrorBound where+  safeConvert x+    | x >= 0 = Right $ ErrorBound $ MPFloat.setPrecisionUp errorBoundPrecision x+    | otherwise = convError "Trying to construct a negative ErrorBound" x++instance Convertible Integer ErrorBound where+  safeConvert x+    | x >= 0 = Right $ ErrorBound $ MPFloat.fromIntegerUp errorBoundPrecision x+    | otherwise = convError "Trying to construct a negative ErrorBound" x++instance Convertible Int ErrorBound where+  safeConvert = safeConvert . integer++{- comparisons -}++instance HasOrderAsymmetric ErrorBound ErrorBound++instance HasOrderAsymmetric ErrorBound MPFloat where+  lessThan = convertFirst lessThan+  leq = convertFirst leq+instance HasOrderAsymmetric MPFloat ErrorBound where+  lessThan = convertSecond lessThan+  leq = convertSecond leq++instance HasEqAsymmetric ErrorBound Rational where+  equalTo = convertFirst equalTo+instance HasEqAsymmetric Rational ErrorBound where+  equalTo = convertSecond equalTo+instance HasOrderAsymmetric ErrorBound Rational where+  lessThan = convertFirst lessThan+  leq = convertFirst leq+instance HasOrderAsymmetric Rational ErrorBound where+  lessThan = convertSecond lessThan+  leq = convertSecond leq++instance HasEqAsymmetric ErrorBound Integer where+  equalTo a b = equalTo (dyadic a) (dyadic b)+instance HasEqAsymmetric Integer ErrorBound where+  equalTo a b = equalTo (dyadic a) (dyadic b)+instance HasOrderAsymmetric ErrorBound Integer where+  lessThan a b = lessThan (dyadic a) (dyadic b)+  leq a b = leq (dyadic a) (dyadic b)+instance HasOrderAsymmetric Integer ErrorBound where+  lessThan a b = lessThan (dyadic a) (dyadic b)+  leq a b = leq (dyadic a) (dyadic b)++instance HasEqAsymmetric ErrorBound Int where+  equalTo a b = equalTo (dyadic a) (dyadic b)+instance HasEqAsymmetric Int ErrorBound where+  equalTo a b = equalTo (dyadic a) (dyadic b)+instance HasOrderAsymmetric ErrorBound Int where+  lessThan a b = lessThan (dyadic a) (dyadic b)+  leq a b = leq (dyadic a) (dyadic b)+instance HasOrderAsymmetric Int ErrorBound where+  lessThan a b = lessThan (dyadic a) (dyadic b)+  leq a b = leq (dyadic a) (dyadic b)++instance CanMinMaxAsymmetric ErrorBound ErrorBound++{- converting operations -}++subMP :: MPFloat -> MPFloat -> ErrorBound+a `subMP` b = errorBound $ a -^ b++absMP :: MPFloat -> ErrorBound+absMP = errorBound . abs++{- up-rounded operations -}++instance CanAddAsymmetric ErrorBound ErrorBound where+    add (ErrorBound a) (ErrorBound b) = ErrorBound $ a +^ b++instance CanAddAsymmetric ErrorBound MPFloat where+  type AddType ErrorBound MPFloat = ErrorBound+  add = convertSecondUsing (\ _ f -> convert f) add+instance CanAddAsymmetric MPFloat ErrorBound where+  type AddType MPFloat ErrorBound = ErrorBound+  add = convertFirstUsing (\ f _ -> convert f) add++instance CanMulAsymmetric ErrorBound ErrorBound where+    mul (ErrorBound a) (ErrorBound b) = ErrorBound $ a *^ b++instance CanMulAsymmetric ErrorBound MPFloat where+  type MulType ErrorBound MPFloat = ErrorBound+  mul = convertSecondUsing (\ _ f -> convert f) mul+instance CanMulAsymmetric MPFloat ErrorBound where+  type MulType MPFloat ErrorBound = ErrorBound+  mul = convertFirstUsing (\ f _ -> convert f) mul++instance CanMulAsymmetric ErrorBound Integer where+    type MulType ErrorBound Integer = ErrorBound+    mul (ErrorBound a) i+        | i >= 0 = ErrorBound $ a *^ (MPFloat.fromIntegerUp errorBoundPrecision i)+        | otherwise = error "trying to multiply ErrorBound by a negative integer"+instance CanMulAsymmetric Integer ErrorBound where+    type MulType Integer ErrorBound = ErrorBound+    mul i (ErrorBound b)+        | i >= 0 = ErrorBound $ (MPFloat.fromIntegerUp errorBoundPrecision i) *^ b+        | otherwise = error "trying to multiply ErrorBound by a negative integer"++instance CanMulAsymmetric ErrorBound Rational where+    type MulType ErrorBound Rational = ErrorBound+    mul (ErrorBound a) r+        | r >= 0.0 = ErrorBound $ a *^ (MPFloat.fromRationalUp errorBoundPrecision r)+        | otherwise = error "trying to multiply ErrorBound by a negative integer"+instance CanMulAsymmetric Rational ErrorBound where+    type MulType Rational ErrorBound = ErrorBound+    mul r (ErrorBound b)+        | r >= 0.0 = ErrorBound $ (MPFloat.fromRationalUp errorBoundPrecision r) *^ b+        | otherwise = error "trying to multiply ErrorBound by a negative integer"++instance CanDiv ErrorBound Integer where+    type DivTypeNoCN ErrorBound Integer = ErrorBound+    type DivType ErrorBound Integer = ErrorBound+    divideNoCN = divide+    divide (ErrorBound a) i+        | i > 0 = ErrorBound $ a /^ (MPFloat.fromIntegerUp errorBoundPrecision i)+        | otherwise = error "trying to multiply ErrorBound by a non-positive integer"++instance Arbitrary ErrorBound where+  arbitrary =+    do+    giveSpecialValue <- frequencyElements [(5, False),(1, True)]+    aux giveSpecialValue+    where+      aux giveSpecialValue+        | giveSpecialValue =+            elements (map convert [0.0,0.0,0.0,10.0,1.0,0.5,0.125])+        | otherwise =+          do+          (s :: Integer) <- arbitrary+          let resultR = ((abs s) `mod` (2^!35))/!(2^!32)+          let result = convert resultR+          return result
+ src/AERN2/MP/UseMPFR/Float.hs view
@@ -0,0 +1,60 @@+{-|+    Module      :  AERN2.MP.UseMPFR.Float+    Description :  Arbitrary precision floating point numbers+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Arbitrary precision floating-point numbers with up/down-rounded operations.++    Currently, we use hmpfr when compiling with ghc 7.10 and higher+    and haskell-mpfr when compiling with ghc 7.8.+-}++module AERN2.MP.UseMPFR.Float+  (+   -- * Precision operations+   module AERN2.MP.Precision+   -- * The type definition and basic operations+   , module AERN2.MP.UseMPFR.Float.Type+   -- * Arithmetic operations+   , module AERN2.MP.UseMPFR.Float.Arithmetic+   , distUp, distDown, avgUp, avgDown+   -- * Conversions, comparisons and norm+   , module AERN2.MP.UseMPFR.Float.Conversions+   -- * Infix operators for up/down-rounded operations+   , module AERN2.MP.UseMPFR.Float.Operators+   -- * Constants such as NaN, infinity+   , module AERN2.MP.UseMPFR.Float.Constants+   -- * Tests+   , module AERN2.MP.UseMPFR.Float.Tests+   )+where++import MixedTypesNumPrelude+-- import qualified Prelude as P++import AERN2.MP.Precision+import AERN2.MP.UseMPFR.Float.Type+import AERN2.MP.UseMPFR.Float.Arithmetic+import AERN2.MP.UseMPFR.Float.Conversions+import AERN2.MP.UseMPFR.Float.Operators+import AERN2.MP.UseMPFR.Float.Constants+import AERN2.MP.UseMPFR.Float.Tests++-- | Computes an upper bound to the distance @|x - y|@ of @x@ and @y@.+distUp :: MPFloat -> MPFloat -> MPFloat+distUp x y = if x >= y then x -^ y else y -^ x++-- | Computes a lower bound to the distance @|x - y|@ of @x@ and @y@.+distDown :: MPFloat -> MPFloat -> MPFloat+distDown x y = if x >= y then x -. y else y -. x++avgUp :: MPFloat -> MPFloat -> MPFloat+avgUp x y = (x +^ y) /^ (mpFloat 2)++avgDown :: MPFloat -> MPFloat -> MPFloat+avgDown x y = (x +. y) /. (mpFloat 2)
+ src/AERN2/MP/UseMPFR/Float/Arithmetic.hs view
@@ -0,0 +1,170 @@+{-# LANGUAGE CPP #-}+{-|+    Module      :  AERN2.MP.UseMPFR.Float.Arithmetic+    Description :  Arbitrary precision floating point numbers+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Arbitrary precision floating-point numbers with up/down-rounded operations.++    Currently, we use hmpfr when compiling with ghc 7.10 and higher+    and haskell-mpfr when compiling with ghc 7.8.+-}++module AERN2.MP.UseMPFR.Float.Arithmetic+  (+   -- * MPFloat basic arithmetic+     addUp, addDown, subUp, subDown+   , mulUp, mulDown, divUp, divDown, recipUp, recipDown+   -- * MPFloat selected constants and operations+   , piUp, piDown+   , cosUp, cosDown, sinUp, sinDown+   , sqrtUp, sqrtDown, expUp, expDown, logUp, logDown+   )+where++import MixedTypesNumPrelude+import qualified Prelude as P++import AERN2.MP.Precision+import AERN2.MP.UseMPFR.Float.Type++#ifdef HaskellMPFR+import qualified Data.Approximate.MPFRLowLevel as MPLow++one :: MPFloat+one = MPLow.fromInt MPLow.Up (P.fromInteger 10) (int 1)++#endif++#ifdef MPFRRounded+import qualified AERN2.MP.UseMPFR.Float.RoundedAdaptor as MPLow++one :: MPFloat+one = MPLow.one++#endif++#ifdef HMPFR+import qualified Data.Number.MPFR as MPLow++one :: MPFloat+one = MPLow.one+#endif++{- common functions -}++instance CanNeg MPFloat where+  negate = unaryUp MPLow.neg++instance CanAbs MPFloat where+  abs x+    | x P.< MPLow.zero = negate x+    | otherwise = x++addUp, addDown :: MPFloat -> MPFloat -> MPFloat+addUp = binaryUp True MPLow.add+addDown = binaryDown True MPLow.add++subUp, subDown :: MPFloat -> MPFloat -> MPFloat+subUp = binaryUp True MPLow.sub+subDown = binaryDown True MPLow.sub++mulUp, mulDown :: MPFloat -> MPFloat -> MPFloat+mulUp = binaryUp True MPLow.mul+mulDown = binaryDown True MPLow.mul++divUp,divDown :: MPFloat -> MPFloat -> MPFloat+divUp = binaryUp False MPLow.div+divDown = binaryDown False MPLow.div++recipUp :: MPFloat -> MPFloat+recipUp x = divUp one x++recipDown :: MPFloat -> MPFloat+recipDown x = divDown one x+++{- special constants and functions -}++piUp :: Precision -> MPFloat+piUp p =+    MPLow.pi MPLow.Up (p2mpfrPrec p)++piDown :: Precision -> MPFloat+piDown p =+    MPLow.pi MPLow.Down (p2mpfrPrec p)++cosUp :: MPFloat -> MPFloat+cosUp = unaryUp MPLow.cos++cosDown :: MPFloat -> MPFloat+cosDown = unaryDown MPLow.cos++sinUp :: MPFloat -> MPFloat+sinUp = unaryUp MPLow.sin++sinDown :: MPFloat -> MPFloat+sinDown = unaryDown MPLow.sin++sqrtUp :: MPFloat -> MPFloat+sqrtUp = unaryUp MPLow.sqrt++sqrtDown :: MPFloat -> MPFloat+sqrtDown = unaryDown MPLow.sqrt++expUp :: MPFloat -> MPFloat+expUp = unaryUp MPLow.exp++expDown :: MPFloat -> MPFloat+expDown = unaryDown MPLow.exp++logUp :: MPFloat -> MPFloat+logUp = unaryUp MPLow.log++logDown :: MPFloat -> MPFloat+logDown = unaryDown MPLow.log++{- auxiliary functions to automatically determine result precision from operand precisions -}++unaryUp ::+    (MPLow.RoundMode -> MPLow.Precision -> MPFloat -> MPFloat) ->+    (MPFloat -> MPFloat)+unaryUp opRP x = opRP MPLow.Up p x+    where+    p = MPLow.getPrec x++unaryDown ::+    (MPLow.RoundMode -> MPLow.Precision -> MPFloat -> MPFloat) ->+    (MPFloat -> MPFloat)+unaryDown opRP x = opRP MPLow.Down p x+    where+    p = MPLow.getPrec x++binaryUp ::+    Bool ->+    (MPLow.RoundMode -> MPLow.Precision -> MPFloat -> MPFloat -> MPFloat) ->+    (MPFloat -> MPFloat -> MPFloat)+binaryUp = binaryApprox True++binaryDown ::+    Bool ->+    (MPLow.RoundMode -> MPLow.Precision -> MPFloat -> MPFloat -> MPFloat) ->+    (MPFloat -> MPFloat -> MPFloat)+binaryDown = binaryApprox False++binaryApprox ::+    Bool -> Bool ->+    (MPLow.RoundMode -> MPLow.Precision -> MPFloat -> MPFloat -> MPFloat) ->+    (MPFloat -> MPFloat -> MPFloat)+binaryApprox isUp _canBeExact opRP x y =+    withPrec pMax+    where+    pMax = (getPrecision x) `max` (getPrecision y)+    withPrec p+        | isUp = opRP MPLow.Up (p2mpfrPrec p) x y+        | otherwise = opRP MPLow.Down (p2mpfrPrec p) x y
+ src/AERN2/MP/UseMPFR/Float/Constants.hs view
@@ -0,0 +1,52 @@+{-# LANGUAGE CPP #-}+{-|+    Module      :  AERN2.MP.UseMPFR.Float.Constants+    Description :  Special constants NaN, infinity etc+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Special constants NaN, infinity etc+-}++module AERN2.MP.UseMPFR.Float.Constants+  (+    zero, one+    , nan, infinity+  )+where++import MixedTypesNumPrelude+import qualified Prelude as P+-- import Data.Ratio++import AERN2.MP.UseMPFR.Float.Type+import AERN2.MP.UseMPFR.Float.Conversions+import AERN2.MP.UseMPFR.Float.Operators++zero, one :: MPFloat+zero = mpFloat 0+one = mpFloat 1++nan, infinity :: MPFloat+nan = zero /. zero+infinity = one /. zero++itisNaN :: MPFloat -> Bool+itisNaN x = x *^ one /= x++itisInfinite :: MPFloat -> Bool+itisInfinite x =+  x *^ (mpFloat 2) P.== x+  &&+  x P./= (mpFloat 0)++instance CanTestFinite MPFloat where+  isInfinite = itisInfinite+  isFinite x = not (itisInfinite x || itisNaN x)++instance CanTestNaN MPFloat where+  isNaN = itisNaN
+ src/AERN2/MP/UseMPFR/Float/Conversions.hs view
@@ -0,0 +1,208 @@+{-# LANGUAGE CPP #-}+{-|+    Module      :  AERN2.MP.UseMPFR.Float.Conversions+    Description :  Conversions and comparisons of arbitrary precision floats+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Conversions and comparisons of arbitrary precision floating point numbers++    Currently, we use hmpfr when compiling with ghc 7.10 and higher+    and haskell-mpfr when compiling with ghc 7.8.+-}++module AERN2.MP.UseMPFR.Float.Conversions+  (+   -- * MPFloat to other types (see also instances)+   toDoubleUp, toDoubleDown+   -- * MPFloat constructors (see also instances)+   , CanBeMPFloat, mpFloat+   , fromIntegerUp, fromIntegerDown+   , fromRationalUp, fromRationalDown+   )+where++import MixedTypesNumPrelude+import qualified Prelude as P++import Data.Ratio+import Data.Convertible++import AERN2.Norm+import AERN2.MP.Precision++import AERN2.MP.UseMPFR.Float.Type+import AERN2.MP.UseMPFR.Float.Arithmetic++#ifdef HaskellMPFR+import qualified Data.Approximate.MPFRLowLevel as MPLow++mpToDouble :: MPLow.RoundMode -> MPFloat -> Double+mpToDouble = MPLow.toDoubleA++mpToRational :: MPFloat -> Rational+mpToRational x+    | x == 0 = 0.0+    | otherwise = MPLow.toRationalA x++mpFromRationalA :: MPLow.RoundMode -> MPLow.Precision -> Rational -> MPFloat+mpFromRationalA = MPLow.fromRationalA++#endif+#ifdef MPFRRounded+import qualified AERN2.MP.UseMPFR.Float.RoundedAdaptor as MPLow++mpToDouble :: MPLow.RoundMode -> MPFloat -> Double+mpToDouble = MPLow.toDoubleA++mpToRational :: MPFloat -> Rational+mpToRational x+  | x == 0 = 0.0+  | otherwise = MPLow.toRationalA x++mpFromRationalA :: MPLow.RoundMode -> MPLow.Precision -> Rational -> MPFloat+mpFromRationalA = MPLow.fromRationalA++#endif+#ifdef HMPFR+import qualified Data.Number.MPFR as MPLow++mpToDouble :: MPLow.RoundMode -> MPLow.MPFR -> Double+mpToDouble = MPLow.toDouble++mpToRational :: MPFloat -> Rational+mpToRational x+    | x == 0 = 0.0+    | otherwise = mantissa * 2.0^!e+    where+    (mantissa, ePre) = MPLow.decompose x+    e = P.toInteger ePre++mpFromRationalA :: MPLow.RoundMode -> MPLow.Precision -> Rational -> MPFloat+mpFromRationalA dir p q+  | q < 0 =+    MPLow.fromIntegerA dir p (numerator q) `divDir` MPLow.fromIntegerA dir p (denominator q)+  | otherwise =+    MPLow.fromIntegerA dir p (numerator q) `divDir` MPLow.fromIntegerA dirOpp p (denominator q)+  where+  (divDir, dirOpp) =+    case dir of+      MPLow.Down -> (divDown, MPLow.Up)+      MPLow.Up -> (divUp, MPLow.Down)+      _ -> error "in mpFromRationalA"++#endif++instance HasNorm MPFloat where+  getNormLog x+    | x == 0 = NormZero+    | otherwise = NormBits (P.toInteger $ MPLow.getExp x)++{- conversions -}++instance CanRound MPFloat where+  properFraction x = (n,f)+    where+      r = rational x+      n = (numerator r) `quot` (denominator r)+      f = x `subUp` (mpFloat n)++instance ConvertibleExactly MPFloat Rational where+  safeConvertExactly = Right . mpToRational++toDoubleUp :: MPFloat -> Double+toDoubleUp = mpToDouble MPLow.Up++toDoubleDown :: MPFloat -> Double+toDoubleDown = mpToDouble MPLow.Down++fromIntegerUp :: Precision -> Integer -> MPFloat+fromIntegerUp p i = MPLow.fromIntegerA MPLow.Up (p2mpfrPrec p) i++fromIntegerDown :: Precision -> Integer -> MPFloat+fromIntegerDown p i = MPLow.fromIntegerA MPLow.Down (p2mpfrPrec p) i++type CanBeMPFloat t = ConvertibleExactly t MPFloat+mpFloat :: (CanBeMPFloat t) => t -> MPFloat+mpFloat = convertExactly++instance ConvertibleExactly Integer MPFloat where+    safeConvertExactly n =+        findExact $ map upDown $ standardPrecisions initPrec+        where+        initPrec =+            case getNormLog n of+              NormBits b -> prec (b + 8)+              _ -> prec 8+        upDown p = (fromIntegerDown p n, fromIntegerUp p n)+        findExact [] =+            convError "integer too high to represent exactly" n+        findExact ((nDown, nUp) : rest)+            | nDown == nUp = Right nUp+            | otherwise = findExact rest++instance ConvertibleExactly Int MPFloat where+    safeConvertExactly = safeConvertExactly . integer++fromRationalUp :: Precision -> Rational -> MPFloat+fromRationalUp p x =+    mpFromRationalA MPLow.Up (p2mpfrPrec p) x++fromRationalDown :: Precision -> Rational -> MPFloat+fromRationalDown p x =+    mpFromRationalA MPLow.Down (p2mpfrPrec p) x++instance Convertible MPFloat Double where+  safeConvert x+    | isFinite dbl = Right dbl+    | otherwise = convError "conversion to double: out of bounds" x+    where+    dbl = toDoubleUp x++{- comparisons -}++instance HasEqAsymmetric MPFloat MPFloat+instance HasEqAsymmetric MPFloat Integer where+  equalTo = convertSecond equalTo+instance HasEqAsymmetric Integer MPFloat where+  equalTo = convertFirst equalTo+instance HasEqAsymmetric MPFloat Int where+  equalTo = convertSecond equalTo+instance HasEqAsymmetric Int MPFloat where+  equalTo = convertFirst equalTo+instance HasEqAsymmetric MPFloat Rational where+  equalTo = convertFirst equalTo+instance HasEqAsymmetric Rational MPFloat where+  equalTo = convertSecond equalTo++instance CanTestZero MPFloat++instance HasOrderAsymmetric MPFloat MPFloat+instance HasOrderAsymmetric MPFloat Integer where+  lessThan = convertSecond lessThan+  leq = convertSecond leq+instance HasOrderAsymmetric Integer MPFloat where+  lessThan = convertFirst lessThan+  leq = convertFirst leq+instance HasOrderAsymmetric MPFloat Int where+  lessThan = convertSecond lessThan+  leq = convertSecond leq+instance HasOrderAsymmetric Int MPFloat where+  lessThan = convertFirst lessThan+  leq = convertFirst leq+instance HasOrderAsymmetric Rational MPFloat where+  lessThan = convertSecond lessThan+  leq = convertSecond leq+instance HasOrderAsymmetric MPFloat Rational where+  lessThan = convertFirst lessThan+  leq = convertFirst leq++instance CanTestPosNeg MPFloat++{- min, max -}++instance CanMinMaxAsymmetric MPFloat MPFloat
+ src/AERN2/MP/UseMPFR/Float/Operators.hs view
@@ -0,0 +1,38 @@+{-|+    Module      :  AERN2.MP.UseMPFR.Float.Operators+    Description :  Infix operators for up/down-rounded floating-point numbers+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Infix operators for up/down-rounded floating-point numbers+-}++module AERN2.MP.UseMPFR.Float.Operators where++import AERN2.MP.UseMPFR.Float.Type+import AERN2.MP.UseMPFR.Float.Arithmetic++infixl 6  +^, -^, +., -.+infixl 7  *^, *., /^, /.++(+^) :: MPFloat -> MPFloat -> MPFloat+(+^) = addUp+(-^) :: MPFloat -> MPFloat -> MPFloat+(-^) = subUp+(*^) :: MPFloat -> MPFloat -> MPFloat+(*^) = mulUp+(/^) :: MPFloat -> MPFloat -> MPFloat+(/^) = divUp++(+.) :: MPFloat -> MPFloat -> MPFloat+(+.) = addDown+(-.) :: MPFloat -> MPFloat -> MPFloat+(-.) = subDown+(*.) :: MPFloat -> MPFloat -> MPFloat+(*.) = mulDown+(/.) :: MPFloat -> MPFloat -> MPFloat+(/.) = divDown
+ src/AERN2/MP/UseMPFR/Float/RoundedAdaptor.hs view
@@ -0,0 +1,73 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE DataKinds, ExistentialQuantification, RankNTypes #-}+-- {-# LANGUAGE DeriveGeneric, DeriveDataTypeable, StandaloneDeriving #-}+{-|+    Module      :  AERN2.MP.UseMPFR.Float.RoundedAdaptor+    Description :  Numeric.Rounded + variable precision+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Numeric.Rounded + variable precision+-}+module AERN2.MP.UseMPFR.Float.RoundedAdaptor+#ifndef MPFRRounded+() where+#else+(+  module AERN2.MP.UseMPFR.Float.RoundedAdaptor+, module Numeric.RoundedSimple+)+where++import Prelude hiding (div, pi)+-- import qualified Prelude as P++import Numeric.RoundedSimple hiding (prec)+-- import qualified Numeric.RoundedSimple as R++type RoundMode = RoundingMode++defaultPrecision :: Precision+defaultPrecision = 10++pi :: RoundMode -> Precision -> Rounded+pi = kPi++fromIntegerA :: RoundMode -> Precision -> Integer -> Rounded+fromIntegerA = fromInteger'++zero, one :: Rounded+zero = fromIntegerA Up defaultPrecision 0+one = fromIntegerA Up defaultPrecision 1++toDoubleA :: RoundMode -> Rounded -> Double+toDoubleA = toDouble++fromRationalA :: RoundMode -> Precision -> Rational -> Rounded+fromRationalA = fromRational'++toRationalA :: Rounded -> Rational+toRationalA = toRational'++add, sub, mul, div, atan2 :: RoundMode -> Precision -> Rounded -> Rounded -> Rounded+add = add_+sub = sub_+mul = mul_+div = div_+atan2 = atan2_++neg, abs, sqrt, exp, log, sin, cos :: RoundMode -> Precision -> Rounded -> Rounded+neg = negate_+abs = abs_+sqrt = sqrt_+exp = exp_+log = log_+sin = sin_+cos = cos_+-- TODO: add more ops++#endif
+ src/AERN2/MP/UseMPFR/Float/Tests.hs view
@@ -0,0 +1,424 @@+{-# LANGUAGE CPP #-}+{-|+    Module      :  AERN2.MP.UseMPFR.Float.Tests+    Description :  Tests for operations on arbitrary precision floats+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Tests for operations on arbitrary precision floats.++    To run the tests using stack, execute:++    @+    stack test aern2-mp --test-arguments "-a 1000 -m MPFloat"+    @+-}++module AERN2.MP.UseMPFR.Float.Tests+  (+    specMPFloat, tMPFloat+    , (=~=), approxEqual, approxEqualWithArgs+    , frequencyElements+  )+where++import MixedTypesNumPrelude+-- import qualified Prelude as P+-- import Data.Ratio+import Text.Printf+import Data.Maybe++import Test.Hspec+import Test.QuickCheck+-- import qualified Test.Hspec.SmallCheck as SC+++import AERN2.Norm+import AERN2.MP.Precision++import AERN2.MP.UseMPFR.Float.Type+import AERN2.MP.UseMPFR.Float.Arithmetic+import AERN2.MP.UseMPFR.Float.Conversions+import AERN2.MP.UseMPFR.Float.Operators+import AERN2.MP.UseMPFR.Float.Constants++instance Arbitrary MPFloat where+  arbitrary =+    do+    giveSpecialValue <- frequencyElements [(9, False),(1, True)]+    aux giveSpecialValue+    where+      aux giveSpecialValue+        | giveSpecialValue =+            elements [nan, infinity, -infinity, zero, one, -one]+        | otherwise =+          do+          (p :: Precision) <- arbitrary+          (s :: Integer) <- arbitrary+          ex <- choose (-20,10)+          let resultR = s * (10.0^!ex)+          let result = fromRationalUp p resultR+          return result++frequencyElements :: ConvertibleExactly t Int => [(t, a)] -> Gen a+frequencyElements elems = frequency [(int n, return e) | (n,e) <- elems]++{- approximate comparison -}++infix 4 =~=++(=~=) :: MPFloat -> MPFloat -> Property+l =~= r =+  approxEqualWithArgs [] l r++{-|+  Assert equality of two MPFloat's with tolerance @1/2^p@.+-}+approxEqual ::+  Integer {-^ @p@ precision to guide tolerance -} ->+  MPFloat {-^ LHS of equation-} ->+  MPFloat {-^ RHS of equation -}->+  Bool+approxEqual e x y+  | isNaN x && isNaN y = True+  | isNaN x && isInfinite y = True+  | isInfinite x && isNaN y = True+  | isNaN x || isNaN y = False+  | isInfinite x || isInfinite y = x == y+  | otherwise =+      abs (x -. y) <= 0.5^!e++{-|+  Assert equality of two MPFloat's with tolerance derived from the size and precision+  of the given intermediate values.+  When the assertion fails, report the given values using the given names.+-}+approxEqualWithArgs ::+  [(MPFloat, String)] {-^ intermediate values from which to determine tolerance, their names to report when the equality fails -} ->+  MPFloat {-^ LHS of equation-} ->+  MPFloat {-^ RHS of equation -}->+  Property+approxEqualWithArgs argsPre l r =+  counterexample description $ approxEqual e l r+  where+    args = argsPre ++ [(l, "L"), (r, "R"), (abs (l-.r),"|L-R|")]+    e =+      (foldl min 1000000 $ catMaybes $ map getNminusP args)+      - (length argsPre)+    getNminusP (x,_) =+      case norm of+        NormZero -> Nothing -- ideally infinity+        NormBits b -> Just (pI-b-1)+      where+      norm = getNormLog x+      pI = integer $ getPrecision x+    description =+      printf "args:\n%s tolerance: <= %s (e=%d)" argsS (show (double (0.5^!e))) e+    argsS =+      unlines+        [printf "    %s = %s (p=%s)" argS (show arg) (show $ getPrecision arg) | (arg, argS) <- args]++{-|+  A runtime representative of type @MPFloat@.+  Used for specialising polymorphic tests to concrete types.+-}+tMPFloat :: T MPFloat+tMPFloat = T "MPFloat"++specMPFloat :: Spec+specMPFloat =+  describe ("MPFloat") $ do+    specCanSetPrecision tMPFloat (printArgsIfFails2 "=~=" (=~=))+    specCanRound tMPFloat+    specCanNegNum tMPFloat+    specCanAbs tMPFloat+    specCanMinMaxNotMixed tMPFloat+    -- specCanMinMax tMPFloat tInteger tMPFloat+    describe "special values" $ do+      it "0 * infinity = NaN" $ do+        isNaN (zero *^ infinity)+        &&+        isNaN (zero *. infinity)+      it "infinity / infinity = NaN" $ do+        isNaN (infinity /^ infinity)+        &&+        isNaN (infinity /. infinity)+      it "infinity - infinity = NaN" $ do+        isNaN (infinity -^ infinity)+        &&+        isNaN (infinity -. infinity)+    describe "approximate addition" $ do+      it "down <= up" $ do+        property $ \ (x :: MPFloat) (y :: MPFloat) ->+          not (isNaN (x +. y))+          ==>+          x +. y <= x +^ y+      it "up ~ down" $ do+        property $ \ (x :: MPFloat) (y :: MPFloat) ->+          x +. y =~= x +^ y+      it "absorbs 0" $ do+        property $ \ (x :: MPFloat) ->+          (not $ isNaN x) ==>+            x +. (mpFloat 0) == x+      it "approximately commutative" $ do+        property $ \ (x :: MPFloat) (y :: MPFloat) ->+          (not $ isNaN $ x +. y) ==>+          x +. y <= y +^ x+          &&+          x +^ y >= y +. x+      it "approximately associative" $ do+        property $ \ (x :: MPFloat) (y :: MPFloat) (z :: MPFloat) ->+          (not $ isNaN $ x +. y +. z) ==>+          (x +. y) +. z <= x +^ (y +^ z)+          &&+          (x +^ y) +^ z >= x +. (y +. z)+    describe "approximate subtraction" $ do+      it "down <= up" $ do+        property $ \ (x :: MPFloat) (y :: MPFloat) ->+          not (isNaN (x -. y))+          ==>+          x -. y <= x -^ y+      it "up ~ down" $ do+        property $ \ (x :: MPFloat) (y :: MPFloat) ->+          x -. y =~= x -^ y+      it "same as negate and add" $ do+        property $ \ (x :: MPFloat) (y :: MPFloat) ->+          (not $ isNaN $ x -. y) ==>+          x -. y <= x +^ (-y)+          &&+          x -^ y >= x +. (-y)+    describe "approximate multiplication" $ do+      it "down <= up" $ do+        property $ \ (x :: MPFloat) (y :: MPFloat) ->+          not (isNaN (x *. y))+          ==>+          x *. y <= x *^ y+      it "up ~ down" $ do+        property $ \ (x :: MPFloat) (y :: MPFloat) ->+          x *. y =~= x *^ y+      it "absorbs 1" $ do+        property $ \ (x :: MPFloat) ->+          (not $ isNaN x) ==>+            x *. (mpFloat 1) == x+      it "approximately commutative" $ do+        property $ \ (x :: MPFloat) (y :: MPFloat) ->+          not (isNaN (x *. y)) ==>+          x *. y <= y *^ x+          &&+          x *^ y >= y *. x+      it "approximately associative" $ do+        property $ \ (x :: MPFloat) (y :: MPFloat) (z :: MPFloat) ->+          (x >= 0 && y >= 0 && z >= 0+           && not (isInfinite x) && not (isInfinite y) && not (isInfinite z)) ==>+          (x *. y) *. z <= x *^ (y *^ z)+          &&+          (x *^ y) *^ z >= x *. (y *. z)+      it "approximately distributes over addition" $ do+        property $ \ (x :: MPFloat) (y :: MPFloat) (z :: MPFloat) ->+          (x >= 0 && y >= 0 && z >= 0+           && not (isInfinite x) && not (isInfinite y) && not (isInfinite z)) ==>+          x *. (y +. z) <= (x *^ y) +^ (x *^ z)+          &&+          x *^ (y +^ z) >= (x *. y) +. (x *. z)+    describe "approximate division" $ do+      it "down <= up" $ do+        property $ \ (x :: MPFloat) (y :: MPFloat) ->+          not (isNaN (x /. y))+          ==>+          x /. y <= x /^ y+      it "up ~ down" $ do+        property $ \ (x :: MPFloat) (y :: MPFloat) ->+          let+            (=~~=) = approxEqualWithArgs [(x /. y,"x/.y")]+            infix 4 =~~=+          in+          not (isNaN (x /. y))+          ==>+          x /. y =~~= x /^ y+      it "recip(recip x) = x" $ do+        property $ \ (x :: MPFloat) ->+          (x > 0 || x < 0) ==>+          one /. (one /^ x) <= x+          &&+          one /^ (one /. x) >= x+      it "x/1 = x" $ do+        property $ \ (x :: MPFloat) ->+          not (isNaN x) ==>+          (x /. one) == x+      it "x/x = 1" $ do+        property $ \ (x :: MPFloat) ->+          (isCertainlyNonZero x && (not $ isNaN $ x /. x)) ==>+            (x /. x) <= one+            &&+            (x /^ x) >= one+      it "x/y = x*(1/y)" $ do+        property $ \ (x :: MPFloat) (y :: MPFloat) ->+          (y > 0 && x >= 0 && x/.y >= 0) ==>+          (x /. y) <= x *^ (one /^ y)+          &&+          (x /^ y) >= x *. (one /. y)+    describe "approximate sqrt" $ do+      it "down <= up" $ do+        property $ \ (x :: MPFloat) ->+          not (isNaN (sqrtDown x))+          ==>+          sqrtDown x <= sqrtUp x+      it "up ~ down" $ do+        property $ \ (x :: MPFloat) ->+          (x >= 0)+          ==>+          sqrtDown x =~= sqrtUp x+      it "sqrt(x) >= 0" $ do+        property $ \ (x :: MPFloat) ->+          (x >= 0)+          ==>+          sqrtUp x >= 0+      it "sqrt(x)^2 ~ x" $ do+        property $ \ (x :: MPFloat) ->+          (x >= 0)+          ==>+          (sqrtDown x) *. (sqrtDown x) <= x+          &&+          (sqrtUp x) *^ (sqrtUp x) >= x+    describe "approximate exp" $ do+      it "down <= up" $ do+        property $ \ (x :: MPFloat) ->+          (abs x < 1000000)+          ==>+          expDown x <= expUp x+      it "up ~ down" $ do+        property $ \ (x :: MPFloat) ->+          (abs x < 1000000)+          ==>+          let+            (=~~=) = approxEqualWithArgs [(x,"x")]+            infix 4 =~~=+          in+          expDown x =~~= expUp x+      it "exp(-x) == 1/(exp x)" $ do+        property $ \ (x :: MPFloat) ->+          (abs x < 1000000)+          ==>+          one /. (expUp x) <= expUp (-x)+          &&+          one /^ (expDown x) >= expDown (-x)+      it "exp(x+y) = exp(x)*exp(y)" $ do+        property $ \ (x :: MPFloat) (y :: MPFloat) ->+          (abs x < 1000000 && abs y < 1000000)+          ==>+          expDown (x +. y) <= (expUp x) *^ (expUp y)+          &&+          expUp (x +^ y) >= (expDown x) *. (expDown y)+    describe "approximate log" $ do+      it "down <= up" $ do+        property $ \ (x :: MPFloat) ->+          (x > 0)+          ==>+          logDown x <= logUp x+      it "up ~ down" $ do+        property $ \ (x :: MPFloat) ->+          (x > 0)+          ==>+          logDown x =~= logUp x+      it "log(1/x) == -(log x)" $ do+        property $ \ (x :: MPFloat) ->+          (x > 0)+          ==>+          logDown (one /. x) <= -(logDown x)+          &&+          logUp (one /^ x) >= -(logUp x)+      it "log(x*y) = log(x)+log(y)" $ do+        property $ \ (x :: MPFloat) (y :: MPFloat) ->+          (x > 0 && y > 0)+          ==>+          logDown (x *. y) <= (logUp x) +^ (logUp y)+          &&+          logUp (x *^ y) >= (logDown x) +. (logDown y)+      it "log(exp x) == x" $ do+        property $ \ (x :: MPFloat) ->+          (abs x < 1000000)+          ==>+          logDown (expDown x) <= x+          &&+          logUp (expUp x) >= x+    describe "approximate sine" $ do+      it "down <= up" $ do+        property $ \ (x :: MPFloat) ->+          (abs x < 1000000)+          ==>+          sinDown x <= sinUp x+      it "up ~ down" $ do+        property $ \ (x :: MPFloat) ->+          (abs x < 1000000)+          ==>+          let+            (=~~=) = approxEqualWithArgs [(x,"x")]+            infix 4 =~~=+          in+          sinDown x =~~= sinUp x+      it "sin(pi)=0" $ do+        property $ \ (p :: Precision) ->+          let+            (=~~=) = approxEqualWithArgs [(piDown p,"pi")]+            infix 4 =~~=+          in+          sinUp(piDown p) =~~= (fromIntegerUp p 0)+      it "in [-1,1]" $ do+        property $ \ (x :: MPFloat) ->+          (abs x < 1000000)+          ==>+          sinDown x <= one+          &&+          sinUp x >= -one+    describe "approximate cosine" $ do+      it "down <= up" $ do+        property $ \ (x :: MPFloat) ->+          (abs x < 1000000)+          ==>+          cosDown x <= cosUp x+      it "up ~ down" $ do+        property $ \ (x :: MPFloat) ->+          (abs x < 1000000)+          ==>+          let+            (=~~=) = approxEqualWithArgs [(x,"x")]+            infix 4 =~~=+          in+          cosDown x =~~= cosUp x+      it "in [-1,1]" $ do+        property $ \ (x :: MPFloat) ->+          (abs x < 1000000)+          ==>+          cosDown x <= one+          &&+          cosUp x >= -one+      it "cos(pi)=-1" $ do+        property $ \ (p :: Precision) ->+          cosUp(piDown p) =~= (fromIntegerUp p (-1))+      it "cos(x)^2 + sin(x)^2 = 1" $ do+        property $ \ (x :: MPFloat) ->+          (abs x < 1000000)+          ==>+          let+            cosxU = cosUp x+            cosxD = cosDown x+            cosx2U = (cosxU *^ cosxU) `max` (cosxD *^ cosxD)+            cosx2D+              | cosxD > 0 = cosxD *. cosxD+              | cosxU < 0 = cosxU *. cosxU+              | otherwise = mpFloat 0+            sinxU = sinUp x+            sinxD = sinDown x+            sinx2U = (sinxU *^ sinxU) `max` (sinxD *^ sinxD)+            sinx2D+              | sinxD > 0 = sinxD *. sinxD+              | sinxU < 0 = sinxU *. sinxU+              | otherwise = mpFloat 0+          in+          (cosx2D +. sinx2D) <= one+          &&+          (cosx2U +^ sinx2U) >= one
+ src/AERN2/MP/UseMPFR/Float/Type.hs view
@@ -0,0 +1,91 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE DeriveGeneric, DeriveDataTypeable, StandaloneDeriving #-}+{-|+    Module      :  AERN2.MP.UseMPFR.Float.Type+    Description :  Arbitrary precision floating point numbers+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Arbitrary precision floating-point numbers++    Currently, we use hmpfr when compiling with ghc 7.10 and higher+    and haskell-mpfr when compiling with ghc 7.8.+-}++module AERN2.MP.UseMPFR.Float.Type+  (+   -- * MPFloat numbers and their basic operations+   MPFloat, setPrecisionUp, setPrecisionDown+   , p2mpfrPrec+   )+where++import MixedTypesNumPrelude+import qualified Prelude as P++import AERN2.MP.Precision++#ifdef IntegerBackend+import qualified AERN2.MP.UseMPFR.Float.Native as MPLow++type MPFloat = MPLow.MPFloat++p2mpfrPrec :: Precision -> Precision+p2mpfrPrec = id++#endif++#ifdef HaskellMPFR+import qualified Data.Approximate.MPFRLowLevel as MPLow+import Data.Typeable++{-| Multiple-precision floating-point type based on MPFR via haskell-mpfr. -}+type MPFloat = MPLow.Rounded++deriving instance (Typeable MPFloat)++p2mpfrPrec :: Precision -> MPLow.Precision+p2mpfrPrec = P.fromInteger . integer++#endif++#ifdef MPFRRounded+import qualified AERN2.MP.UseMPFR.Float.RoundedAdaptor as MPLow+import Data.Typeable++{-| Multiple-precision floating-point type based on MPFR via rounded. -}+type MPFloat = MPLow.Rounded++deriving instance (Typeable MPFloat)++p2mpfrPrec :: Precision -> MPLow.Precision+p2mpfrPrec = P.fromInteger . integer++#endif++#ifdef HMPFR+import qualified Data.Number.MPFR as MPLow++{-| Multiple-precision floating-point type based on MPFR via hmpfr. -}+type MPFloat = MPLow.MPFR++p2mpfrPrec :: Precision -> MPLow.Precision+p2mpfrPrec = P.fromInteger . integer++#endif++instance HasPrecision MPFloat where+  getPrecision x = prec (P.toInteger $ MPLow.getPrec x)++instance CanSetPrecision MPFloat where+  setPrecision = setPrecisionUp++setPrecisionUp :: Precision -> MPFloat -> MPFloat+setPrecisionUp p = MPLow.set MPLow.Up (p2mpfrPrec p)++setPrecisionDown :: Precision -> MPFloat -> MPFloat+setPrecisionDown p = MPLow.set MPLow.Down (p2mpfrPrec p)
+ src/AERN2/Norm.hs view
@@ -0,0 +1,67 @@+{-|+    Module      :  AERN2.Norm+    Description :  Rough logarithmic norm+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable+-}+module AERN2.Norm+(+    HasNorm(..), NormLog(..), invertNormLog+)+where++import MixedTypesNumPrelude+import qualified Prelude as P++import Data.Complex++import Math.NumberTheory.Logarithms (integerLog2)+++class HasNorm a where+    {-|+        For a value @x@, return @NormBits j@ where @j@ is close+        to the smallest @i@ with @|x| <= 2^i@.+        If @x == 0@ then return @NormZero@.+    -}+    getNormLog :: a -> NormLog++data NormLog+    = NormZero -- ^ ie NormBits (-infinity)+    | NormBits Integer+    deriving (P.Eq, P.Ord, Show)++instance HasEqAsymmetric NormLog NormLog+instance HasOrderAsymmetric NormLog NormLog+instance CanMinMaxAsymmetric NormLog NormLog++invertNormLog :: NormLog -> NormLog+invertNormLog NormZero = error "cannot invert NormZero"+invertNormLog (NormBits b) = NormBits (-b)++instance HasNorm Integer where+    getNormLog n+        | n == 0 = NormZero+        | abs n == 1 = NormBits 0+        | otherwise = NormBits $ 1 + (integer $ integerLog2 $ abs n - 1)++instance HasNorm Int where+    getNormLog = getNormLog . integer++instance HasNorm Rational where+    getNormLog x+        | x == 0.0 = NormZero+        | abs x >= 1.0 = getNormLog $ ceiling $ abs x+        | otherwise = NormBits $ negate $ integer $ integerLog2 $ floor $ (1 /! (abs x))++instance+  (HasNorm t)+  =>+  HasNorm (Complex t)+  where+  getNormLog (a :+ i) =+    (getNormLog a) `max` (getNormLog i)
+ src/AERN2/Normalize.hs view
@@ -0,0 +1,21 @@+{-|+    Module      :  AERN2.Normalize+    Description :  Overloadable "normalize" function+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable+-}+module AERN2.Normalize+(+    CanNormalize(..)+)+where++-- import MixedTypesNumPrelude+-- import qualified Prelude as P++class CanNormalize t where+  normalize :: t -> t
+ src/AERN2/Utils/Bench.hs view
@@ -0,0 +1,33 @@+{-|+    Module      :  AERN2.Utils.Bench+    Description :  utilities for benchmarks+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable+-}+module AERN2.Utils.Bench+(+    listFromGen+)+where++import Test.QuickCheck+import Test.QuickCheck.Random (mkQCGen)+import Test.QuickCheck.Gen (Gen(..))++import MixedTypesNumPrelude+-- import qualified Prelude as P++listFromGen :: Gen a -> [a]+listFromGen gen =+    list+    where+    list =+        concat $ map genSome [1..]+        where+        genSome size =+            unGen (sequence $ replicate 10 gen) qcGen (int size)+    qcGen = mkQCGen (int 148548830)
+ test/AERN2/MP/BallSpec.hs view
@@ -0,0 +1,20 @@+{-|+    Module      :  AERN2.MP.BallSpec+    Description :  hspec tests for MPBall+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable+-}++module AERN2.MP.BallSpec (spec) where++-- import MixedTypesNumPrelude+import AERN2.MP.Ball.Tests++import Test.Hspec++spec :: Spec+spec = specMPBall
+ test/AERN2/MP/DyadicSpec.hs view
@@ -0,0 +1,20 @@+{-|+    Module      :  AERN2.MP.DyadicSpec+    Description :  hspec tests for Dyadic+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable+-}++module AERN2.MP.DyadicSpec (spec) where++-- import MixedTypesNumPrelude+import AERN2.MP.Dyadic++import Test.Hspec++spec :: Spec+spec = specDyadic
+ test/AERN2/MP/UseMPFR/FloatSpec.hs view
@@ -0,0 +1,20 @@+{-|+    Module      :  AERN2.MP.FloatSpec+    Description :  hspec tests for MPFloat+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable+-}++module AERN2.MP.UseMPFR.FloatSpec (spec) where++-- import MixedTypesNumPrelude+import AERN2.MP.UseMPFR.Float.Tests++import Test.Hspec++spec :: Spec+spec = specMPFloat
+ test/Spec.hs view
@@ -0,0 +1,1 @@+{-# OPTIONS_GHC -F -pgmF hspec-discover #-}