---
title: "Arbitrary-precision arithmetic"
author: "Mikkel Meyer Andersen and Søren Højsgaard"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Arbitrary-precision arithmetic}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

```{r, message=FALSE}
library(Ryacas)
```

In `R`, the package `Rmpfr` package provide some functions for arbitrary-precision arithmetic. 
The way it is done is to allow for using more memory for storing numbers. 
This is very good for some use-cases, e.g. the example in `Rmpfr`'s a vignette "Accurately Computing $\log(1 - \exp(-|a|))$". 

Here, we demonstrate how to do investigate other aspects of arbitrary-precision arithmetic using `Ryacas`, e.g. solving systems of linear equations and matrix inverses.

# Arbitrary-precision arithmetic

First see the problem when using floating-point arithmetic:

```{r}
(0.3 - 0.2) - 0.1
yac_str("(0.3 - 0.2) - 0.1") # decimal number does not 
                             # always work well; often 
                             # it is better to represent as 
                             # rational number if possible

# (1/3 - 1/5) = 5/15 - 3/15 = 2/15
(1/3 - 1/5) - 2/15
yac_str("(1/3 - 1/5) - 2/15")
```


The `yacas` function [`N()`](https://yacas.readthedocs.io/en/latest/reference_manual/arithmetic.html#N) gives the numeric (floating-point) value for a precision.

```{r}
yac_str("1/3")
yac_str("N(1/3)")
yac_str("N(1/3, 1)")
yac_str("N(1/3, 200)")
```


# Evaluating polynomials

Consider the polynomial $(x-1)^7$ with root 1 (with multiplicity 7). We can consider it both in factorised form and in expanded form:
```{r}
pol <- "(x-1)^7"
p1 <- pol %>% yac_expr()
p1
p2 <- pol %>% y_fn("Expand") %>% yac_expr()
p2
```

Mathematically, these are the same objects, but when it comes to numerical
computations, the picture is different:

```{r}
eval(p1, list(x = 1.001))
eval(p2, list(x = 1.001))
```

We can of course also evaluate the polynomials in the different forms in `yacas` using [`WithValue()`](https://yacas.readthedocs.io/en/latest/reference_manual/controlflow.html#WithValue):

```{r}
# First try with 1.001:
pol_val <- paste0("WithValue(x, 1.001, ", pol, ")")
pol_val
yac_str(pol_val)
#... to get result symbolically, use instead the number as a fraction
pol_val <- paste0("WithValue(x, 1001/1000, ", pol, ")")
pol_val
yac_str(pol_val)
pol_val %>% y_fn("Denom") %>% yac_str()
pol_val %>% y_fn("Denom") %>% y_fn("IntLog", "10") %>% yac_str()
pol_val <- paste0("WithValue(x, 1001/1000, Expand(", pol, "))")
pol_val
yac_str(pol_val)
```

The reason for the difference is the near cancellation problem when using floating-point representation: Subtraction of nearly identical quantities. 
This phenomenon is
illustrated in the plot below.

```{r, fig.height=4, fig.width = 6, fig.pos="H", fig.cap="solid line: $(x-1)^7$; dashed line: factorization of $(x-1)^7$.", echo=FALSE}
xval <- seq(.99, 1.01, by = 0.001)
p1_val <- sapply(xval, function(x) eval(p1, list(x = x)))
p2_val <- sapply(xval, function(x) eval(p2, list(x = x)))
plot(xval, p1_val, type = "l", xlab = "x", ylab = "y")
lines(xval, p2_val, lty = 2)
legend("bottomright",
       legend = c(
         parse(text = paste0("p[1]: y == ", as.character(p1))),
         parse(text = paste0("p[2]: y == ", as.character(p2)))),
       lty = c(1, 2),
       cex = 0.7)
```

This example could of course have been illustrated directly in `R` but
it is convenient that `Ryacas` provides expansion of polynomials
directly so that students can experiment with other polynomials
themselves.


# Inverse of Hilbert matrix

In `R`'s `solve()` help file the [Hilbert matrix](https://en.wikipedia.org/wiki/Hilbert_matrix) is introduced:

```{r}
hilbert <- function(n) { 
  i <- 1:n
  H <- 1 / outer(i - 1, i, "+")
  return(H)
}
```

We will now demonstrate the known fact that it is ill-conditioned (meaning e.g. that it is difficult to determine its inverse numerically).

First we make the Hilbert matrix symbolically:
```{r}
hilbert_sym <- function(n) { 
  mat <- matrix("", nrow = n, ncol = n)
  
  for (i in 1:n) {
    for (j in 1:n) {
      mat[i, j] <- paste0("1 / (", (i-1), " + ", j, ")")
    }
  }
  
  return(mat)
}
```

```{r}
A <- hilbert(8)
A
B1 <- hilbert_sym(8)
B1
B <- ysym(B1) # convert to yacas symbol
B
```

```{r}
Ainv <- solve(A)
Ainv
Binv1 <- solve(B) # result is still yacas symbol
Binv1
Binv <- as_r(Binv1) # convert to R numeric matrix
Binv
Ainv - Binv
max(abs(Ainv - Binv))
max(abs((Ainv - Binv) / Binv))
```

As seen, already for $n=8$ is the numeric solution inaccurate.