diff --git a/arb_hypgeom.h b/arb_hypgeom.h
index 7735eb53..7f4e359e 100644
--- a/arb_hypgeom.h
+++ b/arb_hypgeom.h
@@ -85,6 +85,9 @@ void arb_hypgeom_erfi(arb_t res, const arb_t z, slong prec);
void _arb_hypgeom_erfi_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec);
void arb_hypgeom_erfi_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec);
+void arb_hypgeom_erfinv(arb_t res, const arb_t x, slong prec);
+void arb_hypgeom_erfcinv(arb_t res, const arb_t x, slong prec);
+
void arb_hypgeom_fresnel(arb_t res1, arb_t res2, const arb_t z, int normalized, slong prec);
void _arb_hypgeom_fresnel_series(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen, int normalized, slong len, slong prec);
void arb_hypgeom_fresnel_series(arb_poly_t s, arb_poly_t c, const arb_poly_t h, int normalized, slong len, slong prec);
diff --git a/arb_hypgeom/erfinv.c b/arb_hypgeom/erfinv.c
new file mode 100644
index 00000000..9e04b470
--- /dev/null
+++ b/arb_hypgeom/erfinv.c
@@ -0,0 +1,482 @@
+/*
+ Copyright (C) 2021 Fredrik Johansson
+
+ This file is part of Arb.
+
+ Arb is free software: you can redistribute it and/or modify it under
+ the terms of the GNU Lesser General Public License (LGPL) as published
+ by the Free Software Foundation; either version 2.1 of the License, or
+ (at your option) any later version. See .
+*/
+
+#include "arb_hypgeom.h"
+#include "arb_fmpz_poly.h"
+
+/* Actually an enclosure for |x| < 0.99. */
+static void
+arb_erfinv_approx_tiny(arb_t res, const arb_t x, slong prec)
+{
+ arb_t t;
+ mag_t err;
+ arb_init(t);
+ mag_init(err);
+ arb_get_mag(err, x);
+ mag_pow_ui(err, err, 3);
+ arb_const_sqrt_pi(t, prec);
+ arb_mul_2exp_si(t, t, -1);
+ arb_mul(res, x, t, prec);
+ arb_add_error_mag(res, err);
+ arb_clear(t);
+ mag_clear(err);
+}
+
+/* First terms of asymptotic expansion. */
+/* http://www.ams.org/journals/mcom/1976-30-136/S0025-5718-1976-0421040-7/S0025-5718-1976-0421040-7.pdf */
+static double erfinv_approx_big(double one_sub_x)
+{
+ double eta, l, y;
+
+ eta = -log(one_sub_x*sqrt(3.1415926535897932));
+ l = log(eta);
+ y = sqrt(eta - 0.5*l + (0.25*l - 0.5)/eta + (1/16.*l*l - 3/8.*l + 7./8)/(eta*eta)
+ + (l*l*l/48 - 7*l*l/32 + 17*l/16 - 107./48)/(eta*eta*eta)
+ + (l*l*l*l/128 - 23*l*l*l/192 + 29*l*l/32 - 31*l/8 + 1489/192.)/(eta*eta*eta*eta));
+ return y;
+}
+
+/* First terms of asymptotic expansion, when a double can overflow. */
+static void
+arb_erfinv_approx_huge(arb_t res, const arb_t one_minus_x, slong prec)
+{
+ arb_t eta, l, y;
+ fmpz c[5];
+ arb_ptr poly;
+ mag_t err;
+
+ arb_init(eta);
+ arb_init(l);
+ arb_init(y);
+ mag_init(err);
+ poly = _arb_vec_init(5);
+
+ arb_const_sqrt_pi(eta, prec);
+ arb_mul(eta, eta, one_minus_x, prec);
+ arb_log(eta, eta, prec);
+ arb_neg(eta, eta);
+ arb_log(l, eta, prec);
+ arb_mul_ui(y, eta, 12, prec);
+ arb_inv(y, y, prec);
+ arb_mul_2exp_si(poly + 0, l, -1);
+ arb_neg(poly + 0, poly + 0);
+ c[0] = -2 * 3; c[1] = 3;
+ _arb_fmpz_poly_evaluate_arb(poly + 1, c, 2, l, prec);
+ c[0] = 14 * 9; c[1] = -6 * 9; c[2] = 9;
+ _arb_fmpz_poly_evaluate_arb(poly + 2, c, 3, l, prec);
+ c[0] = -214 * 18; c[1] = 102 * 18; c[2] = -21 * 18; c[3] = 2 * 18;
+ _arb_fmpz_poly_evaluate_arb(poly + 3, c, 4, l, prec);
+ c[0] = 2978 * 54; c[1] = -1488 * 54; c[2] = 348 * 54; c[3] = -46 * 54; c[4] = 3 * 54;
+ _arb_fmpz_poly_evaluate_arb(poly + 4, c, 5, l, prec);
+ _arb_poly_evaluate(res, poly, 5, y, prec);
+ arb_add(res, res, eta, prec);
+ arb_sqrt(res, res, prec);
+
+ /* Should be an enclosure for 1-x < 1e-300 (but not proved). */
+ arb_get_mag(err, res);
+ mag_mul_2exp_si(err, err, -50);
+
+ arb_clear(eta);
+ arb_clear(l);
+ arb_clear(y);
+ mag_clear(err);
+ _arb_vec_clear(poly, 5);
+}
+
+/* Adapted from https://people.maths.ox.ac.uk/gilesm/codes/erfinv/
+ with permission */
+/* Only good up to about 1 - 1e-15 */
+/* Todo: use a good approximation for erfcinv */
+static double erfinv_approx(double x, double one_sub_x)
+{
+ double w, p;
+
+ w = -log(one_sub_x * (1.0 + x));
+
+ if (w < 6.250000)
+ {
+ w = w - 3.125000;
+ p = -3.6444120640178196996e-21;
+ p = -1.685059138182016589e-19 + p*w;
+ p = 1.2858480715256400167e-18 + p*w;
+ p = 1.115787767802518096e-17 + p*w;
+ p = -1.333171662854620906e-16 + p*w;
+ p = 2.0972767875968561637e-17 + p*w;
+ p = 6.6376381343583238325e-15 + p*w;
+ p = -4.0545662729752068639e-14 + p*w;
+ p = -8.1519341976054721522e-14 + p*w;
+ p = 2.6335093153082322977e-12 + p*w;
+ p = -1.2975133253453532498e-11 + p*w;
+ p = -5.4154120542946279317e-11 + p*w;
+ p = 1.051212273321532285e-09 + p*w;
+ p = -4.1126339803469836976e-09 + p*w;
+ p = -2.9070369957882005086e-08 + p*w;
+ p = 4.2347877827932403518e-07 + p*w;
+ p = -1.3654692000834678645e-06 + p*w;
+ p = -1.3882523362786468719e-05 + p*w;
+ p = 0.0001867342080340571352 + p*w;
+ p = -0.00074070253416626697512 + p*w;
+ p = -0.0060336708714301490533 + p*w;
+ p = 0.24015818242558961693 + p*w;
+ p = 1.6536545626831027356 + p*w;
+ }
+ else if (w < 16.000000)
+ {
+ w = sqrt(w) - 3.250000;
+ p = 2.2137376921775787049e-09;
+ p = 9.0756561938885390979e-08 + p*w;
+ p = -2.7517406297064545428e-07 + p*w;
+ p = 1.8239629214389227755e-08 + p*w;
+ p = 1.5027403968909827627e-06 + p*w;
+ p = -4.013867526981545969e-06 + p*w;
+ p = 2.9234449089955446044e-06 + p*w;
+ p = 1.2475304481671778723e-05 + p*w;
+ p = -4.7318229009055733981e-05 + p*w;
+ p = 6.8284851459573175448e-05 + p*w;
+ p = 2.4031110387097893999e-05 + p*w;
+ p = -0.0003550375203628474796 + p*w;
+ p = 0.00095328937973738049703 + p*w;
+ p = -0.0016882755560235047313 + p*w;
+ p = 0.0024914420961078508066 + p*w;
+ p = -0.0037512085075692412107 + p*w;
+ p = 0.005370914553590063617 + p*w;
+ p = 1.0052589676941592334 + p*w;
+ p = 3.0838856104922207635 + p*w;
+ }
+ else
+ {
+ w = sqrt(w) - 5.000000;
+ p = -2.7109920616438573243e-11;
+ p = -2.5556418169965252055e-10 + p*w;
+ p = 1.5076572693500548083e-09 + p*w;
+ p = -3.7894654401267369937e-09 + p*w;
+ p = 7.6157012080783393804e-09 + p*w;
+ p = -1.4960026627149240478e-08 + p*w;
+ p = 2.9147953450901080826e-08 + p*w;
+ p = -6.7711997758452339498e-08 + p*w;
+ p = 2.2900482228026654717e-07 + p*w;
+ p = -9.9298272942317002539e-07 + p*w;
+ p = 4.5260625972231537039e-06 + p*w;
+ p = -1.9681778105531670567e-05 + p*w;
+ p = 7.5995277030017761139e-05 + p*w;
+ p = -0.00021503011930044477347 + p*w;
+ p = -0.00013871931833623122026 + p*w;
+ p = 1.0103004648645343977 + p*w;
+ p = 4.8499064014085844221 + p*w;
+ }
+
+ return p*x;
+}
+
+/* floating-point approximation of erfinv(x), to be validated */
+static void
+arb_hypgeom_erfinv_guess(arb_t res, const arb_t x, const arb_t one_sub_x, slong extraprec)
+{
+ if (arf_cmpabs_2exp_si(arb_midref(x), -30) < 0)
+ {
+ arb_erfinv_approx_tiny(res, x, 128);
+ }
+ else if (arf_cmpabs_2exp_si(arb_midref(one_sub_x), -52) >= 0)
+ {
+ double y;
+
+ y = erfinv_approx(arf_get_d(arb_midref(x), ARF_RND_NEAR),
+ arf_get_d(arb_midref(one_sub_x), ARF_RND_NEAR));
+
+ arf_set_d(arb_midref(res), y);
+ mag_set_d(arb_radref(res), ldexp(y, -50));
+ }
+ else if (arf_cmpabs_2exp_si(arb_midref(one_sub_x), -1000) >= 0)
+ {
+ double t, y;
+
+ t = arf_get_d(arb_midref(one_sub_x), ARF_RND_NEAR);
+ y = erfinv_approx_big(t);
+
+ arf_set_d(arb_midref(res), y);
+ mag_set_d(arb_radref(res), ldexp(y, -26 + 0.1 * log(t)));
+ }
+ else
+ {
+ arb_erfinv_approx_huge(res, one_sub_x, 30 + extraprec);
+ }
+}
+
+void
+arb_hypgeom_erfinv_precise(arb_t res, const arb_t x, const arb_t one_sub_x, int near_one, slong prec)
+{
+ slong wp;
+ arb_t f, fprime, root, mid, t;
+ slong extraprec, goal;
+ int validated;
+
+ if (arb_is_zero(x))
+ {
+ arb_zero(res);
+ return;
+ }
+
+ arb_init(f);
+ arb_init(fprime);
+ arb_init(root);
+ arb_init(mid);
+ arb_init(t);
+
+ goal = prec * 1.001 + 5;
+
+ extraprec = fmpz_bits(ARF_EXPREF(arb_midref(one_sub_x)));
+ extraprec += 15;
+
+ /* Start with a guess */
+ arb_hypgeom_erfinv_guess(root, x, one_sub_x, extraprec);
+ validated = 0;
+
+ while (!validated || arb_rel_accuracy_bits(root) <= goal)
+ {
+ /* We should get double the accuracy. */
+ wp = arb_rel_accuracy_bits(root) * 2 + extraprec;
+ /* But don't set the precision unrealistically high. */
+ wp = FLINT_MIN(wp, 4 * (goal + extraprec));
+
+ /* In case of quadratic convergence, avoid doing the
+ penultimate iteration at higher precision than needed. */
+ if (validated && wp < goal && wp > 0.7 * goal + 2 * extraprec)
+ wp = goal / 2 + 2 * extraprec;
+
+ arb_set(mid, root);
+ mag_zero(arb_radref(mid));
+
+ /* f(y) = erf(y) - x OR (1 - x) - erfc(y) */
+ /* 1/f'(y) = exp(y^2) * sqrt(pi)/2 */
+ if (near_one)
+ {
+ arb_hypgeom_erfc(f, mid, wp);
+ arb_sub(f, one_sub_x, f, wp);
+ }
+ else
+ {
+ arb_hypgeom_erf(f, mid, wp);
+ arb_sub(f, f, x, wp);
+ }
+ arb_sqr(fprime, root, wp);
+ arb_exp(fprime, fprime, wp);
+ arb_const_sqrt_pi(t, wp);
+ arb_mul(fprime, fprime, t, wp);
+ arb_mul_2exp_si(fprime, fprime, -1);
+ arb_mul(t, f, fprime, wp);
+ arb_sub(t, mid, t, wp);
+
+ if (arb_contains_interior(root, t))
+ {
+ /* Interval Newton proves inclusion. */
+ /* printf("newton %d -> 1 %ld: ", validated, wp); arb_printd(t, 50); printf("\n"); */
+ validated = 1;
+ arb_swap(root, t);
+ }
+ else
+ {
+ /* Try to improve the guess with a floating-point Newton step. */
+ /* printf("newton %d -> 0 %ld: ", validated, wp); arb_printd(t, 50); printf("\n"); */
+ arb_sqr(fprime, mid, wp);
+ arb_exp(fprime, fprime, wp);
+ arb_const_sqrt_pi(t, wp);
+ arb_mul(fprime, fprime, t, wp);
+ arb_mul_2exp_si(fprime, fprime, -1);
+ arb_mul(t, f, fprime, wp);
+ /* The Newton correction is a good guess for the error. */
+ arb_get_mag(arb_radref(root), t);
+ mag_mul_2exp_si(arb_radref(root), arb_radref(root), 1);
+ arb_sub(t, mid, t, wp);
+ arf_swap(arb_midref(root), arb_midref(t));
+ /* Use more working precision to get out of any numerical difficulties */
+ extraprec = extraprec * 1.05 + 10;
+ validated = 0;
+ }
+
+ /* Something went wrong. Could fall back to bisection if we
+ want to guarantee convergence? */
+ if (extraprec > 10 * prec + 10000)
+ {
+ arb_indeterminate(root);
+ break;
+ }
+ }
+
+ arb_set_round(res, root, prec);
+
+ arb_clear(f);
+ arb_clear(fprime);
+ arb_clear(root);
+ arb_clear(mid);
+ arb_clear(t);
+}
+
+
+
+
+
+
+
+
+
+static slong
+arb_adjust_precision(slong prec, slong acc)
+{
+ acc = FLINT_MIN(acc, prec);
+ acc = FLINT_MAX(acc, 0);
+ prec = FLINT_MIN(prec, acc + MAG_BITS);
+ prec = FLINT_MAX(prec, 2);
+ return prec;
+}
+
+
+
+void
+arb_hypgeom_erfinv(arb_t res, const arb_t x, slong prec)
+{
+ arb_t x1;
+ int near_one;
+
+ if (arb_is_zero(x))
+ {
+ arb_zero(res);
+ return;
+ }
+
+ if (arf_sgn(arb_midref(x)) < 0)
+ {
+ arb_neg(res, x);
+ arb_hypgeom_erfinv(res, res, prec);
+ arb_neg(res, res);
+ return;
+ }
+
+ if (arb_is_one(x))
+ {
+ arb_pos_inf(res);
+ return;
+ }
+
+ arb_init(x1);
+ near_one = ARF_EXP(arb_midref(x)) == 0;
+
+ if (near_one)
+ {
+ arb_sub_ui(x1, x, 1, ARF_PREC_EXACT);
+ arb_neg(x1, x1);
+ }
+ else
+ {
+ arb_sub_ui(x1, x, 1, prec + 30);
+ arb_neg(x1, x1);
+ }
+
+ if (arb_is_positive(x1))
+ {
+ mag_t err;
+ slong acc;
+ arb_t xm;
+
+ mag_init(err);
+ arb_init(xm);
+
+ /* Propagated error bound based on derivative. */
+ /* erfinv'(x) <= (1/2) sqrt(pi) / (1 - |x|) */
+ arb_get_mag_lower(err, x1);
+ mag_inv(err, err);
+ mag_mul(err, err, arb_radref(x));
+ mag_mul_ui(err, err, 227);
+ mag_mul_2exp_si(err, err, -8);
+
+ acc = arb_rel_accuracy_bits(x);
+ prec = arb_adjust_precision(prec, acc);
+
+ arb_get_mid_arb(xm, x);
+
+ if (near_one)
+ {
+ arb_sub_ui(x1, xm, 1, ARF_PREC_EXACT);
+ arb_neg(x1, x1);
+ }
+ else
+ {
+ arb_sub_ui(x1, xm, 1, prec + 30);
+ arb_neg(x1, x1);
+ }
+
+ arb_hypgeom_erfinv_precise(res, xm, x1, near_one, prec);
+ arb_add_error_mag(res, err);
+
+ mag_clear(err);
+ arb_clear(xm);
+ }
+ else
+ {
+ arb_indeterminate(res);
+ }
+
+ arb_clear(x1);
+}
+
+void
+arb_hypgeom_erfcinv(arb_t res, const arb_t x1, slong prec)
+{
+ arb_t x;
+
+ if (arb_is_one(x1))
+ {
+ arb_zero(res);
+ return;
+ }
+
+ arb_init(x);
+
+ if (arf_cmp_d(arb_midref(x1), 0.01) <= 0 && arb_is_positive(x1))
+ {
+ mag_t err;
+ slong acc;
+ arb_t x1m, xm;
+
+ mag_init(err);
+ arb_init(x1m);
+ arb_init(xm);
+
+ /* Propagated error bound based on derivative. */
+ /* erfinv'(x) <= (1/2) sqrt(pi) / (1 - |x|) */
+ arb_get_mag_lower(err, x1);
+ mag_inv(err, err);
+ mag_mul(err, err, arb_radref(x1));
+ mag_mul_ui(err, err, 227);
+ mag_mul_2exp_si(err, err, -8);
+
+ acc = arb_rel_accuracy_bits(x1);
+ prec = arb_adjust_precision(prec, acc);
+
+ arb_get_mid_arb(x1m, x1);
+ arb_sub_ui(xm, x1m, 1, 2 * prec + 100);
+ arb_neg(xm, xm);
+
+ arb_hypgeom_erfinv_precise(res, xm, x1m, 1, prec);
+ arb_add_error_mag(res, err);
+
+ mag_clear(err);
+ arb_clear(xm);
+ arb_clear(x1m);
+ }
+ else
+ {
+ arb_sub_ui(x, x1, 1, 2 * prec + 100);
+ arb_neg(x, x);
+ arb_hypgeom_erfinv(res, x, prec);
+ }
+
+ arb_clear(x);
+}
diff --git a/arb_hypgeom/test/t-erfinv.c b/arb_hypgeom/test/t-erfinv.c
new file mode 100644
index 00000000..e7e85a45
--- /dev/null
+++ b/arb_hypgeom/test/t-erfinv.c
@@ -0,0 +1,177 @@
+/*
+ Copyright (C) 2021 Fredrik Johansson
+
+ This file is part of Arb.
+
+ Arb is free software: you can redistribute it and/or modify it under
+ the terms of the GNU Lesser General Public License (LGPL) as published
+ by the Free Software Foundation; either version 2.1 of the License, or
+ (at your option) any later version. See .
+*/
+
+#include "arb_hypgeom.h"
+
+int main()
+{
+ slong iter;
+ flint_rand_t state;
+
+ flint_printf("erfinv....");
+ fflush(stdout);
+
+ flint_randinit(state);
+
+ for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
+ {
+ arb_t x, y, z, w;
+ slong prec, prec2;
+
+ arb_init(x);
+ arb_init(y);
+ arb_init(z);
+ arb_init(w);
+
+ prec = 2 + n_randint(state, 500);
+ prec2 = 2 + n_randint(state, 500);
+
+ arb_init(x);
+ arb_init(y);
+
+ arb_randtest(x, state, prec, 10 + n_randint(state, 300));
+ arb_randtest(y, state, prec, 100);
+
+ arb_hypgeom_erfcinv(y, x, prec);
+ arb_hypgeom_erfc(z, y, prec + 10);
+
+ if (!arb_overlaps(x, z))
+ {
+ flint_printf("FAIL: overlap\n\n");
+ flint_printf("prec = %wd\n", prec);
+ flint_printf("x = "); arb_printn(x, 100, 0); flint_printf("\n\n");
+ flint_printf("y = "); arb_printn(y, 100, 0); flint_printf("\n\n");
+ flint_printf("z = "); arb_printn(z, 100, 0); flint_printf("\n\n");
+ flint_abort();
+ }
+
+ if (arb_is_exact(x) && arf_sgn(arb_midref(x)) > 0 &&
+ arf_cmp_2exp_si(arb_midref(x), 1) < 0 &&
+ arb_rel_accuracy_bits(y) < prec - 2)
+ {
+ flint_printf("FAIL: accuracy\n\n");
+ flint_printf("iter = %wd, prec = %wd\n", iter, prec);
+ flint_printf("x = "); arb_printn(x, 100, 0); flint_printf("\n\n");
+ flint_printf("y = "); arb_printn(y, 100, 0); flint_printf("\n\n");
+ flint_printf("z = "); arb_printn(z, 100, 0); flint_printf("\n\n");
+ flint_abort();
+ }
+
+ arb_set(z, x);
+ mag_zero(arb_radref(z));
+
+ if (n_randint(state, 2))
+ {
+ arf_set_mag(arb_midref(w), arb_radref(x));
+
+ if (n_randint(state, 2))
+ arb_neg(w, w);
+
+ arb_add(z, z, w, prec);
+ }
+
+ arb_hypgeom_erfcinv(z, z, prec2);
+
+ if (!arb_overlaps(y, z))
+ {
+ flint_printf("FAIL: overlap 2\n\n");
+ flint_printf("prec = %wd\n", prec);
+ flint_printf("x = "); arb_printn(x, 100, 0); flint_printf("\n\n");
+ flint_printf("y = "); arb_printn(y, 100, 0); flint_printf("\n\n");
+ flint_printf("z = "); arb_printn(z, 100, 0); flint_printf("\n\n");
+ flint_abort();
+ }
+
+ arb_clear(x);
+ arb_clear(y);
+ arb_clear(z);
+ arb_clear(w);
+ }
+
+ for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
+ {
+ arb_t x, y, z, w;
+ slong prec, prec2;
+
+ arb_init(x);
+ arb_init(y);
+ arb_init(z);
+ arb_init(w);
+
+ prec = 2 + n_randint(state, 500);
+ prec2 = 2 + n_randint(state, 500);
+
+ arb_init(x);
+ arb_init(y);
+
+ arb_randtest(x, state, prec, 100);
+ arb_randtest(y, state, prec, 100);
+
+ arb_hypgeom_erfinv(y, x, prec);
+ arb_hypgeom_erf(z, y, prec + 10);
+
+ if (!arb_overlaps(x, z))
+ {
+ flint_printf("FAIL: overlap\n\n");
+ flint_printf("prec = %wd\n", prec);
+ flint_printf("x = "); arb_printn(x, 100, 0); flint_printf("\n\n");
+ flint_printf("y = "); arb_printn(y, 100, 0); flint_printf("\n\n");
+ flint_printf("z = "); arb_printn(z, 100, 0); flint_printf("\n\n");
+ flint_abort();
+ }
+
+ if (arb_is_exact(x) && arf_cmpabs_2exp_si(arb_midref(x), 0) < 0 &&
+ arb_rel_accuracy_bits(y) < prec - 2)
+ {
+ flint_printf("FAIL: accuracy\n\n");
+ flint_printf("prec = %wd\n", prec);
+ flint_printf("x = "); arb_printn(x, 100, 0); flint_printf("\n\n");
+ flint_printf("y = "); arb_printn(y, 100, 0); flint_printf("\n\n");
+ flint_printf("z = "); arb_printn(z, 100, 0); flint_printf("\n\n");
+ flint_abort();
+ }
+
+ arb_set(z, x);
+ mag_zero(arb_radref(z));
+
+ if (n_randint(state, 2))
+ {
+ arf_set_mag(arb_midref(w), arb_radref(x));
+
+ if (n_randint(state, 2))
+ arb_neg(w, w);
+
+ arb_add(z, z, w, prec);
+ }
+
+ arb_hypgeom_erfinv(z, z, prec2);
+
+ if (!arb_overlaps(y, z))
+ {
+ flint_printf("FAIL: overlap 2\n\n");
+ flint_printf("prec = %wd\n", prec);
+ flint_printf("x = "); arb_printn(x, 100, 0); flint_printf("\n\n");
+ flint_printf("y = "); arb_printn(y, 100, 0); flint_printf("\n\n");
+ flint_printf("z = "); arb_printn(z, 100, 0); flint_printf("\n\n");
+ flint_abort();
+ }
+
+ arb_clear(x);
+ arb_clear(y);
+ arb_clear(z);
+ arb_clear(w);
+ }
+
+ flint_randclear(state);
+ flint_cleanup();
+ flint_printf("PASS\n");
+ return EXIT_SUCCESS;
+}
diff --git a/doc/source/arb_hypgeom.rst b/doc/source/arb_hypgeom.rst
index 335ccc81..b8a8bb8d 100644
--- a/doc/source/arb_hypgeom.rst
+++ b/doc/source/arb_hypgeom.rst
@@ -208,6 +208,12 @@ Error functions and Fresnel integrals
Computes the imaginary error function of the power series *z*,
truncated to length *len*.
+.. function:: void arb_hypgeom_erfinv(arb_t res, const arb_t z, slong prec)
+ void arb_hypgeom_erfcinv(arb_t res, const arb_t z, slong prec)
+
+ Computes the inverse error function `\operatorname{erf}^{-1}(z)`
+ or inverse complementary error function `\operatorname{erfc}^{-1}(z)`.
+
.. function:: void arb_hypgeom_fresnel(arb_t res1, arb_t res2, const arb_t z, int normalized, slong prec)
Sets *res1* to the Fresnel sine integral `S(z)` and *res2* to