aern2-real 0.2.1.0 → 0.2.4.0

raw patch · 9 files changed

+830/−158 lines, 9 filesdep +aern2-realdep ~collect-errorsdep ~mixed-types-numPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: aern2-real

Dependency ranges changed: collect-errors, mixed-types-num

API changes (from Hackage documentation)

Files

• README.md +198/−52
• aern2-real.cabal +54/−7
• changelog.md +6/−0
• examples/AERN2/Real/Examples/ClosestPairDist.hs +146/−0
• examples/AERN2/Real/Examples/Introduction.hs +294/−0
• src/AERN2/Real/Tests.hs +94/−94
• src/AERN2/Real/Type.hs +17/−5
• test/AERN2/RealSpec.hs +20/−0
• test/Spec.hs +1/−0

@@ -2,13 +2,15 @@  Exact real arithmetic -API documentation available on the [Hackage page](https://hackage.haskell.org/package/aern2-real).+API documentation is available on the [Hackage page](https://hackage.haskell.org/package/aern2-real). +The remainder of this text is an introductory tutorial.  The code for the examples contained here is also available in file [Introduction.hs](src/AERN2/Real/Introduction.hs).+ ## Table of contents <!-- omit in toc --> -- [1. Numeric data types](#1-numeric-data-types)-- [2. Basic usage with Prelude](#2-basic-usage-with-prelude)-- [3. Basic usage with MixedTypesNumPrelude](#3-basic-usage-with-mixedtypesnumprelude)+- [1. Data types](#1-data-types)+- [2. Usage with Prelude](#2-usage-with-prelude)+- [3. Usage with MixedTypesNumPrelude](#3-usage-with-mixedtypesnumprelude) - [4. Partial functions and error handling](#4-partial-functions-and-error-handling) - [5. Limits](#5-limits) - [6. Multivalued selection](#6-multivalued-selection)@@ -16,7 +18,7 @@   - [6.2. Multi-valued selection](#62-multi-valued-selection) - [7. Specification and tests](#7-specification-and-tests) -## 1. Numeric data types+## 1. Data types  This package provides the following two data types: @@ -27,52 +29,83 @@ The type `CReal` has instances of both [mixed-types-num](https://hackage.haskell.org/package/mixed-types-num) type classes such as `CanAdd`, `CanSqrt` as well as with traditional Prelude type classes such as `Ord`, `Num` and `Floating`. The type `CKleenean` supports the usual Boolean operations. -## 2. Basic usage with Prelude+Real numbers are represented by converging sequences of dyadic intervals: +```Haskell+type CReal = CSequence MPBall+```++A `CSequence` is a list of approximations computed with increasing *precision*.+Precision here does *not* guarantee a certain *accuracy*.+Precision roughly corresponds to the number of *significant digits* used+in all intermediate computations.+With increasing precision the intervals eventually converge to exact values.++The elements of a `CSequence` use the `CN` error-collecting wrapper.+A convergent sequence must be error-free from some point onwards.+A sequence is allowed not to converge, but only if all its elements contain the same error.  +Such a sequence can be thought of as converging to this error.+++## 2. Usage with Prelude+ First, let us load the package with **Prelude** operations:  ```Text-$ stack ghci aern2-real:lib --no-load --ghci-options AERN2.Real-*AERN2.MP> import Prelude hiding (pi)-*AERN2.MP Prelude>+$ stack ghci aern2-real:lib --no-load --ghci-options "AERN2.Real -Wno-type-defaults"+*AERN2.Real> import Prelude hiding (pi)+*AERN2.Real Prelude> ``` -We can obtain approximations of a real number with any **requested accuracies**:+We can obtain approximations of a real number with a chosen *precision*:  ```Text-...> pi ? (bits 1000)-[3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117... ± ~0.0000 ~2^(-1230)]+...> (sin 1 ::CReal) ? (prec 120)+[0.84147098480789650665250232... ± ~4.6644e-35 ~2^(-114)] -...> pi ? (bits 1000000)-[3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117... ± ~0.0000 ~2^(-1028468)]-(4.12 secs, 270,972,152 bytes)+...> (sin 1 ::CReal) ? (prec 10000)+[0.84147098480789650665250232... ± ~0.0000 ~2^(-13530)] ``` -Instead of accuracy, we can request that the computation is performed with a certain **precision**, which roughly corresponds to the number of significant bits.  This usually trades speed with guaranteed accuracy:+Notice that sometimes the accuracy of the interval is lower than the working precision.  Instead of precision, we can request that the computation is performed with a certain *guaranteed accuracy*:  ```Text-...> (sin pi) ? (bits 10000) -- guaranteed accuracy at least 10000-[-0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000... ± ~0.0000 ~2^(-13539)]-(0.27 secs, 196,580,192 bytes)+...> (sin 1 ::CReal) ? (bits 120)+[0.84147098480789650665250232... ± ~2.2431e-55 ~2^(-181)] -...> (sin pi) ? (prec 10000) -- no guaranteed accuracy-[-0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000... ± ~0.0000 ~2^(-13539)]-(0.21 secs, 107,844,784 bytes)+Nevertheless, this sometimes comes with a performance penalty, since internally the computation may need to be restarted with a higher accuracy:++...> sumSines n = sum [sin (creal i) | i <- [1..n::Integer]]++...> sumSines 100 ? (prec 120)+[-0.12717101366042011543675217... ± ~2.8393e-33 ~2^(-108)]+(0.03 secs, 26,203,776 bytes)++...> sumSines 100 ? (bits 120)+[-0.12717101366042011543675217... ± ~1.2220e-53 ~2^(-175)]+(0.05 secs, 60,537,128 bytes) ``` -When formatting a real number, a **default precision** is used:+Which can be obtained faster if directly guessing that we need precision at least 130:  ```Text+...> (sumSines1 100) ? (prec 130)+[-0.12717101366042011543675217... ± ~1.2220e-53 ~2^(-175)]+(0.03 secs, 35,209,088 bytes)+```++When formatting a real number, a *default precision* is used:++```Text ...> pi-{?(prec 36): [3.141592653584666550159454345703125 ± ~1.4552e-11 ~2^(-36)]}+{?(prec 36): [3.14159265358466655015945434... ± ~1.4552e-11 ~2^(-36)]} ``` -The Prelude power operator works only for integral types:+The **Prelude** power operator works only for integral types:  ```Text ...> pi ^ 2-[9.8696044010893586188344909998725639610631902560... ± ~8.1120e-30 ~2^(-96)]-{?(prec 36): [9.8696044009993784129619598388671875 ± ~1.4964e-10 ~2^(-32)]}+{?(prec 36): [9.86960440099937841296195983... ± ~1.4964e-10 ~2^(-32)]}  ...> pi ^ pi <interactive>:18:1: error:@@ -99,32 +132,54 @@ *** Exception: Failed to decide equality of Sequences.  If you switch to MixedTypesNumPrelude instead of Prelude, comparison of Sequences returns CSequence Kleenean or similar instead of Bool. ``` -## 3. Basic usage with MixedTypesNumPrelude+## 3. Usage with MixedTypesNumPrelude -We see that some things do not work with Prelude. Let us use **MixedTypesNumPrelude** operations instead:+We see that some things do not work with Prelude. Let us use [MixedTypesNumPrelude](https://hackage.haskell.org/package/mixed-types-num) operations instead:  ```Text $ stack ghci aern2-real:lib --no-load --ghci-options AERN2.Real-*AERN2.MP> import MixedTypesNumPrelude-*AERN2.MP MixedTypesNumPrelude>+*AERN2.Real> import MixedTypesNumPrelude+*AERN2.Real MixedTypesNumPrelude> ``` -We get a more general power operator:+First, our Prelude expressions +- `(sin 1 :: CReal)`+- `sum [sin (creal i) | i <- [1..n::Integer]]`++can now be simplified as follows:+ ```Text+...> :t sin 1+sin 1 :: CReal++...> sumSines n = sum [sin i | i <- [1..n]]+...> :t sumSines+sumSines :: Integer -> CReal+```++Moreover, we get a more general power operator:++```Text ...> 2^0.5-{?(prec 36): [1.414213562371930730340852514178195642186126256312482171419747717302107387071785637999710161238908... ± ~1.0305e-11 ~2^(-36)]}+{?(prec 36): [1.41421356237193073034085251... ± ~1.0305e-11 ~2^(-36)]}  ...> pi ^ pi-{?(prec 36): [36.462159605538498632520418490483602438877178488347481362195878876490337527904728176508797332644462... ± ~2.7112e-9 ~2^(-28)]}+{?(prec 36): [36.46215960553849863252041849... ± ~2.7112e-9 ~2^(-28)]}  ...> (pi ^ pi) ? (bits 10000)-[36.462159607207911770990826022692123666365508402228818738709335922934074368881699904620079875706774... ± ~0.0000 ~2^(-13532)]-(0.90 secs, 631,865,912 bytes)+[36.46215960720791177099082602... ± ~0.0000 ~2^(-13532)]+(0.83 secs, 631,232,904 bytes) ``` -Real comparison now returns a `CKleenean` instead of `Bool`, supporting undecided comparisons and comparisons with a specified precision:+Real comparison now returns a `CKleenean` instead of `Bool`, where +```Haskell+type CKleenean = CSequence Kleenean+```++As a three-value truth type, `Kleenean` supports undecided comparisons.  Being a sequence, `CKleenean` supports comparisons with a specified precision:+ ```Text ...> pi > 0 {?(prec 36): CertainTrue}@@ -139,34 +194,101 @@ CertainFalse ``` +When the numbers are known exactly, an equality test succeeds:++```Test+...> (creal 0) == 0+{?(prec 36): CertainTrue}+```+ ## 4. Partial functions and error handling -Since comparisons can be only semi-decided, also errors such as division by zero or logarithm of a negative number can be only semi-detected.-Therefore, an invalid input leads to a normal `CReal` value, and the error is demonstrated only when we extract an approximation, and sometimes an error cannot be determined with certainty:+Normally in Haskell, computation such as `1/0` or `sqrt (-1)` result in **NaN** or run-time exceptions.+Since `CReal` uses the [CN wrapper](https://hackage.haskell.org/package/collect-errors), for `CReal` these expressions instead return special values that describe the error. +Since comparisons can be only semi-decided, also such errors can be only semi-detected.+Therefore, an invalid input leads to a normal `CReal` value, and the error is demonstrated only when we extract an approximation:+ ```Text-...> bad1 = pi/0+...> bad1 = sqrt (-1) ...> bad1 ? (prec 100)-{{ERROR: division by 0}}}+{{ERROR: out of domain: negative sqrt argument}}+```+ +and sometimes an error cannot be determined with certainty: -...> bad2 = 1/(pi-pi)+```Text+...> a_third = creal (1/3)++...> bad2 = 1/(a_third-a_third) ...> bad2 ? (prec 100) {{POTENTIAL ERROR: division by 0}}++...> bad2 ? (bits 100)+{{POTENTIAL ERROR: numeric error: failed to find an approximation with sufficient accuracy}} ``` +A query for guaranteed precision may take a long time because before it fails, the computation is attempted iteratively for higher and higher precisions, up to precision around 5,000,000 bits:++```Text+...> bad3 = 1/(pi-pi)+...> bad3 ? (prec 100)+{{POTENTIAL ERROR: division by 0}}++...> bad3 ? (bits 100)+-- TAKES A VERY LONG TIME+```+ When we are sure that potential errors are harmless, we can clear them:  ```Text-...> ok3 = sqrt (pi-pi)-...> ok3 ? (prec 10)-[0.022097086912079610143710452219156792352805496193468570709228515625 ± ~2.2097e-2 ~2^(-5)]{{POTENTIAL ERROR: out of domain: negative sqrt argument}}-...> clearPotentialErrors $ ok3 ? (prec 10)-[0.022097086912079610143710452219156792352805496193468570709228515625 ± ~2.2097e-2 ~2^(-5)]+...> ok4 = sqrt (pi-pi)+...> ok4 ? (prec 100)+[0.00000000000000000061331736... ± ~6.1332e-19 ~2^(-60)]{{POTENTIAL ERROR: out of domain: negative sqrt argument}}}++...> ok5 = clearPotentialErrors $ sqrt (pi-pi)+...> ok5 ? (prec 100)+[0.00000000000000000061331736... ± ~6.1332e-19 ~2^(-60)] ``` +Attempting to clear a certain error is harmless:++```Text+...> bad6 = clearPotentialErrors (sqrt (pi-pi-1))+...> bad6 ? (prec 100)+{{ERROR: out of domain: negative sqrt argument}}+```++But clearing a potential error which is a real error is unsound:++```Text+...> bad7 = clearPotentialErrors (sqrt (pi-pi-2^(-1000)))+...> bad7 ? (prec 100)+[0.00000000000000000061331736... ± ~6.1332e-19 ~2^(-60)]+...> bad7 ? (prec 1000)+{{ERROR: out of domain: negative sqrt argument}}+```++Errors can be investigated, eg as follows:++```Text+...> detectCN r = if not (CN.hasError r) then Just r else Nothing++...> detectCN (sqrt (-1) ? (prec 100))+Nothing++...> detectCN (sqrt 0 ? (prec 100))+Just [0 ± 0]+```++There is also `CN.hasCertainError` which ignores potential errors.+ ## 5. Limits  Computing a limit of a fast converging sequence of numbers or functions is one of the most fundamental operations for real numbers.+A sequence `a_n` is fast converging if each+`a_n` is no more than `0.5^n` distant from the limit.+ For example, we can compute `e` as the limit of the partial sums of terms `1/n!` for `n` ranging from `0` onwards:  ```Text@@ -174,8 +296,18 @@ ... MixedTypesNumPrelude> e_sum n = sum $ map (recip . fact) [0..n] ``` -TODO+The difference between `e` and `e_sum n` is no more than `3/(fact (n+1))` which is less than `0.5^(n-2)`.+Thus the sequence `\n -> e_sum (n+2)` is fast converging and the following limit is valid: +```Text+...> my_e = limit $ \(n :: Integer) -> e_sum (n+2)++...> my_e ? (prec 1000)+[2.71828182845904523536028747... ± ~0.0000 ~2^(-1217)]+```++The type declaration for `n` is required because `limit` is generic and works also for sequences indexed by `Int` or even positive rational numbers.+ ## 6. Multivalued selection  When a comparison is needed for branching, its semi-decidability becomes a challenge.  As an example, consider the task of defining the `abs` function by cases.@@ -184,8 +316,9 @@ ### 6.1. Parallel branching  ```Text-... MixedTypesNumPrelude> abs1 x = if x < 0 then -x else x-... MixedTypesNumPrelude> abs1 (pi - pi)+...> absR1 x = if x < 0 then -x else x++...> absR1 (pi - pi) {?(prec 36): [0 ± ~2.9104e-11 ~2^(-35)]} ``` @@ -193,8 +326,21 @@  ### 6.2. Multi-valued selection -TODO+A more general mechanism for dealing with branching based on semi-decidable conditions such as real-number comparisons is non-deterministic `select`. If given two lazy Kleeneans, `select` will enquire them concurrently with increasing precisions until one of them becomes `CertainTrue`.  By convention `select` returns a `Bool` which is `True` if the first branch succeeds and `False` if the second branch succeeds.   +Here we use `select` to implement a *soft* sign test with some tolerance `eps` and define `absR2` to be the limit of a sequence of approximate implementations of `abs` with different `eps`:++```Text+...> absR2_approx x (q :: Rational) = if select (x > -q) (x < q) then x else -x++...> absR2 x = limit $ absR2_approx x++...> absR2 (pi - pi)+{?(prec 36): [0 ± ~4.3656e-11 ~2^(-34)]}+```+ ## 7. Specification and tests -The approximations obtained using `? (bits n)` or `? (prec p)` are intervals of type `CN MPBall` from package [aern2-mp](../aern2-mp/README.md).  This type is also used internally for all `CReal` arithmetic.  The `MPBall` arithmetic is tested against a fairly complete hspec/QuickCheck specification of algebraic properties.+Most `CReal` operations are simply lifts of the corresponding `CN MPBall` operations, which are tested in package [aern2-mp](../aern2-mp/README.md) against a fairly complete hspec/QuickCheck specification of algebraic properties.++TODO: limit and select tests

@@ -1,13 +1,13 @@ cabal-version: 1.12 --- This file has been generated from package.yaml by hpack version 0.33.0.+-- This file has been generated from package.yaml by hpack version 0.34.4. -- -- see: https://github.com/sol/hpack ----- hash: 2eeae9f84732bd487f03a8ca2f34af6a8da594cc92af23524cd6e594a530ad03+-- hash: dcd52485dd451831f9e0d0903913a199f2dd7764f22b7fc61f3a3b292ba85676  name:           aern2-real-version:        0.2.1.0+version:        0.2.4.0 synopsis:       Real numbers as sequences of MPBalls description:    Please see the README on GitHub at <https://github.com/michalkonecny/aern2/#readme> category:       Math@@ -38,19 +38,66 @@       AERN2.Real.Limit       AERN2.Real.Tests       AERN2.Real.Type+      AERN2.Real.Examples.ClosestPairDist+      AERN2.Real.Examples.Introduction   other-modules:       Paths_aern2_real   hs-source-dirs:       src-  default-extensions: RebindableSyntax, ScopedTypeVariables, TypeFamilies, TypeOperators, ConstraintKinds, DefaultSignatures, MultiParamTypeClasses, FlexibleContexts, FlexibleInstances, UndecidableInstances-  other-extensions: TemplateHaskell+      examples+  default-extensions:+      RebindableSyntax,+      ScopedTypeVariables,+      TypeFamilies,+      TypeOperators,+      ConstraintKinds,+      DefaultSignatures,+      MultiParamTypeClasses,+      FlexibleContexts,+      FlexibleInstances,+      UndecidableInstances+  other-extensions:+      TemplateHaskell   ghc-options: -Wall   build-depends:       QuickCheck     , aern2-mp >=0.2.1     , base ==4.*-    , collect-errors >=0.1+    , collect-errors >=0.1.5     , hspec     , integer-logarithms-    , mixed-types-num >=0.5.1+    , mixed-types-num >=0.5.3+  default-language: Haskell2010++test-suite aern2-real-test+  type: exitcode-stdio-1.0+  main-is: Spec.hs+  other-modules:+      AERN2.RealSpec+      Paths_aern2_real+  hs-source-dirs:+      test+  default-extensions:+      RebindableSyntax,+      ScopedTypeVariables,+      TypeFamilies,+      TypeOperators,+      ConstraintKinds,+      DefaultSignatures,+      MultiParamTypeClasses,+      FlexibleContexts,+      FlexibleInstances,+      UndecidableInstances+  other-extensions:+      TemplateHaskell+  ghc-options: -threaded -rtsopts -with-rtsopts=-N -Wall+  build-depends:+      QuickCheck+    , aern2-mp >=0.2.1+    , aern2-real+    , base ==4.*+    , collect-errors >=0.1.5+    , hspec+    , integer-logarithms+    , mixed-types-num >=0.5.3   default-language: Haskell2010

@@ -1,5 +1,11 @@ # Change log for aern2-real +* v 0.2.4 2021-05-26+  * overhaul README and examples+  * stop "very inaccurate" errors breaking ? (bits n) queries+  * optimisation: ? (bits n) queries start from precision n+  * add tests for accuracy queries, limit and select+  * fix div by 0 during low-accuracy integer powers * v 0.2.1 2021-05-18   * add conversion from WithAnyPrec * v 0.2.0 2021-05-17

@@ -0,0 +1,146 @@+{-# LANGUAGE PartialTypeSignatures #-}+{-# OPTIONS_GHC -Wno-missing-signatures #-}+{-# OPTIONS_GHC -Wno-partial-type-signatures #-}+{-|+    Module      :  AERN2.Real.Introduction+    Description :  aern2-real introductory examples+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    You can run the examples in this file in ghci.+    If you installed AERN2 using the official instructions,+    you can start ghci using the following command in the base+    folder:++    @+    stack repl aern2-real/examples/AERN2/Real/Examples/ClosestPairDist.hs+    @+-}+module AERN2.Real.Examples.ClosestPairDist where++import MixedTypesNumPrelude+-- import qualified Prelude as P+-- import Text.Printf++import Test.QuickCheck+import qualified Data.List as List++import AERN2.MP+import AERN2.Real++-- import Debug.Trace++-- define a short name for the type of real numbers:+type R = CReal++----------------------------------+-- Finding the smallest distance within a set of real numbers+----------------------------------++closestPairDist_naive ::+  _ => [t] -> t+closestPairDist_naive pts+  | length pts < 2 = error "closestPairDist_naive: too few points"+  | otherwise =+      (foldl1 min (map distance (distinctPairs pts)))++distance :: (CanSubSameType t, CanAbsSameType t) => (t, t) -> t+distance (a,b) = abs (a-b)++closestPairDist_run ::+  _ =>+  ([t] -> t) ->+  Integer -> t+closestPairDist_run (closestPairDist :: [t] -> t) n =+  closestPairDist [sin (convertExactly i :: t) | i <- [1..n]]++closestPairDist_run_naive :: Integer -> R+closestPairDist_run_naive =+  closestPairDist_run closestPairDist_naive ++closestPairDist_run_split :: Integer -> R+closestPairDist_run_split =+  closestPairDist_run $ closestPairDist_split compRApprox++{- Example runs:++*AERN2.Real.Examples.ClosestPairDist> closestPairDist_run_naive 1000 ? (prec 1000)+[0.00000013295546744391165086... ± ~0.0000 ~2^(-1221)]+(13.80 secs, 12,017,593,904 bytes)++*AERN2.Real.Examples.ClosestPairDist> closestPairDist_run_split 1000 ? (prec 1000)+[0.00000013295546744391165086... ± ~0.0000 ~2^(-1221)]+(4.95 secs, 9,979,768,504 bytes)++-}++{- specification and randomised tests -}++closestPairDist_spec closestPairDist (getFinite :: r -> t) numbers =+  (length numbers) < 2+  .||.+  (getFinite (closestPairDist numbersR)) ?==?$ (closestPairDist_naive numbers)+  where+  numbersR = map convertExactly numbers :: [r]+  a ?==?$ b = printArgsIfFails2 "?==?" (?==?) a b++closestPairDist_runTests1 =+  quickCheck (closestPairDist_spec (closestPairDist_split compRApprox) (?bits 100) :: [Integer] -> Property)+closestPairDist_runTests2 =+  quickCheck (closestPairDist_spec (closestPairDist_split compMPBall) id :: [Integer] -> Property)++sample_integers = sample' (arbitrary :: Gen [Integer]) >>= mapM_ print+sample_rationals = sample' (arbitrary :: Gen [Rational]) >>= mapM_ print++{- a version that splits, recurses and merges the results -}+closestPairDist_split ::+  _ => (t -> t -> Bool) -> [t] -> t+closestPairDist_split (.<) pts+  | length ptsL < 2 || length ptsR < 2 =+      closestPairDist_naive pts+  | otherwise =+      recurseAndMerge+  where+  (ptsL,ptsR) = List.partition isCertainlyLeft pts+    where+    isCertainlyLeft x = x .< average pts+  recurseAndMerge =+    foldl1 min [dL, dLR, dR]+    where+    dL = closestPairDist_split (.<) ptsL+    dLR = distance (largest ptsL, smallest ptsR)+    dR = closestPairDist_split (.<) ptsR++compRApprox :: R -> R -> Bool+compRApprox a b = (a?ac) !<! (b?ac)+  where+  ac = bits 100++compMPBall :: MPBall -> MPBall -> Bool+compMPBall = (!<!)++{- auxiliary functions -}++-- hull :: MPBall -> MPBall -> MPBall+-- hull = hullMPBall++average :: (HasIntegers t, CanAddSameType t, CanDivBy t Integer) => [t] -> t+average xs = (sum xs) / (length xs)++largest :: (CanMinMaxSameType t) => [t] -> t+largest pts = foldl1 max pts++smallest :: (CanMinMaxSameType t) => [t] -> t+smallest pts = foldl1 min pts++distinctPairs :: [t] -> [(t,t)]+distinctPairs xs = [(x,y) | (x:rest) <- tails1 xs, y <- rest]++{-| non-empty tails -}+tails1 :: [t] -> [[t]]+tails1 list =+  take (length list - 1) $ List.tails list

@@ -0,0 +1,294 @@+{-# LANGUAGE PartialTypeSignatures #-}+{-# OPTIONS_GHC -Wno-missing-signatures #-}+{-# OPTIONS_GHC -Wno-partial-type-signatures #-}+{-|+    Module      :  AERN2.Real.Examples.Introduction+    Description :  aern2-real introductory examples+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Introductory examples for packages aern2-mp and aern2-real.++    Please see aern2-real/README.md for explanations.++    You can run the following examples in ghci.+    If you installed AERN2 using the official instructions,+    you can start ghci using the following command in the base+    folder:++    @+    stack repl aern2-real/examples/AERN2/Real/Examples/Introduction.hs+    @+-}+module AERN2.Real.Examples.Introduction where++import MixedTypesNumPrelude++import qualified Numeric.CollectErrors as CN++import AERN2.MP+import AERN2.Real++-- import Debug.Trace++------------------------------+-- real numbers+------------------------------++-- Start with a simple real number:++sine1 = sin 1++sine1_run1 = sine1 ? (prec 120)+-- result: [0.84147098480789650665250232... ± ~4.6644e-35 ~2^(-114)]+sine1_run2 = sine1 ? (bits 120)+-- result: [0.84147098480789650665250232... ± ~2.2431e-55 ~2^(-181)]++-- Next, do a bit more work:++sumSines1 :: Integer -> CReal+sumSines1 n = sum [sin i | i <- [1..n]]++-- Request the above expression with n = 100 using roughly 100 significant binary digits:+sumSines1_run1 :: CN MPBall+sumSines1_run1 = (sumSines1 100) ? (prec 120)+{- ghci log:++*AERN2.Real.Introduction> sumSines1_run1+[-0.12717101366042011543675217... ± ~2.8393e-33 ~2^(-108)]+(0.03 secs, 26,203,776 bytes)+-}++-- Same as above but request guaranteed 100 bits of accuracy:+sumSines1_run2 = (sumSines1 100) ? (bits 100)+{- ghci log:++*AERN2.Real.Introduction> sumSines1_run2+[-0.12717101366042011543675217... ± ~2.8393e-33 ~2^(-108)]+(0.19 secs, 319,789,600 bytes)++This is considetably slower because there is some backtracking when target accuracy is not reached.  +-}++------------------------------+-- real number comparisons+------------------------------++{-+  First consider comparisons of real number approximations.+  These may be decided or undecided, using a 'Kleenean'.+-}++pi100 :: CN MPBall+pi100 = pi?(bits 100)++compare_run1 :: CN Kleenean+compare_run1 = pi100 > 0+-- returns: CertainTrue++compare_run2 :: CN Kleenean+compare_run2 = pi100 == pi100+-- returns: TrueOrFalse++compare_run3 :: CKleenean+compare_run3 = pi > 0+-- in ghci prints: {?(prec 36): CertainTrue}+-- (evaluated using default precision 36)++compare_run4 = pi == pi + 2^(-100)+-- in ghci prints: {?(prec 36): TrueOrFalse}++compare_run5 = (pi == pi + 2^(-100)) ? (prec 1000)+-- returns: CertainFalse++compare_run6 = (creal 0) == 0+-- in ghci prints: {?(prec 36): CertainTrue}+-- this is decided in finite time because 0 is represented exactly++compare_run7 = pi == pi ? (prec 10000)+-- returns: TrueOrFalse+-- (cannot confirm pi=pi in finite time)++------------------------------+-- checking partial functions+------------------------------++partialfn_bad1 = sqrt (-1)+{- ghci log:++*AERN2.Real.Introduction> partialfn_bad1 ? (bits 100)+{{ERROR: out of domain: negative sqrt argument}}++-}++a_third = creal (1/3)++partialfn_bad2 = 1/(a_third-a_third)+{- ghci log:++*AERN2.Real.Introduction> partialfn_bad2 ? (prec 100)+{{POTENTIAL ERROR: division by 0}}++*AERN2.Real.Introduction> partialfn_bad2 ? (bits 100)+{{POTENTIAL ERROR: numeric error: failed to find an approximation with sufficient accuracy}}++-}++partialfn_bad3 = 1/(pi-pi)+{- ghci log:++*AERN2.Real.Introduction> partialfn_bad3 ? (prec 100)+{{POTENTIAL ERROR: division by 0}}++*AERN2.Real.Introduction> partialfn_bad3 ? (bits 100)+-- TAKES A VERY LONG TIME++-}++{-+ When computing on approximations which do not have enough information+ to check whether an error occurs, we get a *potential* error:+-}++partialfn_ok4 = sqrt (pi-pi)+{- ghci log:++*AERN2.Real.Introduction> partialfn_ok4 ? (prec 100)+[0.00000000000000000061331736... ± ~6.1332e-19 ~2^(-60)]{{POTENTIAL ERROR: out of domain: negative sqrt argument}}+-}++partialfn_ok5 = clearPotentialErrors (sqrt (pi-pi))+{- ghci log:++*AERN2.Real.Introduction> partialfn_ok5 ? (prec 100)+[0.00000000000000000061331736... ± ~6.1332e-19 ~2^(-60)]++-}++partialfn_bad6 = clearPotentialErrors (sqrt (pi-pi-1))+{- ghci log:++*AERN2.Real.Introduction> partialfn_bad6 ? (prec 100) +{{ERROR: out of domain: negative sqrt argument}}++-}++partialfn_bad7 = clearPotentialErrors (sqrt (pi-pi-2^(-1000)))+{- ghci log:++*AERN2.Real.Introduction> partialfn_bad7 ? (prec 100)+[0.00000000000000000061331736... ± ~6.1332e-19 ~2^(-60)]++*AERN2.Real.Introduction> partialfn_bad7 ? (prec 1000)+{{ERROR: out of domain: negative sqrt argument}}++-}++detectCN :: CN.CanTestErrorsPresent a => a -> Maybe a+detectCN r = if not (CN.hasError r) then Just r else Nothing+{- ghci log:++*AERN2.Real.Introduction> detectCN (sqrt (-1) ? (prec 100))+Nothing++*AERN2.Real.Introduction> detectCN (sqrt 0 ? (prec 100))+Just [0 ± 0]++-}++---------------------------------+-- Computing limits+--------------------------------++fact :: Integer -> CReal+fact n = creal $ product [1..n]++e_sum :: Integer -> CReal+e_sum n = sum $ map (recip . fact) [0..n]++my_e :: CReal+my_e = limit $ \(n :: Integer) -> e_sum (n+2)++{- ghci log:++*AERN2.Real.Introduction> my_e ? (prec 1000)+[2.71828182845904523536028747... ± ~0.0000 ~2^(-1217)]++-}++-- a faster version:++e_sum2 :: Integer -> CReal+e_sum2 n = foldl aux (creal 1) $ reverse [1..n]+  where aux x m = 1 + x / m++my_e2 :: CReal+my_e2 = limit $ \(n :: Integer) -> e_sum2 (n+2)++---------------------------------+-- "parallel" branching for real numbers+--------------------------------++absQ :: Rational -> Rational+absQ x = if x < 0 then -x else x++absR1 :: CReal -> CReal+absR1 x = if x < 0 then -x else x+++pif_run1 = absR1 (pi-pi)+{- ghci log:++*AERN2.Real.Introduction> pif_run1+{?(prec 36): [0 ± ~2.9104e-11 ~2^(-35)]}++-}++-- pif_run2 = foldl1 (.) (replicate 100 (absR1 . (100*))) (pi-pi)++absR2_approx x (q :: Rational) = if select (x > -q) (x < q) then x else -x++absR2 :: CReal -> CReal+absR2 x = limit $ absR2_approx x++select_run1 = absR2 (pi-pi)+{- ghci log:++*AERN2.Real.Introduction> select_run1+{?(prec 36): [0 ± ~4.3656e-11 ~2^(-34)]}++-}++-----------------------------------------+-- Cauchy reals vs iRRAM style execution+-----------------------------------------++logistic1 :: _ => Rational -> Integer -> t -> t+logistic1 c n x0 =+  (foldl1 (.) (replicate n lg)) x0+  where+  lg x = c * x * (1-x)++logistic1_CReal_run :: Integer -> CReal+logistic1_CReal_run n = logistic1 3.82 n (creal 0.5)++-- TODO: define logistic1_iter++{-  Example uses:++*AERN2.Real.Examples.Introduction> logistic1_CReal_run 100 ? (bits 100)+[0.95087585116480286419338875... ± ~2.9792e-32 ~2^(-104)]+  +*AERN2.Real.Examples.Introduction> logistic1_CReal_run 10000 ? (bits 100)+[0.20775682944252359241450861... ± ~0.0000 ~2^(-2566)]+(2.06 secs, 2,970,188,704 bytes)++-}++{-+  Recommended further reading:  ClosestPairDist.hs+-}

@@ -14,12 +14,12 @@     To run the tests using stack, execute:      @-    stack test aern2-real --test-arguments "-a 1000 -m Real"+    stack test aern2-real --test-arguments "-a 1000 -m CReal"     @ -} module AERN2.Real.Tests   (-    -- specCauchyReal, tCReal+    specCReal, tCReal   ) where @@ -40,8 +40,12 @@ import AERN2.MP import AERN2.MP.Dyadic +import AERN2.Limit import AERN2.Real.Type+import AERN2.Real.CKleenean import AERN2.Real.Field ()+import AERN2.Real.Elementary ()+import AERN2.Real.Limit ()  instance Arbitrary CReal where   arbitrary =@@ -75,12 +79,12 @@         nextBit _ _ = error "in Arbitrary CReal"  arbitrarySmall :: (Arbitrary a, HasOrderCertainly a Integer) => Integer -> Gen a-arbitrarySmall limit = aux+arbitrarySmall bound = aux   where   aux =     do     x <- arbitrary-    if -limit !<=! x && x !<=! limit+    if -bound !<=! x && x !<=! bound       then return x       else aux @@ -92,9 +96,6 @@ tCReal :: T CReal tCReal = T "CReal" --- tCauchyRealAtAccuracy :: T CauchyRealAtAccuracy--- tCauchyRealAtAccuracy = T "CReal(ac)"- specCRrespectsAccuracy1 ::   String ->   (CReal -> CReal) ->@@ -109,6 +110,21 @@           Right v -> getAccuracy v >=$ ac           _ -> property True +specCRrespectsAccuracy2 ::+  String ->+  (CReal -> CReal -> CReal) ->+  (CReal -> Accuracy -> Bool) ->+  (CReal -> Accuracy -> Bool) ->+  Spec+specCRrespectsAccuracy2 opName op precond1 precond2 =+  it (opName ++ " respects accuracy requests") $ do+    property $+      \ (x :: CReal) (y :: CReal) (ac :: Accuracy) ->+        ac < (bits 1000) && precond1 x ac && precond2 y ac ==>+        case CN.toEither ((op x y) ? ac) of+          Right v -> getAccuracy v >=$ ac+          _ -> property True+ (>=$) :: Accuracy -> Accuracy -> Property (>=$) = printArgsIfFails2 ">=" (>=) @@ -137,50 +153,22 @@ -- specCRrespectsAccuracy2 opName op = --   specCRrespectsAccuracy2CN opName (\ a b -> cn (op a b)) --- specCRrespectsAccuracy2CN ::---   String ->---   (CReal -> CReal -> CauchyRealCN) ->---   (CReal -> Accuracy -> Bool) ->---   (CReal -> Accuracy -> Bool) ->---   Spec--- specCRrespectsAccuracy2CN opName op precond1 precond2 =---   it (opName ++ " respects accuracy requests") $ do---     property $---       \ (x :: CReal) (y :: CReal) (ac :: Accuracy) ->---         let acSG = accuracySG ac in---         ac < (bits 1000) && precond1 x acSG && precond2 y acSG  ==>---         case getMaybeValueCN ((op x y) ? acSG) of---           Just v -> getAccuracy v >=$ ac-          -- _ -> property True---- specCRrespectsAccuracy2T ::---   (Arbitrary t, Show t) =>---   T t ->---   String ->---   (CReal -> t -> CReal) ->---   (CReal -> Accuracy -> Bool) ->---   (t -> Bool) ->---   Spec--- specCRrespectsAccuracy2T tt opName op =---   specCRrespectsAccuracy2TCN tt opName (\ a b -> cn (op a b))---- specCRrespectsAccuracy2TCN ::---   (Arbitrary t, Show t) =>---   T t ->---   String ->---   (CReal -> t -> CauchyRealCN) ->---   (CReal -> Accuracy -> Bool) ->---   (t -> Bool) ->---   Spec--- specCRrespectsAccuracy2TCN (T tName :: T t) opName op precond1 precond2 =---   it (opName ++ " with " ++ tName ++ " respects accuracy requests") $ do---     property $---       \ (x :: CReal) (t :: t) (ac :: Accuracy) ->---         let acSG = accuracySG ac in---         ac < (bits 1000) && precond1 x acSG && precond2 t  ==>---         case getMaybeValueCN ((op x t) ? acSG) of---           Just v -> getAccuracy v >=$ ac---           _ -> property True+specCRrespectsAccuracy2T ::+  (Arbitrary t, Show t) =>+  T t ->+  String ->+  (CReal -> t -> CReal) ->+  (CReal -> Accuracy -> Bool) ->+  (t -> Bool) ->+  Spec+specCRrespectsAccuracy2T  (T tName :: T t) opName op precond1 precond2 =+  it (opName ++ " with " ++ tName ++ " respects accuracy requests") $ do+    property $+      \ (x :: CReal) (t :: t) (ac :: Accuracy) ->+        ac < (bits 1000) && precond1 x ac && precond2 t ==>+        case CN.toEither ((op x t) ? ac) of+          Right v -> getAccuracy v >=$ ac+          _ -> property True  precondAnyT :: t -> Bool precondAnyT _t = True@@ -191,46 +179,58 @@ precondSmallT :: (HasOrderCertainly t Integer) => t -> Bool precondSmallT t = -1000 !<=! t && t !<=! 1000 --- specCauchyReal :: Spec--- specCauchyReal =---   describe ("CReal") $ do---     -- specConversion tInteger tCauchyReal real (fst . integerBounds)---     describe "order" $ do---       specHasEqNotMixed tCReal---       -- specHasEq tInt tCReal tRational---       -- specCanPickNonZero tCReal---       specHasOrderNotMixed tCReal---       -- specHasOrder tInt tCReal tRational-    -- describe "min/max/abs" $ do-    --   specCRrespectsAccuracy1 "abs" abs precondAnyReal-    --   specCRrespectsAccuracy2 "max" max precondAnyReal precondAnyReal-    --   specCRrespectsAccuracy2 "min" min precondAnyReal precondAnyReal-    -- describe "ring" $ do-    --   specCRrespectsAccuracy1 "negate" negate precondAnyReal-    --   specCRrespectsAccuracy2 "+" add precondAnyReal precondAnyReal-    --   specCRrespectsAccuracy2T tInteger "+" add precondAnyReal precondAnyT-    --   specCRrespectsAccuracy2T tRational "+" add precondAnyReal precondAnyT-    --   specCRrespectsAccuracy2T tDyadic "+" add precondAnyReal precondAnyT-    --   specCRrespectsAccuracy2 "a-b" sub precondAnyReal precondAnyReal-    --   specCRrespectsAccuracy2T tInteger "a-b" sub precondAnyReal precondAnyT-    --   specCRrespectsAccuracy2T tRational "a-b" sub precondAnyReal precondAnyT-    --   specCRrespectsAccuracy2T tDyadic "a-b" sub precondAnyReal precondAnyT-    --   specCRrespectsAccuracy2 "*" mul precondAnyReal precondAnyReal-    --   specCRrespectsAccuracy2T tInteger "*" mul precondAnyReal precondAnyT-    --   specCRrespectsAccuracy2T tRational "*" mul precondAnyReal precondAnyT-    --   specCRrespectsAccuracy2T tDyadic "*" mul precondAnyReal precondAnyT-    -- describe "field" $ do-    --   specCRrespectsAccuracy2CN "/" divide precondAnyReal precondNonZeroReal-    --   specCRrespectsAccuracy2TCN tInteger "/" divide precondAnyReal precondNonZeroT-    --   specCRrespectsAccuracy2TCN tRational "/" divide precondAnyReal precondNonZeroT-    --   specCRrespectsAccuracy2TCN tDyadic "/" divide precondAnyReal precondNonZeroT-    -- describe "elementary" $ do-    --   specCRrespectsAccuracy1CN "sqrt" sqrt precondPositiveReal-    --   specCRrespectsAccuracy1 "exp" exp precondSmallReal-    --   specCRrespectsAccuracy1CN "log" log precondPositiveSmallReal-    --   specCRrespectsAccuracy2CN "pow" pow precondPositiveSmallReal precondSmallReal-    --   specCRrespectsAccuracy2TCN tInteger "pow" pow precondNonZeroReal precondSmallT-    --   specCRrespectsAccuracy2TCN tRational "pow" pow precondPositiveSmallReal precondSmallT-    --   specCRrespectsAccuracy2TCN tDyadic "pow" pow precondPositiveSmallReal precondSmallT-    --   specCRrespectsAccuracy1 "cos" cos precondAnyReal-    --   specCRrespectsAccuracy1 "sine" sin precondAnyReal+specCReal :: Spec+specCReal =+  describe ("CReal") $ do+    -- specConversion tInteger tCReal creal (fst . integerBounds)+    -- describe "order" $ do+    --   specHasEqNotMixed tCReal+    --   specHasEq tInt tCReal tRational+    --   specCanPickNonZero tCReal+    --   specHasOrderNotMixed tCReal+    --   specHasOrder tInt tCReal tRational+    describe "min/max/abs" $ do+      specCRrespectsAccuracy1 "abs" abs precondAnyReal+      specCRrespectsAccuracy2 "max" max precondAnyReal precondAnyReal+      specCRrespectsAccuracy2 "min" min precondAnyReal precondAnyReal+    describe "ring" $ do+      specCRrespectsAccuracy1 "negate" negate precondAnyReal+      specCRrespectsAccuracy2 "+" add precondAnyReal precondAnyReal+      specCRrespectsAccuracy2T tInteger "+" add precondAnyReal precondAnyT+      specCRrespectsAccuracy2T tRational "+" add precondAnyReal precondAnyT+      specCRrespectsAccuracy2T tDyadic "+" add precondAnyReal precondAnyT+      specCRrespectsAccuracy2 "a-b" sub precondAnyReal precondAnyReal+      specCRrespectsAccuracy2T tInteger "a-b" sub precondAnyReal precondAnyT+      specCRrespectsAccuracy2T tRational "a-b" sub precondAnyReal precondAnyT+      specCRrespectsAccuracy2T tDyadic "a-b" sub precondAnyReal precondAnyT+      specCRrespectsAccuracy2 "*" mul precondAnyReal precondAnyReal+      specCRrespectsAccuracy2T tInteger "*" mul precondAnyReal precondAnyT+      specCRrespectsAccuracy2T tRational "*" mul precondAnyReal precondAnyT+      specCRrespectsAccuracy2T tDyadic "*" mul precondAnyReal precondAnyT+    describe "field" $ do+      specCRrespectsAccuracy2 "/" divide precondAnyReal precondNonZeroReal+      specCRrespectsAccuracy2T tInteger "/" divide precondAnyReal precondNonZeroT+      specCRrespectsAccuracy2T tRational "/" divide precondAnyReal precondNonZeroT+      specCRrespectsAccuracy2T tDyadic "/" divide precondAnyReal precondNonZeroT+    describe "elementary" $ do+      specCRrespectsAccuracy1 "sqrt" sqrt precondPositiveReal+      specCRrespectsAccuracy1 "exp" exp precondSmallReal+      specCRrespectsAccuracy1 "log" log precondPositiveSmallReal+      specCRrespectsAccuracy2 "pow" pow precondPositiveSmallReal precondSmallReal+      specCRrespectsAccuracy2T tInteger "pow" pow precondNonZeroReal precondSmallT+      specCRrespectsAccuracy2T tRational "pow" pow precondPositiveSmallReal precondSmallT+      -- specCRrespectsAccuracy2T tDyadic "pow" pow precondPositiveSmallReal precondSmallT+      specCRrespectsAccuracy1 "cos" cos precondAnyReal+      specCRrespectsAccuracy1 "sine" sin precondAnyReal+    describe "select" $ do+      it "soft abs via select" $ do+        property $ \ (x :: CReal) (p :: Precision) (q :: Rational) ->+          (1 < q) ==>+          let eps = 1/q in+          (abs (abs x - (if select (x > -eps) (x < eps) then x else -x)) ? p) ?<? 2*eps+    describe "limit" $ do+      it "computing e as a limit of Taylor series" $ do+        property $ \ (p :: Precision) ->+          ((exp (mpBallP p 1.0)) ?==?) $+            (limit $ \(n :: Integer) -> sum $ map (recip . creal) $ take (n+3) $ scanl (*) 1 [1..(n)]) ? p+

@@ -75,9 +75,15 @@   where   withP p = runWithPrec p withCurrentP :: CN b -crealFromWithCurrentPrec :: (forall p. (KnownNat p) => WithCurrentPrec (CN MPBall) p) -> CReal-crealFromWithCurrentPrec = cseqFromWithCurrentPrec+{- Error handling -} +instance CN.CanTakeErrors CN.NumErrors (CSequence t) where+  takeErrors es (CSequence s) = CSequence $ map (CN.takeErrors es) s+  takeErrorsNoValue es = CSequence $ repeat (CN.takeErrorsNoValue es)++instance CN.CanClearPotentialErrors (CSequence t) where+  clearPotentialErrors (CSequence s) = CSequence $ map clearPotentialErrors s+ {- Cauchy real numbers -}  type CReal = CSequence MPBall@@ -106,8 +112,14 @@  instance (HasAccuracy t) => CanExtractApproximation (CSequence t) Accuracy where   type ExtractedApproximation (CSequence t) Accuracy = CN t-  extractApproximation (CSequence s) ac = aux s-    where+  extractApproximation (CSequence s) ac = +    aux $ drop (cseqIndexForPrecision p - 1) s+    where    +    p = +      case ac of +        Exact -> defaultPrecision+        NoInformation -> prec 2+        _ -> ac2prec ac     aux (bCN : rest)        | CN.hasCertainError bCN = bCN       | getAccuracy bCN >= ac = bCN@@ -147,7 +159,7 @@   safeConvertExactly = safeConvertExactly . rational  instance ConvertibleExactly (WithAnyPrec (CN MPBall)) CReal where-  safeConvertExactly (WithAnyPrec wcp) = Right $ crealFromWithCurrentPrec wcp+  safeConvertExactly (WithAnyPrec wcp) = Right $ cseqFromWithCurrentPrec wcp  _example1 :: CReal _example1 = creal 1.0

@@ -0,0 +1,20 @@+{-|+    Module      :  AERN2.RealSpec+    Description :  hspec tests for CauchyReal+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable+-}++module AERN2.RealSpec (spec) where++-- import MixedTypesNumPrelude+import AERN2.Real.Tests++import Test.Hspec++spec :: Spec+spec = specCReal

@@ -0,0 +1,1 @@+{-# OPTIONS_GHC -F -pgmF hspec-discover #-}