From 89b489e64d30224fea9659cb078fdd4195db789d Mon Sep 17 00:00:00 2001 From: postmath Date: Fri, 11 Mar 2022 22:58:19 -0500 Subject: [PATCH 1/2] Positive * pos_inf should be pos_inf, not 0 +- inf. Same for other combinations of signs. This addresses github issue #408. --- arb/addmul.c | 21 +++++++++++++++++++++ arb/fma.c | 21 +++++++++++++++++++++ arb/mul.c | 19 +++++++++++++++++++ arb/submul.c | 21 +++++++++++++++++++++ arb/test/t-addmul.c | 28 ++++++++++++++++++---------- arb/test/t-mul.c | 28 +++++++++++++++++++++------- arb/test/t-submul.c | 28 ++++++++++++++++++---------- 7 files changed, 139 insertions(+), 27 deletions(-) diff --git a/arb/addmul.c b/arb/addmul.c index 4b85a885..d40e91d6 100644 --- a/arb/addmul.c +++ b/arb/addmul.c @@ -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)); diff --git a/arb/fma.c b/arb/fma.c index 5abac1ba..54e8f1d0 100644 --- a/arb/fma.c +++ b/arb/fma.c @@ -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)); diff --git a/arb/mul.c b/arb/mul.c index 12096c5d..cac13db1 100644 --- a/arb/mul.c +++ b/arb/mul.c @@ -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)); diff --git a/arb/submul.c b/arb/submul.c index 001f8d8e..f364206b 100644 --- a/arb/submul.c +++ b/arb/submul.c @@ -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)); diff --git a/arb/test/t-addmul.c b/arb/test/t-addmul.c index 14b0977a..eee1d113 100644 --- a/arb/test/t-addmul.c +++ b/arb/test/t-addmul.c @@ -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"); diff --git a/arb/test/t-mul.c b/arb/test/t-mul.c index 9d342868..aee343d0 100644 --- a/arb/test/t-mul.c +++ b/arb/test/t-mul.c @@ -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"); diff --git a/arb/test/t-submul.c b/arb/test/t-submul.c index bf2e9ff0..25cd02ee 100644 --- a/arb/test/t-submul.c +++ b/arb/test/t-submul.c @@ -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"); From a9f88223619889b74887fefa711371f4d9b89353 Mon Sep 17 00:00:00 2001 From: postmath Date: Mon, 14 Mar 2022 17:05:14 -0400 Subject: [PATCH 2/2] Optimization: LAGOM is finite. Add test. --- arb/addmul.c | 42 ++--- arb/fma.c | 42 ++--- arb/mul.c | 38 ++--- arb/submul.c | 42 ++--- arb/test/t-pos_times_posinf.c | 283 ++++++++++++++++++++++++++++++++++ 5 files changed, 365 insertions(+), 82 deletions(-) create mode 100644 arb/test/t-pos_times_posinf.c diff --git a/arb/addmul.c b/arb/addmul.c index d40e91d6..85007088 100644 --- a/arb/addmul.c +++ b/arb/addmul.c @@ -24,13 +24,6 @@ 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); @@ -40,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); @@ -67,20 +67,6 @@ 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)); @@ -99,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)); diff --git a/arb/fma.c b/arb/fma.c index 54e8f1d0..d00a5e04 100644 --- a/arb/fma.c +++ b/arb/fma.c @@ -26,13 +26,6 @@ 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; @@ -46,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; @@ -80,20 +80,6 @@ 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)); @@ -112,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)); diff --git a/arb/mul.c b/arb/mul.c index cac13db1..4c491fbc 100644 --- a/arb/mul.c +++ b/arb/mul.c @@ -26,15 +26,6 @@ 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); @@ -49,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); @@ -83,16 +83,6 @@ 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)); @@ -110,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)); diff --git a/arb/submul.c b/arb/submul.c index f364206b..ad9adbf9 100644 --- a/arb/submul.c +++ b/arb/submul.c @@ -24,13 +24,6 @@ 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); @@ -40,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); @@ -67,20 +67,6 @@ 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)); @@ -99,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)); diff --git a/arb/test/t-pos_times_posinf.c b/arb/test/t-pos_times_posinf.c new file mode 100644 index 00000000..026d9e30 --- /dev/null +++ b/arb/test/t-pos_times_posinf.c @@ -0,0 +1,283 @@ +#include +#include "arb.h" + +#define PRINT_PRECISION 20 + +typedef enum { + ALL_POSITIVE = 0, + CONTAINS_ZERO = 1, + ALL_NEGATIVE = 2, + ANY_FINITE = 3, /* used only to represent either of the previous three */ + POSITIVE_INFINITY = 4, + NEGATIVE_INFINITY = 5, + WHOLE_LINE = 6, + ANY_INFINITE = 7, /* used only to represent either of the previous three */ + NOT_A_NUMBER = 8, + NUM_VALUE_TYPES = 9 /* used only to count the number of enum values */ +} value_type; + +typedef struct { + const char *value; + value_type type; +} string_with_type; + +const string_with_type test_values_arb[] = { + {"5.3 +/- 1.2", ALL_POSITIVE}, + {"1.2 +/- 5.3", CONTAINS_ZERO}, + {"-1.2 +/- 5.3", CONTAINS_ZERO}, + {"-5.3 +/- 1.2", ALL_NEGATIVE}, + {"inf +/- 17", POSITIVE_INFINITY}, + {"-inf +/- 17", NEGATIVE_INFINITY}, + {"17 +/- inf", WHOLE_LINE}, + {"-inf +/- inf", WHOLE_LINE}, + {"nan +/- 17", NOT_A_NUMBER}, + {"nan +/- inf", NOT_A_NUMBER}, + {NULL, NOT_A_NUMBER}, +}; + +const string_with_type test_values_arf[] = { + {"5.3", ALL_POSITIVE}, + {"0", CONTAINS_ZERO}, + {"-5.3", ALL_NEGATIVE}, + {"inf", POSITIVE_INFINITY}, + {"-inf", NEGATIVE_INFINITY}, + {"nan", NOT_A_NUMBER}, + {NULL, NOT_A_NUMBER}, +}; + +/* Multiplying value_type i with value_type j should yield value_type multiplication[i][j]. */ +const value_type multiplication[NUM_VALUE_TYPES][NUM_VALUE_TYPES] = { + {ALL_POSITIVE, CONTAINS_ZERO, ALL_NEGATIVE, ANY_FINITE, POSITIVE_INFINITY,NEGATIVE_INFINITY,WHOLE_LINE, ANY_INFINITE, NOT_A_NUMBER}, + {CONTAINS_ZERO, CONTAINS_ZERO, CONTAINS_ZERO, ANY_FINITE, WHOLE_LINE, WHOLE_LINE, WHOLE_LINE, ANY_INFINITE, NOT_A_NUMBER}, + {ALL_NEGATIVE, CONTAINS_ZERO, ALL_POSITIVE, ANY_FINITE, NEGATIVE_INFINITY,POSITIVE_INFINITY,WHOLE_LINE, ANY_INFINITE, NOT_A_NUMBER}, + {ANY_FINITE, ANY_FINITE, ANY_FINITE, ANY_FINITE, ANY_INFINITE, ANY_INFINITE, WHOLE_LINE, ANY_INFINITE, NOT_A_NUMBER}, + {POSITIVE_INFINITY,WHOLE_LINE, NEGATIVE_INFINITY,ANY_INFINITE, POSITIVE_INFINITY,NEGATIVE_INFINITY,WHOLE_LINE, ANY_INFINITE, NOT_A_NUMBER}, + {NEGATIVE_INFINITY,WHOLE_LINE, POSITIVE_INFINITY,ANY_INFINITE, NEGATIVE_INFINITY,POSITIVE_INFINITY,WHOLE_LINE, ANY_INFINITE, NOT_A_NUMBER}, + {WHOLE_LINE, WHOLE_LINE, WHOLE_LINE, WHOLE_LINE, WHOLE_LINE, WHOLE_LINE, WHOLE_LINE, WHOLE_LINE, NOT_A_NUMBER}, + {ANY_INFINITE, ANY_INFINITE, ANY_INFINITE, ANY_INFINITE, ANY_INFINITE, ANY_INFINITE, WHOLE_LINE, ANY_INFINITE, NOT_A_NUMBER}, + {NOT_A_NUMBER, NOT_A_NUMBER, NOT_A_NUMBER, NOT_A_NUMBER, NOT_A_NUMBER, NOT_A_NUMBER, NOT_A_NUMBER, NOT_A_NUMBER, NOT_A_NUMBER}, +}; + +/* Adding value_type i to value_type j should yield value_type addition[i][j]. */ +const value_type addition[NUM_VALUE_TYPES][NUM_VALUE_TYPES] = { + {ALL_POSITIVE, ANY_FINITE, ANY_FINITE, ANY_FINITE, POSITIVE_INFINITY,NEGATIVE_INFINITY,WHOLE_LINE, ANY_INFINITE, NOT_A_NUMBER}, + {ANY_FINITE, CONTAINS_ZERO, ANY_FINITE, ANY_FINITE, POSITIVE_INFINITY,NEGATIVE_INFINITY,WHOLE_LINE, ANY_INFINITE, NOT_A_NUMBER}, + {ANY_FINITE, ANY_FINITE, ALL_NEGATIVE, ANY_FINITE, POSITIVE_INFINITY,NEGATIVE_INFINITY,WHOLE_LINE, ANY_INFINITE, NOT_A_NUMBER}, + {ANY_FINITE, ANY_FINITE, ANY_FINITE, ANY_FINITE, POSITIVE_INFINITY,NEGATIVE_INFINITY,WHOLE_LINE, ANY_INFINITE, NOT_A_NUMBER}, + {POSITIVE_INFINITY,POSITIVE_INFINITY,POSITIVE_INFINITY,POSITIVE_INFINITY,POSITIVE_INFINITY,WHOLE_LINE, WHOLE_LINE, ANY_INFINITE, NOT_A_NUMBER}, + {NEGATIVE_INFINITY,NEGATIVE_INFINITY,NEGATIVE_INFINITY,NEGATIVE_INFINITY,WHOLE_LINE, NEGATIVE_INFINITY,WHOLE_LINE, ANY_INFINITE, NOT_A_NUMBER}, + {WHOLE_LINE, WHOLE_LINE, WHOLE_LINE, WHOLE_LINE, WHOLE_LINE, WHOLE_LINE, WHOLE_LINE, WHOLE_LINE, NOT_A_NUMBER}, + {ANY_INFINITE, ANY_INFINITE, ANY_INFINITE, ANY_INFINITE, ANY_INFINITE, ANY_INFINITE, WHOLE_LINE, ANY_INFINITE, NOT_A_NUMBER}, + {NOT_A_NUMBER, NOT_A_NUMBER, NOT_A_NUMBER, NOT_A_NUMBER, NOT_A_NUMBER, NOT_A_NUMBER, NOT_A_NUMBER, NOT_A_NUMBER, NOT_A_NUMBER}, +}; + +int arb_satisfies_value_type(const arb_t x, const value_type t) +{ + switch(t) { + case ALL_POSITIVE: + return arb_is_finite(x) && arb_is_positive(x); + case CONTAINS_ZERO: + return arb_is_finite(x) && arb_contains_zero(x); + case ALL_NEGATIVE: + return arb_is_finite(x) && arb_is_negative(x); + case ANY_FINITE: + return arb_is_finite(x); + case POSITIVE_INFINITY: + return !arb_is_finite(x) && arb_is_positive(x); + case NEGATIVE_INFINITY: + return !arb_is_finite(x) && arb_is_negative(x); + case WHOLE_LINE: + return !arb_is_finite(x) && arb_contains_zero(x); + case ANY_INFINITE: + return !arb_is_finite(x) && !arf_is_nan(arb_midref(x)); + case NOT_A_NUMBER: + return arf_is_nan(arb_midref(x)); + default: + printf("unexpected condition\n"); + flint_abort(); + } +} + +void print_arb_and_type(arb_t x, const char *s, const value_type t) +{ + flint_printf("%s: ", s); + arb_printd(x, PRINT_PRECISION); + flint_printf(" of type: %d\n", t); +} +void print_arf_and_type(arf_t x, const char *s, const value_type t) +{ + flint_printf("%s: ", s); + arf_printd(x, PRINT_PRECISION); + flint_printf(" of type: %d\n", t); +} + +int main() +{ + slong i, j, k, prec; + arb_t t, u, v, w; + arf_t x; + + flint_printf("pos_times_posinf...."); + fflush(stdout); + + arb_init(t); + arb_init(u); + arb_init(v); + arb_init(w); + arf_init(x); + + for(prec = 20; prec <= 100; prec += 80) + for(i = 0; test_values_arb[i].value != NULL; ++i) + { + for(j = 0; test_values_arb[j].value != NULL; ++j) + { + value_type vt = test_values_arb[i].type, vu = test_values_arb[j].type, expected = multiplication[vt][vu]; + arb_set_str(t, test_values_arb[i].value, prec); + arb_set_str(u, test_values_arb[j].value, prec); + + arb_mul(v, t, u, prec); + + if(!arb_satisfies_value_type(v, expected)) + { + flint_printf("FAIL (arb_mul): %s, %s\n", test_values_arb[i].value, test_values_arb[j].value); + print_arb_and_type(t, "t", vt); + print_arb_and_type(u, "u", vu); + flint_printf("v: "); + arb_printd(v, PRINT_PRECISION); + flint_printf("\nexpected type: %d\n", expected); + flint_abort(); + } + + for(k = 0; test_values_arb[k].value != NULL; ++k) + { + value_type vw = test_values_arb[k].type, expected = addition[vw][multiplication[vt][vu]]; + arb_set_str(w, test_values_arb[k].value, prec); + + arb_addmul(w, t, u, prec); + + if(!arb_satisfies_value_type(w, expected)) + { + flint_printf("FAIL (arb_addmul): %s, %s, %s\n", test_values_arb[k].value, test_values_arb[i].value, test_values_arb[j].value); + print_arb_and_type(t, "t", vt); + print_arb_and_type(u, "u", vu); + flint_printf("w (after): "); + arb_printd(w, PRINT_PRECISION); + flint_printf("\nexpected type: %d\n", expected); + flint_abort(); + } + + expected = addition[vw][multiplication[ALL_NEGATIVE][multiplication[vt][vu]]]; + arb_set_str(w, test_values_arb[k].value, prec); + + arb_submul(w, t, u, prec); + + if(!arb_satisfies_value_type(w, expected)) + { + flint_printf("FAIL (arb_submul): %s, %s, %s\n", test_values_arb[k].value, test_values_arb[i].value, test_values_arb[j].value); + print_arb_and_type(t, "t", vt); + print_arb_and_type(u, "u", vu); + flint_printf("w (after): "); + arb_printd(w, PRINT_PRECISION); + flint_printf("\nexpected type: %d\n", expected); + flint_abort(); + } + + expected = addition[vw][multiplication[vt][vu]]; + arb_set_str(w, test_values_arb[k].value, prec); + + arb_fma(v, t, u, w, prec); + if(!arb_satisfies_value_type(v, expected)) + { + flint_printf("FAIL (arb_fma): %s, %s, %s\n", test_values_arb[k].value, test_values_arb[i].value, test_values_arb[j].value); + print_arb_and_type(t, "t", vt); + print_arb_and_type(u, "u", vu); + flint_printf("v (after): "); + arb_printd(v, PRINT_PRECISION); + flint_printf("\nexpected type: %d\n", expected); + flint_abort(); + } + } + } + + for(j = 0; test_values_arf[j].value != NULL; ++j) + { + value_type vt = test_values_arb[i].type, vx = test_values_arf[j].type; + arb_set_str(t, test_values_arb[i].value, prec); + arb_set_str(w, test_values_arf[j].value, prec); + arf_set(x, arb_midref(w)); + + arb_mul_arf(v, t, x, prec); + + if(!arb_satisfies_value_type(v, multiplication[vt][vx])) + { + flint_printf("FAIL (arb_mul_arf): %s, %s\n", test_values_arb[i].value, test_values_arf[j].value); + print_arb_and_type(t, "t", vt); + print_arf_and_type(x, "x", vx); + flint_printf("v: "); + arb_printd(v, PRINT_PRECISION); + flint_printf("\nexpected type: %d\n", multiplication[vt][vx]); + flint_abort(); + } + + for(k = 0; test_values_arb[k].value != NULL; ++k) + { + value_type vw = test_values_arb[k].type, expected = addition[vw][multiplication[vt][vx]]; + arb_set_str(w, test_values_arb[k].value, prec); + + arb_addmul_arf(w, t, x, prec); + + if(!arb_satisfies_value_type(w, expected)) + { + flint_printf("FAIL (arb_addmul_arf): %s, %s, %s\n", test_values_arb[k].value, test_values_arb[i].value, test_values_arf[j].value); + print_arb_and_type(t, "t", vt); + print_arf_and_type(x, "x", vx); + flint_printf("w (after): "); + arb_printd(w, PRINT_PRECISION); + flint_printf("\nexpected type: %d\n", expected); + flint_abort(); + } + + expected = addition[vw][multiplication[ALL_NEGATIVE][multiplication[vt][vx]]]; + arb_set_str(w, test_values_arb[k].value, prec); + + arb_submul_arf(w, t, x, prec); + + if(!arb_satisfies_value_type(w, expected)) + { + flint_printf("FAIL (arb_submul_arf): %s, %s, %s\n", test_values_arb[k].value, test_values_arb[i].value, test_values_arf[j].value); + print_arb_and_type(t, "t", vt); + print_arf_and_type(x, "x", vx); + flint_printf("w (after): "); + arb_printd(w, PRINT_PRECISION); + flint_printf("\nexpected type: %d\n", expected); + flint_abort(); + } + + expected = addition[vw][multiplication[vt][vx]]; + arb_set_str(w, test_values_arb[k].value, prec); + + arb_fma_arf(v, t, x, w, prec); + if(!arb_satisfies_value_type(v, expected)) + { + flint_printf("FAIL (arb_fma): %s, %s, %s\n", test_values_arb[k].value, test_values_arb[i].value, test_values_arf[j].value); + print_arb_and_type(t, "t", vt); + print_arf_and_type(x, "x", vx); + flint_printf("v (after): "); + arb_printd(v, PRINT_PRECISION); + flint_printf("\nexpected type: %d\n", expected); + flint_abort(); + } + } + } + } + + arf_clear(x); + arb_clear(w); + arb_clear(v); + arb_clear(u); + arb_clear(t); + + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +}