From owner-freebsd-standards@FreeBSD.ORG Mon Sep 26 05:41:17 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 6545816A41F for ; Mon, 26 Sep 2005 05:41:17 +0000 (GMT) (envelope-from muhaeocikbxcwqustell@ci.berwyn.il.us) Received: from ci.berwyn.il.us (IGLD-83-130-109-176.inter.net.il [83.130.109.176]) by mx1.FreeBSD.org (Postfix) with SMTP id B149643D4C for ; Mon, 26 Sep 2005 05:41:15 +0000 (GMT) (envelope-from muhaeocikbxcwqustell@ci.berwyn.il.us) Received: from [192.168.204.229] (helo=classconsciousness) by ci.berwyn.il.us with SMTP (Receivingorder ow 3.36 (Awestruck)) id xIkMZn-aNOMwv-Ot for freebsd-standards@freebsd.org; Mon, 26 Sep 2005 00:41:03 -0500 Message-ID: <146625.UKSKBTFPCXN@classconsciousness> From: "Muhammad Stelling" To: "Jimmy Sinquefield" Date: Mon, 26 Sep 2005 00:41:01 -0500 MIME-Version: 1.0 X-Priority: 3 X-Mailer: Receivingorder ow 3.36 Content-Type: text/plain; charset="us-ascii" Content-Transfer-Encoding: quoted-printable X-Content-Filtered-By: Mailman/MimeDel 2.1.5 Subject: No Failure , Pharrmacxy X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list Reply-To: Muhammad Stelling List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Mon, 26 Sep 2005 05:41:17 -0000 CelValProMeCia= ViaUltAmXanLev ebrexiumpeciaridialisgrarambienaxitra $ $$ = 3.75 1.213.33 = http://www.multikopest.= com freedman, imitated many of these fables in Latin iambics about to a = neighboring cedar, The first step has lost us all. If we using them. filching = = from their very altars a part of the sacrifice offered Whatever you do, do = with all your might. that the burlesque style of writing adopted by Scarron = and From owner-freebsd-standards@FreeBSD.ORG Mon Sep 26 05:41:18 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 8E1CB16A423 for ; Mon, 26 Sep 2005 05:41:18 +0000 (GMT) (envelope-from winwa@advancedmail.com) Received: from advancedmail.com (IGLD-83-130-109-176.inter.net.il [83.130.109.176]) by mx1.FreeBSD.org (Postfix) with SMTP id 02C8343D5A for ; Mon, 26 Sep 2005 05:41:16 +0000 (GMT) (envelope-from winwa@advancedmail.com) Received: from [192.168.204.229] (helo=classconsciousness) by advancedmail.com with SMTP (Receivingorder ow 3.36 (Awestruck)) id xIkMZn-aNOMwv-Ot for freebsd-standards@hub.freebsd.org; Mon, 26 Sep 2005 00:41:03 -0500 Message-ID: <146625.UKSKBTFPCXN@classconsciousness> From: "Jacaline Winward" To: "Dwyn Hartman" Date: Mon, 26 Sep 2005 00:41:01 -0500 MIME-Version: 1.0 X-Priority: 3 X-Mailer: Receivingorder ow 3.36 Content-Type: text/plain; charset="us-ascii" Content-Transfer-Encoding: quoted-printable X-Content-Filtered-By: Mailman/MimeDel 2.1.5 Cc: Subject: No Failure , Pharrmacxy X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list Reply-To: Jacaline Winward List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Mon, 26 Sep 2005 05:41:18 -0000 CelValProMeCia= ViaUltAmXanLev ebrexiumpeciaridialisgrarambienaxitra $ $$ = 3.75 1.213.33 = http://www.multikopest.= com freedman, imitated many of these fables in Latin iambics about to a = neighboring cedar, The first step has lost us all. If we using them. filching = = from their very altars a part of the sacrifice offered Whatever you do, do = with all your might. that the burlesque style of writing adopted by Scarron = and From owner-freebsd-standards@FreeBSD.ORG Mon Sep 26 11:02: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 52F6916A430 for ; Mon, 26 Sep 2005 11:02:27 +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 9730A43D64 for ; Mon, 26 Sep 2005 11:02:22 +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 j8QB2M65027238 for ; Mon, 26 Sep 2005 11:02:22 GMT (envelope-from owner-bugmaster@freebsd.org) Received: (from peter@localhost) by freefall.freebsd.org (8.13.3/8.13.1/Submit) id j8QB2Lah027204 for freebsd-standards@freebsd.org; Mon, 26 Sep 2005 11:02:21 GMT (envelope-from owner-bugmaster@freebsd.org) Date: Mon, 26 Sep 2005 11:02:21 GMT Message-Id: <200509261102.j8QB2Lah027204@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, 26 Sep 2005 11:02:27 -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 Thu Sep 29 19:55:52 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 E26EB16A41F for ; Thu, 29 Sep 2005 19:55:52 +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 8DE0543D49 for ; Thu, 29 Sep 2005 19:55:52 +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 j8TJtqki015020 for ; Thu, 29 Sep 2005 12:55:52 -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 j8TJtqai015019 for freebsd-standards@freebsd.org; Thu, 29 Sep 2005 12:55:52 -0700 (PDT) (envelope-from sgk) Date: Thu, 29 Sep 2005 12:55:52 -0700 From: Steve Kargl To: freebsd-standards@freebsd.org Message-ID: <20050929195552.GA14982@troutmask.apl.washington.edu> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline User-Agent: Mutt/1.4.2.1i Subject: 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, 29 Sep 2005 19:55:53 -0000 Is it permissible to implement the complex.h math functions in terms of functions in math.h. For example, if z = x + I * y, then cos(z) = cos(x) * cosh(y) - I sin(x) * sinh(y) This can be (naively?) implemented as double complex ccos(double complex z) { double x, y; x = creal(z); y = cimag(y); return (cosh(y) * (cos(x) - I * sin(x) * tanh(y))); } I don't own a copy of C99, so I can't easily check the wording of the standard. -- Steve From owner-freebsd-standards@FreeBSD.ORG Fri Sep 30 14:03:07 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 B63D916A41F for ; Fri, 30 Sep 2005 14:03:07 +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 1952943D49 for ; Fri, 30 Sep 2005 14:03:06 +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 j8UE34Jn001971; Sat, 1 Oct 2005 00:03:04 +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 j8UE31MA017499; Sat, 1 Oct 2005 00:03:02 +1000 Date: Sat, 1 Oct 2005 00:03:01 +1000 (EST) From: Bruce Evans X-X-Sender: bde@delplex.bde.org To: Steve Kargl In-Reply-To: <20050929195552.GA14982@troutmask.apl.washington.edu> Message-ID: <20050930221846.T34595@delplex.bde.org> References: <20050929195552.GA14982@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: MULTIPART/MIXED; BOUNDARY="0-1213393512-1128088981=:34595" 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, 30 Sep 2005 14:03:07 -0000 This message is in MIME format. The first part should be readable text, while the remaining parts are likely unreadable without MIME-aware tools. --0-1213393512-1128088981=:34595 Content-Type: TEXT/PLAIN; charset=X-UNKNOWN; format=flowed Content-Transfer-Encoding: QUOTED-PRINTABLE On Thu, 29 Sep 2005, Steve Kargl wrote: > Is it permissible to implement the complex.h math > functions in terms of functions in math.h. For There are no namespace problems AFAIK, and the C standard doesn't require any particular quality of implementation. > example, if z =3D x + I * y, then > > cos(z) =3D cos(x) * cosh(y) - I sin(x) * sinh(y) > > This can be (naively?) implemented as > > double complex > ccos(double complex z) > { > double x, y; > x =3D creal(z); > y =3D cimag(y); > return (cosh(y) * (cos(x) - I * sin(x) * tanh(y))); > } The original formula gives a much better implementation. It is missing serious bugs at points where tanh(y) =3D=3D +-Inf; it may be missing bugs at points where y is a NaN (*); tanh(y) takes more work to compute than sinh(y), and then multiplying it by cosh(y) gives 1 extra multultiplication that doesn't take much longer (relatively) but loses precision (relatively more). (*) The best treatment of exceptional args is unclear, but C99 specifies a treatment in detail in at least n869.txt. This spec probably rules out simply using formulae as above -- exceptional args need to be handled specially unless you can be sure that applying the formulae results in the specified NaNs, etc. From n869.txt: % ccos(z) =3D ccosh(iz) % ... Some functions, including ccos() are specified in terms of others like this= =2E % G.5.2.4 The ccosh functions %=20 % [#1] %=20 % -- ccosh(conj(z)) =3D conj(ccosh(z)) and ccosh is even. A strict reading says that we must ensure that if a formula is used, the rounding must always go in such a way that these symmetries are preserved _exactly_. Rounding to nearest may give this, but this is not clear. Rounding towards Inf obviously wouldn't give this. In the following, "I" is spelled "i" and "*" (multiplication) is omitted. %=20 % -- ccosh(+0+i0) returns 1+i0. The formula should give this automatically provided the real functions aren't broken. Handly of -0 and -i0 (in various combinations with +0's) doesn't seem to be specified for this function, although it is for others. The sign of the real 0 doesn't matter because the result is 1. I suppose the sign for ccosh(-i0) should track the sign of the arg in the same way as for similar real functions. Since ccosh(1) =3D=3D 1, the sign of the real 0 can't affect the sign of the imaginary part of the result, unlike for some other functions. I just realized that many of the complications here are for type-generic functions. 0.0 is quite different from (0.0 + I * 0.0) since it has type double but the latter has type Complex. When combined with +0 being slightly different from -0, this gives zillions of cases. %=20 % -- ccosh(+0+i) returns NaN=B1i0 (where the sign of the % imaginary part of the result is unspecified) and raises % the invalid exception. This is apparently a bug in n869.txt. ccosh() on finite values is never a NaN. %=20 % -- ccosh(++i0) returns ++i0. Here ++ is not the C operator, but something like +-. It may be another bug that this says ++ and not +-. Shouldn't the result reverse the sign? %=20 % -- ccosh(++i) returns ++iNaN and raises the invalid % exception. %=20 % -- ccosh(x+i) returns NaN+iNaN and raises the invalid % exception, for finite nonzero x. More bugs in n869.txt. %=20 % -- ccosh(++iy) returns +cis(y), for finite nonzero y. cis(y) is (cos(y) + I * sin(y)). This says that if the arg is pure imaginary then the implicit real part of the arg (=3D 0) must not affect the result, and it says something about signs. %=20 % -- ccosh(+0+iNaN) returns NaN=B1i0 (where the sign of the % imaginary part of the result is unspecified). This requires not propagating a NaN in the imaginary part to the real part if the real part is precisely 0. %=20 % -- ccosh(++iNaN) returns ++iNaN. The case of a pure imaginary NaN may be slightly different. I don't understand the non-ASCII character in the spec for ccosh(+0+iNaN). It is '\xb1'. %=20 % -- ccosh(x+iNaN) returns NaN+iNaN and optionally raises % the invalid exception, for finite nonzero x. Nonzero x causes the NaN to be propagated to the real part. Formulae tend to do this automatically; it is the previous 2 cases that using the formula for ccosh() breaks. %=20 % -- ccosh(NaN+i0) returns NaN=B1i0 (where the sign of the % imaginary part of the result is unspecified). %=20 % -- ccosh(NaN+iy) returns NaN+iNaN and optionally raises % the invalid exception, for all nonzero numbers y. Similarly if only the real part of the arg is a NaN. %=20 % -- ccosh(NaN+iNaN) returns NaN+iNaN. This requirement is obviously right and is hard to avoid satisfy using formulae. About accuratcy of complex functions: It is unreasonable to require or expect the result to be as accurate as possible for the real and imaginary parts independently, since even ordinary operations don't have this propery. E.g., consider the simple polynomial function f(z) =3D z*z. For z =3D x+Iy, f(z) =3D (x*x-y*y)+2Ix*y= =2E Here 2Ix*y is almost as accurate as x and y, but x*x-y*y only has 1 or 2 bits of accuracy when x is as close as possible to y (but different). It might be useful for cpow(z, 2) to be much more accurate than z*z, but this would be hard to implement. Standard polynomial approximation methods just won't give enough accuracy. OTOH, for cexp() and probably for most related functions, equal (relative) accuracy for the real and imaginary parts occurs automatically using standard formulae (since cexp(z) depends on Re(z) and Im(z) independently): cexp(z) =3D exp(x) * (cos(x) + I * sin(x)) The multiplication expands the roundoff errors for the real and imaginary parts independently, so the final relative error for each is < 2 ulps if the errors for exp(), cos() and sin() are always < 1 ulp. I'm not happy with an error of 2 ulps of course, but don't know any morwe accurate way implementing cexp() short of doing everything in more precision. Standard polynomial approximations with complex args might be more accurate=20 relative to |cexp(z)| =3D exp(x), but couldn't be very accurate for the real and imaginary parts independently. This corresponds to the arg reduction requirements for cos(x) and sin(x) being very different. Remez methods also seem to be completely unsuitable for for complex args since according to my limited understanding of them they depend on delicate cancellations which are only likely to occur when their complex arg is in a discrete set of directions. For ccos(z) implemented as: ccos(z) =3D cos(x) * cosh(y) - I sin(x) * sinh(y) the error is again < 2 ulps for each of the real and imaginary parts provided it is < 1 ulp for all the real functions. Bruce --0-1213393512-1128088981=:34595-- From owner-freebsd-standards@FreeBSD.ORG Fri Sep 30 15:27: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 82F6116A41F for ; Fri, 30 Sep 2005 15:27: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 3196F43D48 for ; Fri, 30 Sep 2005 15:27:40 +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 j8UFRdnQ020840; Fri, 30 Sep 2005 08:27: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 j8UFRYPG020839; Fri, 30 Sep 2005 08:27:34 -0700 (PDT) (envelope-from sgk) Date: Fri, 30 Sep 2005 08:27:34 -0700 From: Steve Kargl To: Bruce Evans Message-ID: <20050930152734.GA20655@troutmask.apl.washington.edu> References: <20050929195552.GA14982@troutmask.apl.washington.edu> <20050930221846.T34595@delplex.bde.org> Mime-Version: 1.0 Content-Type: text/plain; charset=unknown-8bit Content-Disposition: inline Content-Transfer-Encoding: 8bit In-Reply-To: <20050930221846.T34595@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, 30 Sep 2005 15:27:40 -0000 On Sat, Oct 01, 2005 at 12:03:01AM +1000, Bruce Evans wrote: > On Thu, 29 Sep 2005, Steve Kargl wrote: > > >Is it permissible to implement the complex.h math > >functions in terms of functions in math.h. For > > There are no namespace problems AFAIK, and the C standard doesn't > require any particular quality of implementation. Good. I found a copy of n869 via google. Hopefully, the actual standard and the draft do not greatly disagree. As you may of guessed, I've put logl() on the back burner for the moment. I was hoping the complex.h functions would be easier to do. > >example, if z = x + I * y, then > > > >cos(z) = cos(x) * cosh(y) - I sin(x) * sinh(y) > > > >This can be (naively?) implemented as > > > >double complex > >ccos(double complex z) > >{ > > double x, y; > > x = creal(z); > > y = cimag(y); > > return (cosh(y) * (cos(x) - I * sin(x) * tanh(y))); > >} > > The original formula gives a much better implementation. It is missing > serious bugs at points where tanh(y) == +-Inf; it may be missing bugs > at points where y is a NaN (*); tanh(y) takes more work to compute than > sinh(y), and then multiplying it by cosh(y) gives 1 extra multultiplication > that doesn't take much longer (relatively) but loses precision (relatively > more). tanh(y) can't be +-Inf. |tanh(y)| <= 1 for all y (except NaN), so cos(), sin(), and tanh() are bounded and only cosh() diverges to +Inf for large argument. 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 can restore the sinh(y). > (*) The best treatment of exceptional args is unclear, but C99 specifies > a treatment in detail in at least n869.txt. This spec probably rules out > simply using formulae as above -- exceptional args need to be handled > specially unless you can be sure that applying the formulae results in > the specified NaNs, etc. From n869.txt: > > % ccos(z) = ccosh(iz) > % ... > > Some functions, including ccos() are specified in terms of others like this. OK, I'll use the implementation details of ccosh to infer the handling of exceptional values. > % G.5.2.4 The ccosh functions > % > % -- ccosh(conj(z)) = conj(ccosh(z)) and ccosh is even. > > A strict reading says that we must ensure that if a formula is used, the > rounding must always go in such a way that these symmetries are preserved > _exactly_. Rounding to nearest may give this, but this is not clear. > Rounding towards Inf obviously wouldn't give this. I haven't looked closely enough at libm, but I suspect that it it at best handles round-to-nearest. The default rounding in FreeBSD is round-to-nearest. To place greater restrictions on the complex.h with respect to proper handling of all rounding modes than we have on libm would seem to be an undue requirement. If a programmer sets another rounding mode, then its her responsibility o test that libm has the expected behavior. > I just realized that many of the complications here are for type-generic > functions. 0.0 is quite different from (0.0 + I * 0.0) since it has > type double but the latter has type Complex. When combined with +0 > being slightly different from -0, this gives zillions of cases. Signed zero are going to be important especially for those complex.h functions with branch cuts (e.g., csqrt). We do not want to return a value on the wrong Riemann sheet. Of course, this is an implementation detail that I'll worry about when I consider a fucntion with a branch cut. > % -- ccosh(+0+i) returns NaNħi0 (where the sign of the > % imaginary part of the result is unspecified) and raises > % the invalid exception. > > This is apparently a bug in n869.txt. ccosh() on finite values is never > a NaN. > > % > % -- ccosh(++i0) returns ++i0. > > Here ++ is not the C operator, but something like +-. It may be another > bug that this says ++ and not +-. Shouldn't the result reverse the sign? I found n869.pdf. I wonder if the above is suppose to be ccosh(+Inf+i0) = +inf+i0 because the pdf has no ++i0 case. > > % > % -- ccosh(++i) returns ++iNaN and raises the invalid > % exception. > % > % -- ccosh(x+i) returns NaN+iNaN and raises the invalid > % exception, for finite nonzero x. > > More bugs in n869.txt. > > % > % -- ccosh(++iy) returns +cis(y), for finite nonzero y. > > cis(y) is (cos(y) + I * sin(y)). This says that if the arg is pure > imaginary then the implicit real part of the arg (= 0) must not affect > the result, and it says something about signs. Again, the pdf file has ccosh(+Inf+iy) = +Inf cis(y) (finite nonzero y). which I infer means the sign on Inf is determined from cis(y). Thanks for visiting the accuracy issue. I came to essentially the same conclusion that you have. If we had reasonable implementations of the long double complex function, we could use these to compute float complex and double complex to achieve better accuracy. Of course, performance would suffer. -- Steve From owner-freebsd-standards@FreeBSD.ORG Sat Oct 1 11:38:53 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 1F87A16A41F for ; Sat, 1 Oct 2005 11:38:53 +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 732A043D46 for ; Sat, 1 Oct 2005 11:38:52 +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 j91BclJ6025302; Sat, 1 Oct 2005 21:38: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 j91BciJ2020242; Sat, 1 Oct 2005 21:38:45 +1000 Date: Sat, 1 Oct 2005 21:38:44 +1000 (EST) From: Bruce Evans X-X-Sender: bde@delplex.bde.org To: Steve Kargl In-Reply-To: <20050930152734.GA20655@troutmask.apl.washington.edu> Message-ID: <20051001204005.N37704@delplex.bde.org> References: <20050929195552.GA14982@troutmask.apl.washington.edu> <20050930221846.T34595@delplex.bde.org> <20050930152734.GA20655@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: MULTIPART/MIXED; BOUNDARY="0-1263416119-1128166724=:37704" 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: Sat, 01 Oct 2005 11:38:53 -0000 This message is in MIME format. The first part should be readable text, while the remaining parts are likely unreadable without MIME-aware tools. --0-1263416119-1128166724=:37704 Content-Type: TEXT/PLAIN; charset=X-UNKNOWN; format=flowed Content-Transfer-Encoding: QUOTED-PRINTABLE On Fri, 30 Sep 2005, Steve Kargl wrote: > On Sat, Oct 01, 2005 at 12:03:01AM +1000, Bruce Evans wrote: >> On Thu, 29 Sep 2005, Steve Kargl wrote: >>> cos(z) =3D cos(x) * cosh(y) - I sin(x) * sinh(y) >>> >>> This can be (naively?) implemented as >>> >>> double complex >>> ccos(double complex z) >>> { >>> double x, y; >>> x =3D creal(z); >>> y =3D cimag(y); >>> return (cosh(y) * (cos(x) - I * sin(x) * tanh(y))); >>> } >> >> The original formula gives a much better implementation. It is missing >> serious bugs at points where tanh(y) =3D=3D +-Inf; it may be missing bug= s >> at points where y is a NaN (*); tanh(y) takes more work to compute than >> sinh(y), and then multiplying it by cosh(y) gives 1 extra multultiplicat= ion >> that doesn't take much longer (relatively) but loses precision (relative= ly >> more). > > tanh(y) can't be +-Inf. |tanh(y)| <=3D 1 for all y (except NaN), > so cos(), sin(), and tanh() are bounded and only cosh() diverges > to +Inf for large argument. Oops. > A slightly better implementation than > the above would be > > if (y =3D=3D 0.) return (cos(x)); > if (x =3D=3D 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 =3D=3D 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. > ... > OK, I'll use the implementation details of ccosh to infer the > handling of exceptional values. > >> % G.5.2.4 The ccosh functions >> % >> % -- ccosh(conj(z)) =3D conj(ccosh(z)) and ccosh is even. >> >> A strict reading says that we must ensure that if a formula is used, the >> rounding must always go in such a way that these symmetries are preserve= d >> _exactly_. Rounding to nearest may give this, but this is not clear. >> Rounding towards Inf obviously wouldn't give this. > > I haven't looked closely enough at libm, but I suspect that it > it at best handles round-to-nearest. The default rounding in > FreeBSD is round-to-nearest. To place greater restrictions on > the complex.h with respect to proper handling of all rounding > modes than we have on libm would seem to be an undue requirement. > If a programmer sets another rounding mode, then its her responsibility > o test that libm has the expected behavior. My main point is that even with round-to-nearest, it is unclear that the rounding goes the right way for standard formulae to work right. >> I just realized that many of the complications here are for type-generic >> functions. 0.0 is quite different from (0.0 + I * 0.0) since it has >> type double but the latter has type Complex. When combined with +0 >> being slightly different from -0, this gives zillions of cases. > > Signed zero are going to be important especially for those complex.h > functions with branch cuts (e.g., csqrt). We do not want to return > a value on the wrong Riemann sheet. Of course, this is an implementation > detail that I'll worry about when I consider a fucntion with a branch > cut. Except signed zeros (and infinities) don't really work for the complex case. Signed zeros give 8 directions (multiples of Pi/4) for approaching the origin, but there are an infinite number of possible directions (all angles). >> % -- ccosh(+0+i) returns NaN=B1i0 (where the sign of the >> % imaginary part of the result is unspecified) and raises >> % the invalid exception. >> >> This is apparently a bug in n869.txt. ccosh() on finite values is never >> a NaN. >> >> % >> % -- ccosh(++i0) returns ++i0. >> >> Here ++ is not the C operator, but something like +-. It may be another >> bug that this says ++ and not +-. Shouldn't the result reverse the sign= ? > > I found n869.pdf. I wonder if the above is suppose to be > > ccosh(+Inf+i0) =3D +inf+i0 > > because the pdf has no ++i0 case. Yes, the infinity sign apparently got lost in translation to ASCII. The previous one, ccosh(+0+i), is actually ccosh(+0+i*Inf), but doesn't look as weird when mistranslated since the infinity is lost at the end. You could also look at the POSIX spec. It is more up to date and has an index. However, it seems to be missing the detail spec for the complex functions in at least the text version of the 2001 draft (this specifies handling of NaNs in detail only for the real functions, and also uses the \xb1 character which I don't understand yet). > Thanks for visiting the accuracy issue. I came to essentially > the same conclusion that you have. If we had reasonable > implementations of the long double complex function, we could > use these to compute float complex and double complex to achieve > better accuracy. Of course, performance would suffer. Float complex can use double precision with tiny or negative loss of performance given our current low performance float functions. Also, long double precision is almost free for the double exp() family (exp/log/trig) on some machines (i386 and amd64) that can do these functions in hardware given our current not very high performance (hardware or not) double functions. I would use long double (real) versions for the initial implementation of the more exotic complex functions if possible. Unfortunately, it doesn't seem to be possible to implement the more exotic complex functions in terms of real functions. Bruce --0-1263416119-1128166724=:37704--