Re: numerical conditioning MAGNITUDE and /
Bradley Lucier 22 Jun 2006 02:52 UTC
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