Re: Fast Complex Division.

[Toon Moene <toon at moene dot indiv dot nluug dot nl>]

Richard Henderson wrote:

[ Oh, BTW, just s/nominator/numerator/g throughout - except for
  "denominator".  Not a native speaker and all that ... ]

> On Sun, Apr 11, 1999 at 02:18:23PM +0200, Toon Moene wrote:
> >       if (abs(c) < abs(d))
> >               den = (1 + (c/d)^2)
> >               nom = (a/d)*(c/d) + (b/d) + i (-(a/d) + (b/d)*(c/d))
> >       else
> >               den = (1 + (d/c)^2)
> >               nom = (a/c) + (b/c)*(d/c) + i (-(a/c)*(d/c) + (b/c))

> These are just close enough in form I would recommend recoding this as
> 
>         tmp1 = c; tmp2 = d; tmp3 = ...
>         if (abs(c) < abs(d))
>                 tmp1 = d; tmp2 = c; tmp3 = ...
>         den = 1 + (tmp2/tmp1)^2
>         nom = ...

Thanks - BTW, I just found a formulation in Numerical Recipes that will
save us a floating point division on each complex divide (by not
dividing numerator and denominator by d^2, but by d).

> > The problematic part is the representation of the if (...) then ... else
> > ...  I do not know how to code that.
> 
>         mode = GET_MODE (abs_c);
>         size = GEN_INT (GET_MODE_BITSIZE (mode));
>         align = GEN_INT (GET_MODE_SIZE (mode));
>         lab = gen_label ();
>         emit_cmp_and_jump_insns(abs_c, abs_d, LT, size, mode, align, 0, lab);
> 
>         // emit then code
> 
>         lab2 = gen_label ();
>         emit_jump_insn (gen_jump (lab2));
>         emit_label (lab);
> 
>         // emit else code
> 
>         emit_label (lab2);

Thanks; that's exactly the guidance I needed.

> > [ BTW, I now see that this hinges on the concept of "gradual underflow",
> >   i.e. it will work on systems where, given that abs(c) < abs(d), c/d
> >   will be zero for abs(c) << abs(d) - this means that it still fails
> >   on the Alpha in some cases ]
> 
> It'll still work if Alpha underflows to zero, won't it?  That's the
> normal case when not -mieee.  Just so long as there aren't denormal
> _inputs_ we should be fine.

Yep - sorry, I'm not *that* familiar with the Alpha floating point rules
(perhaps when I have one on my desk I will).

Besides, reading Numerical Recipes on complex arithmetic, they scared me
with another thing I hadn't given a thought:  One needs to choose a
branch for taking "the" square root of a complex number (goes to show
that I never use complex arithmetic - duh!):

	There are two solutions to:		z1^2 = z2

I'll work further on this in the coming weekend (in time, I hope, for
the feature freeze of the 21st).

[ I'll probably punt on the definition of "square root of complex mode
  integers", because I have no idea what that would *mean* ]

Regards and thanks for your input,

-- 
Toon Moene (toon@moene.indiv.nluug.nl)
Saturnushof 14, 3738 XG  Maartensdijk, The Netherlands
Phone: +31 346 214290; Fax: +31 346 214286
g77 Support: fortran@gnu.org; egcs: egcs-bugs@cygnus.com