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

Positive * pos_inf should be pos_inf, not 0 +- inf. Same for other combinations of signs. This addresses github issue #408.

Browse source
This commit is contained in:
postmath 2022-03-11 22:58:19 -05:00
parent 076249eaee
commit 89b489e64d
7 changed files with 139 additions and 27 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

21
arb/addmul.c
View file

@ -24,6 +24,13 @@ arb_addmul_arf(arb_t z, const arb_t x, const arf_t y, slong prec)
if (inexact)
arf_mag_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 if (ARB_IS_LAGOM(x) && ARF_IS_LAGOM(y) && ARB_IS_LAGOM(z))
{
mag_fast_init_set_arf(ym, y);
@ -60,6 +67,20 @@ arb_addmul(arb_t z, const arb_t x, const arb_t y, slong prec)
{
arb_addmul_arf(z, y, arb_midref(x), prec);
}
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 if (ARB_IS_LAGOM(x) && ARB_IS_LAGOM(y) && ARB_IS_LAGOM(z))
{
mag_fast_init_set_arf(xm, arb_midref(x));

21
arb/fma.c
View file

@ -26,6 +26,13 @@ arb_fma_arf(arb_t res, const arb_t x, const arf_t y, const arb_t z, slong prec)
else
mag_set(arb_radref(res), arb_radref(z));
}
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 if (ARB_IS_LAGOM(res) && ARB_IS_LAGOM(x) && ARF_IS_LAGOM(y) && ARB_IS_LAGOM(z))
{
mag_t tm;
@ -73,6 +80,20 @@ arb_fma(arb_t res, const arb_t x, const arb_t y, const arb_t z, slong prec)
{
arb_fma_arf(res, y, arb_midref(x), z, prec);
}
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 if (ARB_IS_LAGOM(res) && ARB_IS_LAGOM(x) && ARB_IS_LAGOM(y) && ARB_IS_LAGOM(z))
{
mag_fast_init_set_arf(xm, arb_midref(x));

19
arb/mul.c
View file

@ -26,6 +26,15 @@ arb_mul_arf(arb_t z, const arb_t x, const arf_t y, slong prec)
else
mag_zero(arb_radref(z));
}
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 if (ARB_IS_LAGOM(x) && ARF_IS_LAGOM(y) && ARB_IS_LAGOM(z))
{
mag_fast_init_set_arf(ym, y);
@ -74,6 +83,16 @@ arb_mul(arb_t z, const arb_t x, const arb_t y, slong prec)
{
arb_mul_arf(z, x, arb_midref(y), prec);
}
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 if (ARB_IS_LAGOM(x) && ARB_IS_LAGOM(y) && ARB_IS_LAGOM(z))
{
mag_fast_init_set_arf(xm, arb_midref(x));

21
arb/submul.c
View file

@ -24,6 +24,13 @@ arb_submul_arf(arb_t z, const arb_t x, const arf_t y, slong prec)
if (inexact)
arf_mag_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 if (ARB_IS_LAGOM(x) && ARF_IS_LAGOM(y) && ARB_IS_LAGOM(z))
{
mag_fast_init_set_arf(ym, y);
@ -60,6 +67,20 @@ arb_submul(arb_t z, const arb_t x, const arb_t y, slong prec)
{
arb_submul_arf(z, y, arb_midref(x), prec);
}
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 if (ARB_IS_LAGOM(x) && ARB_IS_LAGOM(y) && ARB_IS_LAGOM(z))
{
mag_fast_init_set_arf(xm, arb_midref(x));

28
arb/test/t-addmul.c
View file

@ -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)
{
@ -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");
@ -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");
@ -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");
@ -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");
@ -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");

28
arb/test/t-mul.c
View file

@ -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)
{
@ -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));
@ -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");
@ -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");

28
arb/test/t-submul.c
View file

@ -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_submul_naive(arb_t z, const arb_t x, const arb_t y, slong prec)
{
@ -214,8 +227,7 @@ int main()
arb_submul(z, x, y, prec);
arb_submul_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");
@ -232,8 +244,7 @@ int main()
arb_submul(z, x, y, prec);
arb_submul(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");
@ -248,8 +259,7 @@ int main()
arb_submul(v, x, x, prec);
arb_submul(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");
@ -263,8 +273,7 @@ int main()
arb_submul(v, x, y, prec);
arb_submul(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");
@ -279,8 +288,7 @@ int main()
arb_submul(v, x, y, prec);
arb_submul(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");