FreeBSD numerics - cpow()
brde at optusnet.com.au
Tue Feb 14 10:45:20 UTC 2017
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 <alan.braslau at 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));
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.
More information about the freebsd-numerics