From owner-freebsd-numerics@FreeBSD.ORG Sun Aug 12 23:13:54 2012 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [69.147.83.52]) by hub.freebsd.org (Postfix) with ESMTP id F39C7106564A for ; Sun, 12 Aug 2012 23:13:53 +0000 (UTC) (envelope-from peter@rulingia.com) Received: from vps.rulingia.com (host-122-100-2-194.octopus.com.au [122.100.2.194]) by mx1.freebsd.org (Postfix) with ESMTP id 77D548FC1B for ; Sun, 12 Aug 2012 23:13:53 +0000 (UTC) Received: from server.rulingia.com (c220-239-249-137.belrs5.nsw.optusnet.com.au [220.239.249.137]) by vps.rulingia.com (8.14.5/8.14.5) with ESMTP id q7CNDrPP075983 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Mon, 13 Aug 2012 09:13:53 +1000 (EST) (envelope-from peter@rulingia.com) X-Bogosity: Ham, spamicity=0.000000 Received: from server.rulingia.com (localhost.rulingia.com [127.0.0.1]) by server.rulingia.com (8.14.5/8.14.5) with ESMTP id q7CNDlmw021958 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Mon, 13 Aug 2012 09:13:47 +1000 (EST) (envelope-from peter@server.rulingia.com) Received: (from peter@localhost) by server.rulingia.com (8.14.5/8.14.5/Submit) id q7CNDlNN021957 for freebsd-numerics@freebsd.org; Mon, 13 Aug 2012 09:13:47 +1000 (EST) (envelope-from peter) Resent-From: Peter Jeremy Resent-Date: Mon, 13 Aug 2012 09:13:47 +1000 Resent-Message-ID: <20120812231347.GF20453@server.rulingia.com> Resent-To: freebsd-numerics@freebsd.org Received: from vps.rulingia.com (host-122-100-2-194.octopus.com.au [122.100.2.194]) by server.rulingia.com (8.14.5/8.14.5) with ESMTP id q6JH6je5015859 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Fri, 20 Jul 2012 03:06:46 +1000 (EST) (envelope-from stephen@missouri.edu) Received: from wilberforce.math.missouri.edu (wilberforce.math.missouri.edu [128.206.184.213]) by vps.rulingia.com (8.14.5/8.14.5) with ESMTP id q6JH6g7K081475 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Fri, 20 Jul 2012 03:06:45 +1000 (EST) (envelope-from stephen@missouri.edu) Received: from [127.0.0.1] (wilberforce.math.missouri.edu [128.206.184.213]) by wilberforce.math.missouri.edu (8.14.5/8.14.5) with ESMTP id q6JH6B9x047550; Thu, 19 Jul 2012 12:06:12 -0500 (CDT) (envelope-from stephen@missouri.edu) Message-ID: <50083E83.9090404@missouri.edu> From: Stephen Montgomery-Smith Mail-Followup-To: freebsd-numerics@freebsd.org User-Agent: Mozilla/5.0 (X11; Linux i686; rv:13.0) Gecko/20120615 Thunderbird/13.0.1 MIME-Version: 1.0 To: Bruce Evans References: <20120529045612.GB4445@server.rulingia.com> <20120711223247.GA9964@troutmask.apl.washington.edu> <20120713114100.GB83006@server.rulingia.com> <201207130818.38535.jhb@freebsd.org> <9EB2DA4F-19D7-4BA5-8811-D9451CB1D907@theravensnest.org> <20120713155805.GC81965@zim.MIT.EDU> <20120714120432.GA70706@server.rulingia.com> <20120717084457.U3890@besplex.bde.org> <5004A5C7.1040405@missouri.edu> <5004DEA9.1050001@missouri.edu> <20120717200931.U6624@besplex.bde.org> <5006D13D.2080702@missouri.edu> <20120719144432.N1596@besplex.bde.org> In-Reply-To: <20120719144432.N1596@besplex.bde.org> Content-Type: text/plain; charset=ISO-8859-1; format=flowed Content-Transfer-Encoding: 7bit Cc: Diane Bruce , John Baldwin , David Chisnall , Bruce Evans , Steve Kargl , David Schultz , Peter Jeremy , Warner Losh Subject: Re: Use of C99 extra long double math functions after r236148 X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , Date: Sun, 12 Aug 2012 23:13:54 -0000 X-Original-Date: Thu, 19 Jul 2012 12:06:11 -0500 X-List-Received-Date: Sun, 12 Aug 2012 23:13:54 -0000 On 07/19/2012 01:37 AM, Bruce Evans wrote: > On Wed, 18 Jul 2012, Stephen Montgomery-Smith wrote: > >> I went on a long road trip yesterday, so I didn't get any code >> written, but I did have a lot of thoughts about clog and casinh. >> >> First, the naive formula (here z=x+I*y) >> clog(z) = cpack(log(hypot(x,y)),atan2(x,y)) >> is going to work in a lot more edge cases then one might expect. This >> is because hypot and atan2, especially atan2, already do a rather good >> job getting the edge cases right. I am thinking in particular of when >> x or y are 0 or -0, or one of them is infinity or -infinity. > > Right, clog is deceptively simple. This is because it decomposes perfectly > into 2 real functions of 2 real variables and both of these functions are > standard and already implemented almost as well as possible. ISTR das > saying that it had a complicated case, but I don't see even one. atan2() > is supposed to handle all combinations of +-0. Now I remember a potential > problem. Complex functions should have only poles and zeros, with > projective infinity and "projective zero" (= inverse of projective > infinity). Real functions can and do have affine infinities and zeros > (+-Inf and +-0), with more detailed special cases. It's just impossible > to have useful, detailed special cases for all the ways of approaching > complex (projective) infinity and 0. > I think Kahan wanted projective infinity in IEEE7xx in ~1980. > Intel 8087 had both projective infinity and affine infinities, but > projective infinity didn't make it into the first IEEEE7xx, and > hardly anyone understood it and it was eventually dropped from > Intel FPUs (I think it was in 80287; then in i486 it was reduced > to a bit in the control word that can never be cleared (the bit is > to set affine infinities); then in SSE the bit went away too). > However, C99 tries too hard to make complex functions reduce to real > functions when everything is purely real or purely complex. So most > of the special cases for +-0 and +-Inf affect complex functions (for > other directions of approaching 0 and infinity, not much is specified > but you should try to be as continuous as possible, where continuity > has delicate unclear meanings since it is related to discontinuous > sign functions). Hopefully, the specification of imag(clog()) is > that it has the same sign behaviour as atan2(), so you can just use > atan2(). The sign conventions for both are arbitrary, but they > shouldn't be gratuitously different. You still have to check that > they aren't non-gratuitously different, because different conventions > became established. I checked. Actually the sign conventions are not that arbitrary. But as a mathematician I would say they are a bit useless, e.g. atan(infinity,infinity) = pi/4 = 45 degrees How do you know that the two infinities are the same? One could be double the other. If it had been up to me, there would have been finite numbers, and nan. And none of this -0. > I had forgotten that pari doesn't support -0 at all (AFAIK). I certainly > had to change a sign to get match the pari result for 0+y*I, but it was > the sign of y. Your original code seems to have y where it should have x: Oops. > I don't know what happens with zeros of complex inverse trig functions. > I think they don't have many (like log()), but their real and imaginary > parts do, and they are too general for accurate behaviour of the real > and imaginary parts relative to themselves to fall out. casinh(z) is zero only when z=0, and near that point I could use Taylor's series (but a lot of terms would be needed because the Taylot series converges quite slowly). I can now see that the separate cases of the real part and imaginary parts of casinh being zero is going to be hard. I'll probably end up reading the paper Jeremy suggested, and implementing that. But I always prefer self discovery first.