complex.h math functions

[Steve Kargl]

On Tue, Oct 04, 2005 at 07:31:43PM +1000, Bruce Evans wrote:
> On Sun, 2 Oct 2005, Steve Kargl wrote:
> 
> >One question:  Can we use a #define for complex float or do we
> >need an trivial C program?  That is, in math.h can we do
> >
> >#define ccoshf(z)  ccosh((float complex)(z))
> 
> A program is needed, since (ccoshf)(z) and &(ccoshf) must work even if
> ccoshf() is a macro.

OK.

> I converted to ccoshf() for testing it and comparing with a simpler
> implementation that uses the formula, then converted back.  There are
> lots of special cases for -0, Infs and NaNs, and there were lots of
> bugs, including compiler bugs (gcc is very far from supporting the IEC
> 60559 extension).  It turns out that the version using formula handles
> most of the special cases automatically and is much easier to fix.  I
> fixed both versions and worked around the compiler bugs, and compared
> differences to find the best version of the variant using the formula.

So which implementation do you prefer: your simpler version or
your fixed version of my first cut at ccosh?

Apologies for overlooking -0.  

> % --- s_ccosh.c~	Mon Oct  3 07:10:59 2005
> % +++ s_ccosh.c	Tue Oct  4 16:24:06 2005
> % @@ -40,11 +40,11 @@
> %  #include <math.h>
> % 
> % -#include "math_private.h"
> % +#include "../math_private.h"
> 
> Temporary hack.

I use a hack in my Makefile.  CFLAGS+=-I/usr/src/lib/msun/src. 

> %  double complex
> %  ccosh(double complex z)
> %  {
> % -	int con = 1;
> % -	double x, y;
> % +	double complex f;
> % +	double con, x, y;
> %  	int32_t hx, hy, ix, iy, lx, ly;
> %
> 
> Optimize `con' a little, and fix some style bugs (initialization in
> declaration, and disordered declarations; most other style bugs are not
> fixed or noted).  `f' is used to work around gcc bugs.

You might say that the style violations are the MUTT style.
KNF    == kernel normal format (FreeBSD style)
KNF2   == kargl normal format  (My style)
GNU    == GCC's FSF style
fdlibm == msun style

I'll eventually make it conform to KNF.

> 
> I didn't bother optimizing isnan(x), etc. using hx, etc.
>

For a committable version do you want/prefer the isnan
or the bit twiddling? 

> Handle x == 0, y finite.  There are 2 unobvious details:
> - getting the sign right when y == 0...
> 
> % +		creal(f) = cos(y);
> % +		cimag(f) = x * con;
> % +		return f;

Wow.  I'm used to programming in Fortran and would never have
written "creal(f) = ...".  This looks like your assigning a
value to a function.  (For those in the peanut gallery that
snicker at the mention of Fortran, please see the Fortran 2003
standard.)  

> - working around gcc bugs.

:-)

> Fortunately, we can construct u + I*v by assigning to the real and complex
> parts.  There should be a macro or inline function to do this.

If we go the macro route, do you want it (them?) in math_private.h
or complex.h?  If creal(f) appears on the LHS, is it a generic
reference so that type is forced to match the RHS?  Is

#define CMPLX((z),(x),(y)) {creal((z)) = (x), cimag((z)) = (y)}

acceptable?

> The (cosh(x) * tanh(x)) term should be written as sinh(x), especially now
> that cosh(x) is evaluated twice.  The old way of writing the return value:
> 
>     cosh(x) * (cos(y) + I * con * tanh(x) * sin(y))
> 
> has the interesting property of avoiding some of the compiler bugs.
> The second term is always finite, so gcc doesn't mess up Infs and Nans
> in it, and multiplying it by an infinite real coshx() does less damage
> than multiplying I by an infinite real (sinh(x) * con * siny())..

While I did not know about this property per se.  I realized that
cos(), sin(), and tanh() are all bounded functions for finite
arguments, and cosh() is the only function that could diverge
(for finite floating point arguments).

> y is Inf or Nan, so we can use the simpler formula y - y to generate a
> NaN for it.  ((y - y) / (y - y)) is only needed to generate a NaN from
> a possibly-finite y.

OK.

> %  	}
> % +	abort();
> 
> This should be unreachble, so the previous if statement must always be
> redundant.

Yes, the previous if statement is redundant.  I had filtered the
exceptional cases first in my first cut implementation, and the finite,
nonzero, x and y case were done last.  This is how fdlibm appears to
things for the float and double functions.  With the large number of
exceptional cases, I thought that we should handled those last (as an
optimization).

(2nd version deleted).

> Third, a simple test program that compares implementations of ccosh()
> on a reasonably large set of float complex args (including +0, +-Inf
> and many NaNs):

This is quite valuable, and will be stuffed away in my working
directory for future references.

Once again, thanks for nudging me in the right direction.
One of these days, I'll have some of the long double and
complex.h functions beaten into submission.

-- 
Steve