Date: Mon, 10 Oct 2005 11:51:53 -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: <20051010185153.GA55589@troutmask.apl.washington.edu> In-Reply-To: <20051008052850.S59139@delplex.bde.org> References: <20050930152734.GA20655@troutmask.apl.washington.edu> <20051001204005.N37704@delplex.bde.org> <20051002191841.GA40568@troutmask.apl.washington.edu> <20051004164618.U47262@delplex.bde.org> <20051005032400.GA6736@troutmask.apl.washington.edu> <20051005191936.N51036@delplex.bde.org> <20051006212602.GA40609@troutmask.apl.washington.edu> <20051007231439.F58005@delplex.bde.org> <20051007190239.GA78674@troutmask.apl.washington.edu> <20051008052850.S59139@delplex.bde.org>
next in thread | previous in thread | raw e-mail | index | archive | help
(Attributions are a mess :-)
> >>>>>>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;
> >
> >I think the above is wrong. We need "x * con * sin(y)" because
> >sin(y) is between -1 and 1. For example we may have
> >
> >x * con * sin(y) = (-0) * (1) * (-0.5) = +0
> >
> >but I need to check n1124.pdf to see what the behavior of
> >signed zero arithmetic does.
>
> Right. I found this bug a few days ago and mentioned it in previous
> mail. I have only fixed it for the float versions.
I fixed this in my version of ccosh and used cpack(,) for
constructing the complex values. I tried to use cpack(,)
in your simplified version, but screwed something up in
that the return value is the complex conjugate from that
version.
I've attached the latest version. Hopefully, I caught the
rest of the style(9) problems. The cosh.3 man page has been
updated. I did not hook s_ccosh.c up to the Makefile
because I wasn't sure were to put it in msun/src. If you have
no further comments, can you commit this version (or your
simplified version)? I'll look at the other trigometric and
hyperbolic functions when we converge on s_ccosh.c.
--
Steve
diff -ruN /usr/src/lib/msun/man/cosh.3 freebsd/src/lib/msun/man/cosh.3
--- /usr/src/lib/msun/man/cosh.3 Fri Jan 14 15:28:28 2005
+++ freebsd/src/lib/msun/man/cosh.3 Sat Oct 8 11:46:29 2005
@@ -32,10 +32,12 @@
.\" from: @(#)cosh.3 5.1 (Berkeley) 5/2/91
.\" $FreeBSD: src/lib/msun/man/cosh.3,v 1.12 2005/01/14 23:28:28 das Exp $
.\"
-.Dd January 14, 2005
+.Dd October 8, 2005
.Dt COSH 3
.Os
.Sh NAME
+.Nm ccosh ,
+.Nm ccoshf ,
.Nm cosh ,
.Nm coshf
.Nd hyperbolic cosine functions
@@ -43,6 +45,10 @@
.Lb libm
.Sh SYNOPSIS
.In math.h
+.Ft "double complex"
+.Fn ccosh "double complex x"
+.Ft "float complex"
+.Fn ccosh "float complex x"
.Ft double
.Fn cosh "double x"
.Ft float
diff -ruN /usr/src/lib/msun/src/complex/s_ccosh.c freebsd/src/lib/msun/src/complex/s_ccosh.c
--- /usr/src/lib/msun/src/complex/s_ccosh.c Wed Dec 31 16:00:00 1969
+++ freebsd/src/lib/msun/src/complex/s_ccosh.c Mon Oct 10 11:10:14 2005
@@ -0,0 +1,123 @@
+/*-
+ * Copyright (c) 2005, Bruce D. Evans and 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) cos(y) + i sinh(x) sin(y).
+ *
+ * Exception values are noted in the comments within the source code.
+ * These values and the return value were taken from n1124.pdf.
+ */
+
+#include <complex.h>
+#include <math.h>
+
+#include "../math_private.h"
+
+double complex
+ccosh(double complex z)
+{
+ double con, x, y;
+ int32_t hx, hy, ix, iy, lx, ly;
+
+ x = creal(z);
+ y = cimag(z);
+
+ 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. */
+
+ /* Even function: ccosh(-z) = ccosh(z). */
+ if (hx & 0x80000000 && !isnan(x)) {
+ x = -x;
+ hx = ix;
+ if (!isnan(y)) {
+ y = -y;
+ hy ^= 0x80000000;
+ }
+ }
+
+ /* Required identity: ccosh(conj(z)) = conj(ccosh(z)). */
+ con = 1;
+ if (hy & 0x80000000 && !isnan(y)) {
+ con = -1;
+ y = -y;
+ hy = iy;
+ }
+
+ /* Handle the nearly-non-exceptional cases where x and y are finite. */
+ if (ix < 0x7ff00000 && iy < 0x7ff00000)
+ return (cpack(cosh(x) * cos(y), con * sinh(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.
+ *
+ * Use the sign of (sinh(original_x) * sin(original_y)) for
+ * the optional sign. This expression reduces to (y - y) for
+ * all cases checked (mainly i386's). The fdlibm sin(y) is
+ * certainly (y - y), but in the +-Inf cases we may have
+ * flipped the sign of original_y to get y, and we want to
+ * multiply by the sinh() term. These complications have no
+ * effect with normal NaN semantics (-Inf - -Inf) = same NaN
+ * as (Inf - Inf), and (finite * NaN) = same NaN).
+ */
+ if ((ix | lx) == 0 && iy >= 0x7ff00000)
+ return (cpack(y - y, copysign(0, 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 (cpack(y - y, y - y));
+ /*
+ * 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 (cpack(x * x, con * y));
+ else if (iy >= 0x7ff00000)
+ return (cpack(x * x, y - y));
+ return (cpack((x * x) * cos(y), (x * x) * 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 ((iy | ly) == 0)
+ return (cpack(x - x, copysign(0, x - x)));
+ return (cpack((x - x) + (y - y), (x - x) + (y - y)));
+}
diff -ruN /usr/src/lib/msun/src/complex/s_ccoshf.c freebsd/src/lib/msun/src/complex/s_ccoshf.c
--- /usr/src/lib/msun/src/complex/s_ccoshf.c Wed Dec 31 16:00:00 1969
+++ freebsd/src/lib/msun/src/complex/s_ccoshf.c Fri Oct 7 17:30:36 2005
@@ -0,0 +1,39 @@
+/*-
+ * 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).
+ */
+
+#include <complex.h>
+
+double complex ccosh(double complex);
+
+float complex
+ccoshf(float complex z)
+{
+ return ((float complex) ccosh((double complex) z));
+}
diff -ruN /usr/src/lib/msun/src/math_private.h freebsd/src/lib/msun/src/math_private.h
--- /usr/src/lib/msun/src/math_private.h Fri Feb 4 12:05:39 2005
+++ freebsd/src/lib/msun/src/math_private.h Fri Oct 7 17:30:40 2005
@@ -155,6 +155,36 @@
} while (0)
/*
+ * Inline functions that can be used to construct complex values.
+ */
+static __inline float complex
+cpackf(float __x, float __y)
+{
+ float complex __z;
+ __real__ __z = __x;
+ __imag__ __z = __y;
+ return (__z);
+}
+
+static __inline double complex
+cpack(double __x, double __y)
+{
+ double complex __z;
+ __real__ __z = __x;
+ __imag__ __z = __y;
+ return (__z);
+}
+
+static __inline long double complex
+cpackl(long double __x, long double __y)
+{
+ long double complex __z;
+ __real__ __z = __x;
+ __imag__ __z = __y;
+ return (__z);
+}
+
+/*
* ieee style elementary functions
*
* We rename functions here to improve other sources' diffability
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20051010185153.GA55589>