diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -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 -- [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 :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 diff --git a/aern2-real.cabal b/aern2-real.cabal --- a/aern2-real.cabal +++ b/aern2-real.cabal @@ -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 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 diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -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 diff --git a/examples/AERN2/Real/Examples/ClosestPairDist.hs b/examples/AERN2/Real/Examples/ClosestPairDist.hs new file mode 100644 --- /dev/null +++ b/examples/AERN2/Real/Examples/ClosestPairDist.hs @@ -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) ! MPBall -> Bool +compMPBall = (! 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 diff --git a/examples/AERN2/Real/Examples/Introduction.hs b/examples/AERN2/Real/Examples/Introduction.hs new file mode 100644 --- /dev/null +++ b/examples/AERN2/Real/Examples/Introduction.hs @@ -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 +-} diff --git a/src/AERN2/Real/Tests.hs b/src/AERN2/Real/Tests.hs --- a/src/AERN2/Real/Tests.hs +++ b/src/AERN2/Real/Tests.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) ? + ((exp (mpBallP p 1.0)) ?==?) $ + (limit $ \(n :: Integer) -> sum $ map (recip . creal) $ take (n+3) $ scanl (*) 1 [1..(n)]) ? p + diff --git a/src/AERN2/Real/Type.hs b/src/AERN2/Real/Type.hs --- a/src/AERN2/Real/Type.hs +++ b/src/AERN2/Real/Type.hs @@ -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 diff --git a/test/AERN2/RealSpec.hs b/test/AERN2/RealSpec.hs new file mode 100644 --- /dev/null +++ b/test/AERN2/RealSpec.hs @@ -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 diff --git a/test/Spec.hs b/test/Spec.hs new file mode 100644 --- /dev/null +++ b/test/Spec.hs @@ -0,0 +1,1 @@ +{-# OPTIONS_GHC -F -pgmF hspec-discover #-}