From owner-freebsd-numerics@FreeBSD.ORG Sun Aug 12 23:06:31 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 C7D731065675 for ; Sun, 12 Aug 2012 23:06:31 +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 4C5A88FC0C for ; Sun, 12 Aug 2012 23:06:31 +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 q7CN6V3o075713 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Mon, 13 Aug 2012 09:06:31 +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 q7CN6ON2021362 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Mon, 13 Aug 2012 09:06:24 +1000 (EST) (envelope-from peter@server.rulingia.com) Received: (from peter@localhost) by server.rulingia.com (8.14.5/8.14.5/Submit) id q7CN6O5D021361 for freebsd-numerics@freebsd.org; Mon, 13 Aug 2012 09:06:24 +1000 (EST) (envelope-from peter) Resent-From: Peter Jeremy Resent-Date: Mon, 13 Aug 2012 09:06:24 +1000 Resent-Message-ID: <20120812230624.GY20453@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 q6HD2Otk075222 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Tue, 17 Jul 2012 23:02:25 +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 q6HD2MLE067318 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Tue, 17 Jul 2012 23:02:24 +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 q6HD1xUl056214; Tue, 17 Jul 2012 08:02:00 -0500 (CDT) (envelope-from stephen@missouri.edu) Message-ID: <50056247.2000800@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> In-Reply-To: <20120717200931.U6624@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:06:31 -0000 X-Original-Date: Tue, 17 Jul 2012 08:01:59 -0500 X-List-Received-Date: Sun, 12 Aug 2012 23:06:31 -0000 On 07/17/2012 06:13 AM, Bruce Evans wrote: > On Mon, 16 Jul 2012, Stephen Montgomery-Smith wrote: > >> On 07/16/2012 06:37 PM, Stephen Montgomery-Smith wrote: >>> ... >>> The difficulty is that catanh, cacosh, and cacosh do not have unique >>> answers unless one makes some kind of accepted convention as to what >>> they should be. >>> >>> I did find some definitions at >>> http://publib.boulder.ibm.com/infocenter/zos/v1r11/index.jsp?topic=/com.ibm.zos.r11.bpxbd00/catan.htm. >>> >>> I'm going to guess that there is a typo in that the range of the >>> imaginary part should be [0,i pi], because the usual convention is that >>> acos of a real number is a real number between 0 and pi. > > I think C99 specificies the branch cus and boundary behaviour completely. > >>> We might get lucky, and find that the definitions of csqrt and clog in >>> the C99 standard are already set up so that the naive formulas for >>> cacosh, etc, just work. But whether they do or whether they don't, I >>> think I can do it. (As a first guess, I think that catanh and casinh >>> will work "out of the box" but cacosh is going to take a bit more work.) > > See below what happened for naive formulars for ccosh. > >>> Also casin, catan and cacos are essentially the same as casinh, etc, >>> using formulas like sin(ix) = i sinh(x). (The hardest part is to avoid >>> making a sign error.) >> >> I came up with pseudo code that looks a bit like this. >> >> complex casinh(complex z) { >> double x = z.re, y = z.im; >> >> if (y==0) >> return asinh(x); >> if (x==0) { >> if (fabs(y)<=1) return I*asin(y); >> else return signum(y)* ( >> log(fabs(y)+sqrt(y*y-1)) >> + I*PI/2); >> } >> if (x>0) >> return clog(z+csqrt(z*z+1)); >> else >> return -clog(-z+csqrt(z*z+1)); >> } > > I translated this to pari. There was a sign error for the log() term > (assuming that pari asinh() has the same branch cuts as C99): > > % \p 100 > % % PI = Pi; > % clog(z) = log(z); > % csqrt(z) = sqrt(z); > % fabs(x) = abs(x); > % signum(x) = sign(x); > % % casinh(z) = > % { > % local(x, y); > % % x = real(z); > % y = imag(z); > % if (y == 0, > % return (asinh(x)); > % ); > % if (x == 0, > % if (fabs(y) <= 1, > % return (I * asin(y)); > % , > % return (signum(y) * > % (-log(fabs(y) + sqrt(y * y - 1)) + I * PI / 2)); > % ); > % ); > % if (x > 0, > % return (clog(z + csqrt(z * z + 1))); > % , > % return (-clog(-z + csqrt(z * z + 1))); > % ); > % } > % % { > % forstep (x = -10, 10, 0.1, > % forstep (y = -10, 10, 0.1, > % z = x + I*y; > % r = casinh(z) - asinh(z); > % if (abs(r) > 1e-30, > % print("z = " z); > % print("casinh(z) = " casinh(z)); > % print(" asinh(z) = " asinh(z)); > % print("diff = " r); > % ); > % ); > % ); > % } > > (No differences found after fixing the sign.) > > Pari of course does all the calculations almost perfectly accurately (I > told it to provide 100 decimal digits). Most multi-precision calculators > can do the same. So pari can be used to develop the logic, but it is hard > to use it develop accurate routines for limited precision. > > The most obvious immediate difficulty in translating the above into C is > that y*y and z*z may overflow when the result shouldn't. hypot() and > cabs() handle similar problems, with remarkably large complications. > C99 with IEEE754 FP actually handles some aspects of overflow better > than pari. E.g., exp(10^9) gives oveflow in pari, with no way of handling > it AFAIK, while in C99 + IEEE754 it gives +Inf with FE_OVERFLOW. > > Bruce > > Excellent. I think I will use pari to write the test code to check the ULP of the result. And I'll look into using hypot, or its logic, to compute sqrt(y*y-1).