arb integration question

[Hansjoerg Seybold]

Hi,

Not sure if this is the right place, to post my flint/arb user question, but I could not find a more appropriate channel for asking.

I am trying to integrate a function involving a fractional power law over the interval [0, 1] for real x.

Here an example: \int_0^b x**a = b**(1 + 1/a)/(1 + 1/a), 0 < a <= 1

However, I get (nan + nanj)  +/-  (+inf, +infj) as a result except for a=0.5 and a=1

I understand that the function has a branch-cut along the negative real axis, which may impact the interval arithmetic.

Thus my question here is

(a) How should I approach such a degenerate integration case

(b) is there is a way to circumvent this degenerate case and what would be the best way to do so using the tools provided by flint

Below the source code of a small example program.

Thanks a lot,

                   Hansjoerg

------------------------------------------------------------------------------------------------------------------------------------------------

CODE

#include <flint/acb.h>
#include <flint/acb_calc.h>
#include <math.h>
#include <stdlib.h>
#define NDIGITS 16
#define PRECISION 53
====================================================
/* Parameters for integrating fn */
typedef struct {
    acb_t a;
} params_t;
====================================================
int
fn (acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
    params_t* p = (params_t*) param;

    if (order > 1)
        flint_abort();  /* Would be needed for Taylor method. */

    acb_inv(res, p->a, prec);
    acb_pow(res, z, res, prec);

    return 0;
}

====================================================

 int main (int argc, char **argv)
{
    double a = atof(argv[1]);
    // Integration limits
    double from = 0.0;
    double to = 1.0;
    double expected = pow(to, (1 + 1/a)) / (1 + 1/a);

    acb_t _result, _from, _to, _expected;
    slong prec, goal;
    mag_t tol;
    acb_calc_integrate_opt_t options;
    params_t par;

    acb_init(_from);
    acb_init(_to);
    acb_init(_result);
    acb_init(_expected);
    mag_init(tol);

    acb_set_d(_from, from);
    acb_set_d(_to, to);
    acb_set_d(_expected, expected);

    prec = PRECISION;
    acb_calc_integrate_opt_init(options);
    goal = (slong)prec;
    options->verbose = 2;

    mag_set_ui_2exp_si(tol, 1, -prec);

    printf("\n#========================\n");

    acb_init(par.a);
    acb_set_d(par.a, a);

    printf("# Integrating F(x) = x**a ");
    printf("# from = "); acb_printd(_from, NDIGITS); printf(", ");
       
    printf("# to = "); acb_printd(_to, NDIGITS); printf("\n");
    acb_calc_integrate(_result, fn, &par,
                       _from, _to, goal, tol, options, prec);

    printf("# expected = ");
    acb_printd(_expected, NDIGITS);
    printf("\n");
    printf("# result = ");
    acb_printd(_result, NDIGITS);
    printf("\n");

    acb_clear(par.a);
    acb_clear(_from);
    acb_clear(_to);
    acb_clear(_result);
    acb_clear(_expected);
    mag_clear(tol);

    return 0;
}