improving the complex number division

[Gabriel Dos Reis]

Jerome Houdayer <houdayer@spht.saclay.cea.fr> writes:

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

Jérôme --

  What you're quoting above is the generic implementation which is
never ever used for std::complex<float>, std::complex<double> and
std::complex<long double>.  For these specializations, we use:

  template<typename _Tp>
    inline complex<float>&
    complex<float>::operator/=(const complex<_Tp>& __z)
    {
      _ComplexT __t;
      __real__ __t = __z.real();
      __imag__ __t = __z.imag();
      _M_value /= __t;
      return *this;
    }
 
(same for "double" and "long double").  Here, we rely on the builtin
operations (generated by the compiler) for

   __complex__ float, __complex__ double, __complex__ long double.
  
The above cases cover most uses.

| 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)

Yeah, I'm well aware of that.  In either way one needs to make
assumptions for the generic implementation.  I've choosen the grammar
school implementation (less costly) for the moment.

-- Gaby