From owner-freebsd-standards@FreeBSD.ORG Sun Oct 2 19:19:02 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 E81FA16A420 for ; Sun, 2 Oct 2005 19:19:02 +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 39D8E43D5E for ; Sun, 2 Oct 2005 19:18:49 +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 j92JIl3H040651; Sun, 2 Oct 2005 12:18:47 -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 j92JIfen040650; Sun, 2 Oct 2005 12:18:41 -0700 (PDT) (envelope-from sgk) Date: Sun, 2 Oct 2005 12:18:41 -0700 From: Steve Kargl To: Bruce Evans Message-ID: <20051002191841.GA40568@troutmask.apl.washington.edu> References: <20050929195552.GA14982@troutmask.apl.washington.edu> <20050930221846.T34595@delplex.bde.org> <20050930152734.GA20655@troutmask.apl.washington.edu> <20051001204005.N37704@delplex.bde.org> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20051001204005.N37704@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: Sun, 02 Oct 2005 19:19:03 -0000 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 #include #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; } } From owner-freebsd-standards@FreeBSD.ORG Mon Oct 3 11:02:47 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 B28CA16A421 for ; Mon, 3 Oct 2005 11:02:47 +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 445DC43D58 for ; Mon, 3 Oct 2005 11:02:24 +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 j93B2NBk066428 for ; Mon, 3 Oct 2005 11:02:23 GMT (envelope-from owner-bugmaster@freebsd.org) Received: (from peter@localhost) by freefall.freebsd.org (8.13.3/8.13.1/Submit) id j93B2Mhd066422 for freebsd-standards@freebsd.org; Mon, 3 Oct 2005 11:02:22 GMT (envelope-from owner-bugmaster@freebsd.org) Date: Mon, 3 Oct 2005 11:02:22 GMT Message-Id: <200510031102.j93B2Mhd066422@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, 03 Oct 2005 11:02:48 -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 p [2004/02/21] standards/63173standards Patch to add getopt_long_only(3) to libc 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 32 problems total. From owner-freebsd-standards@FreeBSD.ORG Tue Oct 4 09:31:49 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 A6DA916A41F for ; Tue, 4 Oct 2005 09:31:49 +0000 (GMT) (envelope-from bde@zeta.org.au) Received: from mailout1.pacific.net.au (mailout1.pacific.net.au [61.8.0.84]) by mx1.FreeBSD.org (Postfix) with ESMTP id E7E0743D48 for ; Tue, 4 Oct 2005 09:31:48 +0000 (GMT) (envelope-from bde@zeta.org.au) Received: from mailproxy2.pacific.net.au (mailproxy2.pacific.net.au [61.8.0.87]) by mailout1.pacific.net.au (8.13.4/8.13.4/Debian-3) with ESMTP id j949Vldb020422; Tue, 4 Oct 2005 19:31:47 +1000 Received: from katana.zip.com.au (katana.zip.com.au [61.8.7.246]) by mailproxy2.pacific.net.au (8.13.4/8.13.4/Debian-3) with ESMTP id j949Vhl0012775; Tue, 4 Oct 2005 19:31:43 +1000 Date: Tue, 4 Oct 2005 19:31:43 +1000 (EST) From: Bruce Evans X-X-Sender: bde@delplex.bde.org To: Steve Kargl In-Reply-To: <20051002191841.GA40568@troutmask.apl.washington.edu> Message-ID: <20051004164618.U47262@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> <20051002191841.GA40568@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: Tue, 04 Oct 2005 09:31:49 -0000 On Sun, 2 Oct 2005, Steve Kargl wrote: >>> ... >>> 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)) A program is needed, since (ccoshf)(z) and &(ccoshf) must work even if ccoshf() is a macro. I converted to ccoshf() for testing it and comparing with a simpler implementation that uses the formula, then converted back. There are lots of special cases for -0, Infs and NaNs, and there were lots of bugs, including compiler bugs (gcc is very far from supporting the IEC 60559 extension). It turns out that the version using formula handles most of the special cases automatically and is much easier to fix. I fixed both versions and worked around the compiler bugs, and compared differences to find the best version of the variant using the formula. First, fixes for your version: % --- s_ccosh.c~ Mon Oct 3 07:10:59 2005 % +++ s_ccosh.c Tue Oct 4 16:24:06 2005 % @@ -40,11 +40,11 @@ % #include % % -#include "math_private.h" % +#include "../math_private.h" Temporary hack. % % double complex % ccosh(double complex z) % { % - int con = 1; % - double x, y; % + double complex f; % + double con, x, y; % int32_t hx, hy, ix, iy, lx, ly; % Optimize `con' a little, and fix some style bugs (initialization in declaration, and disordered declarations; most other style bugs are not fixed or noted). `f' is used to work around gcc bugs. % @@ -52,15 +52,4 @@ % 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; % - } % - } % - Move this to after extracting words so that we can look at the sign bit in ix and iy to classify -0. Note that n869.pdf only describes the behaviour explicitly for +0 and +Inf so that it doesn't have to describe many more cases in detail, but ut is clear that the evenness and conjugacy rules to the sign bit in +-0 although this isn't stated explicitly at least near the spec for ccosh(). Evenness is stated explicitly for cos(+-0) and oddness is stated explicitly for sin(+-0) (e.g., sin(-0) is specified to be -0, and this is implemented by fdlibm and i387 hardware sin). Although changing signs so that all the sign bits are 1 simplifies n869.pdf, it mostly complicates things for us since we can depend on multiplication to do the right things with signs. We need complications to avoid flipping the sign bit for NaNs so as to set the sign bits in the correct MD way for the cases where the standard doesn't specify the sign. Most operations on NaNs don't change the sign bit on most (?) machines, but -NaN changes it on i386's at least (0-NaN is different from -NaN on i386's!). % EXTRACT_WORDS(hx, lx, x); % EXTRACT_WORDS(hy, ly, y); % @@ -69,10 +58,37 @@ % iy = 0x7fffffff&hy; /* if iy >= 0x7ff00000, then inf or NaN */ % % - /* Filter out the non-exceptional cases. x and y are finite. */ This was changed to match the code, but didn't really move. % + /* 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; % + } The symmetry for conjugacy is independent of the one for evenness, so flipping the sign for it must be moved out of the if for flipping the sign(s) for evenness. Avoid flipping the sign bit for NaNs. I didn't bother optimizing isnan(x), etc. using hx, etc. % + % + /* Handle the nearly-non-exceptional cases where x and y are finite. */ Change comment to match code. % if (ix < 0x7ff00000 && iy < 0x7ff00000) { % - /* cosh(+0 + I 0) = 1. */ % - if ((ix|lx) == 0 && (iy|ly) == 0) % - return 1.; Adjust this to handle all cases where x is 0 (and y is finite). This fixes handling of cosh(+-0 + I*+-0) and handles nonzero y at negative extra cost. The result of cosh(+-0 + I*+-0) must be (1. + I*(+-0)*(+-0)). We have forgotten the original x and y, else we could return (1. + I*x*y) after working around the compiler bugs. (x * con) would give a correctly-signed 0, but it is simpler to use this in new code that handles all finite y (and zero x). % - return (cosh(x) * (cos(y) + I * con * tanh(x) * sin(y))); This changes but doesn't really move. % + if ((ix | lx) == 0) { 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; % + return f; ... - working around gcc bugs. gcc's support for IEC 60559 is so noexistent that u + I*v is unusable in general to construct a double complex from 2 doubles. Here the main bug is that I*-0 gets corrupted to I*+0. gcc generates code that does extra work to implement this bug. gcc's handling of full complex operations is more reasonable -- it generates code that doesn't do the extra, very expensive work to avoid related bugs. See n869.pdf section G.4.1 for how difficult it is to implement complex multiplication perfectly correctly. The example there handles Infs and and NaNs and some overflows, but not underflows or signed 0's. The complications are related to the ones here and to the ones that I meantioned for x*x. E.g., the real part of (a+ib)(c+id) is ac-bd on real numbers, but for finite floating point numbers the multiplications can overflow without the their difference overflowing, and their difference can underflow; Infs, NaNs and signed 0's give additional complications. Fortunately, we can construct u + I*v by assigning to the real and complex parts. There should be a macro or inline function to do this. % + } % + if ((iy | ly) == 0) { % + creal(f) = cosh(x); % + cimag(f) = con * y; % + return f; % + } This special case for x finite and y zero is much more critical. cosh(x) can be infinite for finite x, and then we must not perform or use the multiplication (cosh(x) * tanh(x) * con * sin(y)), since that would be (Inf * +-0) so it would raise the invalid exception and give result NaN, but we want a result of +-0 for the imaginary part (just like for x infinite and y 0). We work around the compiler bug as usual. % + creal(f) = cosh(x) * cos(y); % + cimag(f) = cosh(x) * tanh(x) * con * sin(y); % + return f; % } % /* The only change for finite nonzero x and y is to avoid the compiler bug. When sinh(x) is Inf and y is nonzero, it is correct for the imaginary part to be NaN. The (cosh(x) * tanh(x)) term should be written as sinh(x), especially now that cosh(x) is evaluated twice. The old way of writing the return value: cosh(x) * (cos(y) + I * con * tanh(x) * sin(y)) has the interesting property of avoiding some of the compiler bugs. The second term is always finite, so gcc doesn't mess up Infs and Nans in it, and multiplying it by an infinite real coshx() does less damage than multiplying I by an infinite real (sinh(x) * con * siny()).. Using cosh(x) * tanh(x) instead of sinh(x) is also wrong because in theory cosh(x) might be Inf when sinh(x) is finite. I think this doesn't happen in practice. Similarly, sinh(x) might be infinite while the infinite precision sinh(x) * sin(y) is finite -- this is a subcase of the problem of multiplying 2 double complexes perfectly correctly. Regression tests showed that ccoshf() implemented using float precision differed significantly from ccosh() implemented using double precision only because the double precision version mostly overflow here (except when the result overflows a float). % @@ -81,6 +97,19 @@ % * cosh(+0 + I NaN) = NaN +- I 0, sign of 0 is unspecified. % */ % - if ((ix|lx) == 0 && iy >= 0x7ff00000) % - return (y - y) / (y - y); % + if ((ix|lx) == 0 && iy >= 0x7ff00000) { % + creal(f) = y - y; Work around compiler bugs as usual. y is Inf or Nan, so we can use the simpler formula y - y to generate a NaN for it. ((y - y) / (y - y)) is only needed to generate a NaN from a possibly-finite y. % + /* % + * 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). % + */ % + cimag(f) = copysign(0, y - y); % + return f; % + } % /* % * cosh(x + I Inf) = NaN + I NaN, finite x != 0, invalid exception. Intricacies to set the sign to the least surprising/most reproducible value. % @@ -88,5 +117,5 @@ % */ % if (ix < 0x7ff00000 && iy >= 0x7ff00000) % - return (x - x) / (x - x) + I * (x - x) / (x - x); % + return y - y + I * (y - y); % /* % * cosh(+Inf + I 0) = +Inf + I 0 Use the right arg in the return value so that NaNs get propagated to the result. Here x is finite (nonzero) and y is +-Inf or NaN. For y = +-Inf, y - y gives a standard fixed NaN "real indefinite" on i386's at least, the same one as the expression involving x, but for y = NaN we want to return the same NaN and not change to "real indefinite". Use a simpler expression for the resulting NaN, as above. Testing showed that working around the compiler bugs is not needed here or in some other places. Here this is because everything ends up as NaN, so the spurious NaNs introduced by gcc have low enough precedence and/or differences that they don't affect the result. % @@ -96,8 +125,10 @@ % */ % 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); % + if ((iy|ly) == 0) { % + creal(f) = x * x; % + cimag(f) = con * y; % + return f; The usual workaround, plus fix the sign for the imaginary part. % + } else if (iy >= 0x7ff00000) % + return x * x + I * (y - y); Use fewer parentheses for the real part and a simpler formula for the imaginary part. % return (x * x) * (cos(y) + I * con * sin(y)); % } % @@ -108,7 +139,11 @@ % */ % if (ix >= 0x7ff00000 && ((hx&0xfffff)|lx) != 0) { % - if ((iy|ly) == 0) % - return x * x; % - return x * x + I * x * x; % + if ((iy|ly) == 0) { % + creal(f) = x - x; % + cimag(f) = copysign(0, x - x); % + return f; The usual workaround, plus fix the sign for the imaginary part. x is NaN, so `con' doesn't apply. % + } % + return (x - x) + (y - y) + I * ((x - x) + (y - y)); Here x and y are both NaNs. Mix both of them into the output. The parentheses are essential, perhaps due to extra precision in intermediate values. On i386's, this results in the largest (with the NaNs considered to be (signed?) 64-bit integers) being propagated. Propagation of NaNs doesn't always preserve them. On i386's the only changes is to convert signaling NaNs to quiet NaNs by setting the quietness bit. We use the formula (y - y) extensively to convert infinities to NaNs and get whatever MD conversions the hardware wants to do to propagate NaNs. % } % + abort(); This should be unreachble, so the previous if statement must always be redundant. % } --- Second, a simple version: % /* % * Hyperbolic cosine of a complex argument. % * % * Most exceptional values are automatically correctly handled by the % * standard formula. See n689.pdf for details. % */ % % #include % #include % % double complex % ccosh1(double complex z) % { % double complex f; % 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. The result must be almost pure % * real or a pure imaginary, except it has type double 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. % * % * In theory, any of sinh(x), cos(y) and sin(y) may be +-0 for a % * nonzero arg, because although the infinitely precise value cannot % * even be rational (Gelfond-Schneider theorem?), the rounded result % * may be 0. For such values, we would have to avoid multiplying % * +-0 by an infinity or a NaN, and cos() and sin() would have to % * be careful to get the sign right. Exhaustive checking shows that % * this never happens for args that are representable as floats. % * % * Note: the expression u + I*v is unusable since gcc introduces % * many overflow, underflow, sign and efficiency bugs by rewriting % * I*v as (0+I)*(v+0*I) and laboriously computing the full complex % * product. In particular, I*Inf is corrupted to NaN+I*Inf, and % * I*-0 is corrupted to 0+I*0. % */ % if (x == 0) { % creal(f) = cos(y); % cimag(f) = isfinite(y) ? x * copysign(1, y) : % copysign(0, y - y); % return (f); % } % if (y == 0) { % creal(f) = cosh(x); % cimag(f) = !isnan(x) ? copysign(1, x) * y : % copysign(0, x - x); % return (f); % } % % creal(f) = cosh(x) * cos(y); % cimag(f) = sinh(x) * sin(y); % return (f); % } I hope this is described in too much detail in its main comment. Third, a simple test program that compares implementations of ccosh() on a reasonably large set of float complex args (including +0, +-Inf and many NaNs): % #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; % vf0 = FUNC(z); % vf1 = __CONCAT(FUNC, 1)(z); % 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 > 2 || eri > 2) { % printf("z = %#8x+I%-#8x %.15g+I%-.15g\n", % xu, yu, x, y); % printf( % __XSTRING(FUNC) "(z) = %#8x+I%-#8x %.15g+I%-.15g\n", % vf0ru, vf0iu, vf0r, vf0i); % printf( % __XSTRING(FUNC) "1(z) = %#8x+I%-#8x %.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 == 0) { % 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" #ifdef nothere % #include "s_ccoshf.c" % #include "s_ccosh1f.c" #endif % % /* 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 This must be built in a subdir of lib/msun/src. I used this mainly with ccosh[1]f() and hacked it to work with ccosh[1](). float complex arg space is too large to check exhaustively, and double complex arg space is much larger and harder to step through and display in the error messages, so I kept using float complex arg space. The coverage with stride 0x80000 should be more than enough for NaNs. Testing with stride 0x80000 showed a maximum error of 2 ulps (for the real and imaginary parts of the result) from using coshf(x) * tanhf(x) instead of sinhf(x) if the non-complex float functions are compiled with -O, but the error is 3-4 ulps if they are compiled with -O0. Bruce From owner-freebsd-standards@FreeBSD.ORG Wed Oct 5 03:24:10 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 C152116A41F for ; Wed, 5 Oct 2005 03:24:10 +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 6827343D45 for ; Wed, 5 Oct 2005 03:24:10 +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 j953O64t006902; Tue, 4 Oct 2005 20:24:06 -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 j953O0Lp006901; Tue, 4 Oct 2005 20:24:00 -0700 (PDT) (envelope-from sgk) Date: Tue, 4 Oct 2005 20:24:00 -0700 From: Steve Kargl To: Bruce Evans Message-ID: <20051005032400.GA6736@troutmask.apl.washington.edu> References: <20050929195552.GA14982@troutmask.apl.washington.edu> <20050930221846.T34595@delplex.bde.org> <20050930152734.GA20655@troutmask.apl.washington.edu> <20051001204005.N37704@delplex.bde.org> <20051002191841.GA40568@troutmask.apl.washington.edu> <20051004164618.U47262@delplex.bde.org> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20051004164618.U47262@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: Wed, 05 Oct 2005 03:24:10 -0000 On Tue, Oct 04, 2005 at 07:31:43PM +1000, Bruce Evans wrote: > On Sun, 2 Oct 2005, Steve Kargl wrote: > > >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)) > > A program is needed, since (ccoshf)(z) and &(ccoshf) must work even if > ccoshf() is a macro. OK. > I converted to ccoshf() for testing it and comparing with a simpler > implementation that uses the formula, then converted back. There are > lots of special cases for -0, Infs and NaNs, and there were lots of > bugs, including compiler bugs (gcc is very far from supporting the IEC > 60559 extension). It turns out that the version using formula handles > most of the special cases automatically and is much easier to fix. I > fixed both versions and worked around the compiler bugs, and compared > differences to find the best version of the variant using the formula. So which implementation do you prefer: your simpler version or your fixed version of my first cut at ccosh? Apologies for overlooking -0. > % --- s_ccosh.c~ Mon Oct 3 07:10:59 2005 > % +++ s_ccosh.c Tue Oct 4 16:24:06 2005 > % @@ -40,11 +40,11 @@ > % #include > % > % -#include "math_private.h" > % +#include "../math_private.h" > > Temporary hack. I use a hack in my Makefile. CFLAGS+=-I/usr/src/lib/msun/src. > % double complex > % ccosh(double complex z) > % { > % - int con = 1; > % - double x, y; > % + double complex f; > % + double con, x, y; > % int32_t hx, hy, ix, iy, lx, ly; > % > > Optimize `con' a little, and fix some style bugs (initialization in > declaration, and disordered declarations; most other style bugs are not > fixed or noted). `f' is used to work around gcc bugs. You might say that the style violations are the MUTT style. KNF == kernel normal format (FreeBSD style) KNF2 == kargl normal format (My style) GNU == GCC's FSF style fdlibm == msun style I'll eventually make it conform to KNF. > > I didn't bother optimizing isnan(x), etc. using hx, etc. > For a committable version do you want/prefer the isnan or the bit twiddling? > 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; > % + return f; Wow. I'm used to programming in Fortran and would never have written "creal(f) = ...". This looks like your assigning a value to a function. (For those in the peanut gallery that snicker at the mention of Fortran, please see the Fortran 2003 standard.) > - working around gcc bugs. :-) > Fortunately, we can construct u + I*v by assigning to the real and complex > parts. There should be a macro or inline function to do this. If we go the macro route, do you want it (them?) in math_private.h or complex.h? If creal(f) appears on the LHS, is it a generic reference so that type is forced to match the RHS? Is #define CMPLX((z),(x),(y)) {creal((z)) = (x), cimag((z)) = (y)} acceptable? > The (cosh(x) * tanh(x)) term should be written as sinh(x), especially now > that cosh(x) is evaluated twice. The old way of writing the return value: > > cosh(x) * (cos(y) + I * con * tanh(x) * sin(y)) > > has the interesting property of avoiding some of the compiler bugs. > The second term is always finite, so gcc doesn't mess up Infs and Nans > in it, and multiplying it by an infinite real coshx() does less damage > than multiplying I by an infinite real (sinh(x) * con * siny()).. While I did not know about this property per se. I realized that cos(), sin(), and tanh() are all bounded functions for finite arguments, and cosh() is the only function that could diverge (for finite floating point arguments). > y is Inf or Nan, so we can use the simpler formula y - y to generate a > NaN for it. ((y - y) / (y - y)) is only needed to generate a NaN from > a possibly-finite y. OK. > % } > % + abort(); > > This should be unreachble, so the previous if statement must always be > redundant. Yes, the previous if statement is redundant. I had filtered the exceptional cases first in my first cut implementation, and the finite, nonzero, x and y case were done last. This is how fdlibm appears to things for the float and double functions. With the large number of exceptional cases, I thought that we should handled those last (as an optimization). (2nd version deleted). > Third, a simple test program that compares implementations of ccosh() > on a reasonably large set of float complex args (including +0, +-Inf > and many NaNs): This is quite valuable, and will be stuffed away in my working directory for future references. Once again, thanks for nudging me in the right direction. One of these days, I'll have some of the long double and complex.h functions beaten into submission. -- Steve From owner-freebsd-standards@FreeBSD.ORG Wed Oct 5 10:16:08 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 160AB16A41F for ; Wed, 5 Oct 2005 10:16:08 +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 6B7F943D46 for ; Wed, 5 Oct 2005 10:16:07 +0000 (GMT) (envelope-from bde@zeta.org.au) Received: from mailproxy2.pacific.net.au (mailproxy2.pacific.net.au [61.8.0.87]) by mailout2.pacific.net.au (8.13.4/8.13.4/Debian-3) with ESMTP id j95AG5GZ027488; Wed, 5 Oct 2005 20:16:05 +1000 Received: from katana.zip.com.au (katana.zip.com.au [61.8.7.246]) by mailproxy2.pacific.net.au (8.13.4/8.13.4/Debian-3) with ESMTP id j95AG35b005146; Wed, 5 Oct 2005 20:16:03 +1000 Date: Wed, 5 Oct 2005 20:16:03 +1000 (EST) From: Bruce Evans X-X-Sender: bde@delplex.bde.org To: Steve Kargl In-Reply-To: <20051005032400.GA6736@troutmask.apl.washington.edu> Message-ID: <20051005191936.N51036@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> <20051002191841.GA40568@troutmask.apl.washington.edu> <20051004164618.U47262@delplex.bde.org> <20051005032400.GA6736@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, 05 Oct 2005 10:16:08 -0000 On Tue, 4 Oct 2005, Steve Kargl wrote: > On Tue, Oct 04, 2005 at 07:31:43PM +1000, Bruce Evans wrote: >> On Sun, 2 Oct 2005, Steve Kargl wrote: > >> I converted to ccoshf() for testing it and comparing with a simpler >> implementation that uses the formula, then converted back. There are > So which implementation do you prefer: your simpler version or > your fixed version of my first cut at ccosh? If course I prefer the simpler version :-). I think version that handles all the exceptional cases explicitly is useful mainly for bootstrapping. Once things work the exceptional cases are better documented by regression checks for them. I looked at ccosh() in glibc and cephes. The generic version in glibc-2.3.2 is similar to your version except, it uses fpclassify() instead of bit fiddling, feraisexcept (FE_INVALID) instead of invalid operations, and many gccisms. It has the same bugs for ccosh(x+Iy) where x is finite, cosh(x) is infinite and y is 0. glibc doesn't have ccosh() in asm except on m68k's, but it has cexp() in asm on i386's; for that it uses far too much code to handle the exceptional case and has the same bug when cosh(x) is infinite, etc. cephes has all c9x complex functions and some others including cgamma() and zetac(), but has almost no support for exceptional values -- for ccosh() it just uses the formula. >> I didn't bother optimizing isnan(x), etc. using hx, etc. >> > > For a committable version do you want/prefer the isnan > or the bit twiddling? Depends if isnan() is significantly inefficient. The classification macros are used a lot in your version so I think it is worth avoiding them, but in my version they are only used in rare cases. >> 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; >> % + return f; > > Wow. I'm used to programming in Fortran and would never have > written "creal(f) = ...". This looks like your assigning a > value to a function. (For those in the peanut gallery that > snicker at the mention of Fortran, please see the Fortran 2003 > standard.) > >> - working around gcc bugs. > > :-) The main gccisms in the glibc version are near these assignments. `creal(f) = x' is spelled `__real__ f = x'. This spelling is also used in rvalues. `I' is never used, even in constants. I found one bug, and it was for the above case in both version. The above case is x = 0, y finite; the sign is that of sin(y) but the above assumes that it is that of y (i.e., 1 since we normalized y). One corrected version of the above is: cimag(f) = x * con * copysign(0, sin(y)); but it is simpler to delete the code for this special case and just use the formula. The x finite, y = 0 case is quite different since cosh(x) can overflow so the formula doesn't work, and the sign of sinh(x) is that of x since sinh() doesn't oscillate like sin() so we don't need a copysign() or an sinh() to determine the sign. >> Fortunately, we can construct u + I*v by assigning to the real and complex >> parts. There should be a macro or inline function to do this. > > If we go the macro route, do you want it (them?) in math_private.h > or complex.h? If creal(f) appears on the LHS, is it a generic > reference so that type is forced to match the RHS? Is > > #define CMPLX((z),(x),(y)) {creal((z)) = (x), cimag((z)) = (y)} > > acceptable? I think it should go in math_private.h only and be an inline function. Here is the handling of the exceptional-exceptional cases (ones which are not handled by the formula) for the simpler version after fixing the (x finite, y finite) case: % if (x == 0 && !isfinite(y)) { % creal(f) = y - y; % cimag(f) = copysign(0, y - y); % return (f); % } % if (y == 0) { % creal(f) = cosh(x); % cimag(f) = isnan(x) ? copysign(0, x - x) : % copysign(0, x) * y; % return (f); % } Bruce From owner-freebsd-standards@FreeBSD.ORG Wed Oct 5 20:25:47 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 749DD16A41F for ; Wed, 5 Oct 2005 20:25:47 +0000 (GMT) (envelope-from des@des.no) Received: from tim.des.no (tim.des.no [194.63.250.121]) by mx1.FreeBSD.org (Postfix) with ESMTP id F2CA943D45 for ; Wed, 5 Oct 2005 20:25:46 +0000 (GMT) (envelope-from des@des.no) Received: from tim.des.no (localhost [127.0.0.1]) by spam.des.no (Postfix) with ESMTP id 464286160; Wed, 5 Oct 2005 22:25:30 +0200 (CEST) Received: from xps.des.no (des.no [80.203.228.37]) by tim.des.no (Postfix) with ESMTP id 361F66156; Wed, 5 Oct 2005 22:25:30 +0200 (CEST) Received: by xps.des.no (Postfix, from userid 1001) id A784033C1D; Wed, 5 Oct 2005 22:25:38 +0200 (CEST) To: Steve Kargl References: <20050929195552.GA14982@troutmask.apl.washington.edu> <20050930221846.T34595@delplex.bde.org> <20050930152734.GA20655@troutmask.apl.washington.edu> From: des@des.no (=?iso-8859-1?q?Dag-Erling_Sm=F8rgrav?=) Date: Wed, 05 Oct 2005 22:25:38 +0200 In-Reply-To: <20050930152734.GA20655@troutmask.apl.washington.edu> (Steve Kargl's message of "Fri, 30 Sep 2005 08:27:34 -0700") Message-ID: <86zmpnk8xp.fsf@xps.des.no> User-Agent: Gnus/5.110002 (No Gnus v0.2) Emacs/21.3 (berkeley-unix) MIME-Version: 1.0 Content-Type: text/plain; charset=iso-8859-1 Content-Transfer-Encoding: quoted-printable X-Spam-Tests: ALL_TRUSTED,AWL,BAYES_00 X-Spam-Learn: ham X-Spam-Score: -5.2/3.0 X-Spam-Checker-Version: SpamAssassin 3.0.4 (2005-06-05) on tim.des.no 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, 05 Oct 2005 20:25:47 -0000 Steve Kargl writes: > Good. I found a copy of n869 via google. Hopefully, the > actual standard and the draft do not greatly disagree. 25 quid will get you a hardcopy of the real thing, including the rationale, with TR1: http://www.amazon.co.uk/exec/obidos/ASIN/0470845732 DES --=20 Dag-Erling Sm=F8rgrav - des@des.no From owner-freebsd-standards@FreeBSD.ORG Thu Oct 6 21:26:15 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 9498C16A41F for ; Thu, 6 Oct 2005 21:26:15 +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 CBDD443D53 for ; Thu, 6 Oct 2005 21:26:13 +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 j96LQ8RO040870; Thu, 6 Oct 2005 14:26:08 -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 j96LQ3ge040869; Thu, 6 Oct 2005 14:26:03 -0700 (PDT) (envelope-from sgk) Date: Thu, 6 Oct 2005 14:26:02 -0700 From: Steve Kargl To: Bruce Evans Message-ID: <20051006212602.GA40609@troutmask.apl.washington.edu> References: <20050929195552.GA14982@troutmask.apl.washington.edu> <20050930221846.T34595@delplex.bde.org> <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> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20051005191936.N51036@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: Thu, 06 Oct 2005 21:26:15 -0000 On Wed, Oct 05, 2005 at 08:16:03PM +1000, Bruce Evans wrote: > On Tue, 4 Oct 2005, Steve Kargl wrote: > > >On Tue, Oct 04, 2005 at 07:31:43PM +1000, Bruce Evans wrote: > >>On Sun, 2 Oct 2005, Steve Kargl wrote: > > > >>I converted to ccoshf() for testing it and comparing with a simpler > >>implementation that uses the formula, then converted back. There are > > >So which implementation do you prefer: your simpler version or > >your fixed version of my first cut at ccosh? > > If course I prefer the simpler version :-). I think version that handles > all the exceptional cases explicitly is useful mainly for bootstrapping. I've implemented some inline functions (see below) for the assignments of the real and imaginary parts. After replacing creal(f) = cosh(x) * cos(y); cimag(f) = sinh(x) * sin(y); return (f); in your simpler function by return (cmplx(cosh(x) * cos(y), sinh(x) * sin(y)); your test program is generating a large number of z = 0x7f800000+I0x7fa00000 inf+Inan ccosh(z) = 0x7f800000+I0x7fe00000 inf+Inan ccosh1(z) = 0x7fe00000+I0x7fe00000 nan+Inan err = +0x600000+I+.0 where ccosh1 is your simpler function and ccosh is your fixed version of my ccosh (both have been converted to use the inline functions). n869.pdf says ccosh(inf+Inan) = inf+Inan > >>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; > >>% + return f; > > > >Wow. I'm used to programming in Fortran and would never have > >written "creal(f) = ...". This looks like your assigning a > >value to a function. (For those in the peanut gallery that > >snicker at the mention of Fortran, please see the Fortran 2003 > >standard.) Apparently, my version of gcc does not like the above code. How did you compile your s_ccosh1.c? > >If we go the macro route, do you want it (them?) in math_private.h > >or complex.h? If creal(f) appears on the LHS, is it a generic > >reference so that type is forced to match the RHS? Is > > > >#define CMPLX((z),(x),(y)) {creal((z)) = (x), cimag((z)) = (y)} > > > >acceptable? > > I think it should go in math_private.h only and be an inline function. --- /mnt1/src/lib/msun/src/math_private.h Fri Feb 4 12:05:39 2005 +++ ../math_private.h Thu Oct 6 14:06:37 2005 @@ -155,6 +155,43 @@ } while (0) /* + * Inline functions that can be used to construct complex values. + */ + +static __inline float complex +cmplxf(float __x, float __y) +{ + float *__ptr; + float complex __z; + __ptr = (float *) &__z; + *__ptr++ = __x; + *__ptr = __y; + return (__z); +} + +static __inline double complex +cmplx(double __x, double __y) +{ + double *__ptr; + double complex __z; + __ptr = (double *) &__z; + *__ptr++ = __x; + *__ptr = __y; + return (__z); +} + +static __inline long double complex +cmplxl(long double __x, long double __y) +{ + long double *__ptr; + long double complex __z; + __ptr = (long double *) &__z; + *__ptr++ = __x; + *__ptr = __y; + return (__z); +} + +/* * ieee style elementary functions * * We rename functions here to improve other sources' diffability -- Steve From owner-freebsd-standards@FreeBSD.ORG Fri Oct 7 14:51:27 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 56BB716A424 for ; Fri, 7 Oct 2005 14:51:27 +0000 (GMT) (envelope-from bde@zeta.org.au) Received: from mailout1.pacific.net.au (mailout1.pacific.net.au [61.8.0.84]) by mx1.FreeBSD.org (Postfix) with ESMTP id 5CB3843D48 for ; Fri, 7 Oct 2005 14:51:26 +0000 (GMT) (envelope-from bde@zeta.org.au) Received: from mailproxy2.pacific.net.au (mailproxy2.pacific.net.au [61.8.0.87]) by mailout1.pacific.net.au (8.13.4/8.13.4/Debian-3) with ESMTP id j97EpOmi031653; Sat, 8 Oct 2005 00:51:24 +1000 Received: from katana.zip.com.au (katana.zip.com.au [61.8.7.246]) by mailproxy2.pacific.net.au (8.13.4/8.13.4/Debian-3) with ESMTP id j97EpMWX028252; Sat, 8 Oct 2005 00:51:22 +1000 Date: Sat, 8 Oct 2005 00:51:22 +1000 (EST) From: Bruce Evans X-X-Sender: bde@delplex.bde.org To: Steve Kargl In-Reply-To: <20051006212602.GA40609@troutmask.apl.washington.edu> Message-ID: <20051007231439.F58005@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> <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> 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: Fri, 07 Oct 2005 14:51:27 -0000 On Thu, 6 Oct 2005, Steve Kargl wrote: > On Wed, Oct 05, 2005 at 08:16:03PM +1000, Bruce Evans wrote: >> If course I prefer the simpler version :-). I think version that handles >> all the exceptional cases explicitly is useful mainly for bootstrapping. > > I've implemented some inline functions (see below) for the assignments of > the real and imaginary parts. After replacing > > creal(f) = cosh(x) * cos(y); > cimag(f) = sinh(x) * sin(y); > return (f); > > in your simpler function by > > return (cmplx(cosh(x) * cos(y), sinh(x) * sin(y)); > > your test program is generating a large number of > > z = 0x7f800000+I0x7fa00000 inf+Inan > ccosh(z) = 0x7f800000+I0x7fe00000 inf+Inan > ccosh1(z) = 0x7fe00000+I0x7fe00000 nan+Inan > err = +0x600000+I+.0 > > where ccosh1 is your simpler function and ccosh is your fixed > version of my ccosh (both have been converted to use the inline > functions). n869.pdf says > > ccosh(inf+Inan) = inf+Inan This is because: (1) n869.pdf is broken here. For +inf*nan, if we knew that the nan meant an indeterminate but positive number (possibly including inf), then +inf*nan should be +inf, but in both of the expressions cosh(+inf) * cos(nan) and sinh(+inf) * sin(nan) we don't know the signs of the trig terms so the result should be nan for both the real and the imaginary parts. (2) s_ccosh.c still uses essentially your version for the return value in this case. It is "x * x + I * (y - y)". x * x is supposed to give inf, but the gcc bug propagates the nan y * y into the real part. (3) s_ccosh1.c is bug for bug compatible with s_ccosh.c here. Passing the regression test forced it to be incompatible. (4) using cmplx(u, v) instead of u + I*v fixed s_ccosh.c, and the regression test shows the difference. Unless this is fixed in C99[+TC], s_ccosh1.c needs to be broken to be compatible. The C99 behaviour has implications for cproj(). cproj(inf+Inan) is inf+Icopysign(0,nan), but cproj(nan+Inan) is nan+Inan. I think it is a larger bug to completely lose the nan-ness here since the real part is not really inf. BTW, many spurious infs and maybe some spurious nans occur when the terms overflow but the result doesn't. E.g., cosh(x) might be inf so cosh(x) * cos(x) is +-inf for all nonzero x, but if cos(x) is say 0.5 then overflow shouldn't occur until slightly larger x. I recently noticed the great lengths to which fdlibm goes to avoid similar overflows for the real functions. E.g., for cosh(), cosh(x) ~= exp(x)/2 and the RHS of this would overflow for x slightly less large than for the LHS. So instead, for large x fdlibm evaluates the RHS as exp(x/2)/2.*exp(x/2). This is imperfect due to extra roundoff error from the second multiplication. In practice, it gives errors of up to ~1.5 ulps in cosh() and then errors of up ~3 ulps in ccosh1(). So an extra-precision exp() is apparently needed to get 1-ulp accuracy even for the real cosh(). ccosh() could creep up on oveflow similarly. It would need about 53+[1 or 2] bits of extra precision to handle zero points of the trig terms where cosh() only nneds 1 or 2 bits extra to handle the 1./2 term. I don't think this will be implemented soon. >>>> 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; >>>> % + return f; >>> >>> Wow. I'm used to programming in Fortran and would never have >>> written "creal(f) = ...". This looks like your assigning a >>> value to a function. (For those in the peanut gallery that >>> snicker at the mention of Fortran, please see the Fortran 2003 >>> standard.) > > Apparently, my version of gcc does not like the above > code. How did you compile your s_ccosh1.c? I used gcc-3.3.3. The problem is that creal(f) never was an lvalue since it is declared as a function in , and gcc has finally fixed its its default of violating the C standard's requirements for lvalues. The uglier syntax used in glibc (e.g., __real__ f = cos(y)) still works. > --- /mnt1/src/lib/msun/src/math_private.h Fri Feb 4 12:05:39 2005 > +++ ../math_private.h Thu Oct 6 14:06:37 2005 > @@ -155,6 +155,43 @@ > } while (0) > > /* > + * Inline functions that can be used to construct complex values. > + */ > + > +static __inline float complex > +cmplxf(float __x, float __y) > +{ > + float *__ptr; > + float complex __z; > + __ptr = (float *) &__z; > + *__ptr++ = __x; > + *__ptr = __y; > + return (__z); > +} Use __real__/__imag or something less ugly if possible. These and creal()/cimag() are mentioned in gcc.info but neither is really documented there. Neither is mentioned to work as an lvalue. But the layout of "float complex" needed for the pointer hacks in the above to work is documented to not always work there -- it is documented that the the floats may noncontiguous, even with 1 in memory and the other in a register. Please think of better names. I think "complex" should always be abbreviated to "c" in libm. Maybe cpackf()? Bruce From owner-freebsd-standards@FreeBSD.ORG Fri Oct 7 19:02:40 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 5DA2E16A41F for ; Fri, 7 Oct 2005 19:02:40 +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 D3B1843D48 for ; Fri, 7 Oct 2005 19:02:39 +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 j97J2dJi056751; Fri, 7 Oct 2005 12:02:39 -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 j97J2d0V056750; Fri, 7 Oct 2005 12:02:39 -0700 (PDT) (envelope-from sgk) Date: Fri, 7 Oct 2005 12:02:39 -0700 From: Steve Kargl To: Bruce Evans Message-ID: <20051007190239.GA78674@troutmask.apl.washington.edu> References: <20050929195552.GA14982@troutmask.apl.washington.edu> <20050930221846.T34595@delplex.bde.org> <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> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20051007231439.F58005@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: Fri, 07 Oct 2005 19:02:40 -0000 On Sat, Oct 08, 2005 at 12:51:22AM +1000, Bruce Evans wrote: > On Thu, 6 Oct 2005, Steve Kargl wrote: > > >your test program is generating a large number of > > > >z = 0x7f800000+I0x7fa00000 inf+Inan > >ccosh(z) = 0x7f800000+I0x7fe00000 inf+Inan > >ccosh1(z) = 0x7fe00000+I0x7fe00000 nan+Inan > >err = +0x600000+I+.0 > > > >where ccosh1 is your simpler function and ccosh is your fixed > >version of my ccosh (both have been converted to use the inline > >functions). n869.pdf says > > > >ccosh(inf+Inan) = inf+Inan > > This is because: > (1) n869.pdf is broken here. For +inf*nan, if we knew that the nan > meant an indeterminate but positive number (possibly including inf), > then +inf*nan should be +inf, but in both of the expressions cosh(+inf) > * cos(nan) and sinh(+inf) * sin(nan) we don't know the signs of the > trig terms so the result should be nan for both the real and the > imaginary parts. It's still broken in n1124.pdf. > (2) s_ccosh.c still uses essentially your version for the return value > in this case. It is "x * x + I * (y - y)". x * x is supposed to give > inf, but the gcc bug propagates the nan y * y into the real part. I replaced the "x * x + I * (y - y)" with cmplx(x * x, y - y). > (4) using cmplx(u, v) instead of u + I*v fixed s_ccosh.c, and the regression > test shows the difference. Ah, this is the key. (ULP discussion of cosh(x) deleted). > >>>>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. > >>>>% + return f; > >>> > >>>Wow. I'm used to programming in Fortran and would never have > >>>written "creal(f) = ...". This looks like your assigning a > >>>value to a function. (For those in the peanut gallery that > >>>snicker at the mention of Fortran, please see the Fortran 2003 > >>>standard.) > > > >Apparently, my version of gcc does not like the above > >code. How did you compile your s_ccosh1.c? > > I used gcc-3.3.3. The problem is that creal(f) never was an lvalue since > it is declared as a function in , and gcc has finally fixed its > its default of violating the C standard's requirements for lvalues. The > uglier syntax used in glibc (e.g., __real__ f = cos(y)) still works. > > >+static __inline float complex > >+cmplxf(float __x, float __y) > >+{ > >+ float *__ptr; > >+ float complex __z; > >+ __ptr = (float *) &__z; > >+ *__ptr++ = __x; > >+ *__ptr = __y; > >+ return (__z); > >+} > > Use __real__/__imag or something less ugly if possible. These and > creal()/cimag() are mentioned in gcc.info but neither is really > documented there. Neither is mentioned to work as an lvalue. But > the layout of "float complex" needed for the pointer hacks in the > above to work is documented to not always work there -- it is documented > that the the floats may noncontiguous, even with 1 in memory and the > other in a register. Ugh. If I read Section 6.2.5(13), 6.2.5(20), and 6.2.6.1(2) correctly, then z in the above needs to occupy contiguous memory. I'll use the __real__/__imag__ gccism. > Please think of better names. I think "complex" should always be > abbreviated to "c" in libm. Maybe cpackf()? You don't do Fortran, do you? :-) program a complex z real :: x = 1., y = 2. z = cmplx(x,y) end program a OK, I'll use cpackf, cpack, and cpackl. -- Steve From owner-freebsd-standards@FreeBSD.ORG Fri Oct 7 20:35:30 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 37F6E16A41F for ; Fri, 7 Oct 2005 20:35:30 +0000 (GMT) (envelope-from bde@zeta.org.au) Received: from mailout1.pacific.net.au (mailout1.pacific.net.au [61.8.0.84]) by mx1.FreeBSD.org (Postfix) with ESMTP id 40EC043D48 for ; Fri, 7 Oct 2005 20:35:29 +0000 (GMT) (envelope-from bde@zeta.org.au) Received: from mailproxy1.pacific.net.au (mailproxy1.pacific.net.au [61.8.0.86]) by mailout1.pacific.net.au (8.13.4/8.13.4/Debian-3) with ESMTP id j97KZRiJ018770; Sat, 8 Oct 2005 06:35:27 +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 j97KZPoV006856; Sat, 8 Oct 2005 06:35:26 +1000 Date: Sat, 8 Oct 2005 06:35:26 +1000 (EST) From: Bruce Evans X-X-Sender: bde@delplex.bde.org To: Steve Kargl In-Reply-To: <20051007190239.GA78674@troutmask.apl.washington.edu> Message-ID: <20051008052850.S59139@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> <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> 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: Fri, 07 Oct 2005 20:35:30 -0000 On Fri, 7 Oct 2005, Steve Kargl wrote: > On Sat, Oct 08, 2005 at 12:51:22AM +1000, Bruce Evans wrote: >> On Thu, 6 Oct 2005, Steve Kargl wrote: >>> ... n869.pdf says >>> >>> ccosh(inf+Inan) = inf+Inan >> >> This is because: >> (1) n869.pdf is broken here. For +inf*nan, if we knew that the nan >> meant an indeterminate but positive number (possibly including inf), >> then +inf*nan should be +inf, but in both of the expressions cosh(+inf) >> * cos(nan) and sinh(+inf) * sin(nan) we don't know the signs of the >> trig terms so the result should be nan for both the real and the >> imaginary parts. > > It's still broken in n1124.pdf. I just checked that too. The rationale (both the 1998 n850.pdf one and the new one near n1124.pdf) says a bit about this. Infinities are supposed to be preserved as much as possible. In particular, it is almost right for Inf * NaN to be Inf, not NaN, since NaN is best thought of as an indeterminate number and Inf * non_NaN is +-Inf except when non_NaN is 0. There is a problem with the sign if the resulting Inf being indeterminate, but that vanishes if we use projective infinity. However the problem with Inf * 0 is too large, so both Inf * 0 and Inf * NaN are NaN for IEEE floating point. The rationale for ccosh(inf+Inan) = inf+Inan is partly that after projectivizing the value is unambigously the unique projective infinity. I think this is true for the finite precision case because cos(y) and sin(y) are never (?) zero for real rational x != 0. Thus the cosh(x) and sinh(x) terms eventually grow large enough to give +-Inf when multiplied by the cos(y) and sin(y) terms (calculating everything in infinite precision until the final step of rounding to double precision gives +-Inf). There is a similar rationale for pow(-2.0, Inf) being +Inf (its sign is uniquely determined since the domain is limited to even values if the base is 2). I like this rationale for pow() but not for ccosh(). Where pow(-2.0, y) increases monotonically to +Inf as y increases to +Inf, ccosh(z) has an essentially singularity as z -> projective infinity, so in infinite precision it takes on all values infinitely many times and in finite precision it oscilates wildly. >>>>>> 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. >> Please think of better names. I think "complex" should always be >> abbreviated to "c" in libm. Maybe cpackf()? > > You don't do Fortran, do you? :-) > > program a > complex z > real :: x = 1., y = 2. > z = cmplx(x,y) > end program a > > OK, I'll use cpackf, cpack, and cpackl. I'm used to C, and forget what Fortran does. I think of (x, y) -> x+I*y as a constructor but couldn't think of a good (short) name with "constructor" in it. Other ideas that don't quite seem to work: - only provide an imaginifying constructor from y to I*y. I think this would give pessimizations from extra additions for addition x to it. - use C++ and overload + and * in x+I*y with non-broken operators so as to keep the sources clean and not have to clean them up when x+I*y is fixed in gcc. I would have to learn C++. Bruce From owner-freebsd-standards@FreeBSD.ORG Sat Oct 8 23:42:07 2005 Return-Path: X-Original-To: freebsd-standards@hub.freebsd.org Delivered-To: freebsd-standards@hub.freebsd.org Received: from mx1.FreeBSD.org (mx1.freebsd.org [216.136.204.125]) by hub.freebsd.org (Postfix) with ESMTP id DE41F16A41F; Sat, 8 Oct 2005 23:42:07 +0000 (GMT) (envelope-from rodrigc@FreeBSD.org) Received: from freefall.freebsd.org (freefall.freebsd.org [216.136.204.21]) by mx1.FreeBSD.org (Postfix) with ESMTP id 9A76E43D45; Sat, 8 Oct 2005 23:42:07 +0000 (GMT) (envelope-from rodrigc@FreeBSD.org) Received: from freefall.freebsd.org (rodrigc@localhost [127.0.0.1]) by freefall.freebsd.org (8.13.3/8.13.3) with ESMTP id j98Ng7CO038379; Sat, 8 Oct 2005 23:42:07 GMT (envelope-from rodrigc@freefall.freebsd.org) Received: (from rodrigc@localhost) by freefall.freebsd.org (8.13.3/8.13.1/Submit) id j98Ng73D038375; Sat, 8 Oct 2005 23:42:07 GMT (envelope-from rodrigc) Date: Sat, 8 Oct 2005 23:42:07 GMT From: Craig Rodrigues Message-Id: <200510082342.j98Ng73D038375@freefall.freebsd.org> To: marius@alchemy.franken.de, rodrigc@FreeBSD.org, freebsd-standards@FreeBSD.org Cc: Subject: Re: standards/63173: Patch to add getopt_long_only(3) to libc 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: Sat, 08 Oct 2005 23:42:08 -0000 Synopsis: Patch to add getopt_long_only(3) to libc State-Changed-From-To: patched->closed State-Changed-By: rodrigc State-Changed-When: Sat Oct 8 23:39:47 GMT 2005 State-Changed-Why: This patch is in RELENG_5 http://www.freebsd.org/cgi/query-pr.cgi?pr=63173