Re: guile & r5rs

[Radey Shouman <Radey_Shouman@splashtech.com>]

Harvey J. Stein writes:
 >  > (define *float-expt* 1024)  ;; this depends on float representation.
 >  >                             ;; probably there should be some way of
 >  >                             ;; determining this from the Scheme level.
 >  > (define (big-i/ num den)
 >  >   (let* ((nf (max 0
 >  >                   (- (max (integer-length num)
 >  >                           (integer-length den))
 >  >                      *float-expt*)))
 >  >          (fact/2 (ash 1 (- nf 1))))
 >  >     (/ (exact->inexact (ash (+ num fact/2) (- nf)))
 >  >        (exact->inexact (ash (+ den fact/2) (- nf))))))
 > 
 > That's very cool.  It took me a while to understand your algorithm.  I
 > didn't realize the fcns needed (integer-length & ash) were in fact
 > available in guile.

They have been inherited from SCM.  It is also possible to write them
in portable Scheme, (I believe they are in SLIB), but I'm sure that
would cost quite a bit in Guile/SCM.

 > is that the former is equal to (truncate (/ num (expt 2 nf))) and the
 > latter is (round (/ num (expt 2 nf))).  Very nice.  Now what do you do
 > to make it round towards even? :).

I actually did think about this, but couldn't think of a way to do it
without allocating more bignums than it seemed to be worth.  A
scale-and-round-bignum sort of primitive might be useful if we can
figure out exactly what it should do.

 > However, there seems to be a bug:
 > 
 > guile> (big-i/ (expt 10 2048) (expt 10 2048))
 > 1.0
 > guile> (big-i/ (expt 2 2048) (expt 2 2048))
 > 1.0
 > guile> (big-i/ (- (expt 2 2048) 1) (- (expt 2 2048) 1))
 > #.#
 > 
 > Either *float-expt* is too big or the computation for nf becomes wrong
 > after adding fact/2 to effect rounding.  This happens when all the
 > high order bits are 1 (hence it fails when the numerator is 2^2048-1).

You're right, this happens when both arguments have the same
integer-length and one or both are rounded up -- all high order bits
don't have to be 1.  This can be fixed by making *float-expt* 1023 for
IEEE doubles.  It could also be fixed by adding a test for the
integer-length of the rounded bignums before converting.  This test
gains nothing when the numerator and denominator are about the same
size, since we throw away all but 52 bits of mantissa, but would help
a little when one is greater than the other by a factor near 2^1024.


 > Btw, couldn't you initialize *float-expt* by doing a binary search for
 > the largest n s.t. (exact->inexact (expt 2 n)) is not +#.# give the
 > appropriate value?

Sure, back when there was a little more variety in floating point
formats this was the done thing -- eg machar.f on netlib.