# The lazy numbers in R

Posted on November 12, 2022 by Stéphane Laurent
Tags: R, C++, Rcpp, maths

My new package lazyNumbers is a port of the lazy numbers implemented in the C++ library CGAL. The lazy numbers represent the real numbers and arithmetic on these numbers is exact.

The ordinary floating-point arithmetic is not exact. Consider for example the simple operation $$1 - 7 \times 0.1$$. In double precision, it is not equal to $$0.3$$:

x <- 1 - 7 * 0.1
x == 0.3
##  FALSE

With the lazy numbers, it is equal to $$0.3$$:

library(lazyNumbers)
x <- lazynb(1) - lazynb(7) * lazynb(0.1)
as.double(x) == 0.3
##  TRUE

Correction: this is wrong!

In fact, when a binary operation involves a lazy number, the other number is automatically converted to a lazy number, so you can shortly enter this operation as follows:

x <- 1 - lazynb(7) * 0.1
as.double(x) == 0.3
##  TRUE

Let’s see a more dramatic example now. Consider the sequence $$(u_n)$$ recursively defined by $$u_1 = 1/7$$ and $$u_{n+1} = 8 u_n - 1$$. You can easily check that $$u_2 = 1/7$$, therefore $$u_n = 1/7$$ for every $$n \geqslant 1$$. However, when implemented in double precision, this sequence quickly goes crazy ($$1/7 \approx 0.1428571$$):

u <- function(n) {
if(n == 1) {
return(1 / 7)
}
8 * u(n-1) - 1
}
u(15)
##  0.1428223
u(18)
##  0.125
u(20)
##  -1
u(30)
##  -1227133513

With the lazy numbers, this sequence never moves from $$1/7$$:

u_lazy <- function(n) {
if(n == 1) {
return(1 / lazynb(7))
}
8 * u_lazy(n-1) - 1
}
as.double(u_lazy(100))
##  0.1428571

Let’s compare with the multiple precision numbers provided by the Rmpfr package. One has to set the precision of these numbers. Let’s try with $$256$$ bits (the double precision corresponds to $$53$$ bits):

library(Rmpfr)
u_mpfr <- function(n) {
if(n == 1) {
return(1 / mpfr(7, prec = 256L))
}
8 * u_mpfr(n-1) - 1
}
asNumeric(u_mpfr(30))
##  0.1428571
asNumeric(u_mpfr(85))
##  0.140625
asNumeric(u_mpfr(100))
##  -78536544841

The sequence goes crazy before the $$100^{\textrm{th}}$$ term. Of course we could increase the precision. With the lazy numbers, there’s no precision to set. Moreover they are faster (at least for this example):

library(microbenchmark)
options(digits = 4L)
microbenchmark(
lazy = u_lazy(200),
mpfr = u_mpfr(200),
times = 20L
)
## Unit: milliseconds
##  expr   min    lq  mean median    uq   max neval cld
##  lazy 38.42 39.57 40.30  40.04 40.68 43.47    20  a
##  mpfr 59.22 60.09 61.01  60.60 61.94 64.89    20   b

Vectors of lazy numbers and matrices of lazy numbers are available in the lazyNumbers package. One can get the inverse and the determinant of a square lazy matrix.

A thing to note is that the usual mathematical functions such as $$\exp$$, $$\cos$$ or $$\sqrt{}$$, are not implemented for lazy numbers. Only the addition, the subtraction, the multiplication, the division and the absolute value are implemented.