arbmpfr.c

#include "arb.h"
#include "acb_hypgeom.h"

void arb_root_ui_algebraic(arb_t res, const arb_t x, ulong k, slong prec);

void
_arb_pow(arb_t res, const arb_t x, const arb_t y, slong prec)
{
    slong rootlim;

    rootlim = FLINT_BIT_COUNT(prec);

    /* slow detection of exact cases */
    if (!arb_is_zero(x) && !arb_is_int(y) && arb_bits(x) < 0.1 * prec && arb_bits(y) < 0.1 * prec &&
        arb_is_exact(y) && arf_is_int_2exp_si(arb_midref(y), -rootlim))
    {
        arb_t t, u;
        fmpz_t m, e;

        arb_init(t);
        arb_init(u);

        fmpz_init(m);
        fmpz_init(e);

        arf_get_fmpz_2exp(m, e, arb_midref(y));
        arb_root_ui_algebraic(t, x, UWORD(1) << fmpz_get_ui(e), prec);
        arb_pow_fmpz(t, t, m, prec);
        arb_set(res, t);

        arb_clear(t);
        arb_clear(u);
    }
    else
    {
        arb_pow(res, x, y, prec);
    }
}

void
_arb_agm(arb_t res, const arb_t x, const arb_t y, slong prec)
{
    if (arb_is_zero(x) || arb_is_zero(y))
    {
        arb_zero(res);
    }
    else
    {
        arb_agm(res, x, y, prec);
    }
}

void
_arb_csch(arb_t res, const arb_t x, slong prec)
{
    arb_sinh(res, x, prec);
    arb_inv(res, res, prec);
}

void
_arb_sech(arb_t res, const arb_t x, slong prec)
{
    arb_cosh(res, x, prec);
    arb_inv(res, res, prec);
}

void
_arb_csc(arb_t res, const arb_t x, slong prec)
{
    arb_sin(res, x, prec);
    arb_inv(res, res, prec);
}

void
_arb_sec(arb_t res, const arb_t x, slong prec)
{
    arb_cos(res, x, prec);
    arb_inv(res, res, prec);
}

void
_arb_log10(arb_t res, const arb_t x, slong prec)
{
    arb_t t;

    /* slow detection of exact cases */
    if (arb_is_positive(x) && arb_is_int(x) && arf_cmp_2exp_si(arb_midref(x), prec) < 0)
    {
        slong flog, clog;
        fmpz_t t;
        fmpz_init(t);
        arf_get_fmpz(t, arb_midref(x), ARF_RND_DOWN);
        flog = fmpz_flog_ui(t, 10);
        clog = fmpz_clog_ui(t, 10);

        if (flog == clog)
        {
            arb_set_si(res, flog);
            fmpz_clear(t);
            return;
        }

        fmpz_clear(t);
    }

    arb_init(t);
    arb_log(res, x, prec);
    arb_const_log10(t, prec);
    arb_div(res, res, t, prec);
    arb_clear(t);
}

void
_arb_log2(arb_t res, const arb_t x, slong prec)
{
    if (arb_is_positive(x) && arb_bits(x) == 1)
    {
        arb_set_fmpz(res, ARF_EXPREF(arb_midref(x)));
        arb_sub_ui(res, res, 1, prec);
    }
    else
    {
        arb_t t;
        arb_init(t);
        arb_log(res, x, prec);
        arb_const_log2(t, prec);
        arb_div(res, res, t, prec);
        arb_clear(t);
    }
}

void
_arb_exp2(arb_t res, const arb_t x, slong prec)
{
    arb_t t;
    arb_init(t);
    arb_set_ui(t, 2);
    arb_pow(res, t, x, prec);
    arb_clear(t);
}

void
_arb_exp10(arb_t res, const arb_t x, slong prec)
{
    arb_t t;
    arb_init(t);
    arb_set_ui(t, 10);
    arb_pow(res, t, x, prec);
    arb_clear(t);
}

void
_arb_cbrt(arb_t res, const arb_t x, slong prec)
{
    if (arf_sgn(arb_midref(x)) < 0)
    {
        arb_neg(res, x);
        arb_root_ui(res, res, 3, prec);
        arb_neg(res, res);
    }
    else
    {
        arb_root_ui(res, x, 3, prec);
    }
}

void
_arb_jn(arb_t res, ulong n, const arb_t x, slong prec)
{
    acb_t t, u;
    acb_init(t);
    acb_init(u);
    acb_set_ui(t, n);
    acb_set_arb(u, x);
    acb_hypgeom_bessel_j(t, t, u, prec);
    arb_set(res, acb_realref(t));
    acb_clear(t);
    acb_clear(u);
}

void
_arb_yn(arb_t res, ulong n, const arb_t x, slong prec)
{
    acb_t t, u;
    acb_init(t);
    acb_init(u);
    acb_set_ui(t, n);
    acb_set_arb(u, x);
    acb_hypgeom_bessel_y(t, t, u, prec);
    arb_set(res, acb_realref(t));
    acb_clear(t);
    acb_clear(u);
}

void
_arb_j0(arb_t res, const arb_t x, slong prec)
{
    _arb_jn(res, 0, x, prec);
}

void
_arb_j1(arb_t res, const arb_t x, slong prec)
{
    _arb_jn(res, 1, x, prec);
}

void
_arb_j2(arb_t res, const arb_t x, slong prec)
{
    _arb_jn(res, 2, x, prec);
}

void
_arb_j100(arb_t res, const arb_t x, slong prec)
{
    _arb_jn(res, 100, x, prec);
}

void
_arb_y0(arb_t res, const arb_t x, slong prec)
{
    _arb_yn(res, 0, x, prec);
}

void
_arb_y1(arb_t res, const arb_t x, slong prec)
{
    _arb_yn(res, 1, x, prec);
}

void
_arb_y2(arb_t res, const arb_t x, slong prec)
{
    _arb_yn(res, 2, x, prec);
}

void
_arb_y100(arb_t res, const arb_t x, slong prec)
{
    _arb_yn(res, 100, x, prec);
}

void
_arb_li2(arb_t res, const arb_t x, slong prec)
{
    acb_t t, u;
    acb_init(t);
    acb_init(u);
    acb_set_ui(t, 2);
    acb_set_arb(u, x);
    acb_polylog(t, t, u, prec);
    arb_set(res, acb_realref(t));
    acb_clear(t);
    acb_clear(u);
}

void
_arb_ai(arb_t res, const arb_t x, slong prec)
{
    acb_t t;
    acb_init(t);
    acb_set_arb(t, x);
    acb_hypgeom_airy(t, NULL, NULL, NULL, t, prec);
    arb_set(res, acb_realref(t));
    acb_clear(t);
}

void
_arb_eint(arb_t res, const arb_t x, slong prec)
{
    acb_t t;
    acb_init(t);
    acb_set_arb(t, x);
    acb_hypgeom_ei(t, t, prec);
    arb_set(res, acb_realref(t));
    acb_clear(t);
}

void
_arb_erf(arb_t res, const arb_t x, slong prec)
{
    acb_t t;
    acb_init(t);
    acb_set_arb(t, x);
    acb_hypgeom_erf(t, t, prec);
    arb_set(res, acb_realref(t));
    acb_clear(t);
}

void
_arb_erfc(arb_t res, const arb_t x, slong prec)
{
    acb_t t;
    acb_init(t);
    acb_set_arb(t, x);
    acb_hypgeom_erfc(t, t, prec);
    arb_set(res, acb_realref(t));
    acb_clear(t);
}

void
_arb_lgamma(arb_t res, const arb_t x, slong prec)
{
    acb_t t;
    acb_init(t);
    acb_set_arb(t, x);
    acb_lgamma(t, t, prec);
    arb_set(res, acb_realref(t));
    acb_clear(t);
}

int
_mpfr_lgamma(mpfr_t res, const mpfr_t x, mpfr_rnd_t rnd)
{
    int tmp;
    return mpfr_lgamma(res, &tmp, x, rnd);
}

int
_mpfr_j2(mpfr_t res, const mpfr_t x, mpfr_rnd_t rnd)
{
    return mpfr_jn(res, 2, x, rnd);
}

int
_mpfr_j100(mpfr_t res, const mpfr_t x, mpfr_rnd_t rnd)
{
    return mpfr_jn(res, 100, x, rnd);
}

int
_mpfr_y2(mpfr_t res, const mpfr_t x, mpfr_rnd_t rnd)
{
    return mpfr_yn(res, 2, x, rnd);
}

int
_mpfr_y100(mpfr_t res, const mpfr_t x, mpfr_rnd_t rnd)
{
    return mpfr_yn(res, 100, x, rnd);
}


#define DEFINE_FUNCTION(fname, funcall) \
int \
fname(mpfr_t res, const mpfr_t x, mpfr_rnd_t rnd) \
{ \
    arb_t t, u; \
    slong prec, wp; \
    int ret; \
 \
    arb_init(t); \
    arb_init(u); \
 \
    prec = mpfr_get_prec(res); \
 \
    arf_set_mpfr(arb_midref(t), x); \
 \
    for (wp = prec + 20; ; wp *= 2) \
    { \
        funcall; \
 \
        if (arb_can_round_mpfr(u, prec, rnd)) \
        { \
            break; \
        } \
 \
        if (wp > 1000 * prec + 1000000) \
        { \
            printf("high precision... (prec = %ld, wp = %ld)\n", prec, wp); \
            arb_printd(t, 200); printf("\n\n"); \
            arb_printd(u, 200); printf("\n\n"); \
            abort(); \
        } \
    } \
 \
    ret = arf_get_mpfr(res, arb_midref(u), rnd); \
 \
    arb_clear(t); \
    arb_clear(u); \
 \
    return ret; \
}

#define DEFINE_FUNCTION2(fname, funcall) \
int \
fname(mpfr_t res, const mpfr_t x, const mpfr_t y, mpfr_rnd_t rnd) \
{ \
    arb_t t, u, v; \
    slong prec, wp; \
    int ret; \
 \
    arb_init(t); \
    arb_init(u); \
    arb_init(v); \
 \
    prec = mpfr_get_prec(res); \
 \
    arf_set_mpfr(arb_midref(t), x); \
    arf_set_mpfr(arb_midref(v), y); \
 \
    for (wp = prec + 20; ; wp *= 2) \
    { \
        funcall; \
 \
        if (arb_can_round_mpfr(u, prec, rnd)) \
        { \
            break; \
        } \
 \
        if (wp > 1000 * prec + 1000000) \
        { \
            printf("high precision... (prec = %ld, wp = %ld)\n", prec, wp); \
            arb_printd(t, 200); printf("\n\n"); \
            arb_printd(v, 200); printf("\n\n"); \
            arb_printd(u, 200); printf("\n\n"); \
            abort(); \
        } \
    } \
 \
    ret = arf_get_mpfr(res, arb_midref(u), rnd); \
 \
    arb_clear(t); \
    arb_clear(u); \
    arb_clear(v); \
 \
    return ret; \
}

DEFINE_FUNCTION2(arb_mpfr_pow, _arb_pow(u, t, v, wp))
DEFINE_FUNCTION2(arb_mpfr_atan2, arb_atan2(u, t, v, wp))
DEFINE_FUNCTION2(arb_mpfr_hypot, arb_hypot(u, t, v, wp))
DEFINE_FUNCTION2(arb_mpfr_agm, _arb_agm(u, t, v, wp))
DEFINE_FUNCTION(arb_mpfr_ai, _arb_ai(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_j0, _arb_j0(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_j1, _arb_j1(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_j2, _arb_j2(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_j100, _arb_j100(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_y0, _arb_y0(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_y1, _arb_y1(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_y2, _arb_y2(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_y100, _arb_y100(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_erfc, _arb_erfc(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_erf, _arb_erf(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_zeta, arb_zeta(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_digamma, arb_digamma(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_lgamma, _arb_lgamma(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_gamma, arb_gamma(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_li2, _arb_li2(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_eint, _arb_eint(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_log1p, arb_log1p(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_expm1, arb_expm1(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_asinh, arb_asinh(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_acosh, arb_acosh(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_atanh, arb_atanh(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_sech, _arb_sech(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_csch, _arb_csch(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_coth, arb_coth(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_sinh, arb_sinh(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_cosh, arb_cosh(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_tanh, arb_tanh(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_asin, arb_asin(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_acos, arb_acos(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_atan, arb_atan(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_csc, _arb_csc(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_sec, _arb_sec(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_cot, arb_cot(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_log10, _arb_log10(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_log2, _arb_log2(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_log, arb_log(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_exp10, _arb_exp10(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_exp2, _arb_exp2(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_exp, arb_exp(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_rec_sqrt, arb_rsqrt(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_cbrt, _arb_cbrt(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_sqrt, arb_sqrt(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_sin, arb_sin(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_cos, arb_cos(u, t, wp))
DEFINE_FUNCTION(arb_mpfr_tan, arb_tan(u, t, wp))

void
mpfr_randinit(mpfr_t x, flint_rand_t state, slong prec)
{
    arf_t t;
    arf_init(t);
    mpfr_init2(x, prec);
    arf_randtest(t, state, prec, 0);
    arf_get_mpfr(x, t, MPFR_RNDN);
    arf_clear(t);
}

typedef int (*mpfrfunc1)(mpfr_t, const mpfr_t, mpfr_rnd_t);
typedef int (*mpfrfunc2)(mpfr_t, const mpfr_t, const mpfr_t, mpfr_rnd_t);

void testfunc1(const char * msg, mpfrfunc1 mpfrf, mpfrfunc1 arbf,
    slong reps, slong maxprec, slong minexp, slong maxexp)
{
    mpfr_t x, y1, y2;
    int r1, r2, i;
    mpfr_rnd_t rnd;
    slong iter, xprec, yprec;
    flint_rand_t state;

    printf("%s ", msg);
    fflush(stdout);

    flint_randinit(state);

    mpfr_set_emin(MPFR_EMIN_MIN);
    mpfr_set_emax(MPFR_EMAX_MAX);

    for (iter = 0; iter < reps; iter++)
    {
        if (iter % (reps / 100) == 0)
        {
            printf(".");
            fflush(stdout);
        }

        if (iter < reps / 2)
        {
            xprec = n_randtest(state) % 100;
            xprec = FLINT_MAX(xprec, 2);
            yprec = n_randtest(state) % 100;
            yprec = FLINT_MAX(yprec, 2);
        }
        else
        {
            xprec = n_randtest(state) % maxprec;
            xprec = FLINT_MAX(xprec, 2);
            yprec = n_randtest(state) % maxprec;
            yprec = FLINT_MAX(yprec, 2);
        }

        mpfr_randinit(x, state, xprec);
        mpfr_randinit(y1, state, yprec);
        mpfr_randinit(y2, state, yprec);

        if (iter < reps / 2)
        {
            mpfr_mul_2ui(x, x, n_randtest(state) % maxexp, MPFR_RNDN);
            mpfr_div_2ui(x, x, n_randtest(state) % minexp, MPFR_RNDN);
        }
        else
        {
            mpfr_mul_2ui(x, x, n_randtest(state) % 5, MPFR_RNDN);
            mpfr_div_2ui(x, x, n_randtest(state) % 5, MPFR_RNDN);
        }

        switch (n_randint(state, 5))
        {
            case 0: rnd = MPFR_RNDD; break;
            case 1: rnd = MPFR_RNDU; break;
            case 2: rnd = MPFR_RNDZ; break;
            case 3: rnd = MPFR_RNDA; break;
            case 4: rnd = MPFR_RNDN; break;
            default: rnd = MPFR_RNDN;
        }

        r1 = mpfrf(y1, x, rnd);

        if (r1 > 0) r1 = 1;
        if (r1 < 0) r1 = -1;

        if (mpfr_number_p(y1))
        {
            r2 = arbf(y2, x, rnd);

            if (r1 != r2 || mpfr_cmp(y1, y2))
            {
                printf("FAIL! (iter %ld)\n", iter);
                printf("xprec = %ld, yprec = %ld, rnd = %s\n", xprec, yprec, mpfr_print_rnd_mode(rnd));
                printf("x = "); mpfr_out_str(stdout, 10, 0, x, MPFR_RNDN); printf("\n");
                printf("y1 = "); mpfr_out_str(stdout, 10, 0, y1, MPFR_RNDN); printf("\n");
                printf("y2 = "); mpfr_out_str(stdout, 10, 0, y2, MPFR_RNDN); printf("\n");
                printf("r1 = %d\n", r1);
                printf("r2 = %d\n", r2);
                abort();
            }
        }

        mpfr_clear(x);
        mpfr_clear(y1);
        mpfr_clear(y2);
    }

    flint_randclear(state);
    flint_cleanup();
    printf("PASS\n");
}

void testfunc2(const char * msg, mpfrfunc2 mpfrf, mpfrfunc2 arbf,
    slong reps, slong maxprec, slong minexp, slong maxexp)
{
    mpfr_t x, w, y1, y2;
    int r1, r2, i;
    mpfr_rnd_t rnd;
    slong iter, xprec, wprec, yprec;
    flint_rand_t state;

    printf("%s ", msg);
    fflush(stdout);

    flint_randinit(state);

    mpfr_set_emin(MPFR_EMIN_MIN);
    mpfr_set_emax(MPFR_EMAX_MAX);

    for (iter = 0; iter < reps; iter++)
    {
        if (iter % (reps / 100) == 0)
        {
            printf(".");
            fflush(stdout);
        }

        if (iter < reps / 2)
        {
            xprec = n_randtest(state) % 100;
            xprec = FLINT_MAX(xprec, 2);
            wprec = n_randtest(state) % 100;
            wprec = FLINT_MAX(xprec, 2);
            yprec = n_randtest(state) % 100;
            yprec = FLINT_MAX(yprec, 2);
        }
        else
        {
            xprec = n_randtest(state) % maxprec;
            xprec = FLINT_MAX(xprec, 2);
            wprec = n_randtest(state) % maxprec;
            wprec = FLINT_MAX(wprec, 2);
            yprec = n_randtest(state) % maxprec;
            yprec = FLINT_MAX(yprec, 2);
        }

        mpfr_randinit(x, state, xprec);
        mpfr_randinit(w, state, wprec);
        mpfr_randinit(y1, state, yprec);
        mpfr_randinit(y2, state, yprec);

        if (iter < reps / 2)
        {
            mpfr_mul_2ui(x, x, n_randtest(state) % maxexp, MPFR_RNDN);
            mpfr_div_2ui(x, x, n_randtest(state) % minexp, MPFR_RNDN);
            mpfr_mul_2ui(w, w, n_randtest(state) % maxexp, MPFR_RNDN);
            mpfr_div_2ui(w, w, n_randtest(state) % minexp, MPFR_RNDN);
        }
        else
        {
            mpfr_mul_2ui(x, x, n_randtest(state) % 5, MPFR_RNDN);
            mpfr_div_2ui(x, x, n_randtest(state) % 5, MPFR_RNDN);
            mpfr_mul_2ui(w, w, n_randtest(state) % 5, MPFR_RNDN);
            mpfr_div_2ui(w, w, n_randtest(state) % 5, MPFR_RNDN);
        }

        switch (n_randint(state, 5))
        {
            case 0: rnd = MPFR_RNDD; break;
            case 1: rnd = MPFR_RNDU; break;
            case 2: rnd = MPFR_RNDZ; break;
            case 3: rnd = MPFR_RNDA; break;
            case 4: rnd = MPFR_RNDN; break;
            default: rnd = MPFR_RNDN;
        }

        r1 = mpfrf(y1, x, w, rnd);

        if (r1 > 0) r1 = 1;
        if (r1 < 0) r1 = -1;

        if (mpfr_number_p(y1))
        {
            r2 = arbf(y2, x, w, rnd);

            if (r1 != r2 || mpfr_cmp(y1, y2))
            {
                printf("FAIL! (iter %ld)\n", iter);
                printf("xprec = %ld, wprec = %ld, yprec = %ld, rnd = %s\n", xprec, wprec, yprec, mpfr_print_rnd_mode(rnd));
                printf("x = "); mpfr_out_str(stdout, 10, 0, x, MPFR_RNDN); printf("\n");
                printf("w = "); mpfr_out_str(stdout, 10, 0, w, MPFR_RNDN); printf("\n");
                printf("y1 = "); mpfr_out_str(stdout, 10, 0, y1, MPFR_RNDN); printf("\n");
                printf("y2 = "); mpfr_out_str(stdout, 10, 0, y2, MPFR_RNDN); printf("\n");
                printf("r1 = %d\n", r1);
                printf("r2 = %d\n", r2);
                abort();
            }
        }

        mpfr_clear(x);
        mpfr_clear(w);
        mpfr_clear(y1);
        mpfr_clear(y2);
    }

    flint_randclear(state);
    flint_cleanup();
    printf("PASS\n");
}

int main()
{

    /* these fail with mpfr 3.1.3 */
    /*
    testfunc1("sqrt", mpfr_sqrt, arb_mpfr_sqrt, 100000, 2000, 1000, 10000);
    testfunc1("zeta", mpfr_zeta, arb_mpfr_zeta, 100000, 500, 10, 10);
    testfunc1("j0", mpfr_j0, arb_mpfr_j0, 1000000, 1000, 1000, 1000);
    testfunc1("j1", mpfr_j1, arb_mpfr_j1, 1000000, 1000, 1000, 1000);
    testfunc1("j2", _mpfr_j2, arb_mpfr_j2, 100000, 1000, 1000, 1000);
    testfunc1("j100", _mpfr_j100, arb_mpfr_j100, 100000, 1000, 1000, 1000);
    testfunc1("y0", mpfr_y0, arb_mpfr_y0, 100000, 1000, 1000, 1000);
    testfunc1("y1", mpfr_y1, arb_mpfr_y1, 100000, 1000, 1000, 1000);
    testfunc1("y2", _mpfr_y2, arb_mpfr_y2, 100000, 1000, 1000, 1000);
    testfunc1("y100", _mpfr_y100, arb_mpfr_y100, 100000, 1000, 1000, 1000);
    */

    /* reps, prec, minexp, maxexp */
    testfunc2("pow", mpfr_pow, arb_mpfr_pow, 1000000, 5000, 1000, 20);
    testfunc2("atan2", mpfr_atan2, arb_mpfr_atan2, 100000, 2000, 10000, 10000);
    testfunc2("hypot", mpfr_hypot, arb_mpfr_hypot, 100000, 10000, 10000, 10000);
    testfunc2("agm", mpfr_agm, arb_mpfr_agm, 100000, 10000, 10000, 10000);
    testfunc1("ai", mpfr_ai, arb_mpfr_ai, 100000, 1000, 100, 8);
    testfunc1("erfc", mpfr_erfc, arb_mpfr_erfc, 100000, 1000, 100, 7);
    testfunc1("erf", mpfr_erf, arb_mpfr_erf, 100000, 1000, 100, 7);
    testfunc1("digamma", mpfr_digamma, arb_mpfr_digamma, 100000, 500, 1000, 1000);
    testfunc1("gamma", mpfr_gamma, arb_mpfr_gamma, 100000, 500, 100, 50);
    testfunc1("lgamma", _mpfr_lgamma, arb_mpfr_lgamma, 100000, 500, 10000, 10000);
    testfunc1("li2", mpfr_li2, arb_mpfr_li2, 100000, 500, 1000, 1000);
    testfunc1("eint", mpfr_eint, arb_mpfr_eint, 100000, 2000, 1000, 20);
    testfunc1("expm1", mpfr_expm1, arb_mpfr_expm1, 100000, 2000, 1000, 10);
    testfunc1("log1p", mpfr_log1p, arb_mpfr_log1p, 100000, 2000, 10000, 10000);
    testfunc1("atan", mpfr_atan, arb_mpfr_atan, 100000, 2000, 1000, 1000);
    testfunc1("acosh", mpfr_acosh, arb_mpfr_acosh, 100000, 2000, 1000, 10000);
    testfunc1("asinh", mpfr_asinh, arb_mpfr_asinh, 100000, 2000, 1000, 10000);
    testfunc1("atanh", mpfr_atanh, arb_mpfr_atanh, 100000, 2000, 1000, 10000);
    testfunc1("tanh", mpfr_tanh, arb_mpfr_tanh, 100000, 2000, 1000, 15);
    testfunc1("coth", mpfr_coth, arb_mpfr_coth, 100000, 2000, 1000, 15);
    testfunc1("csch", mpfr_csch, arb_mpfr_csch, 100000, 2000, 1000, 60);
    testfunc1("sech", mpfr_sech, arb_mpfr_sech, 100000, 2000, 1000, 60);
    testfunc1("cosh", mpfr_cosh, arb_mpfr_cosh, 100000, 2000, 1000, 60);
    testfunc1("sinh", mpfr_sinh, arb_mpfr_sinh, 100000, 2000, 1000, 60);
    testfunc1("acos", mpfr_acos, arb_mpfr_acos, 100000, 1000, 1000, 10);
    testfunc1("asin", mpfr_asin, arb_mpfr_asin, 100000, 1000, 1000, 10);
    testfunc1("cot", mpfr_cot, arb_mpfr_cot, 100000, 2000, 1000, 10000);
    testfunc1("csc", mpfr_csc, arb_mpfr_csc, 100000, 2000, 1000, 10000);
    testfunc1("sec", mpfr_sec, arb_mpfr_sec, 100000, 2000, 1000, 10000);
    testfunc1("log10", mpfr_log10, arb_mpfr_log10, 100000, 5000, 1000, 10000);
    testfunc1("log2", mpfr_log2, arb_mpfr_log2, 100000, 1000, 1000, 10000);
    testfunc1("log", mpfr_log, arb_mpfr_log, 100000, 5000, 10000, 10000);
    testfunc1("exp10", mpfr_exp10, arb_mpfr_exp10, 100000, 5000, 1000, 60);
    testfunc1("exp2", mpfr_exp2, arb_mpfr_exp2, 100000, 5000, 1000, 60);
    testfunc1("exp", mpfr_exp, arb_mpfr_exp, 100000, 10000, 1000, 60);
    testfunc1("rec_sqrt", mpfr_rec_sqrt, arb_mpfr_rec_sqrt, 1000000, 1000, 10000, 10000);
    testfunc1("cbrt", mpfr_cbrt, arb_mpfr_cbrt, 1000000, 1000, 10000, 10000);
    testfunc1("cos", mpfr_cos, arb_mpfr_cos, 100000, 2000, 1000, 10000);
    testfunc1("tan", mpfr_tan, arb_mpfr_tan, 100000, 2000, 1000, 10000);
    testfunc1("sin", mpfr_sin, arb_mpfr_sin, 100000, 2000, 1000, 10000);
}