Skip to content

Commit

Permalink
Merge pull request #409 from postmath/issue408
Browse files Browse the repository at this point in the history 
Positive * pos_inf should be pos_inf, not 0 +- inf. Same for other combinations of signs.
  • Loading branch information
@fredrik-johansson
fredrik-johansson authored Mar 25, 2022
2 parents 076249e + a9f8822 commit ab5e81c
Show file tree
Hide file tree
Showing 8 changed files with 422 additions and 27 deletions.
21 changes: 21 additions & 0 deletions arb/addmul.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,13 @@ arb_addmul_arf(arb_t z, const arb_t x, const arf_t y, slong prec)
if (inexact)
arf_mag_fast_add_ulp(arb_radref(z), arb_radref(z), arb_midref(z), prec);
}
else if (arf_is_inf(y) && arb_is_nonzero(x))
{
if (arf_sgn(arb_midref(x)) > 0)
arb_add_arf(z, z, y, prec);
else
arb_sub_arf(z, z, y, prec);
}
else
{
mag_init_set_arf(ym, y);
Expand Down Expand Up @@ -78,6 +85,20 @@ arb_addmul(arb_t z, const arb_t x, const arb_t y, slong prec)

*arb_radref(z) = *zr;
}
else if (arf_is_inf(arb_midref(x)) && mag_is_finite(arb_radref(x)) && arb_is_nonzero(y))
{
if (arf_sgn(arb_midref(y)) > 0)
arb_add_arf(z, z, arb_midref(x), prec);
else
arb_sub_arf(z, z, arb_midref(x), prec);
}
else if (arf_is_inf(arb_midref(y)) && mag_is_finite(arb_radref(y)) && arb_is_nonzero(x))
{
if (arf_sgn(arb_midref(x)) > 0)
arb_add_arf(z, z, arb_midref(y), prec);
else
arb_sub_arf(z, z, arb_midref(y), prec);
}
else
{
mag_init_set_arf(xm, arb_midref(x));
Expand Down
21 changes: 21 additions & 0 deletions arb/fma.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,13 @@ arb_fma_arf(arb_t res, const arb_t x, const arf_t y, const arb_t z, slong prec)
if (inexact)
arf_mag_fast_add_ulp(arb_radref(res), arb_radref(res), arb_midref(res), prec);
}
else if (arf_is_inf(y) && arb_is_nonzero(x))
{
if (arf_sgn(arb_midref(x)) > 0)
arb_add_arf(res, z, y, prec);
else
arb_sub_arf(res, z, y, prec);
}
else
{
mag_t tm;
Expand Down Expand Up @@ -91,6 +98,20 @@ arb_fma(arb_t res, const arb_t x, const arb_t y, const arb_t z, slong prec)

*arb_radref(res) = *zr;
}
else if (arf_is_inf(arb_midref(x)) && mag_is_finite(arb_radref(x)) && arb_is_nonzero(y))
{
if (arf_sgn(arb_midref(y)) > 0)
arb_add_arf(res, z, arb_midref(x), prec);
else
arb_sub_arf(res, z, arb_midref(x), prec);
}
else if (arf_is_inf(arb_midref(y)) && mag_is_finite(arb_radref(y)) && arb_is_nonzero(x))
{
if (arf_sgn(arb_midref(x)) > 0)
arb_add_arf(res, z, arb_midref(y), prec);
else
arb_sub_arf(res, z, arb_midref(y), prec);
}
else
{
mag_init_set_arf(xm, arb_midref(x));
Expand Down
19 changes: 19 additions & 0 deletions arb/mul.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,15 @@ arb_mul_arf(arb_t z, const arb_t x, const arf_t y, slong prec)

*arb_radref(z) = *zr;
}
else if (arf_is_inf(y) && arb_is_nonzero(x))
{
mag_zero(arb_radref(z));

if (arf_sgn(arb_midref(x)) * arf_sgn(y) > 0)
arf_pos_inf(arb_midref(z));
else
arf_neg_inf(arb_midref(z));
}
else
{
mag_init_set_arf(ym, y);
Expand Down Expand Up @@ -91,6 +100,16 @@ arb_mul(arb_t z, const arb_t x, const arb_t y, slong prec)

*arb_radref(z) = *zr;
}
else if (arf_is_inf(arb_midref(x)) && mag_is_finite(arb_radref(x)) && arb_is_nonzero(y) ||
arf_is_inf(arb_midref(y)) && mag_is_finite(arb_radref(y)) && arb_is_nonzero(x))
{
mag_zero(arb_radref(z));

if (arf_sgn(arb_midref(x)) * arf_sgn(arb_midref(y)) > 0)
arf_pos_inf(arb_midref(z));
else
arf_neg_inf(arb_midref(z));
}
else
{
mag_init_set_arf(xm, arb_midref(x));
Expand Down
21 changes: 21 additions & 0 deletions arb/submul.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,13 @@ arb_submul_arf(arb_t z, const arb_t x, const arf_t y, slong prec)
if (inexact)
arf_mag_fast_add_ulp(arb_radref(z), arb_radref(z), arb_midref(z), prec);
}
else if (arf_is_inf(y) && arb_is_nonzero(x))
{
if (arf_sgn(arb_midref(x)) > 0)
arb_sub_arf(z, z, y, prec);
else
arb_add_arf(z, z, y, prec);
}
else
{
mag_init_set_arf(ym, y);
Expand Down Expand Up @@ -78,6 +85,20 @@ arb_submul(arb_t z, const arb_t x, const arb_t y, slong prec)

*arb_radref(z) = *zr;
}
else if (arf_is_inf(arb_midref(x)) && mag_is_finite(arb_radref(x)) && arb_is_nonzero(y))
{
if (arf_sgn(arb_midref(y)) > 0)
arb_sub_arf(z, z, arb_midref(x), prec);
else
arb_add_arf(z, z, arb_midref(x), prec);
}
else if (arf_is_inf(arb_midref(y)) && mag_is_finite(arb_radref(y)) && arb_is_nonzero(x))
{
if (arf_sgn(arb_midref(x)) > 0)
arb_sub_arf(z, z, arb_midref(y), prec);
else
arb_add_arf(z, z, arb_midref(y), prec);
}
else
{
mag_init_set_arf(xm, arb_midref(x));
Expand Down
28 changes: 18 additions & 10 deletions arb/test/t-addmul.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,19 @@ mag_close(const mag_t am, const mag_t bm)
return res1 && res2;
}

int
arb_equal_mid_close_mag(const arb_t a, const arb_t b)
{
return arf_equal(arb_midref(a), arb_midref(b)) &&
(mag_close(arb_radref(a), arb_radref(b)) ||
/* If a's and b's centres are infinite but their radii are finite, the radii don't need
to be close: they represent signed infinity regardless. If their centres are NaN,
then we should ignore their radii. */
(arf_is_inf(arb_midref(a)) && arf_is_inf(arb_midref(b)) &&
mag_is_finite(arb_radref(a)) && mag_is_finite(arb_radref(b))) ||
(arf_is_nan(arb_midref(a)) && arf_is_nan(arb_midref(b))));
}

void
arb_addmul_naive(arb_t z, const arb_t x, const arb_t y, slong prec)
{
Expand Down Expand Up @@ -214,8 +227,7 @@ int main()
arb_addmul(z, x, y, prec);
arb_addmul_naive(v, x, y, prec);

if (!arf_equal(arb_midref(z), arb_midref(v))
|| !mag_close(arb_radref(z), arb_radref(v)))
if (!arb_equal_mid_close_mag(z, v))
{
flint_printf("FAIL!\n");
flint_printf("x = "); arb_print(x); flint_printf("\n\n");
Expand All @@ -232,8 +244,7 @@ int main()
arb_addmul(z, x, y, prec);
arb_addmul(v, x, x, prec);

if (!arf_equal(arb_midref(z), arb_midref(v))
|| !mag_close(arb_radref(z), arb_radref(v)))
if (!arb_equal_mid_close_mag(z, v))
{
flint_printf("FAIL (aliasing 1)!\n");
flint_printf("x = "); arb_print(x); flint_printf("\n\n");
Expand All @@ -248,8 +259,7 @@ int main()
arb_addmul(v, x, x, prec);
arb_addmul(x, x, x, prec);

if (!arf_equal(arb_midref(x), arb_midref(v))
|| !mag_close(arb_radref(x), arb_radref(v)))
if (!arb_equal_mid_close_mag(x, v))
{
flint_printf("FAIL (aliasing 2)!\n");
flint_printf("x = "); arb_print(x); flint_printf("\n\n");
Expand All @@ -263,8 +273,7 @@ int main()
arb_addmul(v, x, y, prec);
arb_addmul(x, x, y, prec);

if (!arf_equal(arb_midref(x), arb_midref(v))
|| !mag_close(arb_radref(x), arb_radref(v)))
if (!arb_equal_mid_close_mag(x, v))
{
flint_printf("FAIL (aliasing 3)!\n");
flint_printf("x = "); arb_print(x); flint_printf("\n\n");
Expand All @@ -279,8 +288,7 @@ int main()
arb_addmul(v, x, y, prec);
arb_addmul(x, y, x, prec);

if (!arf_equal(arb_midref(x), arb_midref(v))
|| !mag_close(arb_radref(x), arb_radref(v)))
if (!arb_equal_mid_close_mag(x, v))
{
flint_printf("FAIL (aliasing 4)!\n");
flint_printf("x = "); arb_print(x); flint_printf("\n\n");
Expand Down
28 changes: 21 additions & 7 deletions arb/test/t-mul.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,19 @@ mag_close(const mag_t am, const mag_t bm)
return res1 && res2;
}

int
arb_equal_mid_close_mag(const arb_t a, const arb_t b)
{
return arf_equal(arb_midref(a), arb_midref(b)) &&
(mag_close(arb_radref(a), arb_radref(b)) ||
/* If a's and b's centres are infinite but their radii are finite, the radii don't need
to be close: they represent signed infinity regardless. If their centres are NaN,
then we should ignore their radii. */
(arf_is_inf(arb_midref(a)) && arf_is_inf(arb_midref(b)) &&
mag_is_finite(arb_radref(a)) && mag_is_finite(arb_radref(b))) ||
(arf_is_nan(arb_midref(a)) && arf_is_nan(arb_midref(b))));
}

void
arb_mul_naive(arb_t z, const arb_t x, const arb_t y, slong prec)
{
Expand Down Expand Up @@ -73,15 +86,18 @@ arb_mul_naive(arb_t z, const arb_t x, const arb_t y, slong prec)
fmpz_clear(e);
}

/* propagated error */
if (!arb_is_exact(x))
/* propagated error - note that (signed infinity) * nonzero should be a signed
infinity, meaning the error should *not* propagate */
if (!arb_is_exact(x) && !(
arb_is_nonzero(x) && arf_is_inf(arb_midref(y)) && mag_is_finite(arb_radref(y))))
{
arf_set_mag(t, arb_radref(x));
arf_abs(u, arb_midref(y));
arf_addmul(zr, t, u, MAG_BITS, ARF_RND_UP);
}

if (!arb_is_exact(y))
if (!arb_is_exact(y) && !(
arb_is_nonzero(y) && arf_is_inf(arb_midref(x)) && mag_is_finite(arb_radref(x))))
{
arf_set_mag(t, arb_radref(y));
arf_abs(u, arb_midref(x));
Expand Down Expand Up @@ -265,8 +281,7 @@ int main()
arb_mul(z, x, y, prec);
arb_mul_naive(v, x, y, prec);

if (!arf_equal(arb_midref(z), arb_midref(v))
|| !mag_close(arb_radref(z), arb_radref(v)))
if (!arb_equal_mid_close_mag(z, v))
{
flint_printf("FAIL!\n");
flint_printf("x = "); arb_print(x); flint_printf("\n\n");
Expand All @@ -282,8 +297,7 @@ int main()
arb_mul(z, x, y, prec);
arb_mul(v, x, x, prec);

if (!arf_equal(arb_midref(z), arb_midref(v))
|| !mag_close(arb_radref(z), arb_radref(v)))
if (!arb_equal_mid_close_mag(z, v))
{
flint_printf("FAIL (aliasing 1)!\n");
flint_printf("x = "); arb_print(x); flint_printf("\n\n");
Expand Down
Loading

0 comments on commit ab5e81c

Please sign in to comment.