5

I have R version 4.1.2 (2021-11-01).

It seems trunc() function is not consistent when the input number has a large number of decimal values.

trunc(3.99999999999999977799999999999999999999900)
[1] 4

trunc(3.999999999999999777999999999999999999999000)
[1] 3

or

trunc(3.9999999999999997778888888888880)
[1] 4
trunc(3.999999999999999777888888888888)
[1] 3

I'm not sure what caused this inconsistency?

Bằng Rikimaru
  • 1,512
  • 2
  • 24
  • 50
  • 3
    Floating point numbers in R (and most other languages) simply do not support that level of precision. Not every decimal number you can type can be accurately represented in binary form. There are errors when storing the value in a form the computer can understand. Though I cannot replicate these results using R 4.2.1 – MrFlick Apr 04 '23 at 14:34
  • 3
    Voting to reopen because this post asks a question I do not see answered in the purported [original](https://stackoverflow.com/questions/9508518/why-are-these-numbers-not-equal). In each example here, the two numerals differ by a single trailing 0 and represent the same number. With any good software, they ought to be converted to the same floating-point number. This is more a question of how R parses input numerals than about how floating-point works, and it deserves an answer addressing that. – Eric Postpischil Apr 04 '23 at 17:15
  • 2
    I can reproduce with R 4.2.3 / macOS Ventura / Intel. Using `options(digits=22)` to display the values with greater precision confirms that it's the numeric-literal-to-float conversion that differs (so the issue is not specific to `trunc`). – Mark Dickinson Apr 04 '23 at 18:30

3 Answers3

4

There are two questions here: what's the right answer, and why is R getting a different answer under different circumstances?

That number 3.999999999999999777999999999999999999999 is obviously very close to 4. In fact, it's closer to 4 than it is to any other IEEE-754 double-precision floating-point number. The next-lower representable number is about 3.9999999999999995, which is a little bit farther away. So, strictly speaking, trunc(3.999999999999999777999999999999999999999) should be trunc(4.0) which is obviously 4. That is, when R takes the input 3.999999999999999777999999999999999999999, it should immediately convert it to an internal value of 4, even before trying to truncate it. This looks "wrong", because you and I can plainly see that truncating 3.999… should give 3, but the fact that not every real number is representable in a finite-precision floating-point representation does occasionally lead to anomalies like this. (See also these three questions, which gather SO's canonical answers for these sorts of binary floating-point anomalies.)

For the rest of this answer, we leave the realm of using R and delve into the world of implementing R. (I'm a C programmer, not an R user, and this answer is likely to betray that bias, by being ignorant of any nuances of R. Apologies for that.) But, at any rate, R is written in C, and uses C's double type for much of its arithmetic. And on the vast majority of popular general-purpose computers today, a C double is implemented using IEEE-754 double precision, which is why I led off this answer by mentioning that standard.

But why is R getting a different answer depending on how many trailing 0's there are? The answer lies in a function deep down in the R interpreter which is doing the actual conversion of characters typed by the user into internal R data structures.

How might we convert a string like "123.456" into its internal floating-point representation? One way to do it is to temporarily ignore the decimal point and convert it to an integer, resulting in the number 123456, then count the number of digits after the decimal point, and divide by ten to that power. And, indeed, 123456 ÷ 10³ is 123.456.

But using that strategy, converting 3.99999999999999977799999999999999999999900 is going to involve taking a 42-digit number and dividing it by 1041, while converting 3.999999999999999777999999999999999999999000 is going to involve taking a 43-digit number and dividing it by 1042.

And none of these numbers are going to be exactly representable in binary floating point. They're going to be a little bit off, and that's sometimes going to lead to discrepancies. And in particular, when the numbers are this big, there's no guarantee that a ÷ b will give you exactly the same answer as 10a ÷ 10b.

For the current example, the discrepancy is that one division leads to a number that's closer to 4, and one division leads to a number that's closer to 3.9999999999999995. (And, remember, I'm talking about divisions that are happening in C code deep in the R interpreter, not any divisions you thought you were doing in R.)

There are several additional factors involved here. (In particular, R is using "binary exponentiation" to compute 10N, and that ends up making a difference, too.) I don't have time to write those details up just now; maybe later. Interested readers can consult the file src/main/util.c in the R source distribution, specifically the function R_strtod5.

But the take-home lesson is that accurately converting back and forth between binary floating-point numbers and human-readable base-10 representations is hard. Among other things, getting properly-rounded results usually requires doing your calculations in some higher-precision representation, so that you'll have something just barely accurate enough to round off at the end (that is, to yield the desired "properly-rounded result"). Ironically, R's implementation is trying to do the right thing in that respect, computing both numbers (that is, the two numbers to be divided) using C's long double type. I would have thought that would have been sufficient to avoid anomalies like this, but evidently not.

Also this would be worth reporting as a bug in R. A truly high-quality strtod implementation won't have anomalies like this, and having gone the route of implementing its own, R is (I would say) on the hook to reinvent whatever wheels are necessary in order to get a properly-rounded result in all cases.

Steve Summit
  • 45,437
  • 7
  • 70
  • 103
  • Interesting, so for now the most reliable way to trunc high-precision numbers would be to convert it to a char and then chop everything off at the decimal point? – DuckPyjamas Apr 05 '23 at 17:03
  • Hmmm. Python3 does better (converts both to 4), haven't tried the base C-library `strtod` yet ... – Ben Bolker Apr 05 '23 at 17:10
  • @DuckPyjamas: yes. I think very few languages can handle [binary-coded decimal representation](https://en.wikipedia.org/wiki/Binary-coded_decimal) ... – Ben Bolker Apr 05 '23 at 17:13
  • 1
    @BenBolker You don't need BCD to do decimal floating-point. (IBM definitely favors it, it's true.) [IEEE-754-2008](https://en.wikipedia.org/wiki/IEEE_754-2008_revision) actually specifies some [true-decimal variants](https://en.wikipedia.org/wiki/Decimal64_floating-point_format), where the exponent is of a power of 10, and where the significand can be implemented, at your option, in either BCD ("DPD"), or straight binary ("BID"). – Steve Summit Apr 05 '23 at 17:18
  • 2
    @DuckPyjamas It probably depends on where the numbers are coming from. If you really do have numbers like 3.99999999999999977799999999999999999999900 to work with, you may have to go through contortions. But I'm guessing that number arose for Bằng Rikimaru while doing some (unrealistic) testing. More normal numbers would probably work fine. Or, an explicit rounding step might be appropriate. But, really, it depends on, like I said, where the numbers are coming from, and what the purpose (expected result) of the truncation step is. – Steve Summit Apr 05 '23 at 20:56
  • @BenBolker in this case they should not round to 4, since trunc() is supposed to floor. – DuckPyjamas Apr 05 '23 at 21:25
  • 1
    @DuckPyjamas I know what you mean, but the question comes back to, where did that number 3.99999999999999977799999999999999999999900 come from? It may be that it was actually *supposed* to be 4, that it only departs from 4 due to errors that crept in earlier, in which case a rounding step might be perfectly appropriate. (Or not.) – Steve Summit Apr 05 '23 at 21:28
  • @DuckPyjamas: what I meant to say is that Python converts *consistently* regardless of the numbers of trailing zeros. – Ben Bolker Apr 05 '23 at 21:38
2

To complement @SteveSummit's excellent answer, let's store both of these numbers in its own variable and see what they look before we apply trunc(), printing them to the maximum available precision:

x1 <- 3.99999999999999977799999999999999999999900
x2 <- 3.999999999999999777999999999999999999999000
print(x1, digits = 22)
## [1] 4
print(x2, digits = 22)
## [1] 3.999999999999999555911

If you want to see the exact code that @SteveSummit is referring to (find the n-digit number by successively adding 10* the previous value to the next digit, then divide by 10 the appropriate number of times), it's here ...

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
1

I'm not sure what caused this inconsistency?

Implementation weakness.

What's special about 3.99999999999999977799999999999999999999900(0)?

With common encoding of floating point numbers, 4.0 and its subsequent and preceding values are exactly:

// As decimal                                           As hexadecimal
4.00000000000000088817841970012523233890533447265625    0x4.0000000000002
4.0                                                     0x4.0 
3.999999999999999555910790149937383830547332763671875   0x3.ffffffffffffe

The value half-way in between values, along with OP's constants, using wider math, are:

// As decimal                                           As hexadecimal
4.000000000000000444089209850062616169452667236328125   0x4.0000000000001
3.9999999999999997779553950749686919152736663818359375  0x3.fffffffffffff
3.99999999999999977799999999999999999999900             Let us call this C1
3.999999999999999777999999999999999999999000            Let us call this C2  
1 23456789 123456789 123456789 123456789 123            Significant digit count

It is evident, OP's two constants (C1 42 digits and C2 43 digits) are selected to test the text-to-floating-point value conversion of OP's R.

In a perfect text-to-floating-point value conversion, both text C1, C2 would convert to the closer 4.0. As C1 instead converts to the smaller 3.9999999999999995559... simply reflects a weakness in the quality of the R implementation.

Any text "3.999999999999999555910790149937383830547332763671875" or more (and ≤ "4.000000000000000444089209850062616169452667236328125") should become 4.0.

chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256