Skip to content

arb_hypgeom_2f1 output nan #451

@xgluo

Description

@xgluo

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

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type
    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions