---
title: "The lazy numbers in R"
author: "Stéphane Laurent"
date: '2022-11-12'
tags: R, C++, Rcpp, maths
rbloggers: yes
output:
  md_document:
    variant: markdown
    preserve_yaml: true
  html_document:
    highlight: kate
    keep_md: no
highlighter: pandoc-solarized
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE)
```

My new package [**lazyNumbers**](https://github.com/stla/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$:

```{r}
x <- 1 - 7 * 0.1
x == 0.3
```

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

```{r, message=FALSE}
library(lazyNumbers)
x <- lazynb(1) - lazynb(7) * lazynb(0.1)
as.double(x) == 0.3
```

***Correction:*** *this is* [*wrong*](https://laustep.github.io/stlahblog/posts/lazyNumbers2.html)*!*

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:

```{r}
x <- 1 - lazynb(7) * 0.1
as.double(x) == 0.3
```

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$):

```{r}
u <- function(n) {
  if(n == 1) {
    return(1 / 7)
  }
  8 * u(n-1) - 1
}
u(15)
u(18)
u(20)
u(30)
```

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

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

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):

```{r, message=FALSE}
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))
asNumeric(u_mpfr(85))
asNumeric(u_mpfr(100))
```

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):

```{r}
library(microbenchmark)
options(digits = 4L)
microbenchmark(
  lazy = u_lazy(200),
  mpfr = u_mpfr(200),
  times = 20L
)
```

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.