clang is almost useless for complex arithmetic

[Bruce Evans]

On Sat, 31 Jan 2015, Peter Jeremy wrote:

> To resurrect an old thread...
>
> On 2014-Mar-25 17:22:05 -0700, Steve Kargl <sgk at troutmask.apl.washington.edu> wrote:
>> It appears that clang developers have chosen the naive
>> complex division algorithm, and it does not matter whether
>> one turns CX_LIMITED_RANGE on or off.  This means that
>> if one uses clang with complex types, one must be careful
>> with the range of values allowed in complex division.  In
>> other words, implementation of complex libm routines cannot
>> use complex data types and must fallback to a decomposition
>> into real and imaginary components.
>
> Whilst the fixes don't seem to have made it into FreeBSD yet, it seems
> that this has been improved in recent Clang/LLVM - see
> http://reviews.llvm.org/D5698, committed as rL219557.

The libcalls (as in gcc) allow replacing the method easily.

> I re-implemented their division algorithm in C and checked it against
> the algorithms listed in http://arxiv.org/pdf/1210.4539v2.pdf.  Whilst
> it fails on same of the difficult values listed in that paper, when
> faced with a range of random arguments, it seems to perform better[*]
> than the "robust" algorithm.
>
> I didn't look closely at the multiplication algorithm but, based on a
> quick look, it seems that it's still susceptable to spurious overflows.
>
> [*] or at least closer to the result given by the naive algorithm using
>    x87 long doubles.

Using long doubles would fix many cases without even trying.  But clang
breaks their automatic use even on i386 with -march optimizations, by
using SSE instructions if the arch supports them.  With gcc, it takes
more magic compiler options (-mfpmath=sse+387) to get this bug and this
is still documented as experimental.  The libcalls could use long doubles
even more intentionally, but then in FreeBSD they would have to fight with
the default rounding precision being double (so long doubles are much
like doubles and using them only helps for float complex).

I checked some newer gcc (4.8) support for extra precision.  There is
a -fextra-precision=style option, except it is spelled
-fexcess-precision.  -fexcess-precision=fast gives the historical
broken behaviour.  -fexcess-precision=standard is claimed to give the
C99 behaviour.  It may give it, but it also gives the C11 bug of
destroying extra precision on function return (including for inline
functions).  -std=c99 implies -fexcess-precision=standard.  This is
good for purity, but mostly bad.  It is too surprising that -std=c99,
gives a large pessimization.  Fortunately, -std=c99 is often unusuable
for other reasons.  FreeBSD uses -std=gnu99 and that gives the historical
behaviour.

-fexcess-precision=standard is several years old now, but not supported
by clang.  Destruction of extra precision on return is not done by
clang even with -std=c11.  This goes with not supporting C99 yet.

Bruce