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 <profjcn...@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 <profjcn...@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@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.
______________________________________________
R-help@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.
