[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

numerical conditioning MAGNITUDE and /

This page is part of the web mail archives of SRFI 77 from before July 7th, 2015. The new archives for SRFI 77 contain all messages, not just those from before July 7th, 2015.

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)))
		    (/ (- (* (imag-part z1) (real-part z2))
			  (* (real-part z1) (imag-part z2)))

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"

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",