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

Re: numerical conditioning MAGNITUDE and /




On Jun 19, 2006, at 2:47 PM, Aubrey Jaffer wrote:

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

This gives the wrong answer on, e.g.,

> (mag +inf.+inf.i)
+nan.

And the code's not symmetric in d and c, which it should be, so it's immediately suspect.


(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)))))

This code has similar problems with IEEE 754 infinities:

> (div 1.0+0.0i +inf.+inf.i)
+nan.+nan.i