arb_hypgeom_2f1 output nan
xgluo opened this issue · 1 comments
xgluo commented
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
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.