numerical conditioning MAGNITUDE and /

[Aubrey Jaffer <agj@xxxxxxxxxxxx>]

Writing arithmetic routines in Scheme does not eliminate the need to
do the proper numerical analysis.  For example, the naive
implementation of complex magnitude and / are:

(define (square z) (* z z))

(define (mag z)
  (sqrt (+ (square (real-part z)) (square (imag-part z)))))

(define (div z1 z2)
  (define den (+ (square (real-part z2)) (square (imag-part z2))))
  (make-rectangular (/ (+ (* (real-part z1) (real-part z2))
			  (* (imag-part z1) (imag-part z2)))
		       den)
		    (/ (- (* (imag-part z1) (real-part z2))
			  (* (real-part z1) (imag-part z2)))
		       den)))

But these routines both generate intermediate results larger or
smaller than their inputs, overflowing on:

  (/ 1e300+1e300i (* 4 1e300+1e300i))
  (magnitude 1e300+1e300i)

and underflowing on:

  (/ 1e-300+1e-300i (* 4 1e-300+1e-300i))
  (magnitude 1e-300+1e-300i)

MzScheme is broken for both magnitude and /; MIT-Scheme only for /.
Tests for both cases are in "r4rstest.scm"
(http://cvs.savannah.gnu.org/viewcvs/*checkout*/scm/scm/r4rstest.scm?rev=HEAD).

Does SRFI-77 (or R5RS) mandate that floating-point procedures work
for arguments generating the full range of possible outputs?

			      -=-=-=-=-

A correct implementation is:

(define (mag z)
  (define c (abs (real-part z)))
  (define d (abs (imag-part z)))
  (if (< d c)
      (* c (sqrt (+ 1 (square (/ d c)))))
      (if (zero? d) d (* d (sqrt (+ 1 (square (/ c d))))))))

(define (div z1 z2)
  (define a (real-part z1))
  (define b (imag-part z1))
  (define c (real-part z2))
  (define d (imag-part z2))
  (if (< (abs d) (abs c))
      (let ((r (/ d c)))
	(define den (+ c (* d r)))
	(make-rectangular (/ (+ a (* b r)) den)
			  (/ (- b (* a r)) den)))
      (let ((r (/ c d)))
	(define den (+ d (* c r)))
	(make-rectangular (/ (+ b (* a r)) den)
			  (/ (- a (* b r)) den)))))

			      -=-=-=-=-

[1] David Goldberg, "What Every Computer Scientist Should Know About
Floating-Point Arithmetic",
http://cch.loria.fr/documentation/IEEE754/ACM/goldberg.pdf