Re: guile & r5rs

[hjstein@bfr.co.il (Harvey J. Stein)]

Jim Blandy <jimb@red-bean.com> writes:

 > >  > Guile could do better.  To compute B / C, where B and C are
 > >  > bignums, as a floating-point value with a mantissa at most p bits
 > >  > long, Guile should find the closest integer to (B * 2^p) / C ==
 > >  > (B/C) * 2^p, using only exact bignum arithmetic, and then use the
 > >  > IEEE 754 scalb function to simultaneously divide by 2^p and convert
 > >  > to double, without losing bits.
 > > 
 > > The integer you refer to is (quotient (* B (expt 2 mantissa_bits)) C)
 > > or this number + 1 (depending on whether the remainder is > or <
 > > C/2).
 > 
 > Right.
 > 
 > > But this fails when this number is > the largest double, which is when the
 > > result is within 2^p of max_double.
 > 
 > You mean, it fails when the result doesn't fit into a double?  That's
 > fine; this is only to be used when the remainder is non-zero, in which
 > case R5RS leaves us no choice (given that we don't have rationals, or
 > inexact bignums) but to use a double.  Maybe I'm missing something.

The result (number/2^p) might fit in a double, but number might not.

 > > This also fails when C is large relative to B.  In particular, if C >
 > > B*2^p, then the closest integer is zero, even though the quotient is
 > > easily representable in a double.
 > 
 > You're right.  I guess we need to multiply B by 2 until it is between
 > C*2^(p-1) and C*2^p, and then divide by the appropriate amount.

Right.  I haven't thought carefully about all the border cases, but I
think this mostly fixes it.  It doesn't have to do with p being the
number of bits in a double's mantissa, but with choosing q such that
(quotient (* B (expt 2 q)) c) has p bits of precision, or
equivalently, if Bp = the index of the high order bit in b, and Cp =
the index of the high order bit in c, then shift Bp to the right
Cp+p-Bp places.

-- 
Harvey J. Stein
BFM Financial Research
hjstein@bfr.co.il