complex.h math functions

[Steve Kargl]

On Sat, Oct 01, 2005 at 09:38:44PM +1000, Bruce Evans wrote:
> >A slightly better implementation than
> >the above would be
> >
> >       if (y == 0.) return (cos(x));
> >       if (x == 0.) return (cosh(y));
> >       return (cosh(y) * (cos(x) - I * sin(x) * tanh(y)));
> 
> I'm not sure if this is better.  The special cases may be so rare that
> checking for them slows down the average case.  I think special handling
> is unnecessary since cos(0) is exactly 1.0, etc, so the checks are just
> an optimization of the special args.  OTOH, the type-generic version
> should have special checks using something like
> "if (__builtin_const_p(y) && y == 0.)" to classify the pure real case,
> etc.  Inline function args are never constant, so the type-generic
> version must be implemented as a macro.

Type-generic stuff is way above my C competency level.  I'll
leave that to bde, das, of stefanf. :-)

> 
> >...
> >OK, I'll use the implementation details of ccosh to infer the
> >handling of exceptional values.
> >

At the risk of embarassing myself once again on freebsd-standards,
here's an implementation of ccosh.  If this is acceptable, I'll
proceed with implementations for csinh, ctanh, ccos, csin, and
ctan.

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))

-- 
Steve

/*-
 * Copyright (c) 2005, Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice unmodified, this list of conditions, and the following
 *    disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

/*
 *  Hyperbolic cosine of a complex argument z = x + i y such that i = sqrt(-1).
 *
 *  cosh(z) = cosh(x+iy)
 *          = cosh(x)                              if y == 0.
 *          = cos(y)                               if x == 0.
 *          = cosh(x) [cos(y) + i tanh(x) sin(y)]  otherwise.
 *
 *  Exception values are noted in the comments within the source code.
 *  These values and the return result were taken from n869.pdf.
 */

#include <complex.h>
#include <math.h>

#include "math_private.h"

double complex
ccosh(double complex z)
{
	int con = 1;
	double x, y;
	int32_t hx, hy, ix, iy, lx, ly;

	x = creal(z);
	y = cimag(z);

	if (x < 0.) {
	    /*  Even function: ccosh(-z) = ccosh(z). */
	    x = -x;
	    y = -y;
	    /* Required identity: ccosh(conj(z)) = conj(ccosh(z)). */
	    if (y < 0.) {
		con = -1;
		y = -y;
	    }
	}

	EXTRACT_WORDS(hx, lx, x);
	EXTRACT_WORDS(hy, ly, y);

	ix = 0x7fffffff&hx;  /* if ix >= 0x7ff00000, then inf or NaN */
        iy = 0x7fffffff&hy;  /* if iy >= 0x7ff00000, then inf or NaN */

	/* Filter out the non-exceptional cases.  x and y are finite. */
	if (ix < 0x7ff00000 && iy < 0x7ff00000) {
	   /*  cosh(+0 + I 0) = 1.  */
	   if ((ix|lx) == 0 && (iy|ly) == 0)
		return 1.;
	   return (cosh(x) * (cos(y) + I * con * tanh(x) * sin(y)));
	}
 	/*
	 *  cosh(+0 + I Inf) = NaN +- I 0, sign of 0 is unspecified,
	 *				   invalid exception.
	 *  cosh(+0 + I NaN) = NaN +- I 0, sign of 0 is unspecified.
	 */
	if ((ix|lx) == 0 && iy >= 0x7ff00000)
		return (y - y) / (y - y);
	/*
	 *  cosh(x + I Inf) = NaN + I NaN, finite x != 0, invalid exception.
	 *  cosh(x + I NaN) = NaN + I NaN, finite x != 0, opt. inval. except.
	 */
	if (ix < 0x7ff00000 && iy >= 0x7ff00000)
	    return (x - x) / (x - x) + I * (x - x) / (x - x);
	/*
	 *  cosh(+Inf + I 0)   = +Inf + I 0
	 *  cosh(+Inf + I Inf) = +Inf + I NaN, invalid exception.
	 *  cosh(+Inf + I y)   = +Inf [cos(y) + I sin(y)], y != 0 and finite.
	 *  cosh(+Inf + I NaN) = +Inf + I NaN.
	 */
	if (ix >= 0x7ff00000 && ((hx&0xfffff)|lx) == 0) {
	    if ((iy|ly) == 0)
		return (x * x);
	    else if (iy >= 0x7ff00000)
		return (x * x) + I * (y - y) / (y - y);
	    return (x * x) * (cos(y) + I * con * sin(y));
	}
	/*  
	 *  cosh(NaN + I 0)   = NaN +- I 0, sign of 0 unspecified.
	 *  cosh(NaN + I y)   = NaN + I NaN, nonzero y, opt. inval. except.
	 *  cosh(NaN + I NaN) = NaN + I NaN.
	 */
	if (ix >= 0x7ff00000 && ((hx&0xfffff)|lx) != 0) {
	    if ((iy|ly) == 0)
		return x * x;
	    return x * x + I * x * x;
	}
}