Re: guile & r5rs

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

hjstein@bfr.co.il (Harvey J. Stein) writes:

 > Jim Blandy <jimb@red-bean.com> writes:
 >  > 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:

Actually, it's not conformant even in this sense:

   guile> (define a (/ 1 0))
   guile> a
   +#.#
   guile> (inexact? a)
   #f
   guile> (inexact->exact a)
   ERROR: In procedure inexact->exact in expression (inexact->exact a):
   ERROR: Wrong type argument in position 1: +#.#
   ABORT: (wrong-type-arg)

   Type "(backtrace)" to get more information.

+#.# isn't inexact, so / isn't actually coercing to inexact when
necessary, it's returning another type of object.


Here's the kind of thing I was thinking of.  This is just off the top
of my head, so there are probably much better ways of doing these
things that maybe a numerical analyst on this list could suggest, but
here's a start:

(define (bb->i-quotient n d)
;;; Bignum bignum quotient whose result is inexact, assuming it fits
;;; in an exact.
;;; n/d = (quotient n d) + (/ (remainder n d) d)
;;;     = (quotient n d) + (/ (* 2 (remainder n d)) d 2)
;;;     = (quotient n d) + .5*(/ (* 2 (remainder n d)) d)
  (define (eiq n d iters)
    (if (= iters 0)
	0
	(+ (quotient n d)
	   (* .5
	      (eiq (* 2 (remainder n d)) d (- iters 1))))))
  (eiq n d 60))

(define (dumb-rationalize n)
;;; Convert the inexact n to a pair (numerator . denomenator) of
;;; exacts such that (/ numerator denomenator) = n
  (let* ((integer-part (inexact->exact (floor n)))
	 (fractional-part (- n integer-part)))
    (let loop ((num 0)
	       (denom 1))
      (let ((q (bb->i-quotient num denom)))
	(cond ((= fractional-part q) (cons (+ (* denom integer-part) num)
					   denom))
	      (else (if (<= (bb->i-quotient (+ (* 2 num) 1) (* 2 denom))
			    fractional-part)
			(loop (+ (* 2 num) 1)
			      (* 2 denom))
			(loop (* 2 num) (* 2 denom)))))))))

(define (bi->i-quotient n d)
;;; Bignum inexact quotient whose result is inexact, assuming it fits
;;; in an inexact.
  (let ((r (dumb-rationalize d)))
    (bb->i-quotient (* n (cdr r)) (car r))))


(define (bi->i-product big inexact)
;;; Bignum inexact product whose result is inexact, assuming it fits
;;; in an inexact.
  (let ((r (dumb-rationalize inexact)))
    (bb->i-quotient (* big (car r)) (cdr r))))



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