Date: Sun, 12 Aug 2012 23:06:31 -0000 From: Stephen Montgomery-Smith <stephen@missouri.edu> To: Bruce Evans <brde@optusnet.com.au> Cc: Diane Bruce <db@db.net>, John Baldwin <jhb@freebsd.org>, David Chisnall <theraven@freebsd.org>, Bruce Evans <bde@freebsd.org>, Steve Kargl <sgk@troutmask.apl.washington.edu>, David Schultz <das@freebsd.org>, Peter Jeremy <peter@rulingia.com>, Warner Losh <imp@bsdimp.com> Subject: Re: Use of C99 extra long double math functions after r236148 Message-ID: <50056247.2000800@missouri.edu> Resent-Message-ID: <20120812230624.GY20453@server.rulingia.com> In-Reply-To: <20120717200931.U6624@besplex.bde.org> 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> <C527B388-3537-406F-BA6D-2FA45B9EAA3B@FreeBSD.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>
next in thread | previous in thread | raw e-mail | index | archive | help
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).
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?50056247.2000800>