From owner-freebsd-numerics@freebsd.org Tue Feb 14 10:45:20 2017 Return-Path: Delivered-To: freebsd-numerics@mailman.ysv.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) by mailman.ysv.freebsd.org (Postfix) with ESMTP id 988F3CDF915 for ; Tue, 14 Feb 2017 10:45:20 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail110.syd.optusnet.com.au (mail110.syd.optusnet.com.au [211.29.132.97]) by mx1.freebsd.org (Postfix) with ESMTP id 47DC31628 for ; Tue, 14 Feb 2017 10:45:19 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from besplex.bde.org (c122-106-153-191.carlnfd1.nsw.optusnet.com.au [122.106.153.191]) by mail110.syd.optusnet.com.au (Postfix) with ESMTPS id A9D27104236; Tue, 14 Feb 2017 21:26:34 +1100 (AEDT) Date: Tue, 14 Feb 2017 21:26:33 +1100 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: "Montgomery-Smith, Stephen" cc: Peter Jeremy , Alan Braslau , "freebsd-numerics@freebsd.org" Subject: Re: FreeBSD numerics - cpow() In-Reply-To: Message-ID: <20170214185050.R941@besplex.bde.org> References: <20170213091051.600f91e1@zoo.hsd1.co.comcast.net> <20170214053712.GD84013@server.rulingia.com> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed X-Optus-CM-Score: 0 X-Optus-CM-Analysis: v=2.2 cv=H7qr+6Qi c=1 sm=1 tr=0 a=Tj3pCpwHnMupdyZSltBt7Q==:117 a=Tj3pCpwHnMupdyZSltBt7Q==:17 a=kj9zAlcOel0A:10 a=C_IRinGWAAAA:8 a=6I5d2MoRAAAA:8 a=b4LDLZbEAAAA:8 a=3jIIAKEMzJUsS6AJfKIA:9 a=CjuIK1q_8ugA:10 a=VOW5rJRXBamZ5z19bD7L:22 a=IjZwj45LgO3ly-622nXo:22 a=20T61YgZp4ItGotXEy2O:22 X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.23 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Tue, 14 Feb 2017 10:45:20 -0000 On Tue, 14 Feb 2017, Montgomery-Smith, Stephen wrote: > On 02/13/2017 11:37 PM, Peter Jeremy wrote: >> On 2017-Feb-13 09:10:51 -0700, Alan Braslau 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)); >> } It's actually: double complex cpow(double complex x, double complex y) { return (cexp(y * clog(x)); } >> >> 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. This suggests that the official C99 standard is unreasonably bad. 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. 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). 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. 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. 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). > So I say yes, let's do it. FreeBSD doesn't have clog() either, but I have an implementation. Bruce