Re: guile & r5rs

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

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

 > >    R4RS says:
 > > 
 > > 	Implementations are encouraged, but not required, to support exact
 > >       integers and exact rationals of practically unlimited size and
 > >       precision, and to implement the above procedures and the `/'
 > >       procedure in such a way that they always return exact results when
 > >       given exact arguments.  If one of these procedures is unable to
 > >       deliver an exact result when given exact arguments, then it may
 > >       either report a violation of an implementation restriction or it
 > >       may silently coerce its result to an inexact number.  Such a
 > >       coercion may cause an error later.
 > > 
 > >    so I guess this isn't completely allowed behavior for /.  It's not
 > >    coercing the result to a double & it's not reporting an implementation
 > >    restriction violation.
 > 
 > Ahh, this is subtle.  The result fits in a double, but neither of the
 > operands does, so you can't convert to a double and then do the
 > division.  This behavior is unfortunate.

Exactly.  You have to be serious about doing your computations on
bignums.

 > I think Guile's behavior is legal, though.  The "encouraged, but not
 > required" clause applies only to the rest of the sentence.  Since
 > Guile does not support exact rationals, it silently coerces its result
 > to an inexact number.  +#.# is an inexact positive infinity.  R5RS
 > doesn't specify the precision with which the result is computed, only
 > the precision with which it is represented.

Correct, but that sentence is about bignums & rationals.  It's the
next sentence that's at issue:

       If one of these procedures is unable to deliver an exact result
       when given exact arguments, then it may either report a
       violation of an implementation restriction or it may silently
       coerce its result to an inexact number.

So, yes, I guess +#.# is an inexact number.  I guess it's an issue of
what "coerce" is supposed to mean.

 > 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.
 > 
 > My copy of the IEEE floating point standard arrived yesterday, and I'm
 > getting so much satisfaction from knowing what's actually going on.  I
 > feel so geeky.  Unfortunately, they misprinted it(!!!), so order it
 > electronically if you can.

That's the sort of thing I was thinking of.  Or, just do long
division until you get the right number of bits.

 > >    It also says:
 > > 
 > > 	In particular, implementations that use flonum representations must
 > >       follow these rules: A flonum result must be represented with at
 > >       least as much precision as is used to express any of the inexact
 > >       arguments to that operation.  It is desirable (but not required)
 > >       for potentially inexact operations such as `sqrt', when applied to
 > >       exact arguments, to produce exact answers whenever possible (for
 > >       example the square root of an exact 4 ought to be an exact 2).  If,
 > >       however, an exact number is operated upon so as to produce an
 > >       inexact result (as by `sqrt'), and if the result is represented as
 > >       a flonum, then the most precise flonum format available must be
 > >       used; but if the result is represented in some other way then the
 > >       representation must have at least as much precision as the most
 > >       precise flonum format available.
 > > 
 > >    so this is definitely *not* allowed for sqrt.
 > 
 > Yeah, this is basically the same thing; he's converting to double
 > before doing the sqrt.  The result fits in a double, but the argument
 > does not, so the value dies on the way in.  sqrt gets passed an
 > infinity, and returns it.
 > 
 > I don't see that Guile is doing anything illegal, though.  Could you
 > explain exactly why you think Guile is breaking the rules?

The above paragraph first says that the following rules must be
followed in implementations using flonums.  One of the rules is:

       If, however, an exact number is operated upon so as to produce
       an inexact result (as by `sqrt'), and if the result is
       represented as a flonum, then the most precise flonum format
       available must be used;

My reading of this is:

   When an arithmetic operation on exacts produces inexacts and is
   represented internally as a flonum, the result should be as close
   to the true mathematical result of the operation as can be
   represented in the flonum.

Actually, now that I'm looking at it again, it's possible that they're
just saying that, for example, if you're internally representing
inexacts as floats & doubles, then when an operation on exacts
produces an inexact, you must use a double for the result.  It doesn't
say that the value in the double must have any connection with
reality.

But, why specify that you have to use the moset precise format
available but not that the result should be meaningful?  Maybe an R5RS
author can comment on the issue (if anyone has one readily
available).  However, I did find a different paragraph (in R5RS)
relating to the issue, which seems to say that guile isn't conformant:

   If two implementations produce exact results for a computation that
   did not involve inexact intermediate results, the two ultimate
   results will be mathematically equivalent. This is generally not
   true of computations involving inexact numbers since approximate
   methods such as floating point arithmetic may be used, but it is
   the duty of each implementation to make the result as close as
   practical to the mathematically ideal result.

I'm refering the the second half of the last sentence - when creating
inexacts, each implementation must make the inexact as close the the
matematically ideal result as possible.  This is what I was
interpreting the other paragraph to say.

So, taking everything together, I'd say that when they say "coerce",
they mean with as much precision as practical.

 > To compute the square root of a bignum B as a floating-point value
 > with a mantissa p bits long, Guile should use a similar trick: find
 > the closest integer to sqrt (B * 2^(2p)) == sqrt (B) * sqrt (2^(2p))
 > == sqrt (B) * 2^p, using only exact bignum arithmetic, and then use
 > scalb to simultaneously divide by 2^p and convert to double, without
 > losing bits.

Or directly apply a sqrt algorithm.  I'll write up a little code &
post it shortly.

 > IHNW.  IJLS "SDB2^pACTD,WLB".

Are you ok?

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