Re: improving the complex number division

Re: improving the complex number division

[Robert Dewar]

It is indeed surprising to find this entirely naive and quite incorrect
implementation of complex division in gcc. This is certainly a computation
for which absolutely standard algorithms exist. I am not clear whether
Jerome is quoting one of these standard algorithms (in which case a reference
would be helpful), or inventing from scratch (in which case I would still be
a little skeptical, though to my not-really-expert-enough eye, the proposed
replacement algorithm looks correct -- I just don't know if it is the best
choice.

Re: improving the complex number division

[Gabriel Dos Reis]

dewar@gnat.com (Robert Dewar) writes:

| It is indeed surprising to find this entirely naive and quite incorrect

What actually surprised me is that neither Jérôme nor you  noticed
that that implementation isn't being used for anything we have defined
semantics for.

| implementation of complex division in gcc. This is certainly a computation
| for which absolutely standard algorithms exist.

True.  But, as noted in the reply I made to Jérôme's message, that
implementation wasn't never ever used for the cases the C++ standard
defined semantics for (i.e. specializations of std::complex<T> with 
T=float, T=double, T=long double).

| I am not clear whether
| Jerome is quoting one of these standard algorithms

That is standard algorithm.

| (in which case a reference would be helpful),

I believe Goldberg's paper or G. Stewart's book on matrix computations
covers it.

-- Gaby

Re: improving the complex number division

[Gabriel Dos Reis]

dewar@gnat.com (Robert Dewar) writes:

| > What actually surprised me is that neither J=E9r=F4me nor you  noticed
| > that that implementation isn't being used for anything we have defined
| > semantics for.
| 
| But if it is used at all, it is a mistake. 

It is a mistake only if you define precisely under which circumstances
it is a mistake.  It is not a mistake in absolute.

| After all, Only recently has
| C tried to define semantics for fpt at all, but that is not an excuse for
| doing something completely wrong :-)

True, but irrelevant here.

-- Gaby

Re: improving the complex number division

[Robert Dewar]

> What actually surprised me is that neither J=E9r=F4me nor you  noticed
> that that implementation isn't being used for anything we have defined
> semantics for.

But if it is used at all, it is a mistake. After all, Only recently has
C tried to define semantics for fpt at all, but that is not an excuse for
doing something completely wrong :-)

Re: 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

improving the complex number division

[Jerome Houdayer]

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.houdayer@polytechnique.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
===========================================================