as.integer(8952) = 8951?

user_2298323 0
Viewed 98

I accidentally discovered a weird bug with either the as.integer or the det function in R base. Does anyone know what is going on here and how to prevent it?

I was computing the determinant of the following 3-by-3 matrix:

mat <- matrix(c(15, 6, 116, 10, 13, 16, 14, 23, 56), ncol = 3)

Which looks like this:

     [,1] [,2] [,3]
[1,]   15   10   14
[2,]    6   13   23
[3,]  116   16   56

Two things are easy to see: all entries are integers, and each of the six sets of three-entries-no-two-of-which-lie-in-the-same-row-or-column contains at least one even integer. Hence the determinant must be an even integer.

Asking R the actual value of this determinant by typing det(mat) it returns something that looks like an even integer: 8952. But lo and behold: deep down inside R's mind it actually is either a non-integer or an odd integer, since when typing as.integer(det(mat)) we get 8951.

What is going on here? The 8951 is obviously wrong. Also, less obviously, the value 8952 is correct as can be seen with pen and paper.

So my questions are:

  1. what is going on here?

  2. How do I force R to give me the correct integer values when asked to compute the determinant of an integer matrix?

determinants r integer
1 Answers

Underlying causes: Truncation by is.integer rather than rounding and floating point math on logged intermediate values to explain the second result, combined with the default digit level of display in the print portion of the console REPL to explain the initial result of det(mat):

print( det(mat), digits =16)
[1] 8951.999999999993

It's quite possible that the theoretic answer is 8952 but R is not a symbolic math engine.

enter image description here

You can use the Rmpfr package (as suggested by @BenBolker) to increase the level of precision:

 library(Rmpfr)
 mat <- mpfr(mat, 64)
 as.integer( det(mat) )
[1] 8952

as.integer truncates rather than rounds. See ?as.integer. R would be able to handle addition or multiplication of integers without loss of precision but as soon as division occurs there is likely to be floating point error. (Actually the problem arises because the default for det is to use determinant with log=TRUE and then exponentiate the modulus of the complex result.) From the Values section of the help page:

Non-integral numeric values are truncated towards zero (i.e., as.integer(x) equals trunc(x) there), and imaginary parts of complex numbers are discarded (with a warning).

user_18556770
Related