arb_pow could do much better for <a +- a> ^ <b +- c> with b >= c > 0 and a > 0

[postmath]

Consider arb_pow(z, x, y, prec) with x = <a +- a> and y nonnegative. This works as I would expect for x=0, and also for exact half integers y, including y = 0. But if a > 0 and y is not an exact half integer, then arb_pow returns a non-finite answer (because it takes the logarithm of x, then multiplies by y, then takes the exponential; but we can't represent the half-infinite interval that the logarithm requires), whereas there is a finite answer. I propose to do roughly this, after all existing special cases have been done:

    if(arf_cmp_si(arb_midref(x), 0) > 0
        && arf_cmpabs_mag(arb_midref(x), arb_radref(x)) == 0
        && arb_is_nonnegative(y)) {
        /* x is an interval of the form <a +- a>, y is an interval of the form <b +- c> with b >=
          c. 

          (x, y) -> x^y is nondecreasing in x for y >= 0, and it is 0 for x=0 and y>0.  The lower
          bound of the interval for x is 0, and we have points with y>0 (because the case
          arb_is_zero(y) is excluded above). So the lowerbound of the result is 0, and the
          upperbound can be found somewhere at (2*a)^y.

          We compute this upperbound by computing z^y with z = < 3*a/2 +- a/2 >, which has the same
          upperbound (2*a) as x, then keeping the upper bound of this result and setting the lower
          bound to 0. This necessarily drops the precision to MAG_BITS, so we might as well do that
          from the start. */
        prec = (prec > MAG_BITS) ? MAG_BITS : prec;
        arf_mul_ui(arb_midref(z), arb_midref(x), 3, prec, ARF_RND_UP);
        arf_mul_2exp_si(arb_midref(z), arb_midref(z), -1);
        mag_mul_2exp_si(arb_radref(z), arb_radref(x), -1);

        arb_pow(z, z, y);

        /* Now we keep the upper bound of z (we may need to round up) and set the lower bound to
           zero. That is, if currently, z = < zc +/- zr >, we set it to < (zc + zr)/2 +/- (zc +
           zr)/2 >. */
        arb_get_ubound_arf(arb_midref(z), z, prec);
        arf_mul_2exp_si(arb_midref(z), arb_midref(z), -1);
        arf_get_mag(arb_radref(z), arb_midref(z));
        /* Note the above is inexact (rounding up), so need to update arb_midref(z) to match again */
        arf_set_mag(arb_midref(z), arb_radref(z));
    }

[fredrik-johansson]

Agreed, but arb_pow ought to have such a special case for any x containing 0, not necessarily of the form [a +- a].

[postmath]

Hmm, but if there are negative numbers in x, and y is not an exact integer, the result isn't real, is it? So in that case I think the current answer of nan is correct.

[fredrik-johansson]

Ah yes, you're right of course. So this really is a special special case.

The more general case would be relevant in acb_pow and acb_pow_arb.

[postmath]

Agreed, acb_pow and acb_pow_arb are a bit more complicated. I guess for acb_pow_arb you could reason that, with 0 \in x and y real:

x^y = exp(y * ln(x))
    = exp(y * ([-infinity, max(abs(x))] + i * acb_arg(x)))
    = [0, exp(max(y) * max(abs(x)))] * exp(i * y * acb_arg(x))

I suppose we could at least do a case distinction there between x being contained in any quadrants, vs containing 0 in its interior, to compute that second factor.

The acb_pow case would be another step messier. Would you be opposed to merging just this change to arb_pow for now? If not I can prepare a pull request next week.

[fredrik-johansson]

Patching arb_pow is fine.