improving the complex number division

[Jerome Houdayer <houdayer at spht dot saclay dot cea dot fr>]

Hi,

In gcc-3.2.2/libstdc++-v3/include/std/std_complex.h

the complex division is implemented like this:

// 26.2.5/15
// XXX: This is a grammar school implementation.
template<typename _Tp>
template<typename _Up>
complex<_Tp>&
complex<_Tp>::operator/=(const complex<_Up>& __z)
{
  const _Tp __r =  _M_real * __z.real() + _M_imag * __z.imag();
  const _Tp __n = norm(__z);
  _M_imag = (_M_imag * __z.real() - _M_real * __z.imag()) / __n;
  _M_real = __r / __n;
  return *this;
}

The problem is that the function norm() may overflow whereas the whole
division operation would not (typically when abs(__z) is larger than the
sqrt of the maximum value for _Tp)

If the arguments are z1 = a + i * b and z2 = c + i * d  (with i*i = -1)
the current implementation computes:
n = c * c + d * d   (that's norm(z2) and that may overflow)
z1 / z2 = (a * c + b * d) / n + i * (b * c - a * d) / n

it would be better to procede as follows:

if abs(c) >= abs(d)
x = d / c
n = c + d * x
z1 / z2 = (a + b * x) / n + i * (b - a * x) / n

if abs(c) < abs(d)
x = c / d
n = d + c * x
z1 / z2 = (a * x + b) / n + i * (b * x - a) / n

this would gives:

template<typename _Tp>
template<typename _Up>
complex<_Tp>&
complex<_Tp>::operator/=(const complex<_Up>& __z)
{
  if (abs(__z.real())>=abs(__z.imag())) {
    const _Tp __x = __z.imag() / __z.real();
    const _Tp __n = __z.real() + __z.imag() * x;
    const _Tp __r = _M_real + _M_imag * x;
    _M_imag = (_M_imag - _M_real * x) / __n;
    _M_real = __r / __n;
  }
  else {
    const _Tp __x = __z.real() / __z.imag();
    const _Tp __n = __z.imag() + __z.real() * x;
    const _Tp __r = _M_real * x + _M_imag;
    _M_imag = (_M_imag * x - _M_real) / __n;
    _M_real = __r / __n;
  }
  return *this;
}

Note that I didn't test it...

Best regards,

===========================================================
Jérôme HOUDAYER           jerome dot houdayer at polytechnique dot org

Service de physique théorique
CEA/Saclay - Orme des Merisiers      tel: +33 1 69 08 50 33
F-91191 Gif-sur-Yvette Cedex         fax: +33 1 69 08 81 20
FRANCE
===========================================================