Interval enclosure of X**Y for negative x in X
David G Hough at validgh
validgh
Mon Mar 30 19:30:04 PST 1998
As I observed a while back,
> Even more interesting are cases like
> [-2.0,-2.0] ** [0.3,2.7]
> that really bring the exception issues in point x**y into sharp focus. Is it
> really right to evaluate that interval expression as the convex hull of
> ([-2.0,-2.0]) ** 1 U ([-2.0,-2.0]) ** 2 ?
> I'm beginning to think that floating-point x**y with negative x should just
> be an invalid exception, regardless of the integralness of y,
> and that programs that intend to raise negative numbers to
> integer powers should do so explicitly, as is easily
> expressible in Fortran x**n though more problematic in C.
It's one thing to raise positive rational numbers to rational powers,
and another to raise real numbers to real powers.
Considered in the complex plane,
x**y for negative real x and varying y is almost always complex (and
perhaps multivalued at that), and is only
by accident real on a subset of measure zero in y.
Many programmers might prefer the result of x**(1.0/3) to be the same as the
result of x**(2.0/6), and so defining x**(n/m) to be the m'th root of (x**n)
would be problematic for negative x. In contrast, defining x**y to
be exp(y*log(x)) reduces most of the problem of defining x**y to the problem of
defining log(x). In real variables, log(x) is not defined for negative
x, and so neither is x**y. And real floating point arithmetic should work
that way as well.
Now interval folks have argued that real floating-point sqrt([-1,+eps]) should
ignore the negative part of the operand interval, and base the computed result
solely on the operand for which the function is defined, namely [0,+eps].
That is at least a measurable subset as long as eps > 0.
It is a significant further step to assert that interval X**Y, for X < 0,
should be based on that measure-zero subset of the values of Y for which X**Y
happens to have a definable real value. That seems
to me far less kind to programs and to programmers than forcing them to
specify explicitly what they meant in the first place. If this case is
so important in practice, perhaps what
languages for points and intervals need instead is
powfrac(x,n,m)
defined to be
(x**n)**(1/m)
for integers n and m, that
can confidently be applied to floating-point arguments, even negative x**n,
returning then negative answers for odd m and exceptions for even m.
And languages like
C and Java that omit x**n and attempt to overload it onto x**y need to
rectify that mistake as well.
As my colleagues have discovered, the extension of IEEE 754 default
nonstop exception semantics to interval methods is fraught with unexpected
difficulties. The variety of responses that have been proposed or even
implemented can't all be correct in all situations.
The IEEE 754 design goal to optimize the ability of numerical
software experts to write efficient
bullet-resistant code for easy use by non-experts is still a controversial one,
involving, as it does, judgments about expertise, efficiency, robustness,
and ease of use.
More information about the Numeric-interest
mailing list