From owner-freebsd-standards@FreeBSD.ORG Mon Oct 10 18:51:58 2005 Return-Path: X-Original-To: freebsd-standards@freebsd.org Delivered-To: freebsd-standards@freebsd.org Received: from mx1.FreeBSD.org (mx1.freebsd.org [216.136.204.125]) by hub.freebsd.org (Postfix) with ESMTP id BA92B16A41F for ; Mon, 10 Oct 2005 18:51:58 +0000 (GMT) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.208.78.105]) by mx1.FreeBSD.org (Postfix) with ESMTP id 4C83443D45 for ; Mon, 10 Oct 2005 18:51:58 +0000 (GMT) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (localhost [127.0.0.1]) by troutmask.apl.washington.edu (8.13.4/8.13.4) with ESMTP id j9AIpv4S077138; Mon, 10 Oct 2005 11:51:57 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.13.4/8.13.1/Submit) id j9AIprs2077137; Mon, 10 Oct 2005 11:51:53 -0700 (PDT) (envelope-from sgk) Date: Mon, 10 Oct 2005 11:51:53 -0700 From: Steve Kargl To: Bruce Evans Message-ID: <20051010185153.GA55589@troutmask.apl.washington.edu> 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> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20051008052850.S59139@delplex.bde.org> User-Agent: Mutt/1.4.2.1i Cc: freebsd-standards@freebsd.org Subject: Re: complex.h math functions X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Mon, 10 Oct 2005 18:51:58 -0000 (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 +#include + +#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 + +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