Complex Numerical Analysis

	http://people.csail.mit.edu/jaffer/III/ComplexMagnitude
	Complex Numerical Analysis

Although the definitions of complex arithmetic operations in terms of reals are not complicated, they are not a substitute for 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)))
		       den)
		    (/ (- (* (imag-part z1) (real-part z2))
			  (* (real-part z1) (imag-part z2)))
		       den)))

But these routines both generate intermediate results larger or smaller than their inputs, overflowing on:

  (mag 1e155+1e155i)
  (div 1e155+1e155i 4e155+4e155i)

and underflowing on:

  (mag 1e-170+1e-170i)
  (div 1e-170+1e-170i 4e-170+4e-170i)

[The underflow in div causes a subsequent overflow.]

A full-range 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)))))

The number of (real) flonum operations for computing complex magnitude and / are:

mag	+ or -	*	/	abs	sqrt	<
naive	1	2	0	0	1	0
full-range	1	2	1	2	1	1
div	+ or -	*	/	abs	sqrt	<
naive	3	6	2	0	0	0
full-range	3	3	3	2	0	1

Floating-point abs involves just flipping one bit; it should be very fast. The full-range versions are unlikely to be much slower than naive ones; full-range division is perhaps faster. The full-range division algorithm is Smith's formula given in:

David Goldberg,
What Every Computer Scientist Should Know About Floating-Point Arithmetic,
ACM Computing Surveys, Vol. 23 No. 1, 1991

I am a guest and not a member of the MIT Computer Science and Artificial Intelligence Laboratory. My actions and comments do not reflect in any way on MIT.
		CS and Numerical Analysis
	agj @ alum.mit.edu	Go Figure!