I should correct some mistakes in the algorithms I
sketched, in case someone actually uses them.
For number->string with a specified precision < 53,
my second while loop was infinite. This is better:
s := (number->string x);
p := precision;
#| this loop terminates with p <= 53 |#
while (not (eqv? (string->number
(string-append s "|" (number->string p)))
x)) {
p := p + 1;
}
s2 := trimFinalDigitRoundingDown (s);
s3 := trimFinalDigitRoundingUp (s);
#| now trim as many digits from s as you can |#
while (or (eqv? (string->number
(string-append s2 "|" (number->string p)))
x)
(eqv? (string->number
(string-append s3 "|" (number->string p)))
x))
if (eqv? (string->number
(string-append s2 "|" (number->string p))))
then s := s2
else s := s3;
s2 := trimFinalDigitRoundingDown (s);
s3 := trimFinalDigitRoundingUp (s);
}
return (string-append s "|" (number->string p));
That algorithm isn't quite as bad as it looks. If
string->number performs correct rounding, then the
first while loop won't increment p at all. The second
while loop can be sped up by precomputing the maximum
number of digits that might be required for precisions
less than 53, so its body should execute at most once.
In practice, the most likely precisions are going to
be 24, 53, and 64, so the above algorithm will almost
always be used with precision = 24, and you can turn
it into straight-line code for that case.
I made a mistake when describing my quick-and-easy
way to perform correct rounding from 53 to 24 bits.
You just look at the 29 = 53 - 24 bits that you're
going to drop. The only hard case is when those bits
start with a 1 and the remaining bits are 0. You can
use AlgorithmM [1] to perform correct rounding for that
case, or you can cheat by returning the 53-bit
approximation (which will also be the best 25-bit
approximation).
When correctly rounding random inputs from 53 to 24
bits, the above algorithm avoids the hard case about
99.9999998% of the time, so the cost of using
AlgorithmM in the hard cases is unmeasurable. When
the inputs are printed results of a single precision
computation, which is where they will come from in
practice, and denormalized numbers are printed with
the correct mantissa width, then the hard case will
*never* arise.
Let's do Alan Watson's example:
> Think of (number->string 1.1 10 1) on an implementation
> that uses IEEE doubles.
Note that Alan is asking us to round to 1 bit, which
will be the hidden bit, so we are going to discard 52
bits of precision here.
Let's assume a big-endian implementation. The binary
representation of double precision 1.1 is
#x3ff199999999999a
so the 52 low-order bits are #x199999999999a. That is
not #x8000000000000, so this is not a hard case. We
just round down to #x3ff0000000000000, which is (as
expected) the double precision representation of 1.0.
The result of (number->string 1.1 10 1) is therefore
"1.|1".
> I am not sure what problem the editors are trying to solve with |p, but
> despite this ignorance I am willing, nay, eager to put money on |p not
> being the solution.
How much?
> Here is the problem I want to solve. I want to be able to read IEEE
> singles and doubles and have them represented by the nearest inexact
> real. I want to be able to write the IEEE single and double that is
> nearest to a given inexact real. I am prepared to accept a solution that
> works on many but not all implementations and I am prepared to do some
> of the work myself.
If the R6RS were to require implementations to read x|p
as the best p-bit binary floating point approximation to
x, then the R6RS would solve your problem.
The current draft R6RS does not require that, and I don't
see how it could be revised to require that without also
requiring implementations to represent inexact numbers in
binary floating point. In practice, however, I think the
current draft R6RS would accomplish exactly what you want,
assuming implementors of the R6RS pay attention to this.
> (a) We keep |p but give is some kind of sane semantics.
The meaning of x|p is mathematically well-defined, and is
pretty clear from the current draft: x|p "represents the
best binary floating-point approximation of x using a
p-bit significand." The word "represents" might be
changed to "indicates". The word "best" is ambiguous,
but I think it's fairly clear that it should be defined
as in my paper and IEEE-754R, as a nearest approximation,
rounding to the even significand in case of ties [2].
(In IEEE-754R, the best approximation depends on the
rounding mode, but I think we should assume the default
round-to-even mode for Scheme.)
The real questions about the semantics of x|p in the
current draft concern implementations' responsibilities.
The current draft gives implementations a lot of leeway
with respect to inexact arithmetic; for example, the
current draft does not require inexact arithmetic to
be implemented using any kind of floating point, sane
or otherwise.
In my opinion, standard practice should be to read x|p
with no loss of precision or accuracy when p <= 53, and
to represent x|p as accurately as possible if p exceeds
the highest precision provided by a implementation.
For output, I think the best practice would be to omit
the mantissa width when printing IEEE double precision
unless the number is denormalized, and to provide an
explicit mantissa width for denormalized numbers and
for all numbers whose precision is not 53.
Will
[1]
http://larceny.ccs.neu.edu/larceny-trac/browser/trunk/larceny_src/src/Lib/Common/belle.sch
[2]
ftp://ftp.ccs.neu.edu/pub/people/will/howtoread.ps
Received on Fri Jun 29 2007 - 19:11:18 UTC