flintlib / arb

Arb has been merged into FLINT -- use https://github.com/flintlib/flint/ instead
http://arblib.org/
GNU Lesser General Public License v2.1
455 stars 137 forks source link

arb_hypgeom_2f1 output nan #451

Open xgluo opened 1 year ago

xgluo commented 1 year ago

I tried to run the following script. With valid inputs, the function "arb_hypgeom_2f1" does not compute properly but output nan. Any ideas why this is the case?

Code:

#include "arb_hypgeom.h"

int main()
{
    int prec = 200;

    /* compute paramaters  */
    arb_t mu_1, b_1, d_1, mu_2, b_2, d_2;
    arb_init(mu_1);
    arb_init(b_1);
    arb_init(d_1);
    arb_init(mu_2);
    arb_init(b_2);
    arb_init(d_2);
    arb_set_d(mu_1, 3.1);
    arb_set_d(b_1, 9.2);
    arb_set_d(d_1, 8.8);
    arb_set_d(mu_2, 1e-5);
    arb_set_d(b_2, 9);
    arb_set_d(d_2, 8.8);

    arb_t r_1, r_2, gamma_1, gamma_2, omega;
    arb_init(r_1);
    arb_init(r_2);
    arb_init(gamma_1);
    arb_init(gamma_2);
    arb_init(omega);

    arb_div(r_1, mu_1, b_1, prec);
    arb_div(r_2, mu_2, b_2, prec);
    arb_sub(gamma_1, b_1, d_1, prec);
    arb_sub(gamma_2, b_2, d_2, prec);

    arb_t a, b, c;
    arb_init(a);
    arb_init(b);
    arb_init(c);

    arb_t one, two, s_2;
    arb_init(one);
    arb_init(two);
    arb_init(s_2);
    arb_set_d(one, 1);
    arb_set_d(two, 2);
    arb_set_d(s_2, 1 - 1e-12);

    arb_t temp1, temp2;
    arb_init(temp1);
    arb_init(temp2);

    arb_mul(temp1, r_2, gamma_2, prec);
    arb_sub(temp1, gamma_1, temp1, prec);
    arb_div(temp1, temp1, two, prec);
    arb_mul(temp2, r_2, gamma_2, prec);
    arb_mul(temp2, temp2, b_1, prec);
    arb_pow(omega, temp1, two, prec);
    arb_add(omega, omega, temp2, prec);
    arb_sqrt(omega, omega, prec);
    arb_sub(omega, omega, temp1, prec);

    // this->a = omega / gamma_2;
    arb_div(a, omega, gamma_2, prec);

    // this->b = (omega + gamma_1) / gamma_2;
    arb_add(temp1, omega, gamma_1, prec);
    arb_div(b, temp1, gamma_2, prec);

    // this->c = a + b + 1 - r_2;
    arb_add(temp1, a, b, prec);
    arb_add(temp1, temp1, one, prec);
    arb_sub(c, temp1, r_2, prec);

    // clear
    arb_clear(temp1);
    arb_clear(temp2);

    arb_t F, z;
    arb_init(F);
    arb_init(z);
    arb_sub(z, s_2, one, prec);
    arb_mul(z, z, b_2, prec);
    arb_div(z, gamma_2, z, prec);
    arb_add(z, z, one, prec);

    /* compute paramaters ends  */

    printf("a = ");
    arb_printn(a, prec, 0);
    printf("\n");
    printf("b = ");
    arb_printn(b, prec, 0);
    printf("\n");
    printf("c = ");
    arb_printn(c, prec, 0);
    printf("\n");
    printf("z = ");
    arb_printn(z, prec, 0);
    printf("\n");

    arb_hypgeom_2f1(F, a, b, c, z, 0, prec); // this function does not evaluate properly
    printf("F = ");
    arb_printn(F, prec, 0);
    printf("\n");

    arb_clear(F);
    arb_clear(z);
    arb_clear(a);
    arb_clear(b);
    arb_clear(c);
    arb_clear(mu_1);
    arb_clear(b_1);
    arb_clear(d_1);
    arb_clear(mu_2);
    arb_clear(b_2);
    arb_clear(d_2);
    arb_clear(r_1);
    arb_clear(r_2);
    arb_clear(gamma_1);
    arb_clear(gamma_2);
    arb_clear(omega);
    arb_clear(one);
    arb_clear(two);
    arb_clear(s_2);

}

The output of the code is

a = [2.555524321768503210385577445633261842075283570379494615e-5 +/- 7.05e-60]
b = [2.0000255552432176850321038557744563326184207528357037949461 +/- 5.13e-59]
c = [3.0000499993753242589530057081556748526706794364338235486677 +/- 1.72e-59]
z = [-22222713825.8777617408682136116353946930165426890196513822582 +/- 4.83e-50]
F = nan
xgluo commented 1 year ago

For now, the workaround is to convert the parameters a, b, c back into double via

arb_set_d(a, arf_get_d(arb_midref(a), ARF_RND_NEAR));
arb_set_d(b, arf_get_d(arb_midref(b), ARF_RND_NEAR));
arb_set_d(c, arf_get_d(arb_midref(c), ARF_RND_NEAR));

which temporarily solves the problem and prints

a = 2.555524321768503121670747246785282413839013315737247467041015625e-5
b = 2.0000255552432175676358383498154580593109130859375
c = 3.000049999375324460970659856684505939483642578125
z = [-22222713825.8777617408682136116353946930165426890196513822582 +/- 4.83e-50]
F = [0.9994041175506549726202688333505541814603498257699103208789 +/- 5.20e-59]

The original code is quite complicated for me (without software engineering background) to pinpoint where exactly the problem is. How could the additional digits in the parameters result in NaNs? I would highly appreciate it if this could be investigated further, as the high precision and the speed of the arb library (arb_hypgeom in particular) help my research a lot.