hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-05 17:31:38 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

optimize arb_set_interval_mag

Browse source
This commit is contained in:
Fredrik Johansson 2017-12-30 21:49:06 +01:00
parent 3ef95feb3f
commit b5cb696bf2
2 changed files with 186 additions and 41 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

107
arb/set_interval_mag.c
View file

@ -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);
}

120
arb/set_interval_neg_pos_mag.c
View file

@ -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);
}
}