http://people.csail.mit.edu/jaffer/III/ComplexMagnitude | |

## Complex Numerical Analysis |

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