Integration produces a nan

[kellino]

Hi, I'm trying to integrate the function f(x) = x / (1 - x). This integrates just fine in, eg mpmath, with the output 9.2935262951108832. However, it produces a nan in arb.

I'm not a C programmer, so it's perfectly possible I've made some foolish/obvious mistake. Arb version 2.23.0 Operation System: MacOS Monterey version 12.5

Minimum working example

#include <acb.h>
#include <acb_calc.h>
#include <mag.h>

int f(acb_ptr res, const acb_t z, void *param, slong order, slong prec) {

  /* z / (1 - z) */
  acb_t one, one_min_z;
  acb_init(one);
  acb_set_d(one, 1);    // one = 1

  acb_init(one_min_z);
  acb_sub(one_min_z, one, z, prec); // one_min_z = 1 - z

  acb_div(res, z, one_min_z, prec); // res = z / (1-z)

  acb_clear(one);
  acb_clear(one_min_z);

  return 0;
}

int main(int argc, char *argv[]) {

  slong prec = 50;
  acb_t s, a, b;
  acb_init(a);
  acb_init(b);
  acb_init(s);

  acb_set_d(a, 0);
  acb_set_d(b, 1);

  mag_t tol;
  mag_init(tol);
  mag_set_ui_2exp_si(tol, 1, -prec);

  acb_calc_integrate(s, f, NULL, a, b, prec, tol, NULL, prec);

  acb_printd(s, prec);

  acb_clear(a);
  acb_clear(b);
  acb_clear(s);
  mag_clear(tol);

  return 0;
}

[fredrik-johansson]

This integral is infinite.

Even if it were finite, e.g. if you take x/sqrt(1-x) instead, Arb would not be able to handle the blow-up singularity directly.

[kellino]

Thanks for the quick response. What is it that mpmath does differently that is produces an answer? Even if I set the options to limit the number of evaluations to 5 and insert an acb_printd statement on res, then it produces nans immediately, which is (to me) unexpected.

acb_div(res, z, one_min_z, prec); // res = z / (1-z)
acb_printd(res, prec);
printf("\n");

(nan + 0j) +/- (+inf, 0j) (0.50000000279396772384643554687500000000000000000000 + 0j) +/- (0.500, 0j) (nan + 0j) +/- (+inf, 0j) (0.33333333333333303727386009995825588703155517578125 + 0j) +/- (16.7, 13.3j) (nan + nanj) +/- (+inf, +infj) (0.98945680897518251128985866671428084373474121093750 + 0j) +/- (5.27e-15, 0j)

[fredrik-johansson]

mpmath is not rigorous and just returns some random meaningless number when the integral is infinite.

[kellino]

That's good to know! Thank you.