Date: Sun, 2 Oct 2005 12:18:41 -0700 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: Bruce Evans <bde@zeta.org.au> Cc: freebsd-standards@freebsd.org Subject: Re: complex.h math functions Message-ID: <20051002191841.GA40568@troutmask.apl.washington.edu> In-Reply-To: <20051001204005.N37704@delplex.bde.org> References: <20050929195552.GA14982@troutmask.apl.washington.edu> <20050930221846.T34595@delplex.bde.org> <20050930152734.GA20655@troutmask.apl.washington.edu> <20051001204005.N37704@delplex.bde.org>
next in thread | previous in thread | raw e-mail | index | archive | help
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;
}
}
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20051002191841.GA40568>