[R] Interesting quirk with fractions and rounding / using == for floating point

[J C Nash]

Thanks Richard.

I've some stuff too, but I need to look it up. A few years ago I built
a small test spreadsheet for Gnumeric when working with Jody Goldberg.
In the early 2000s, Jody contacted R (I think Duncan Murdoch) to ask if
it was OK for Gnumeric to use R's distribution function approximations (Yes).
Someone can correct me if wrong, but apparently a few weeks later Jody
sent some improvements, and I believe these were incorporated. Good
open-source improvements like that deserve to be recorded in a better way
than I have done here. But in any event, I can dig out that material, as
well as some of Kahan's PARANOIA materials.

Suggest taking the discussion off-line from here on unless we come up
with some important bugs or improvements or ...

Best, JN


On 2017-04-23 11:42 AM, Richard M. Heiberger wrote:
> John,
>
> I would be happy to participate in designing the test suite you suggest.
>
> About a year ago I revised FAQ 7.31, based on my talk at the Aalberg R
> conference.  It now points, in addition to the Goldberg paper that has
> been referenced there for a long time, to my appendix on precision.
> Here is the link from the FAQ
>
>       A discussion with many easily followed examples is in Appendix G
>    "Computational Precision and Floating Point Arithmetic", pages 753-771
>    of _Statistical Analysis and Data Display: An Intermediate Course with
>    Examples in R_, Richard M. Heiberger and Burt Holland (Springer 2015,
>    second edition).  This appendix is a free download from
>    <http://link.springer.com/content/pdf/bbm%3A978-1-4939-2122-5%2F1.pdf>.
>
> The appendix is based on the paper I gave at the Aalberg R
> Conference. It uses the Rmpfr package to illustrate exactly what
> is happening at the level of the machine arithmetic.  The
> investigation and discussion of the effect of optimization on
> floating point arithmetic should use the Rmpfr tools.
>
>
> Rich
>
> On Sun, Apr 23, 2017 at 10:06 AM, J C Nash <profjcnash at gmail.com> wrote:
>> Yes. I should have mentioned "optimizing" compilers, and I can agree with
>> "never
>> trusting exact equality", though I consider conscious use of equality tests
>> useful.
>> Optimizing compilers have bitten me once or twice. Unfortunately, a lot of
>> floating-point work requires attention to detail. In the situation of
>> testing
>> for convergence, the alternatives to the test I propose require quite a lot
>> of
>> code, with variants for different levels of precision e.g., single, double,
>> quad.
>>
>> There's no universal answer and we do have to look "under the hood". A
>> particularly
>> nasty instance (now fortunately long gone) was the Tektronix 4051 graphics
>> station,
>> where the comparisons automatically included a "fuzz". There was a FUZZ
>> command to
>> set this. Sometimes the "good old days" weren't! Today's equivalent is when
>> there
>> is an upstream change to an "optimization" that changes the manner of
>> computation,
>> as in Peter's examples. If we specify x <- y * (a / b), then we should not
>> get
>> x <- (y * a) / b.
>>
>> A slightly related case concerns eigenvectors / singular vectors when there
>> are
>> degenerate values (i.e., two or more equal eigenvalues). The vectors are
>> then
>> determined only to lie in a (hyper)plane. A large computer contracting firm
>> spent
>> two weeks of high-priced but non-numerical help trying to find the "error"
>> in either
>> an IBM or Univac program because they gave very different eigenvectors.
>>
>> And in my own field of function minimization / nonlinear least squares, it
>> is quite
>> common to have multiple minima or a plateau.
>>
>> Does anyone know if R has a test suite to check some of these situations? If
>> not,
>> I'll be happy to participate in generating some. They would not need to be
>> run
>> very often, and could be useful as a didactic tool as well as checking for
>> differences in platforms.
>>
>> JN
>>
>> On 2017-04-23 09:37 AM, peter dalgaard wrote:
>>>
>>>
>>>> On 23 Apr 2017, at 14:49 , J C Nash <profjcnash at gmail.com> wrote:
>>>>
>>>>
>>>> So equality in floating point is not always "wrong", though it should be
>>>> used
>>>> with some attention to what is going on.
>>>>
>>>> Apologies to those (e.g., Peter D.) who have heard this all before. I
>>>> suspect
>>>> there are many to whom it is new.
>>>
>>>
>>> Peter D. still insists on never trusting exact equality, though. There was
>>> at least one case in the R sources where age-old code got itself into a
>>> condition where a residual terme that provably should decrease on every
>>> iteration oscillated between two values of 1-2 ulp in magnitude without ever
>>> reaching 0. The main thing is that you cannot trust optimising compilers
>>> these days. There is, e.g.,  no guarantee that a compiler will not transform
>>>
>>> (x_new + offset) == (x_old + offset)
>>>
>>> to
>>>
>>> (x_new + offset) - (x_old + offset) == 0
>>>
>>> to
>>>
>>> (x_new - x_old) + (offset - offset) == 0
>>>
>>> to.... well, you get the point.
>>>
>>> -pd
>>>
>>
>> ______________________________________________
>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.