Date: Sun, 12 Aug 2012 23:13:54 -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: <50083E83.9090404@missouri.edu> Resent-Message-ID: <20120812231347.GF20453@server.rulingia.com> In-Reply-To: <20120719144432.N1596@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> <5006D13D.2080702@missouri.edu> <20120719144432.N1596@besplex.bde.org>
next in thread | previous in thread | raw e-mail | index | archive | help
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.
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?50083E83.9090404>