aern2-real
Real numbers as convergent sequences of intervals
https://github.com/michalkonecny/aern2#readme
LTS Haskell 22.33: | 0.2.15 |
Stackage Nightly 2024-09-07: | 0.2.15 |
Latest on Hackage: | 0.2.15 |
aern2-real-0.2.15@sha256:98e05c4c1f21e5a2904d999d7916d688bb8d88664b5d278f898ad317f4cff69f,2531
Module documentation for 0.2.15
aern2-real
Exact real arithmetic
API documentation is available on the Hackage page.
The remainder of this text is an introductory tutorial. The code for the examples contained here is also available in file Introduction.hs.
Table of contents
- 1. Data types
- 2. Usage with Prelude
- 3. Usage with MixedTypesNumPrelude
- 4. Partial functions and error handling
- 5. Limits
- 6. Multivalued selection
- 7. Specification and tests
1. Data types
This package provides the following two data types:
-
CReal
: Exact real numbers via lazy sequences of interval approximations -
CKleenean
: Lazy Kleeneans, naturally arising from comparisons ofCReal
s
The type CReal
has instances of both 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.
Real numbers are represented by converging sequences of dyadic intervals:
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:
$ 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 a chosen precision:
...> (sin 1 ::CReal) ? (prec 120)
[0.84147098480789650665250232... ± ~4.6644e-35 ~2^(-114)]
...> (sin 1 ::CReal) ? (prec 10000)
[0.84147098480789650665250232... ± ~0.0000 ~2^(-13530)]
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:
...> (sin 1 ::CReal) ? (bits 120)
[0.84147098480789650665250232... ± ~2.2431e-55 ~2^(-181)]
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)
Which can be obtained faster if directly guessing that we need precision at least 130:
...> (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:
...> pi
{?(prec 36): [3.14159265358466655015945434... ± ~1.4552e-11 ~2^(-36)]}
The Prelude power operator works only for integral types:
...> pi ^ 2
{?(prec 36): [9.86960440099937841296195983... ± ~1.4964e-10 ~2^(-32)]}
...> pi ^ pi
<interactive>:18:1: error:
• No instance for (Integral CReal) arising from a use of ‘^’
Numerical order cannot be decided when the two numbers are equal:
...> pi > 0
True
...> pi == pi
*** 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.
Prelude comparison fails to determine also inequality when the two numbers are very close:
...> pi == pi + 0.00000000000000000000000000000000001
False
...> pi == pi + 0.00000000000000000000000000000000000001
*** 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. Usage with MixedTypesNumPrelude
We see that some things do not work with Prelude. Let us use MixedTypesNumPrelude operations instead:
$ stack ghci aern2-real:lib --no-load --ghci-options AERN2.Real
*AERN2.Real> import MixedTypesNumPrelude
*AERN2.Real MixedTypesNumPrelude>
First, our Prelude expressions
(sin 1 :: CReal)
sum [sin (creal i) | i <- [1..n::Integer]]
can now be simplified as follows:
...> :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:
...> 2^0.5
{?(prec 36): [1.41421356237193073034085251... ± ~1.0305e-11 ~2^(-36)]}
...> pi ^ pi
{?(prec 36): [36.46215960553849863252041849... ± ~2.7112e-9 ~2^(-28)]}
...> (pi ^ pi) ? (bits 10000)
[36.46215960720791177099082602... ± ~0.0000 ~2^(-13532)]
(0.83 secs, 631,232,904 bytes)
Real comparison now returns a CKleenean
instead of Bool
, where
type CKleenean = CSequence Kleenean
As a three-value truth type, Kleenean
supports undecided comparisons. Being a sequence, CKleenean
supports comparisons with a specified precision:
...> pi > 0
{?(prec 36): CertainTrue}
...> pi == pi
{?(prec 36): TrueOrFalse}
...> pi == pi + 2^(-100)
{?(prec 36): TrueOrFalse}
...> (pi == pi + 2^(-100)) ? (prec 1000)
CertainFalse
When the numbers are known exactly, an equality test succeeds:
...> (creal 0) == 0
{?(prec 36): CertainTrue}
4. Partial functions and error handling
Normally in Haskell, computation such as 1/0
or sqrt (-1)
result in NaN or run-time exceptions.
Since CReal
uses the CN wrapper, 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:
...> bad1 = sqrt (-1)
...> bad1 ? (prec 100)
{{ERROR: out of domain: negative sqrt argument}}
and sometimes an error cannot be determined with certainty:
...> 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:
...> 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:
...> 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:
...> 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:
...> 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:
...> 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:
... MixedTypesNumPrelude> fact n = creal $ product [1..n]
... MixedTypesNumPrelude> e_sum n = sum $ map (recip . fact) [0..n]
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:
...> 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.
We have two ways to overcome the challenge:
6.1. Parallel branching
...> absR1 x = if x < 0 then -x else x
...> absR1 (pi - pi)
{?(prec 36): [0 ± ~2.9104e-11 ~2^(-35)]}
This simple definition works even when x = 0 because AERN2 has redefined the if-then-else operator for a CKleenean
condition and real number branches in such a way that in situations where the condition is inconclusive, both branches are computed and the results safely merged. This is convenient, but can lead to inefficient code if the number of branches that need to be considered grows large.
6.2. Multi-valued selection
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
:
...> 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
Most CReal
operations are simply lifts of the corresponding CN MPBall
operations, which are tested in package aern2-mp against a fairly complete hspec/QuickCheck specification of algebraic properties.
There are also tests specific to CReal
, mostly checking that ? (bits n)
queries return sufficiently accurate approximations. There is are tests for select
and limit
.
Changes
Change log for aern2-real
- current
- v 0.2.15 2023-04-11
- CanSelectCountable instance for CKleenean
- v 0.2.14.1 2023-04-10
- fix compilation errors
- v 0.2.13 2023-04-07
- unsafeApproximationExtension
- v 0.2.12 2023-04-07
- more instances for if-then-else with CKleenean
- or_countable and and_countable for CKleenean
- v 0.2.11 2022-08-25
- left-first and/or for CKleeneans
- v 0.2.10 2022-08-20
- add to RealNumber: HasLimitsSameType, CanSelectCNBool, same compare types
- SelectType CKleenean = CN Bool
- v 0.2.9.2 2022-08-19
- add to RealNumber: HasIfThenElse for its SelectType
- v 0.2.9 2022-07-13 (update in aern2-mp)
- v 0.2.8 2021-08-04
- compatibility with ghc 9.0.1
- v 0.2.7 2021-06-02
- add RealNumber type class
- v 0.2.6 2021-05-29
- adapt to new ppow operations
- v 0.2.5 2021-05-27
- CanSelect moved to aern2-mp
- 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
- moving Arrow-based functionality to package aern2-net
- replacing Arrow-based sequences by list-based sequences
- switch to new simplified collect-errors, mixed-types-num 0.5.0
- got rid of EnsureCE etc.
- not introducing CN wrapper unless at least one parameter is already CN
- v 0.1.2 2019-03-19
- adapts to mixed-types-num 0.3.2 (new divI, mod)
- v 0.1.1.0 2017-12-06
- disable aern2-generate-netlog-elm for now
- v 0.1.0.3 2017-12-06
- remove further upper bounds
- v 0.1.0.2 2017-11-14
- remove most upper bounds, building with ghc 8.2
- v 0.1.0.1 2017-09-12
- first release on Hackage
- fast convergent sequences indexed by AccuracySG
- arrow-based networks of nodes communicating via query-answer protocols
- networks executable with cached and parallel strategies
- network execution can be visualised in browser using an Elm frontend