aern2real
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
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 of CReal
s
The type CReal
has instances of both mixedtypesnum 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
errorcollecting wrapper.
A convergent sequence must be errorfree 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 aern2real:lib noload ghcioptions "AERN2.Real Wnotypedefaults"
*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.6644e35 ~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.2431e55 ~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.8393e33 ~2^(108)]
(0.03 secs, 26,203,776 bytes)
...> sumSines 100 ? (bits 120)
[0.12717101366042011543675217... ± ~1.2220e53 ~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.2220e53 ~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.4552e11 ~2^(36)]}
The Prelude power operator works only for integral types:
...> pi ^ 2
{?(prec 36): [9.86960440099937841296195983... ± ~1.4964e10 ~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 aern2real:lib noload ghcioptions 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.0305e11 ~2^(36)]}
...> pi ^ pi
{?(prec 36): [36.46215960553849863252041849... ± ~2.7112e9 ~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 threevalue 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 runtime exceptions.
Since CReal
uses the CN wrapper, for CReal
these expressions instead return special values that describe the error.
Since comparisons can be only semidecided, also such errors can be only semidetected.
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_thirda_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/(pipi)
...> 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 (pipi)
...> ok4 ? (prec 100)
[0.00000000000000000061331736... ± ~6.1332e19 ~2^(60)]{{POTENTIAL ERROR: out of domain: negative sqrt argument}}}
...> ok5 = clearPotentialErrors $ sqrt (pipi)
...> ok5 ? (prec 100)
[0.00000000000000000061331736... ± ~6.1332e19 ~2^(60)]
Attempting to clear a certain error is harmless:
...> bad6 = clearPotentialErrors (sqrt (pipi1))
...> 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 (pipi2^(1000)))
...> bad7 ? (prec 100)
[0.00000000000000000061331736... ± ~6.1332e19 ~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^(n2)
.
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 semidecidability 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.9104e11 ~2^(35)]}
This simple definition works even when x = 0 because AERN2 has redefined the ifthenelse 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. Multivalued selection
A more general mechanism for dealing with branching based on semidecidable conditions such as realnumber comparisons is nondeterministic 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.3656e11 ~2^(34)]}
7. Specification and tests
Most CReal
operations are simply lifts of the corresponding CN MPBall
operations, which are tested in package aern2mp 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
.