>>>>> "Stavros" == Stavros Macrakis <macrakis at alum.mit.edu> writes:
Stavros> Hi, Ray,
Stavros> I saw your comment on ration1 (in rat3e.lisp) wondering about the accuracy
Stavros> of the float-to-rational conversion there.
I tried to dig up some of the bugs. The only one I could find was
that CMUCL's old algorithm failed for 5.964664720239993d-15 (division
by zero), but maxima doesn't seem to have that problem.
Stavros> Out of curiosity, I compared its results to results from
Stavros> cf(rationalize(...)), which are the exact answers. A little experimentation
Note that maxima's rationalize doesn't call ration1. It calls
Barton's $rationalize in maxmin.lisp. (I tried tracing ration1 and
found that it was never called from rationalize.) But
maxima-rationalize (which does call ration1) is still called from
several other places.
Stavros> with random arguments shows that it is +/- 1 ULP off about 5% of the time,
Stavros> and +/- 2 ULP off about 1% of the time. I didn't find examples of larger
Stavros> errors, and I'm not sure how I'd construct cases which would cause them
Stavros> (ideas?). Of course, this does not give us any real upper bound. But I'm
It seems we shouldn't ever be more than 1 ulp off.
Stavros> pretty confident that there is no problem with the default ratepsilon of
Stavros> 2e-8 (which was based on the PDP-10's 27-bit mantissa!).
Stavros> That said, I thought we'd agreed in February 2005 that the default
Stavros> ratepsilon should be much smaller, say 2.0e-15 (for IEEE doubles, which I
Stavros> believe all Maxima implementations use).
Stavros> And I think it would be logical that for ratepsilon = 0,
Stavros> rat(x)=rationalize(x).
Could we not just set ratepsilon to 0 instead of 2e-15 or so?
But given that rationalize doesn't seem to call ration1 it's no longer
clear to me what we should do. Should maxima-rationalize be changed?
Should rationalize actually call maxima-rationalize again? What about
ratepsilon if rationalize doesn't use ratepsilon? What does it mean,
then, when the docs say ratepsilon is used in the conversion of
floating-point numbers to rationals?
Ray