Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

arb_hypgeom_2f1 output nan #451

Open
xgluo opened this issue Jul 13, 2023 · 1 comment
Open

arb_hypgeom_2f1 output nan #451

xgluo opened this issue Jul 13, 2023 · 1 comment

Comments

@xgluo
Copy link

xgluo commented Jul 13, 2023

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
Copy link
Author

xgluo commented Jul 14, 2023

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant