[r6rs-discuss] Confusion over |p

From: Alan Watson <alan>
Date: Sat, 30 Jun 2007 18:59:32 -0500

William D Clinger wrote:

> s := (number->string x);
> p := 53;
> y := x;
> #| this loop terminates with p <= 53 |#
> while (and (< precision p)
> (eqv? y x)) {
> p := p - 1;
> #| compute the nearest p-bit approximation to x |#
> y := roundToNearest (x, p);
> }
> if (< precision p)
> then p := p + 1;
> return (string-append s "|" (number->string p));

This is an algorithm for the output problem:

   (string->number x 10 precision).

If internally we use binary floating-point numbers for inexact reals, I
think calculating p is actually quite easy: we find the smallest
integral m that satisfies

   (= x (* m (expt 2 k)))

in which k is also integral. We then note that the number of bits of
precision required to represent x is just (bitwise-length m). Therefore,
p is

   (max precision (bitwise-length m))

I think this is identical to your explicit loop; I have an implicit loop
when I find the smallest integral m above.

However, I don't think we can then simply write

   (string-append (number->string x) "|" (number->string p))

as this need not be the shortest decimal representation that gives x
when read back. For example, consider the IEEE single-precision number
whose bytes are as follows on an i386:

   #xa4 #x85 #xf2 #x46

If I use the algorithms given by Burger & Dybvig to write this assuming
it will be read back in single precision, I get 0.3104282e5. However, if
I assume it will be read back in double precision, I get
0.310428203125e5. I think we have to use an appropriate algorithm from,
e.g., Steele & White, Gay, or Burger & Dybvig to find the shortest
decimal representation that is equal to x when rounded to p bits.

Now the input problem:
   (number->string "x|p")

We could use your Bellerophon or something similar directly.
Alternatively, we could follow D'Argo and convert x to the full internal
precision using the standard string->number. We could then round to p
bits, which will give either the correct result or a number adjacent to
the correct result. Finally, we use something like your AlgorithmR to
find the correct result.

I think handling subnormals and overflow will be easier in the second
approach. I'll see if I can come up with something.

Regards,

Alan
Received on Sat Jun 30 2007 - 19:59:32 UTC