complex division alogrithm?

[Steve Kargl]

On Thu, May 06, 2004 at 03:42:37PM -0700, Richard Henderson wrote:
> On Wed, May 05, 2004 at 08:37:53PM -0700, Steve Kargl wrote:
> > I haven't had time to implement a f2c-like algorithm, but I
> > should be able to get it this weekend.
> 
> Implement at what level?  This is already done at both the
> gimple and rtl level.  See expand_complex_div_wide for
> instance, for the gimple implementation.
> 
> r~

At whatever level gfc_arith_divide() in gcc/fortran/arith.c is
invoked for complex division.  Lines 1439 through 1452 are

      mpf_mul (x, op2->value.complex.r, op2->value.complex.r);
      mpf_mul (y, op2->value.complex.i, op2->value.complex.i);
      mpf_add (div, x, y);

      mpf_mul (x, op1->value.complex.r, op2->value.complex.r);
      mpf_mul (y, op1->value.complex.i, op2->value.complex.i);
      mpf_add (result->value.complex.r, x, y);
      mpf_div (result->value.complex.r, result->value.complex.r, div);

      mpf_mul (x, op1->value.complex.i, op2->value.complex.r);
      mpf_mul (y, op1->value.complex.r, op2->value.complex.i);
      mpf_sub (result->value.complex.i, x, y);
      mpf_div (result->value.complex.i, result->value.complex.i, div);

Unless GMP is called with a sufficient number of bits, the computation
of div can overflow.  My original question was whether GMP and gfortran
are doing the right thing with the above naive implementation of complex
division.

-- 
Steve