Date: Tue, 14 Feb 2017 13:12:50 +0000 From: "Montgomery-Smith, Stephen" <stephen@missouri.edu> To: Bruce Evans <brde@optusnet.com.au> Cc: Peter Jeremy <peter@rulingia.com>, Alan Braslau <alan.braslau@comcast.net>, "freebsd-numerics@freebsd.org" <freebsd-numerics@freebsd.org> Subject: Re: FreeBSD numerics - cpow() Message-ID: <37cf4f6e-6b72-6732-82c7-9c2f45bd7e6b@missouri.edu> In-Reply-To: <20170214185050.R941@besplex.bde.org> References: <20170213091051.600f91e1@zoo.hsd1.co.comcast.net> <20170214053712.GD84013@server.rulingia.com> <c9754b4d-1a2b-39ce-f20c-29bbcf8eafdf@missouri.edu> <20170214185050.R941@besplex.bde.org>
next in thread | previous in thread | raw e-mail | index | archive | help
On 02/14/2017 04:26 AM, Bruce Evans wrote: > On Tue, 14 Feb 2017, Montgomery-Smith, Stephen wrote: >=20 >> On 02/13/2017 11:37 PM, Peter Jeremy wrote: >>> On 2017-Feb-13 09:10:51 -0700, Alan Braslau >>> <alan.braslau@comcast.net> wrote: >>>> What is the current status of getting cpow() implemented in FreeBSD? >>> >>> There's a WIP in https://svnweb.freebsd.org/base/user/peterj/ >>> but I got caught up trying to work out how to perfectly multiply >>> two doubles and am not currently working on it. >>> >>> I wonder if we should implement something like >>> >>> double cpow(double x, double y) >>> { >>> return cexp(y * clog(x)); >>> } >=20 > It's actually: >=20 > double complex cpow(double complex x, double complex y) > { > return (cexp(y * clog(x)); > } >=20 >>> >>> just to have something to resolve symbols. >> >> If you look at page 478 of the document >> http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1124.pdf (which is, I >> think, the official C99 standard), the footnote suggests that this would >> be a reasonable thing to do. >=20 > This suggests that the official C99 standard is unreasonably bad. >=20 > The footnote is wrong anyway. cpow() is specially weakened by allowing > it to return spurious exceptions, but the same sentence that opens > this hole strengthens it by requiring it to return "appropriate > exceptions". Detecting when it is appropriate to return an overflow > exception requires precise classification exception boundaries and > doesn't seem to be required for any other function, but the boundary > is especially complicated for cpow() since it is 3-dimensional. For > cexp(), it is merely 1-dimensional, but still too hard to classify > precisely -- FreeBSD implementation has complications mainly to get > this boundary almost right. >=20 > Actually, C99 is too strong in its specification of not returning > spurious exceptions for other functions. This requires precisely > classifying everything on the non-exception side of the boundary > (you can reasonably misclassify by a few ulps on the other side, > but this is almost as hard as determining the precise boudary). >=20 > The specification should be that if an exception is detected (or not), > then the exception flags are set appropriately (it is not required to > detect exceptional args appropriately (precisely) except for a few > special cases like +-Inf). The implementation can be as bad as it can > get away with for both normal values and exceptional values, and the > rule for exceptions should be that there are no inconsistent ones > (except for the bad C99 rule for cpow(), and where it is not considered > inconsistent that the normal values are inconsistent -- e.g., exp(x) > should be monotonic in x, but might return DBL_MAX for a certain x, > then Inf for x + 1 ulp, then back to DBL_MAX for x + 2 ulps, then > stable at Inf for larger x. Low, quality like that is easy to avoid > for exp(), but near a 3-dimensional boundary it would take a > higher-precision calculation to even see if adding 0 or 1 ulps to the > 4 real components of 2 complex args crosses the boundary. >=20 > A sloppy implementation like the above is going to have very large, very > inappropriate classification of the boundary, but would meet my weakened > rule for classification of exceptions, with perhaps no special case for > cpow(). The result is whatever it is, with overflow set if the result > is Inf. cpow() is now allowed to set divison-by-zero spuriously, but > I think the sloppy implementation avoids that without really trying. >=20 > The imprecise cpowl() also meets my rule but not C99's. It is not only > imprecise, but misclassifies overflow boundaries by a lot. E.g., for > powl(2.0L, y), the overflow boundary should be at about log2(LDBL_MAX), > but is at log2(DBL_MAX). >=20 >> So I say yes, let's do it. >=20 > FreeBSD doesn't have clog() either, but I have an implementation. One of my big frustrations with this group is that there are many implementations that haven't been committed. Why hasn't clog been committed yet? I've seen the code - I've even tested it, and it seems to work remarkably well. What are we waiting for? Stephen
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?37cf4f6e-6b72-6732-82c7-9c2f45bd7e6b>