Topic: complex non-member operations --- glitch in standard?


Author: Jan van Dijk <jan@epgmod.phys.tue.nl>
Date: Thu, 3 Jun 2010 14:54:55 CST
Raw View
       Good day,

IMHO, the complex non-member operations are not well-specified in the draft
(and previous standard document).

Let's look at 26.4.6, item 7:

 "template<class T>
 complex<T> operator*(const complex<T>& lhs, const complex<T>& rhs);
 template<class T> complex<T> operator*(const complex<T>& lhs, const T& rhs);
 template<class T> complex<T> operator*(const T& lhs, const complex<T>& rhs);

 Returns: complex<T>(lhs) *= rhs."

The first two operations are OK, but left-multiplication of a complex with a
scalar should definitely return

 complex<T>(rhs) *= lhs

instead of the presently mandated value.

This makes a difference in the presence of infinite numbers (and makes
multiplication non-commutative). As an example, consider the result of
inf*(inf,inf). One expects (inf,inf) but the standard mandates
(inf,0)*(inf,inf) = (inf*inf -0*inf,0*inf+inf*inf)=(nan,nan)...

... while (inf,inf)*inf gives (inf,inf).

(There is also a performance penalty: an IEEE-observing implementation cannot
optimise the dummy 0*i terms away, since nans may be generated --- see above.)

This is not as theoretical as it may seem. I ran into this issue when working
on an implementation of complex Bessel functions: I discovered that I relied
on non-conforming behaviour for almost-zero arguments.

This looks like an editorial oversight, and a duly note that gcc's libstdc++
implements the result I expected.

I am looking forward to hearing your opinion.

With kind regards,

       Jan van Dijk.

--
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@netlab.cs.rpi.edu]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.comeaucomputing.com/csc/faq.html                      ]





Author: =3D?ISO-8859-1?Q?Daniel_Kr=3DFCgler?=3D <daniel.kruegler@googlemail.c=.om>
Date: Fri, 4 Jun 2010 11:50:59 CST
Raw View
On 3 Jun., 22:54, Jan van Dijk <j...@epgmod.phys.tue.nl> wrote:
> IMHO, the complex non-member operations are not well-specified in the
draft
> (and previous standard document).
>
> Let's look at 26.4.6, item 7:
>
>  "template<class T>
>  complex<T> operator*(const complex<T>& lhs, const complex<T>& rhs);
>  template<class T> complex<T> operator*(const complex<T>& lhs, const T&
rhs);
>  template<class T> complex<T> operator*(const T& lhs, const complex<T>&
rhs);
>
>  Returns: complex<T>(lhs) *= rhs."
>
> The first two operations are OK, but left-multiplication of a complex wit=
h
a
> scalar should definitely return
>
>  complex<T>(rhs) *= lhs
>
> instead of the presently mandated value.
>
> This makes a difference in the presence of infinite numbers (and makes
> multiplication non-commutative). As an example, consider the result of
> inf*(inf,inf). One expects (inf,inf) but the standard mandates
> (inf,0)*(inf,inf) = (inf*inf -0*inf,0*inf+inf*inf)=(nan,nan)...
>
> ... while (inf,inf)*inf gives (inf,inf).
>
> (There is also a performance penalty: an IEEE-observing implementation
cannot
> optimise the dummy 0*i terms away, since nans may be generated --- see
above.)
>
> This is not as theoretical as it may seem. I ran into this issue when
working
> on an implementation of complex Bessel functions: I discovered that I
relied
> on non-conforming behaviour for almost-zero arguments.
>
> This looks like an editorial oversight, and a duly note that gcc's
libstdc++
> implements the result I expected.
>
> I am looking forward to hearing your opinion.

Nice catch - I completely agree with your analysis.

Greetings from Bremen,

Daniel Kr=FCgler


--
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use
mailto:std-c++@netlab.cs.rpi.edu<std-c%2B%2B@netlab.cs.rpi.edu>
]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.comeaucomputing.com/csc/faq.html                      ]





Author: Howard Hinnant <howard.hinnant@gmail.com>
Date: Sun, 6 Jun 2010 15:43:29 CST
Raw View
On Jun 3, 4:54 pm, Jan van Dijk <j...@epgmod.phys.tue.nl> wrote:
>        Good day,
>
> IMHO, the complex non-member operations are not well-specified in the draft
> (and previous standard document).

<snip>

> As an example, consider the result of
> inf*(inf,inf). One expects (inf,inf) but the standard mandates
> (inf,0)*(inf,inf) = (inf*inf -0*inf,0*inf+inf*inf)=(nan,nan)...

<snip>

> This looks like an editorial oversight, and a duly note that gcc's libstdc++
> implements the result I expected.
>
> I am looking forward to hearing your opinion.

The standard doesn't mandate that complex multiplication be carried
out by the simplistic algorithm:

  (a c - b d) + (b c + a d)i

(though that is conforming).  Nor does the standard even mandate that
the floating point arithmetic types even have an infinity
representation.  So how complex arithmetic handles nans and infinities
is QOI (quality of implementation).  The C standard Annex G does a
nice job of recommending (not demanding) how complex arithmetic might
deal with nans and infinities.  Every implementation I'm aware of
follows this advice on platforms where IEEE floating point is
available (I know for sure gcc does).  libc++ does as well.

Glancing at the libc++ implementation, it does indeed implement
operator*(const T& lhs, const complex<T>& rhs);  just as you fear:

http://llvm.org/svn/llvm-project/libcxx/trunk/include/complex

template<class _Tp>
inline _LIBCPP_INLINE_VISIBILITY
complex<_Tp>
operator+(const _Tp& __x, const complex<_Tp>& __y)
{
   complex<_Tp> __t(__y);
   __t += __x;
   return __t;
}

However this implementation prints out:

(inf, inf)

as the answer of inf*(inf,inf).  This happens because (inf,
0)*(inf,inf) is not computed as (inf*inf -0*inf,0*inf+inf*inf).  See
the C standard, Annex G.

If you want to mandate that C++0X correctly handle nans and
infinities, that job is much, much larger than just mandating the
algorithm for this one signature.  At the very least you'll want to
take into account all operator*() and operator/().

As a practical matter, the C++ <complex> vendors appear to be handling
this QOI issue pretty well for the past decade with the standard as it
is.

-Howard


--
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@netlab.cs.rpi.edu]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.comeaucomputing.com/csc/faq.html                      ]