Re: guile & r5rs

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

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

 > > Here's the kind of thing I was thinking of.  This is just off the top
 > > of my head, so there are probably much better ways of doing these
 > > things that maybe a numerical analyst on this list could suggest, but
 > > here's a start:
 > 
 > Yeesh.  The algorithm I posted runs in a constant number of bignum
 > ops, and lets you specify an exact lower bound on the number of bits
 > of precision you get.

At least it works and is runnable...

But, I didn't mean to make a proposal, just to give an example of how
these things can be done.

BTW, the algorithms I posted also operate in a constant number of
bignum ops.  It's just a large constant :).

You said:

 > 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).

But this fails when this number is > the largest double, which is when the
result is within 2^p of max_double.

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.

Also, I can't find scalb on my system (linux) either in the header
files or in the libs.  How does it differ from ldexp?

I guess you could get around using scalb & get somewhat better
accuracy by doing as you suggest, but in 2 pieces - right-shift the
bignum p bits to get the integer part, convert it to a double, then
add in the bits that were dropped/2^p (as a float).  The GNU MP
package has code for right-shifting bignums, but however guile
implements them this shouldn't be too hard.

 > What I could really use is a good bignum sqrt algorithm.

There's also one in GNU MP.

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