real128 cmplx ATAN, ATANH - quality of implementation?

[Steve Kargl]

On Fri, Feb 03, 2017 at 05:45:15AM -0800, Anton Shterenlikht wrote:
> Behaviour of complex intrinsics such as ATAN or ATANH
> at branch points is not defined by the standard (I think).
> However, for complex of kind real32 and real64
> gfortran 6 and 7 give very helpful values at branch points.
> Not so for real128.
> 
> The full details: http://cmplx.sourceforge.net/
> 
> It would be very helpful if the return values of ATAN and
> ATANH at branch points for real128 matched those for
> real64 and real32 kinds.
> 

You need to do two things.  First, recalling that you are a FreeBSD,
you should file a number of bug reports with FreeBSD for the missing
long double complex functions.  You probably want to have one bug
report for each missing function to hopefully highlight the current
state of FreeBSD's libm.  For example, on FreeBSD atan for REAL(4)
and REAL(8) map to libm's catanf and catan, which according to the
manpage are designed to deal with the branchcut.  For REAL(10) (aka
REAL(REAL128)), atan is implemented in libgfortran's c99_functions.c.
This leads to the second item, you'll need to do: audit c99_functions.c.
A quick inspection of c99_functions.c suggests that there is very little
chance that the fallback implementations of many of the long double
complex functions are correct under a number of conditions.

Here's catanl from c99_functions.c

complex long double
catanl (complex long double z)
{
  return I*clogl ((I+z)/(I-z))/2.0L;
}

'I' here is intended to be the mathmatical i = sqrt(-1), but
gcc (at least older versions as I have not checked trunk) treat
'I' as I = 0 + i1.  If you allow me to write z = x + iy, then 
the above converts to 

   return (0+i1)*clogl((0+i1+x+iy)/(0+i1-x-iy))/2;

which becomes

   return (0+i1)*clogl((x+i(y+1))/(-x-i(y-1)))/2;

inserting, for example, 0-i1, you then get

   return (0+i1)*clogl((0+i0))/(-0-i2))/2;

Now, look at src/lib/msun/math_private.h on FreeBSD, you'll
find this comment:

 * The C99 standard intends x+I*y to be used for this, but x+I*y is
 * currently unusable in general since gcc introduces many overflow,
 * underflow, sign and efficiency bugs by rewriting I*y as
 * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.
 * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted
 * to -0.0+I*0.0.

In short, the fallback functions in c99_functions.c were not designed to
deal with special cases, i.e., branchcuts.

-- 
Steve
20161221 https://www.youtube.com/watch?v=IbCHE-hONow