Menu
▾
▴
Re: [Numpy-discussion] complex division
|
From: Charles R H. <cha...@gm...> - 2006-02-20 00:12:13
|
Hmm...
The new algorithm does look better with respect to overflow and
underflow, but I wonder if it is not a bit of overkill. It seems to me
that the same underflow/overflow problems attend complex
multiplication, which is pretty much all that goes on in the original
algorithm.
One thing I do know is that division is expensive. I wonder if one
division and two multiplications might be cheaper than two divisions.
I'll have to check that out.
Chuck
On 2/19/06, Tim Hochberg <tim...@co...> wrote:
>
> While rummaging around Python's complexobject.c looking for code to
> steal for complex power, I came across the following comment relating to
> complex division:
>
> /****************************************************************=
**
> This was the original algorithm. It's grossly prone to spurious
> overflow and underflow errors. It also merrily divides by 0 desp=
ite
> checking for that(!). The code still serves a doc purpose here, =
as
> the algorithm following is a simple by-cases transformation of th=
is
> one:
>
> Py_complex r;
> double d =3D b.real*b.real + b.imag*b.imag;
> if (d =3D=3D 0.)
> errno =3D EDOM;
> r.real =3D (a.real*b.real + a.imag*b.imag)/d;
> r.imag =3D (a.imag*b.real - a.real*b.imag)/d;
> return r;
> *****************************************************************=
*/
>
> /* This algorithm is better, and is pretty obvious: first
> divide the
> * numerators and denominator by whichever of {b.real, b.imag} ha=
s
> * larger magnitude. The earliest reference I found was to CACM
> * Algorithm 116 (Complex Division, Robert L. Smith, Stanford
> * University). As usual, though, we're still ignoring all IEEE
> * endcases.
> */
>
> The algorithm shown, and maligned, in this comment is pretty much
> exactly what is done in numpy at present. The function goes on to use
> the improved algorithm, which I will include at the bottom of the post.
> It seems nearly certain that using this algorithm will result in some
> speed hit, although I'm not certain how much. I will probably try this
> out at some point and see what the speed hit, but in case I drop the
> ball I thought I'd throw this out there as something we should at least
> look at. In most cases, I'll take accuracy over raw speed (within reason)=
.
>
> -tim
>
>
>
>
>
>
>
> Py_complex r; /* the result */
> const double abs_breal =3D b.real < 0 ? -b.real : b.real;
> const double abs_bimag =3D b.imag < 0 ? -b.imag : b.imag;
>
> if (abs_breal >=3D abs_bimag) {
> /* divide tops and bottom by b.real */
> if (abs_breal =3D=3D 0.0) {
> errno =3D EDOM;
> r.real =3D r.imag =3D 0.0;
> }
> else {
> const double ratio =3D b.imag / b.real;
> const double denom =3D b.real + b.imag * ratio;
> r.real =3D (a.real + a.imag * ratio) / denom;
> r.imag =3D (a.imag - a.real * ratio) / denom;
> }
> }
> else {
> /* divide tops and bottom by b.imag */
> const double ratio =3D b.real / b.imag;
> const double denom =3D b.real * ratio + b.imag;
> assert(b.imag !=3D 0.0);
> r.real =3D (a.real * ratio + a.imag) / denom;
> r.imag =3D (a.imag * ratio - a.real) / denom;
> }
> return r;
>
>
>
>
> -------------------------------------------------------
> This SF.net email is sponsored by: Splunk Inc. Do you grep through log fi=
les
> for problems? Stop! Download the new AJAX search engine that makes
> searching your log files as easy as surfing the web. DOWNLOAD SPLUNK!
> http://sel.as-us.falkag.net/sel?cmd=3Dlnk&kid=3D103432&bid=3D230486&dat=
=3D121642
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>
|