Re: guile & r5rs

[Jim Blandy <jimb@red-bean.com>]

>  > I don't have any plans to implement the full numeric tower at the
>  > moment.  This is explicitly permitted by R5RS.  Guile does support
>  > exact and inexact numbers.  So, as far as I know, Guile does follow
>  > the numeric standard; if you feel this is not so, please report it as
>  > a bug.
> 
> Actually, I don't think it does, and I did report it as a bug (on Oct
> 29th).  I don't have a copy, but it's in the guile archive (I sent it
> to guile & bug-guile):

Right; I remember reading this.  It's on my list now.


>    Consider:
> 
>    guile> (define a 2706660098140645489185923553970566016893099856360336323402524418695428905878508075836626775246052131230477353094376980271701893510989004673972199336498786760105507295770706655889359381824318115244951819800829175510325405549778256791513205371984083396827897605501338582435065123577807680901404556091719054173009876284055841904995815079252887561888975757510847638427607806239796367139527811680401979800)
>    guile> (define b 10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000)
>    guile> (/ a b)
>    +#.#
>    guile> (sqrt a)
>    +#.#
>    guile> (sqrt b)
>    +#.#
> 
>    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.

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.

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.

>    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?

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.

IHNW.  IJLS "SDB2^pACTD,WLB".