I'm happy to take this to r-devel if you think that is the best way of
resolving the question, but I cannot agree with your assessment of the
issue. This is not a philosophical complaint, this is a genuine bug
in R. I am, deliberately and with full knowledge of the issues,
working right down at the limit of the precision of an IEEE double. I
need to be able to see the contribution to the value made by every
last bit of the mantissa. Also, you may not have realized that this
affects saving and reloading of data - because of this bug it is not
possible to save numeric data in textual format and later reload
exactly what you had in memory, which can have serious consequences
for some algorithms.
I shall try to be more concrete about what I am doing. I am working
with statistical methods that are sensitive to ties in the data. In
order to eliminate ties, I am applying "jitter" to the problem data
points; to minimize the effect on all calculations, the jitter is
supposed to be just one ulp either way. To do this, I've written
myself an itty bitty extension that invokes the C math library
function nextafter(). This works, except that (as I said in the
original bug report) R won't print enough decimal places.
#include <math.h>
void R_nextafter(double *val, double *dir)
{ *val = nextafter(*val, *dir); }
dyn.load('nextafter.so')
nextafter <- function(val, dir)
+ .C("R_nextafter", val=as.double(val), dir=as.double(dir))$val
format(nextafter(1,2), digits=16) # doesn't work, grrr
[1] "1"
format(nextafter(1,2)-1, digits=16) # does show the difference
[1] "2.220446049250313e-16"
Contrast what you get from C (i.e. what I consider to be correct behavior):
$ cat test.c
#include <stdio.h>
#include <math.h>
int main(void)
{
double up1 = nextafter(1., 2.);
printf("nextafter(1,2) = %.17g [%a]\n", up1, up1);
printf("nextafter(1,2)-1 = %.17g [%a]\n", up1-1, up1-1);
return 0;
}
$ gcc test.c -lm && ./a.out
nextafter(1,2) = 1.0000000000000002 [0x1.0000000000001p+0]
nextafter(1,2)-1 = 2.2204460492503131e-16 [0x1p-52]
It is obviously much harder to eyeball one's data set for correct
handling of ties if one cannot simply apply format() with an
appropriate number of digits to the entire vector. And if I dump out
the data set to a text file after doing this, and then reload it, the
corrections will be lost.
zw
--
To UNSUBSCRIBE, email to [EMAIL PROTECTED]
with a subject of "unsubscribe". Trouble? Contact [EMAIL PROTECTED]
