Re: numerical conditioning MAGNITUDE and /

[Bradley Lucier <lucier@xxxxxxxxxxxxxxx>]

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