trying to divide complex numbers, division by zero

Question

I'm trying the program below to divide complex numbers, it works for complex numbers but not when the denominator is real (i.e, the complex part is zero). Division by zero occurs in this line ratio = b->r / b->i ;, when the complex part b->i is zero (in the case of a real denominator).

How do I get around this? and why did the programmer do this, instead of the more straightforward rule for complex division

The wikipedia rule seems to be better, and no division by zero error would occur here. Did I miss something? Why did the programmer not use the wikipedia formula??

Thanks

/*! @file dcomplex.c
 * \brief Common arithmetic for complex type
 *
 * <pre>
 * -- SuperLU routine (version 2.0) --
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
 * and Lawrence Berkeley National Lab.
 * November 15, 1997
 *
 * This file defines common arithmetic operations for complex type.
 * </pre>
 */

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include "slu_dcomplex.h"


/*! \brief Complex Division c = a/b */
void z_div(doublecomplex *c, doublecomplex *a, doublecomplex *b)
{
    double ratio, den;
    double abr, abi, cr, ci;

    if( (abr = b->r) < 0.)
        abr = - abr;
    if( (abi = b->i) < 0.)
         abi = - abi;
    if( abr <= abi ) {
        if (abi == 0) {
            fprintf(stderr, "z_div.c: division by zero\n");
            exit(-1);
        }   
        ratio = b->r / b->i ;
        den = b->i * (1 + ratio*ratio);
        cr = (a->r*ratio + a->i) / den;
        ci = (a->i*ratio - a->r) / den;
    } else {
        ratio = b->i / b->r ;
        den = b->r * (1 + ratio*ratio);
        cr = (a->r + a->i*ratio) / den;
        ci = (a->i - a->r*ratio) / den;
    }
    c->r = cr;
    c->i = ci;
}

@Gman. Its a C program. See the heading of the comment, its dcomplex.c. Op should remove the c++ tag — Tom
I believe there could be a edge case bug here. Doesn't make sense to be honest seems the programmer wasn't really a mathematician. PS I agree with Tom. — langerra.com
@Tom: Oh, I swear it was C++... (Oh, it was. Read the tags, what was I thinking.) — GManNickG
It's the SuperLU matrix solver from Berkeley, developed in 1997 I think. I would imagine a bug wouldn't have lasted that long. — user553619

Michael Burr Michael Burr · Accepted Answer · 2011-01-05T02:57:36

Looks to me like the original programmer was ensuring that ratio would end up between 0 and 1 (inclusive) probably to ensure that adequate precision was maintained.

A straight forward coding of the algorithm on Wikipedia looks like would have danger of getting some very large divisors.

Also, as far as I can tell, the following test:

   if (abi == 0) {
        fprintf(stderr, "z_div.c: division by zero\n");
        exit(-1);
   }

can only be true if both abi and abr are zero, so you would be in a real divide by zero situation.

(At that point abi and abr have both been forced to be non-negative, and that code path only gets hit when (abr <= abi)).

Perhaps you can post a small test case that shows the divide-by-zero being hit when r is non-zero.

trying to divide complex numbers, division by zero

3 Answers