Interval enclosure of X**Y for negative x in X

[David G Hough at validgh]

Prof. Neumaier responded to reliable-computing as follows to
my note about interval X**Y:

> Date: Tue, 31 Mar 1998 10:55:52 +0200
> From: Arnold Neumaier <neumacma.univie.ac.at>
> Subject: Re: Interval enclosure of X**Y for negative x in X
> To: reliable_computingainterval.usl.edu
> 
> David Hough wrote,
> 
> >>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.<<
> 
> This is precisely what is needed in practice for interval applications.
> Fractional powers of negative numbers are very useful (though one could
> get by without them, at the expense of programming effort). 
> 
> The main use
> in our global optimization context is to draw conclusions about ranges
> from constraints such as sum(x_i^p) in S, given x_i in X_i.
> In the simplest case we have sums consisting of a single term only, for
> example 
> x^3 in S,
> and I want to be able to conclude that
> x in S^(1/3)
> without having to split S into its positive and negative part.
> And of course I do not want to lose solutions of x^3=-1 just 
> because (-1)^(1/3) results in NaN or NaI or (even worse) empty!
> 
> There are also applications in integer programming via constraint
> satisfaction techniques. Consider the equation
> 
> (-1)^x*y = z
> 
> with x in [1:50] integral, y in [0,1], z in [1,2].
> 
> Here x=even, y=z=1 are the only solutions. I think this is
> representative for more complicated integer constraints that 
> enough users are likely to employ in their models when nonlinear 
> mixed integer programming software becomes more widely used.
> 
> The correct interval evaluation is (-1)^[1,50]=[-1,1], and this
> gives enough information to deduce from the equation that y=z=1,
> as is done in constraint propagation techniques.
> 
> Defining (-1)^[1,50]=empty gives a containment failure. And
> defining (-1)^[1,50]=NaI gives useless information; one is
> forced to split the integer variable into 50 pieces until a numerical
> enclosure is found. This is inefficient, and it becomes prohibitive 
> if there are a number of such constraints for different integer 
> variables, since the work then grows exponentially. 
> 
> Failure to define x**y properly when x<0 makes the power considerably
> more difficult to use in such cases, while a proper definition doesn't 
> harm users that only need x**y for x>=0.
> 
> >>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.<<
> The right definition should be x**q = (-1)**n * |x|**q if x is negative
> and q is a rational number that can be written as quotient of n/d with
> odd d. This respects the wish of what >>Many programmers might prefer<<,
> and the needs of users of interval arithmetic.
> 
> 
> Arnold Neumaier
>

When he writes

> The right definition should be x**q = (-1)**n * |x|**q if x is negative
> and q is a rational number that can be written as quotient of n/d with
> odd d.

I think the underlying issue is that we are thinking of different functions.
I think of x**y as corresponding to the multiple-valued complex function,
analytic in places if we choose a principal branch of log, defined by
exp(y**log(x)).    With the usual choice of log, then e.g.

(-1)**(2/3) is exp(+-(2/3)*pi*i), 

the sign perhaps depending on the sign attributed to the zero imaginary part.
In contrast Prof. Neumaier's definition is based upon a different choice
of log branch to yield 1.   Since the log branch to choose depends upon 
n and d, it doesn't make much sense as a restriction of a branch of a
complex analytic function, 
and has to be understood as another kind of function defined for rational
numbers, which is just the function one needs for mixed integer programming.

I think it's a mistake (in language design) to use one notation ** to
represent two rather different functions.   I decided to check Sun's f77
and also g77 with a test program

        print *,(-1.0)**(1.0/3.0)
        print *,cmplx(-1.0,0.0)**cmplx(1.0/3.0,0.0)
        print *,(-1.0)**(2.0/3.0)
        print *,cmplx(-1.0,0.0)**cmplx(2.0/3.0,0.0)
        end

and got

  NaN           
   (   0.500000,   0.866025)
  NaN           
   (  -0.500000,   0.866025)
 Note: IEEE floating-point exception flags raised: 
    Inexact;  Invalid Operation; 
 See the Numerical Computation Guide, ieee_flags(3M) 

as I'd expect.   Probably both complex pow functions were ultimately
derived from, or at least influenced by,
the same Bell Labs f77 implementation also enshrined in the
f2c run-time libraries, but I think they represent the semantics expected
by most programmers who have any expectation at all.   

So as part of the shift in thinking that interval arithmetic entails, the
languages or libraries need to define new functions that incorporate
rational-power semantics, rather than overload the exponential notation
with different explanations in different contexts.