From a1efd15f0185c243063263d01969b6f5c89da855 Mon Sep 17 00:00:00 2001 From: fredrik Date: Mon, 9 Jul 2018 13:25:19 +0200 Subject: [PATCH] slightly better arb_exp_invexp --- arb/exp.c | 111 ++++++++++++++++++++++++++++++---------- arb/test/t-exp_invexp.c | 34 +++++++++--- 2 files changed, 112 insertions(+), 33 deletions(-) diff --git a/arb/exp.c b/arb/exp.c index 64227e8e..239afd77 100644 --- a/arb/exp.c +++ b/arb/exp.c @@ -562,46 +562,105 @@ arb_expm1(arb_t z, const arb_t x, slong prec) } } -void -arb_exp_invexp(arb_t z, arb_t w, const arb_t x, slong prec) +void arb_exp_invexp(arb_t res, arb_t res2, const arb_t x, slong prec) { - if (arb_is_exact(x)) + slong maglim = MAGLIM(prec); + + if (arf_is_special(arb_midref(x)) || mag_is_special(arb_radref(x))) { - arb_exp_arf(z, arb_midref(x), prec, 0, MAGLIM(prec)); - arb_inv(w, z, prec); - } - else - { - /* exp(a+b) - exp(a) = exp(a) * (exp(b)-1) */ - if (mag_cmp_2exp_si(arb_radref(x), 20) < 0 || !arb_is_finite(x)) + /* [c +/- 0] */ + if (arf_is_finite(arb_midref(x)) && mag_is_zero(arb_radref(x))) { + arb_exp_arf(res, arb_midref(x), prec, 0, maglim); + arb_inv(res2, res, prec); + } /* [nan +/- ?] */ + else if (arf_is_nan(arb_midref(x))) + { + arb_indeterminate(res); + arb_indeterminate(res2); + } /* [c +/- inf] */ + else if (mag_is_inf(arb_radref(x))) + { + arb_zero_pm_inf(res); + arb_zero_pm_inf(res2); + } /* [+inf +/- c] */ + else if (arf_is_pos_inf(arb_midref(x))) + { + arb_pos_inf(res); + arb_zero(res2); + } /* [-inf +/- c] */ + else if (arf_is_neg_inf(arb_midref(x))) + { + arb_zero(res); + arb_pos_inf(res2); + } + else + { + arb_t t; + arb_init(t); + arb_neg(t, x); + arb_exp_wide(res, x, prec, maglim); + arb_exp_wide(res2, t, prec, maglim); + arb_clear(t); + } + } + else /* both finite, non-special */ + { + slong acc, mexp, rexp; + + mexp = ARF_EXP(arb_midref(x)); + rexp = MAG_EXP(arb_radref(x)); + + if (COEFF_IS_MPZ(rexp)) + rexp = (fmpz_sgn(ARF_EXPREF(arb_radref(x))) < 0) ? COEFF_MIN : COEFF_MAX; + if (COEFF_IS_MPZ(mexp)) + mexp = (fmpz_sgn(MAG_EXPREF(arb_midref(x))) < 0) ? COEFF_MIN : COEFF_MAX; + + if (mexp < -prec && rexp < -prec) + { + arb_get_mag(arb_radref(res), x); + mag_expm1(arb_radref(res), arb_radref(res)); + arf_one(arb_midref(res)); + arb_set(res2, res); + return; + } + + acc = -rexp; + acc = FLINT_MAX(acc, 0); + acc = FLINT_MIN(acc, prec); + prec = FLINT_MIN(prec, acc + MAG_BITS); + prec = FLINT_MAX(prec, 2); + + if (acc < 20 && (rexp >= 0 || mexp <= 10)) + { + /* may evaluate at endpoints */ + arb_t t; + arb_init(t); + arb_neg(t, x); + arb_exp_wide(res, x, prec, maglim); + arb_exp_wide(res2, t, prec, maglim); + arb_clear(t); + } + else + { + /* exp(a+b) - exp(a) = exp(a) * (exp(b)-1) */ mag_t t, u; mag_init_set(t, arb_radref(x)); mag_init(u); - arb_exp_arf(z, arb_midref(x), prec, 0, MAGLIM(prec)); - arb_inv(w, z, prec); + arb_exp_arf(res, arb_midref(x), prec, 0, maglim); + arb_inv(res2, res, prec); mag_expm1(t, t); - arb_get_mag(u, z); - mag_addmul(arb_radref(z), t, u); - arb_get_mag(u, w); - mag_addmul(arb_radref(w), t, u); + arb_get_mag(u, res); + mag_addmul(arb_radref(res), t, u); + arb_get_mag(u, res2); + mag_addmul(arb_radref(res2), t, u); mag_clear(t); mag_clear(u); } - else - { - arb_t v; - arb_init(v); - arb_neg(v, x); - arb_exp(z, x, prec); - arb_exp(w, v, prec); - arb_clear(v); - } } } - diff --git a/arb/test/t-exp_invexp.c b/arb/test/t-exp_invexp.c index a5f90904..6429f770 100644 --- a/arb/test/t-exp_invexp.c +++ b/arb/test/t-exp_invexp.c @@ -24,7 +24,7 @@ int main() for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++) { arb_t a, b, c, d; - slong prec; + slong prec, acc1, acc2; if (iter % 10 == 0) prec = 10000; @@ -42,30 +42,51 @@ int main() arb_randtest_special(b, state, 1 + n_randint(state, prec), 1 + n_randint(state, 100)); arb_randtest_special(c, state, 1 + n_randint(state, prec), 1 + n_randint(state, 100)); - arb_exp_invexp(b, c, a, prec); + if (n_randint(state, 2)) + { + arb_exp_invexp(b, c, a, prec); + } + else if (n_randint(state, 2)) + { + arb_set(b, a); + arb_exp_invexp(b, c, b, prec); + } + else + { + arb_set(c, a); + arb_exp_invexp(b, c, c, prec); + } arb_exp(d, a, prec); - if (!arb_overlaps(b, d)) + acc1 = arb_rel_accuracy_bits(b); + acc2 = arb_rel_accuracy_bits(d); + + if (!arb_overlaps(b, d) || ((acc1 > 0 || acc2 > 0) && acc1 < FLINT_MIN(acc2, prec) - 3)) { flint_printf("FAIL: overlap 1\n\n"); + flint_printf("prec = %wd, acc1 = %wd, acc2 = %wd\n\n", prec, acc1, acc2); flint_printf("a = "); arb_print(a); flint_printf("\n\n"); flint_printf("b = "); arb_print(b); flint_printf("\n\n"); flint_printf("c = "); arb_print(c); flint_printf("\n\n"); - flint_printf("d = "); arb_print(c); flint_printf("\n\n"); + flint_printf("d = "); arb_print(d); flint_printf("\n\n"); flint_abort(); } arb_neg(d, a); arb_exp(d, d, prec); - if (!arb_overlaps(c, d)) + acc1 = arb_rel_accuracy_bits(c); + acc2 = arb_rel_accuracy_bits(d); + + if (!arb_overlaps(c, d) || ((acc1 > 0 || acc2 > 0) && acc1 < FLINT_MIN(acc2, prec) - 3)) { flint_printf("FAIL: overlap 2\n\n"); + flint_printf("prec = %wd, acc1 = %wd, acc2 = %wd\n\n", prec, acc1, acc2); flint_printf("a = "); arb_print(a); flint_printf("\n\n"); flint_printf("b = "); arb_print(b); flint_printf("\n\n"); flint_printf("c = "); arb_print(c); flint_printf("\n\n"); - flint_printf("d = "); arb_print(c); flint_printf("\n\n"); + flint_printf("d = "); arb_print(d); flint_printf("\n\n"); flint_abort(); } @@ -80,4 +101,3 @@ int main() flint_printf("PASS\n"); return EXIT_SUCCESS; } -