From owner-freebsd-standards@FreeBSD.ORG Mon Oct 10 11:02:14 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 73ABB16A41F for ; Mon, 10 Oct 2005 11:02:14 +0000 (GMT) (envelope-from owner-bugmaster@freebsd.org) Received: from freefall.freebsd.org (freefall.freebsd.org [216.136.204.21]) by mx1.FreeBSD.org (Postfix) with ESMTP id 3997843D46 for ; Mon, 10 Oct 2005 11:02:14 +0000 (GMT) (envelope-from owner-bugmaster@freebsd.org) Received: from freefall.freebsd.org (peter@localhost [127.0.0.1]) by freefall.freebsd.org (8.13.3/8.13.3) with ESMTP id j9AB2EhV051505 for ; Mon, 10 Oct 2005 11:02:14 GMT (envelope-from owner-bugmaster@freebsd.org) Received: (from peter@localhost) by freefall.freebsd.org (8.13.3/8.13.1/Submit) id j9AB2CXW051499 for freebsd-standards@freebsd.org; Mon, 10 Oct 2005 11:02:12 GMT (envelope-from owner-bugmaster@freebsd.org) Date: Mon, 10 Oct 2005 11:02:12 GMT Message-Id: <200510101102.j9AB2CXW051499@freefall.freebsd.org> X-Authentication-Warning: freefall.freebsd.org: peter set sender to owner-bugmaster@freebsd.org using -f From: FreeBSD bugmaster To: freebsd-standards@FreeBSD.org Cc: Subject: Current problem reports assigned to you 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 11:02:14 -0000 Current FreeBSD problem reports Critical problems Serious problems S Submitted Tracker Resp. Description ------------------------------------------------------------------------------- o [2001/03/05] bin/25542 standards /bin/sh: null char in quoted string o [2002/12/13] kern/46239 standards posix semaphore implementation errors o [2003/04/21] standards/51209standards [PATCH] add a64l()/l64a/l64a_r functions o [2003/07/12] standards/54410standards one-true-awk not POSIX compliant (no exte o [2005/06/25] standards/82654standards C99 long double math functions are missin 5 problems total. Non-critical problems S Submitted Tracker Resp. Description ------------------------------------------------------------------------------- o [2000/09/24] bin/21519 standards sys/dir.h should be deprecated some more o [2001/01/16] bin/24390 standards Replacing old dir-symlinks when using /bi s [2001/01/24] standards/24590standards timezone function not compatible witn Sin s [2001/06/18] kern/28260 standards UIO_MAXIOV needs to be made public p [2001/11/20] standards/32126standards getopt(3) not Unix-98 conformant s [2002/03/19] standards/36076standards Implementation of POSIX fuser command o [2002/06/14] standards/39256standards snprintf/vsnprintf aren't POSIX-conforman p [2002/08/12] standards/41576standards POSIX compliance of ln(1) o [2002/10/23] standards/44425standards getcwd() succeeds even if current dir has o [2002/12/09] standards/46119standards Priority problems for SCHED_OTHER using p o [2002/12/21] standards/46441standards /bin/sh does not do parameter expansion i o [2003/07/25] standards/54833standards [pcvt] more pcvt deficits o [2003/07/25] standards/54839standards [pcvt] pcvt deficits o [2003/07/31] standards/55112standards glob.h, glob_t's gl_pathc should be "size o [2003/09/05] standards/56476standards cd9660 unicode support simple hack o [2003/10/29] standards/58676standards grantpt(3) alters storage used by ptsname s [2004/02/14] standards/62858standards malloc(0) not C99 compliant o [2004/03/29] kern/64875 standards [patch] add a system call: fdatasync() o [2004/05/07] standards/66357standards make POSIX conformance problem ('sh -e' & o [2004/05/11] standards/66531standards _gettemp uses a far smaller set of filena o [2004/08/22] standards/70813standards [PATCH] ls not Posix compliant o [2004/08/26] docs/70985 standards [patch] sh(1): incomplete documentation o o [2004/09/22] standards/72006standards floating point formating in non-C locales o [2005/03/20] standards/79055standards Add an IFS regression test for shells o [2005/03/20] standards/79056standards regex(3) regression tests o [2005/03/21] standards/79067standards /bin/sh should be more intelligent about a [2005/04/23] standards/80293standards sysconf() does not support well-defined u o [2005/05/20] standards/81287standards [PATCH]: fingerd(8) might send a line not o [2005/07/21] standards/83845standards [libm] [patch] add log2() and log2f() sup o [2005/08/18] standards/85080standards output of long double subnormals (with pr o [2005/08/18] standards/85090standards [patch] add memalign() and posix_memalign 31 problems total. 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 From owner-freebsd-standards@FreeBSD.ORG Wed Oct 12 07:23:35 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 24EB716A41F for ; Wed, 12 Oct 2005 07:23:35 +0000 (GMT) (envelope-from bde@zeta.org.au) Received: from mailout2.pacific.net.au (mailout2.pacific.net.au [61.8.0.115]) by mx1.FreeBSD.org (Postfix) with ESMTP id 5B60443D45 for ; Wed, 12 Oct 2005 07:23:34 +0000 (GMT) (envelope-from bde@zeta.org.au) Received: from mailproxy1.pacific.net.au (mailproxy1.pacific.net.au [61.8.0.86]) by mailout2.pacific.net.au (8.13.4/8.13.4/Debian-3) with ESMTP id j9C7NWZp028238; Wed, 12 Oct 2005 17:23:32 +1000 Received: from katana.zip.com.au (katana.zip.com.au [61.8.7.246]) by mailproxy1.pacific.net.au (8.13.4/8.13.4/Debian-3) with ESMTP id j9C7NRbd031939; Wed, 12 Oct 2005 17:23:29 +1000 Date: Wed, 12 Oct 2005 17:23:28 +1000 (EST) From: Bruce Evans X-X-Sender: bde@delplex.bde.org To: Steve Kargl In-Reply-To: <20051010185153.GA55589@troutmask.apl.washington.edu> Message-ID: <20051012160109.I73531@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> <20051010185153.GA55589@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed 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: Wed, 12 Oct 2005 07:23:35 -0000 On Mon, 10 Oct 2005, Steve Kargl wrote: > (Attributions are a mess :-) Now accumulating them again :-). > 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 No :-). > 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. I'm beginning to prefer the unsimplified version, but it needs some more editing after the following changes. The simplified version needed another case for ccosh(+-Inf + I NaN) and that already gives it about half as many cases and the other version. It saves code mainly by falling through to the code that works for finite args, but eventually the case of finite args will need to be more special to avoid overflows and then the non-finite cases won't be able to just use it. Listing all the special cases also serves as documentation. I think the current order of special cases is not quite the best, however. Patch relative to your last version: % --- s_ccosh.c~ Wed Oct 12 05:11:13 2005 % +++ s_ccosh.c Wed Oct 12 09:05:44 2005 % @@ -43,5 +43,5 @@ % ccosh(double complex z) % { % - double con, x, y; % + double x, y; % int32_t hx, hy, ix, iy, lx, ly; % % @@ -55,25 +55,10 @@ % 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; % - } % - Delete the explicit sign flipping since it mostly just complicates things. % /* 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))); % + if (ix < 0x7ff00000 && iy < 0x7ff00000) { % + if ((iy | ly) == 0) % + return (cpack(cosh(x), x * y)); This special case was lost. Without this, cosh(x) * +-0 gives NaN if cosh(x) is +-Inf. x * y for the imaginary part is an efficient way to get +-0 with the right sign. % + return (cpack(cosh(x) * cos(y), sinh(x) * sin(y))); % + } % % /* % @@ -81,16 +66,7 @@ % * 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). % */ Delete this since it became out of date (x and y haven't changed) and getting the sign right has fine details not just here. I don't want comments in the code about all of them. % if ((ix | lx) == 0 && iy >= 0x7ff00000) % - return (cpack(y - y, copysign(0, y - y))); % + return (cpack(y - y, copysign(0, x * (y - y)))); % % /* This change is a no-op on all machines that I know of. It is to get as close as possible to what the hardware would do for sinh(x) * sinh(y). x is +-0 so sinh(x) is the same as x. sinh(y) is (y - y) (a NaN). We can't just muliply these since C99 specifies a result of +-0, so we get as close as possible by taking the sign of the multiple. Note that there is no additional invalid exception for this -- there is already one for the real part. (Actually, there probably is a second one since evaluating (y - y) only once would be an invalid optimization, but these exceptions coalesce.) Multiplying by x has no effect on all machines that I know of. Normally, if y is Inf then (y - y) is some NaN and multiplying by x preserves this, and if y is NaN then (y - y) is the same NaN and that is preseved. % @@ -99,5 +75,5 @@ % */ % if (ix < 0x7ff00000 && iy >= 0x7ff00000) % - return (cpack(y - y, y - y)); % + return (cpack(y - y, x * (y - y))); % /* % * cosh(+Inf + I 0) = +Inf + I 0 Similarly, except here we want a NaN imaginary part so we don't need a copysign(). Cases like this can currently all use the formula. % @@ -108,8 +84,8 @@ % 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))); % + return (cpack(x * x, copysign(0, x) * y)); This is the only place that deleting the sign flipping made more complicated later. The sign that was in "con" now needs to be determined. % + if (iy >= 0x7ff00000) Minor style fix (don't use "return; else"). % + return (cpack(x * x, x * (y - y))); Let the sign of x possibly affect the NaN (y - y) as usual. % + return (cpack((x * x) * cos(y), x * sin(y))); "con" just goes away (it is in y), provided we use the right sign for the x-related term here. (x * x) was logically wrong -- x is +-Inf and the term is sinh(x) = +-Inf = x. (x * x) only worked because we normalized x so it was +Inf. % } % /* % @@ -119,5 +95,5 @@ % */ % if ((iy | ly) == 0) % - return (cpack(x - x, copysign(0, x - x))); % - return (cpack((x - x) + (y - y), (x - x) + (y - y))); % + return (cpack(x * x, copysign(0, (x + x) * y))); Now x is NaN. (x - x) was logically wrong for cosh(x) -- cosh(x) is (x * x) for both Infs and NaNs -- it's a square so that it works for +-Inf, and this gives a choice of sign and value for NaNs too, depending on what the hardware does for NaN * NaN. Just use the same value here. Similarly for the imaginary part: - sinh(x) is (x + x), not (x - x), for Infs and NaNs, so as to preserve the sign for Infs. Just use the same value. - let the sign of y possibly affect the result like the sign of x did in other cases. % + return (cpack((x * x) * (y - y), (x + x) * (y - y))); % } Similarly. We want to let the sign of sin(y) possibly affect the real part. sin(y) is (y - y), so just use that. For the imaginary part, multiply instead of add to do the right thing with signs; just use sinh(x) = (x + x) not (x - x); sin(y) = (y - y) was already correct. % 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 % ... OK; could be more detailed. % 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 Patch was relative to this. % + % +/* % + * Hyperbolic cosine of a complex argument z = x + i y such that i = sqrt(-1). Other changes: don't define "i" here... % + ... % + 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; Still has lots of strange indents -- corrupt tab above, 4-chars only for most second indents. % + /* % + * 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. % + */ Not quite enough space to describe it all on 1 line; the inputs and outputs get mixed cryptically; might use more formal abbreviations. % + * cosh(+Inf + I y) = +Inf [cos(y) + I sin(y)], y != 0 and finite. Still have the strange brackets here. % 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 % ... % +float complex % +ccoshf(float complex z) % +{ % + return ((float complex) ccosh((double complex) z)); % +} I prefer to use implicit conversions. % 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. % + */ I moved my comment about the gcc bug -- the reason for existence of these functions -- to here. % +static __inline float complex % +cpackf(float __x, float __y) I needed a namespace hack to make this compile. Clients that don't use this shouldn't need to include , and this file shouldn't include it. I used "#ifdef I". I think the underscores shouldn't be used here (not even for __inline). This is not a public interface so we don't need to be very careful with the namespace. Otherwise good. These interfaces seem to work well. Here is my current "simpler" version. I plan to merge it into the other version when the other version's indentationis fixed. %%% /* * Hyperbolic cosine of a double complex argument. * * Most exceptional values are automatically correctly handled by the * standard formula. See n1124.pdf for details. */ #include #include #include "../math_private.h" double complex ccosh1(double complex z) { double x, y; x = creal(z); y = cimag(z); /* * This is subtler than it looks. We handle the cases x == 0 and * y == 0 directly not for efficiency, but to avoid multiplications * that don't work like we need. In these cases, the result must * be almost pure real or a pure imaginary, except it has type * float complex and its impure part may be -0. Multiplication of * +-0 by an infinity or a NaN would give a NaN for the impure part, * and would raise an unwanted exception. * * We depend on cos(y) and sin(y) never being precisely +-0 except * when y is +-0 to prevent getting NaNs from other cases of * +-Inf*+-0. This is true in infinite precision (since pi is * irrational), and checking shows that it is also true after * rounding to float precision. */ if (x == 0 && !isfinite(y)) return (cpack(y - y, copysign(0, x * (y - y)))); if (y == 0) return (cpack(cosh(x), isnan(x) ? copysign(0, (x + x) * y) : copysign(0, x) * y)); if (isinf(x) && !isfinite(y)) return (cpack(x * x, x * (y - y))); return (cpack(cosh(x) * cos(y), sinh(x) * sin(y))); } Here is my program for comparing complex functions. It has been cleaned up a bit (it now uses fixed-width types but not GET_FLOAT_WORD() etc. to access them properly). Compile with -DTESTD for simplistic double precision comparison. Oops, it needs to use cpack*(). %%% #include #include #include #include #include #ifndef FUNC #define FUNC ccosh #endif double complex ccosh(double complex z); double complex ccosh1(double complex z); float complex ccoshf(float complex z); float complex ccosh1f(float complex z); int main(void) { float complex vf0, vf1; float complex z; volatile float vf0i, vf0r, vf1i, vf1r; float x, y; uintmax_t na, tot_er; uint32_t eri, err, max_er, stride, vf0iu, vf0ru, vf1iu, vf1ru, xu, yu; max_er = 0; na = 0; stride = 0x80000; tot_er = 0; for (xu = 0;; xu += stride) { for (yu = 0;; yu += stride) { x = *(float *)&xu; y = *(float *)&yu; /* `z = x + I * y;', unbroken: */ crealf(z) = x; cimagf(z) = y; #ifdef TESTD /* * Hack for double precision functions. Testing them only * on float args is OK, but we should report errors in * double precision for them. */ vf0 = __CONCAT(FUNC, )(z); vf1 = __CONCAT(FUNC, 1)(z); #else vf0 = __CONCAT(FUNC, f)(z); vf1 = __CONCAT(FUNC, 1f)(z); #endif vf0r = crealf(vf0); vf0ru = *(volatile uint32_t *)&vf0r; vf0i = cimagf(vf0); vf0iu = *(volatile uint32_t *)&vf0i; vf1r = crealf(vf1); vf1ru = *(volatile uint32_t *)&vf1r; vf1i = cimagf(vf1); vf1iu = *(volatile uint32_t *)&vf1i; na++; err = abs(vf1ru - vf0ru); eri = abs(vf1iu - vf0iu); tot_er += err + eri; if (max_er < err) max_er = err; if (max_er < eri) max_er = eri; if (err > 4 || eri > 4) { printf("z = %#10x+I*%-#10x %.15g+I*%-.15g\n", xu, yu, x, y); printf( __XSTRING(FUNC) "f(z) = %#10x+I*%-#10x %.15g+I*%-.15g\n", vf0ru, vf0iu, vf0r, vf0i); printf( __XSTRING(FUNC) "1f(z) = %#10x+I*%-#10x %.15g+I*%-.15g\n", vf1ru, vf1iu, vf1r, vf1i); printf("err = %s%#x+I*%s.%#x\n", vf0ru <= vf1ru ? "+" : "-", err, vf0iu <= vf1iu ? "+" : "-", eri); fflush(stdout); } if (na % 10000000 == 1) { printf( "%ju: %#x+I*%#x: %g+I*%g: max_error = %#x, avg_error = %.3f\n", na, xu, yu, x, y, max_er, (double)tot_er / na); fflush(stdout); } if (yu == -stride) break; } if (xu == -stride) break; } printf("%ju cases: max_error = %#x, avg_error = %.3f\n", na, max_er, (double)tot_er / na); } #include "s_ccosh.c" #include "s_ccosh1.c" #include "s_ccoshf.c" #include "s_ccosh1f.c" /* For ccosh(), must avoid using the broken hardware cos() and sin(): */ #undef rcsid #define rcsid rcsid6 #include "../s_cos.c" #undef rcsid #define rcsid rcsid7 #include "../s_sin.c" #ifdef DEBUG #undef rcsid #define rcsid rcsid0 #include "../e_coshf.c" #undef one #define one one0 #undef rcsid #define rcsid rcsid1 #include "../e_sinhf.c" #undef one #define one one1 #undef rcsid #define rcsid rcsid2 #include "../s_tanhf.c" #endif %%% Bruce