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;
}
