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

Re: 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.




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