Re: guile & r5rs

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

Radey Shouman <Radey_Shouman@splashtech.com> writes:

 > I think the following function does a reasonable
 > job of dividing two bignums and returning a floating
 > point result.  It doesn't look as though it
 > should be too hard to render into C, or to extend
 > to sqrt and log.
 > 
 > A conservative test for whether something like this is necessary
 > should also be fairly quick, just look at the number of big digits in
 > each of the arguments.
 > 
 > (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.

You're trying to scale down the numerator & denomenator the minimum
necessary so that they they both fit into an inexact.  This should be
shifting them NF to the right.  What took a while is understanding
what the + fact/2 is for?  I finally realized that the difference
between:

   (ash num (- nf))

and:

   (ash (+ num fact/2) (- nf)))

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

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

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?  Of course, doing it at compile time in the C code
is preferable.  At least on my system, /usr/include/values.h includes
DMAXEXP (amongst other useful constants).

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