From b5cb696bf2f4b913ce5293f63d22ad7f8f480fe0 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Sat, 30 Dec 2017 21:49:06 +0100 Subject: [PATCH] optimize arb_set_interval_mag --- arb/set_interval_mag.c | 107 +++++++++++++++++++++-------- arb/set_interval_neg_pos_mag.c | 120 +++++++++++++++++++++++++++++---- 2 files changed, 186 insertions(+), 41 deletions(-) diff --git a/arb/set_interval_mag.c b/arb/set_interval_mag.c index a4f69677..f6aa86ab 100644 --- a/arb/set_interval_mag.c +++ b/arb/set_interval_mag.c @@ -14,37 +14,88 @@ void arb_set_interval_mag(arb_t res, const mag_t a, const mag_t b, slong prec) { - arf_t aa, bb; - int inexact; - - if (mag_cmp(a, b) > 0) + if (MAG_IS_LAGOM(a) && MAG_IS_LAGOM(b)) { - flint_printf("exception: arb_set_interval_mag: endpoints not ordered\n"); - flint_abort(); - } + slong aexp, bexp; + mp_limb_t aman, bman, mman, rman, tmp; - if (mag_is_inf(a)) + aman = MAG_MAN(a); + bman = MAG_MAN(b); + + aexp = MAG_EXP(a); + bexp = MAG_EXP(b); + + if (aman == 0 && bman == 0) + { + arb_zero(res); + return; + } + + if (bman == 0 || (aman != 0 && + (aexp > bexp || (aexp == bexp && aman > bman)))) + { + flint_printf("exception: arb_set_interval_mag: endpoints not ordered\n"); + flint_abort(); + } + + /* now a = 0 or bexp >= aexp */ + if (aman == 0 || bexp - aexp > MAG_BITS) + { + mman = bman; /* midpoint a+b */ + rman = bman + (aman != 0); /* radius b-a */ + } + else + { + tmp = (aman >> (bexp - aexp)); + mman = bman + tmp; /* midpoint a+b */ + rman = bman - tmp; /* radius b-a */ + rman += ((tmp << (bexp - aexp)) != aman); /* rounding error */ + } + + arf_set_ui(arb_midref(res), mman); + /* m can't be zero */ + ARF_EXP(arb_midref(res)) += bexp - MAG_BITS - 1; + + mag_set_ui(arb_radref(res), rman); + if (rman != 0) /* r can be zero */ + MAG_EXP(arb_radref(res)) += bexp - MAG_BITS - 1; + + arb_set_round(res, res, prec); + } + else { - arb_pos_inf(res); - return; + int inexact; + arf_t aa, bb; + + if (mag_cmp(a, b) > 0) + { + flint_printf("exception: arb_set_interval_mag: endpoints not ordered\n"); + flint_abort(); + } + + if (mag_is_inf(a)) + { + arb_pos_inf(res); + return; + } + + if (mag_is_inf(b)) + { + arb_zero_pm_inf(res); + return; + } + + arf_init_set_mag_shallow(aa, a); + arf_init_set_mag_shallow(bb, b); + + inexact = arf_add(arb_midref(res), aa, bb, prec, ARB_RND); + + mag_sub(arb_radref(res), b, a); + + if (inexact) + arf_mag_add_ulp(arb_radref(res), arb_radref(res), arb_midref(res), prec); + + arb_mul_2exp_si(res, res, -1); } - - if (mag_is_inf(b)) - { - arb_zero_pm_inf(res); - return; - } - - arf_init_set_mag_shallow(aa, a); - arf_init_set_mag_shallow(bb, b); - - inexact = arf_add(arb_midref(res), aa, bb, prec, ARB_RND); - - mag_sub(arb_radref(res), b, a); - - if (inexact) - arf_mag_add_ulp(arb_radref(res), arb_radref(res), arb_midref(res), prec); - - arb_mul_2exp_si(res, res, -1); } diff --git a/arb/set_interval_neg_pos_mag.c b/arb/set_interval_neg_pos_mag.c index 9148d40b..3a031736 100644 --- a/arb/set_interval_neg_pos_mag.c +++ b/arb/set_interval_neg_pos_mag.c @@ -14,25 +14,119 @@ void arb_set_interval_neg_pos_mag(arb_t res, const mag_t a, const mag_t b, slong prec) { - arf_t aa, bb; - int inexact; - - if (mag_is_inf(a) || mag_is_inf(b)) + if (MAG_IS_LAGOM(a) && MAG_IS_LAGOM(b)) { - arb_zero_pm_inf(res); - return; + slong aexp, bexp, mexp, shift; + mp_limb_t aman, bman, mman, rman, tmp; + int negative; + + aman = MAG_MAN(a); + bman = MAG_MAN(b); + + aexp = MAG_EXP(a); + bexp = MAG_EXP(b); + + if (aman == 0) + { + if (bman == 0) + { + arb_zero(res); + return; + } + + negative = 0; + mexp = bexp; + mman = bman; + rman = bman; + } + else if (bman == 0) + { + negative = 1; + mexp = aexp; + mman = aman; + rman = aman; + } + else if (aexp == bexp) + { + mexp = aexp; + negative = aman >= bman; + if (negative) + mman = aman - bman; + else + mman = bman - aman; + rman = aman + bman; + } + else if (aexp > bexp) + { + negative = 1; + mexp = aexp; + shift = aexp - bexp; + if (shift > MAG_BITS) + { + mman = aman; + rman = aman + 2; + } + else + { + tmp = bman >> shift; + mman = aman - tmp; + rman = aman + tmp; + rman += 2 * ((tmp << shift) != bman); + } + } + else + { + negative = 0; + mexp = bexp; + shift = bexp - aexp; + if (shift > MAG_BITS) + { + mman = bman; + rman = bman + 2; + } + else + { + tmp = aman >> shift; + mman = bman - tmp; + rman = bman + tmp; + rman += 2 * ((tmp << shift) != aman); + } + } + + arf_set_ui(arb_midref(res), mman); + if (negative) + arf_neg(arb_midref(res), arb_midref(res)); + if (mman != 0) + ARF_EXP(arb_midref(res)) += mexp - MAG_BITS - 1; + + mag_set_ui(arb_radref(res), rman); + /* r can't be zero */ + MAG_EXP(arb_radref(res)) += mexp - MAG_BITS - 1; + + arb_set_round(res, res, prec); } + else + { + arf_t aa, bb; + int inexact; - arf_init_set_mag_shallow(aa, a); - arf_init_set_mag_shallow(bb, b); + if (mag_is_inf(a) || mag_is_inf(b)) + { + arb_zero_pm_inf(res); + return; + } - inexact = arf_sub(arb_midref(res), bb, aa, prec, ARB_RND); + arf_init_set_mag_shallow(aa, a); + arf_init_set_mag_shallow(bb, b); - mag_add(arb_radref(res), b, a); + inexact = arf_sub(arb_midref(res), bb, aa, prec, ARB_RND); - if (inexact) - arf_mag_add_ulp(arb_radref(res), arb_radref(res), arb_midref(res), prec); + mag_add(arb_radref(res), b, a); - arb_mul_2exp_si(res, res, -1); + if (inexact) + arf_mag_add_ulp(arb_radref(res), arb_radref(res), arb_midref(res), prec); + + arb_mul_2exp_si(res, res, -1); + } }