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)⁷ 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)⁷; dashed line: factorization of (x − 1)⁷.

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   -11642401    15567553
## [3,]   20160   -952560   11430720  -58212002   149688006  -204324129
## [4,]  -92400   4656960  -58212002  304920013  -800415033  1109908845
## [5,]  221760 -11642401  149688006 -800415033  2134440086 -2996753877
## [6,] -288288  15567553 -204324129 1109908846 -2996753878  4249941857
## [7,]  192192 -10594585  141261126 -776936192  2118916882 -3030051135
## [8,]  -51480   2882880  -38918882  216216009  -594594023   856215391
##             [,7]       [,8]
## [1,]      192192     -51480
## [2,]   -10594584    2882880
## [3,]   141261126  -38918882
## [4,]  -776936191  216216009
## [5,]  2118916881 -594594022
## [6,] -3030051135  856215390
## [7,]  2175421325 -618377781
## [8,]  -618377781  176679366

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,]  0.0000012595 -6.280053e-05  0.0007762981  -0.004023341    0.01045849
## [2,] -0.0000634834  3.150850e-03 -0.0388042819   0.200490239   -0.51978479
## [3,]  0.0007913917 -3.913434e-02  0.4805039167  -2.476326779    6.40613982
## [4,] -0.0041300914  2.036124e-01 -2.4937464967  12.824575186  -33.11633968
## [5,]  0.0107986356 -5.310077e-01  6.4897316694 -33.314724684   85.89412045
## [6,] -0.0149156042  7.318508e-01 -8.9280686677  45.760925531 -117.82621908
## [7,]  0.0104020056 -5.094212e-01  6.2047709525 -31.759819388   81.68062663
## [8,] -0.0028846884  1.410395e-01 -1.7154899687   8.770542115  -22.53314412
##               [,6]          [,7]          [,8]
## [1,]   -0.01437405    0.00998076  -0.002757407
## [2,]    0.71275051   -0.49393205   0.136219052
## [3,]   -8.76787683    6.06612584  -1.670521952
## [4,]   45.25378895  -31.26587522   8.599618167
## [5,] -117.21701336   80.88934231 -22.225065231
## [6,]  160.60632944 -110.71761990  30.392856956
## [7,] -111.22357941   76.60545301 -21.011935592
## [8,]   30.65552807  -21.09724677   5.782607675

max(abs(Ainv - Binv))

## [1] 160.6063

max(abs((Ainv - Binv) / Binv))

## [1] 5.603513e-08

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