Arbitrary-precision arithmetic

Mikkel Meyer Andersen and Søren Højsgaard

2023-01-16

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:

(0.3 - 0.2) - 0.1
## [1] -2.775558e-17
yac_str("(0.3 - 0.2) - 0.1") # decimal number does not 
## [1] "0"
                             # 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
## [1] -2.775558e-17
yac_str("(1/3 - 1/5) - 2/15")
## [1] "0"

The yacas function N() gives the numeric (floating-point) value for a precision.

yac_str("1/3")
## [1] "1/3"
yac_str("N(1/3)")
## [1] "0.3333333333"
yac_str("N(1/3, 1)")
## [1] "0.3"
yac_str("N(1/3, 200)")
## [1] "0.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333"

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:

pol <- "(x-1)^7"
p1 <- pol %>% yac_expr()
p1
## expression((x - 1)^7)
p2 <- pol %>% y_fn("Expand") %>% yac_expr()
p2
## expression(x^7 - 7 * x^6 + 21 * x^5 - 35 * x^4 + 35 * x^3 - 21 * 
##     x^2 + 7 * x - 1)

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

eval(p1, list(x = 1.001))
## [1] 1e-21
eval(p2, list(x = 1.001))
## [1] 8.881784e-15

We can of course also evaluate the polynomials in the different forms in yacas using WithValue():

# First try with 1.001:
pol_val <- paste0("WithValue(x, 1.001, ", pol, ")")
pol_val
## [1] "WithValue(x, 1.001, (x-1)^7)"
yac_str(pol_val)
## [1] "0"
#... to get result symbolically, use instead the number as a fraction
pol_val <- paste0("WithValue(x, 1001/1000, ", pol, ")")
pol_val
## [1] "WithValue(x, 1001/1000, (x-1)^7)"
yac_str(pol_val)
## [1] "1/1000000000000000000000"
pol_val %>% y_fn("Denom") %>% yac_str()
## [1] "1000000000000000000000"
pol_val %>% y_fn("Denom") %>% y_fn("IntLog", "10") %>% yac_str()
## [1] "21"
pol_val <- paste0("WithValue(x, 1001/1000, Expand(", pol, "))")
pol_val
## [1] "WithValue(x, 1001/1000, Expand((x-1)^7))"
yac_str(pol_val)
## [1] "1/1000000000000000000000"

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.

solid line: \((x-1)^7\); dashed line: factorization of \((x-1)^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 is introduced:

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:

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)
}
A <- hilbert(8)
A
##           [,1]      [,2]      [,3]       [,4]       [,5]       [,6]       [,7]
## [1,] 1.0000000 0.5000000 0.3333333 0.25000000 0.20000000 0.16666667 0.14285714
## [2,] 0.5000000 0.3333333 0.2500000 0.20000000 0.16666667 0.14285714 0.12500000
## [3,] 0.3333333 0.2500000 0.2000000 0.16666667 0.14285714 0.12500000 0.11111111
## [4,] 0.2500000 0.2000000 0.1666667 0.14285714 0.12500000 0.11111111 0.10000000
## [5,] 0.2000000 0.1666667 0.1428571 0.12500000 0.11111111 0.10000000 0.09090909
## [6,] 0.1666667 0.1428571 0.1250000 0.11111111 0.10000000 0.09090909 0.08333333
## [7,] 0.1428571 0.1250000 0.1111111 0.10000000 0.09090909 0.08333333 0.07692308
## [8,] 0.1250000 0.1111111 0.1000000 0.09090909 0.08333333 0.07692308 0.07142857
##            [,8]
## [1,] 0.12500000
## [2,] 0.11111111
## [3,] 0.10000000
## [4,] 0.09090909
## [5,] 0.08333333
## [6,] 0.07692308
## [7,] 0.07142857
## [8,] 0.06666667
B1 <- hilbert_sym(8)
B1
##      [,1]          [,2]          [,3]          [,4]          [,5]         
## [1,] "1 / (0 + 1)" "1 / (0 + 2)" "1 / (0 + 3)" "1 / (0 + 4)" "1 / (0 + 5)"
## [2,] "1 / (1 + 1)" "1 / (1 + 2)" "1 / (1 + 3)" "1 / (1 + 4)" "1 / (1 + 5)"
## [3,] "1 / (2 + 1)" "1 / (2 + 2)" "1 / (2 + 3)" "1 / (2 + 4)" "1 / (2 + 5)"
## [4,] "1 / (3 + 1)" "1 / (3 + 2)" "1 / (3 + 3)" "1 / (3 + 4)" "1 / (3 + 5)"
## [5,] "1 / (4 + 1)" "1 / (4 + 2)" "1 / (4 + 3)" "1 / (4 + 4)" "1 / (4 + 5)"
## [6,] "1 / (5 + 1)" "1 / (5 + 2)" "1 / (5 + 3)" "1 / (5 + 4)" "1 / (5 + 5)"
## [7,] "1 / (6 + 1)" "1 / (6 + 2)" "1 / (6 + 3)" "1 / (6 + 4)" "1 / (6 + 5)"
## [8,] "1 / (7 + 1)" "1 / (7 + 2)" "1 / (7 + 3)" "1 / (7 + 4)" "1 / (7 + 5)"
##      [,6]          [,7]          [,8]         
## [1,] "1 / (0 + 6)" "1 / (0 + 7)" "1 / (0 + 8)"
## [2,] "1 / (1 + 6)" "1 / (1 + 7)" "1 / (1 + 8)"
## [3,] "1 / (2 + 6)" "1 / (2 + 7)" "1 / (2 + 8)"
## [4,] "1 / (3 + 6)" "1 / (3 + 7)" "1 / (3 + 8)"
## [5,] "1 / (4 + 6)" "1 / (4 + 7)" "1 / (4 + 8)"
## [6,] "1 / (5 + 6)" "1 / (5 + 7)" "1 / (5 + 8)"
## [7,] "1 / (6 + 6)" "1 / (6 + 7)" "1 / (6 + 8)"
## [8,] "1 / (7 + 6)" "1 / (7 + 7)" "1 / (7 + 8)"
B <- ysym(B1) # convert to yacas symbol
B
## {{   1,  1/2,  1/3,  1/4,  1/5,  1/6,  1/7,  1/8},
##  { 1/2,  1/3,  1/4,  1/5,  1/6,  1/7,  1/8,  1/9},
##  { 1/3,  1/4,  1/5,  1/6,  1/7,  1/8,  1/9, 1/10},
##  { 1/4,  1/5,  1/6,  1/7,  1/8,  1/9, 1/10, 1/11},
##  { 1/5,  1/6,  1/7,  1/8,  1/9, 1/10, 1/11, 1/12},
##  { 1/6,  1/7,  1/8,  1/9, 1/10, 1/11, 1/12, 1/13},
##  { 1/7,  1/8,  1/9, 1/10, 1/11, 1/12, 1/13, 1/14},
##  { 1/8,  1/9, 1/10, 1/11, 1/12, 1/13, 1/14, 1/15}}
Ainv <- solve(A)
Ainv
##         [,1]      [,2]       [,3]       [,4]        [,5]        [,6]
## [1,]      64     -2016      20160     -92400      221760     -288288
## [2,]   -2016     84672    -952560    4656960   -11642400    15567552
## [3,]   20160   -952560   11430720  -58212000   149688000  -204324119
## [4,]  -92400   4656960  -58212000  304919999  -800414996  1109908794
## [5,]  221760 -11642400  149688000 -800414996  2134439987 -2996753738
## [6,] -288288  15567552 -204324119 1109908793 -2996753738  4249941661
## [7,]  192192 -10594584  141261119 -776936154  2118916782 -3030050996
## [8,]  -51480   2882880  -38918880  216215998  -594593995   856215351
##             [,7]       [,8]
## [1,]      192192     -51480
## [2,]   -10594584    2882880
## [3,]   141261119  -38918880
## [4,]  -776936155  216215998
## [5,]  2118916783 -594593995
## [6,] -3030050996  856215352
## [7,]  2175421226 -618377753
## [8,]  -618377753  176679358
Binv1 <- solve(B) # result is still yacas symbol
Binv1
## {{         64,       -2016,       20160,      -92400,      221760,     -288288,      192192,      -51480},
##  {      -2016,       84672,     -952560,     4656960,   -11642400,    15567552,   -10594584,     2882880},
##  {      20160,     -952560,    11430720,   -58212000,   149688000,  -204324120,   141261120,   -38918880},
##  {     -92400,     4656960,   -58212000,   304920000,  -800415000,  1109908800,  -776936160,   216216000},
##  {     221760,   -11642400,   149688000,  -800415000,  2134440000, -2996753760,  2118916800,  -594594000},
##  {    -288288,    15567552,  -204324120,  1109908800, -2996753760,  4249941696, -3030051024,   856215360},
##  {     192192,   -10594584,   141261120,  -776936160,  2118916800, -3030051024,  2175421248,  -618377760},
##  {     -51480,     2882880,   -38918880,   216216000,  -594594000,   856215360,  -618377760,   176679360}}
Binv <- as_r(Binv1) # convert to R numeric matrix
Binv
##         [,1]      [,2]       [,3]       [,4]        [,5]        [,6]
## [1,]      64     -2016      20160     -92400      221760     -288288
## [2,]   -2016     84672    -952560    4656960   -11642400    15567552
## [3,]   20160   -952560   11430720  -58212000   149688000  -204324120
## [4,]  -92400   4656960  -58212000  304920000  -800415000  1109908800
## [5,]  221760 -11642400  149688000 -800415000  2134440000 -2996753760
## [6,] -288288  15567552 -204324120 1109908800 -2996753760  4249941696
## [7,]  192192 -10594584  141261120 -776936160  2118916800 -3030051024
## [8,]  -51480   2882880  -38918880  216216000  -594594000   856215360
##             [,7]       [,8]
## [1,]      192192     -51480
## [2,]   -10594584    2882880
## [3,]   141261120  -38918880
## [4,]  -776936160  216216000
## [5,]  2118916800 -594594000
## [6,] -3030051024  856215360
## [7,]  2175421248 -618377760
## [8,]  -618377760  176679360
Ainv - Binv
##               [,1]          [,2]          [,3]          [,4]          [,5]
## [1,]  2.022656e-07 -7.897346e-06  0.0000742203 -0.0002794876  4.892563e-04
## [2,] -7.516163e-06  2.443442e-04 -0.0016348482  0.0022857217  8.007960e-03
## [3,]  6.596705e-05 -1.459678e-03 -0.0013010856  0.0928798467 -4.353278e-01
## [4,] -2.222202e-04  6.508809e-04  0.1018696651 -1.0164485574  3.740381e+00
## [5,]  3.084455e-04  1.383389e-02 -0.4781522453  3.8412284851 -1.304656e+01
## [6,] -1.008651e-04 -3.762038e-02  0.9126323462 -6.7141628265  2.190641e+01
## [7,] -1.274265e-04  3.748520e-02 -0.7896100283  5.5373635292 -1.763533e+01
## [8,]  8.357911e-05 -1.313744e-02  0.2562877908 -1.7438534796  5.464779e+00
##               [,6]          [,7]          [,8]
## [1,] -3.859913e-04  9.308662e-05  1.726509e-05
## [2,] -2.777486e-02  2.954033e-02 -1.067452e-02
## [3,]  8.310004e-01 -7.192242e-01  2.335241e-01
## [4,] -6.457403e+00  5.287978e+00 -1.657728e+00
## [5,]  2.159914e+01 -1.723701e+01  5.309898e+00
## [6,] -3.544080e+01  2.786332e+01 -8.493719e+00
## [7,]  2.811761e+01 -2.189091e+01  6.626538e+00
## [8,] -8.625638e+00  6.669370e+00 -2.008779e+00
max(abs(Ainv - Binv))
## [1] 35.4408
max(abs((Ainv - Binv) / Binv))
## [1] 1.136963e-08

As seen, already for \(n=8\) is the numeric solution inaccurate.