Subject: Bfloat conversion of rationals can be inaccurate
From: Raymond Toy
Date: Sat, 12 Dec 2009 22:42:44 -0500
Raymond Toy wrote:
> Maxima currently converts rationals by bfloating the numerator and
> denominator and then performing a bfloat division. This is probably ok,
> but consider what happens if the numerator and/or denominator has more
> than fpprec digits. We then get two rounding operations, and then one
> more for the division.
>
> We can do better than that. We should left shift the numerator by
> enough bits so that the numerator divided by the denominator produces
> fpprec digits. The remainder tells us how to round the result. This
> gives us the nearest bfloat to the rational, with just the one rounding
> needed.
>
> I noticed that some of the differences between the current and fast
> bfloat input routines appear to be actually incorrect values for the
> current algorithm. I think this is caused by the inaccurate conversion
> of rationals.
>
Here is one example of an inaccurate conversion of a rational to a bfloat.
(%i151) fpprec:5;
(%o151) 5
(%i152) (2^20+1)/(2^20-1);
(%o152) 1048577/1048575
(%i153) bfloat(%);
(%o153) 1.0b0
(%i154) :lisp $%o153
((BIGFLOAT SIMP 19) 262144 1)
That is, maxima thinks the answer is exactly 1. But
(%i154) %o152-1;
(%o154) 2/1048575
(%i157) is(%o154>= 1/262144/2);
(%o157) true
that is, the difference between the ratio and 1 is greater than 1/2
unit, so maxima should have rounded the answer up to ((bigfloat 19)
262145 1).
Borrowing some code from cmucl, I have implemented a better conversion
from rationals to bfloats. With this change, maxima produces the
expected answer.
A nice side-effect of this change is that the fast bfloat input routine
now appears to produce answers that differs from the slow routine less
than 1% of the time.
Ray