From b17c5a700cabad28771df8f4d22ae92d1ed59743 Mon Sep 17 00:00:00 2001 From: p15-git-acc <37548430+p15-git-acc@users.noreply.github.com> Date: Mon, 21 Sep 2020 18:26:25 -0500 Subject: [PATCH] inf vs nan --- arb/cot_pi.c | 25 ++++++-- arb/log_hypot.c | 16 ++++- arb/sinc_pi.c | 11 ++++ arb/tan_pi.c | 24 ++++++-- arb/test/t-cot_pi.c | 131 +++++++++++++++++++++++++++++++++++++++++ arb/test/t-log_hypot.c | 19 ++++++ arb/test/t-sinc.c | 47 +++++++++++++++ arb/test/t-sinc_pi.c | 47 +++++++++++++++ arb/test/t-tan_pi.c | 126 +++++++++++++++++++++++++++++++++++++++ 9 files changed, 434 insertions(+), 12 deletions(-) create mode 100644 arb/test/t-cot_pi.c create mode 100644 arb/test/t-tan_pi.c diff --git a/arb/cot_pi.c b/arb/cot_pi.c index 0f504e32..0ccac5e7 100644 --- a/arb/cot_pi.c +++ b/arb/cot_pi.c @@ -14,10 +14,25 @@ void arb_cot_pi(arb_t y, const arb_t x, slong prec) { - arb_t u; - arb_init(u); - arb_sin_cos_pi(y, u, x, prec + 4); - arb_div(y, u, y, prec); - arb_clear(u); + + if (!arb_is_finite(x)) + { + arb_indeterminate(y); + } + else if (arb_is_int_2exp_si(x, -1)) + { + if (arb_is_int(x)) + arb_indeterminate(y); + else + arb_zero(y); + } + else + { + arb_t u; + arb_init(u); + arb_sin_cos_pi(y, u, x, prec + 4); + arb_div(y, u, y, prec); + arb_clear(u); + } } diff --git a/arb/log_hypot.c b/arb/log_hypot.c index aa965c2a..93671878 100644 --- a/arb/log_hypot.c +++ b/arb/log_hypot.c @@ -68,11 +68,23 @@ arb_log_hypot(arb_t res, const arb_t a, const arb_t b, slong prec) if (!arb_is_finite(a) || !arb_is_finite(b)) { - arb_indeterminate(res); + if (arf_is_nan(arb_midref(a)) || arf_is_nan(arb_midref(b))) + { + arb_indeterminate(res); + } + else if ((!arb_is_finite(a) && !arb_contains_zero(a)) || + (!arb_is_finite(b) && !arb_contains_zero(b))) + { + arb_pos_inf(res); + } + else + { + arb_indeterminate(res); + } return; } - /* a close to 1 -- for accurate acb_log1p */ + /* a close to 1 -- for accurate arb_log1p */ if (mag_cmp_2exp_si(arb_radref(a), -3) < 0 && mag_cmp_2exp_si(arb_radref(b), -3) < 0 && arf_cmpabs_2exp_si(arb_midref(b), -3) < 0 && diff --git a/arb/sinc_pi.c b/arb/sinc_pi.c index b9b87621..dc7ae847 100644 --- a/arb/sinc_pi.c +++ b/arb/sinc_pi.c @@ -17,6 +17,17 @@ arb_sinc_pi(arb_t res, const arb_t x, slong prec) mag_t m; arb_t t; + if (!arb_is_finite(x)) + { + if (arf_is_nan(arb_midref(x))) + arb_indeterminate(res); + else if (arb_contains_zero(x)) + arb_zero_pm_one(res); + else + arb_zero(res); + return; + } + if (arb_is_zero(x)) { arb_one(res); diff --git a/arb/tan_pi.c b/arb/tan_pi.c index 48669c28..4f3f5b4e 100644 --- a/arb/tan_pi.c +++ b/arb/tan_pi.c @@ -14,10 +14,24 @@ void arb_tan_pi(arb_t y, const arb_t x, slong prec) { - arb_t u; - arb_init(u); - arb_sin_cos_pi(y, u, x, prec + 4); - arb_div(y, y, u, prec); - arb_clear(u); + if (!arb_is_finite(x)) + { + arb_indeterminate(y); + } + else if (arb_is_int_2exp_si(x, -1)) + { + if (arb_is_int(x)) + arb_zero(y); + else + arb_indeterminate(y); + } + else + { + arb_t u; + arb_init(u); + arb_sin_cos_pi(y, u, x, prec + 4); + arb_div(y, y, u, prec); + arb_clear(u); + } } diff --git a/arb/test/t-cot_pi.c b/arb/test/t-cot_pi.c new file mode 100644 index 00000000..a37b5483 --- /dev/null +++ b/arb/test/t-cot_pi.c @@ -0,0 +1,131 @@ +/* + Copyright (C) 2013 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "arb.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("cot_pi...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++) + { + arb_t x, y, a, b, c; + slong prec1, prec2; + + prec1 = 2 + n_randint(state, 1000); + prec2 = prec1 + 30; + + arb_init(x); + arb_init(y); + arb_init(a); + arb_init(b); + arb_init(c); + + arb_randtest(x, state, 1 + n_randint(state, 1000), 2 + n_randint(state, 100)); + + arb_cot_pi(a, x, prec1); + arb_cot_pi(b, x, prec2); + + /* check consistency */ + if (!arb_overlaps(a, b)) + { + flint_printf("FAIL: overlap\n\n"); + flint_printf("x = "); arb_print(x); flint_printf("\n\n"); + flint_printf("a = "); arb_print(a); flint_printf("\n\n"); + flint_printf("b = "); arb_print(b); flint_printf("\n\n"); + flint_abort(); + } + + /* compare with cot */ + arb_const_pi(c, prec1); + arb_mul(y, x, c, prec1); + arb_cot(c, y, prec1); + + if (!arb_overlaps(a, c)) + { + flint_printf("FAIL: functional equation\n\n"); + flint_printf("x = "); arb_print(x); flint_printf("\n\n"); + flint_printf("y = "); arb_print(y); flint_printf("\n\n"); + flint_printf("a = "); arb_print(a); flint_printf("\n\n"); + flint_printf("c = "); arb_print(c); flint_printf("\n\n"); + flint_abort(); + } + + arb_cot_pi(x, x, prec1); + + if (!arb_overlaps(a, x)) + { + flint_printf("FAIL: aliasing\n\n"); + flint_printf("a = "); arb_print(a); flint_printf("\n\n"); + flint_printf("x = "); arb_print(x); flint_printf("\n\n"); + flint_abort(); + } + + arb_clear(x); + arb_clear(y); + arb_clear(a); + arb_clear(b); + arb_clear(c); + } + + /* Check a few special points. */ + { + arb_t x, y; + slong i, prec; + + arb_init(x); + arb_init(y); + + prec = 32; + + for (i = -1; i <= 1; i++) + { + int arbitrary_integer = 7; + + /* integer */ + arb_set_d(x, arbitrary_integer*i); + arb_cot_pi(y, x, prec); + if (arb_is_finite(y)) + { + flint_printf("FAIL: (integer)\n"); + flint_printf("x = "); arb_printd(x, 30); flint_printf("\n\n"); + flint_printf("y = "); arb_printd(y, 30); flint_printf("\n\n"); + flint_abort(); + } + + /* integer + 1/2 */ + arb_set_d(x, arbitrary_integer*i + 0.5); + arb_cot_pi(y, x, prec); + if (!arb_is_zero(y)) + { + flint_printf("FAIL: (integer + 1/2)\n"); + flint_printf("x = "); arb_printd(x, 30); flint_printf("\n\n"); + flint_printf("y = "); arb_printd(y, 30); flint_printf("\n\n"); + flint_abort(); + } + } + + arb_clear(x); + arb_clear(y); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/arb/test/t-log_hypot.c b/arb/test/t-log_hypot.c index 3b93e913..46624181 100644 --- a/arb/test/t-log_hypot.c +++ b/arb/test/t-log_hypot.c @@ -94,6 +94,25 @@ int main() flint_abort(); } + if (!arf_is_nan(arb_midref(a)) && !arf_is_nan(arb_midref(b))) + { + if ((!arb_is_finite(a) && !arb_contains_zero(a)) || + (!arb_is_finite(b) && !arb_contains_zero(b))) + { + if (arf_is_nan(arb_midref(r)) || arb_is_finite(r) || + !arb_is_positive(r)) + { + flint_printf("FAIL: infinite\n\n"); + flint_printf("prec1 = %wd, acc1 = %wd, acc2 = %wd\n\n", prec1, acc1, acc2); + flint_printf("a = "); arb_printn(a, 50, 0); flint_printf("\n\n"); + flint_printf("b = "); arb_printn(b, 50, 0); flint_printf("\n\n"); + flint_printf("r = "); arb_printn(r, 50, 0); flint_printf("\n\n"); + flint_printf("s = "); arb_printn(s, 50, 0); flint_printf("\n\n"); + flint_abort(); + } + } + } + arb_clear(a); arb_clear(b); arb_clear(r); diff --git a/arb/test/t-sinc.c b/arb/test/t-sinc.c index 78b7026b..c68f56f5 100644 --- a/arb/test/t-sinc.c +++ b/arb/test/t-sinc.c @@ -95,6 +95,53 @@ int main() arb_clear(a); } + /* Check a few special intervals. */ + { + arb_t x, y; + slong prec; + + arb_init(x); + arb_init(y); + + prec = 32; + + arb_neg_inf(x); + arb_sinc(y, x, prec); + if (!arb_is_zero(y)) + { + flint_printf("FAIL: negative infinity\n"); + flint_printf("x = "); arb_printd(x, 30); flint_printf("\n\n"); + flint_printf("y = "); arb_printd(y, 30); flint_printf("\n\n"); + flint_abort(); + } + + arb_pos_inf(x); + arb_sinc(y, x, prec); + if (!arb_is_zero(y)) + { + flint_printf("FAIL: positive infinity\n"); + flint_printf("x = "); arb_printd(x, 30); flint_printf("\n\n"); + flint_printf("y = "); arb_printd(y, 30); flint_printf("\n\n"); + flint_abort(); + } + + arb_zero_pm_inf(x); + arb_sinc(y, x, prec); + if (!(arb_is_finite(y) && + arb_contains_negative(y) && + arb_contains_positive(y) && + arb_contains_zero(y))) + { + flint_printf("FAIL: the whole extended real line\n"); + flint_printf("x = "); arb_printd(x, 30); flint_printf("\n\n"); + flint_printf("y = "); arb_printd(y, 30); flint_printf("\n\n"); + flint_abort(); + } + + arb_clear(x); + arb_clear(y); + } + flint_randclear(state); flint_cleanup(); flint_printf("PASS\n"); diff --git a/arb/test/t-sinc_pi.c b/arb/test/t-sinc_pi.c index 871b8213..4bd92086 100644 --- a/arb/test/t-sinc_pi.c +++ b/arb/test/t-sinc_pi.c @@ -61,6 +61,53 @@ int main() arb_clear(z); } + /* Check a few special intervals. */ + { + arb_t x, y; + slong prec; + + arb_init(x); + arb_init(y); + + prec = 32; + + arb_neg_inf(x); + arb_sinc_pi(y, x, prec); + if (!arb_is_zero(y)) + { + flint_printf("FAIL: negative infinity\n"); + flint_printf("x = "); arb_printd(x, 30); flint_printf("\n\n"); + flint_printf("y = "); arb_printd(y, 30); flint_printf("\n\n"); + flint_abort(); + } + + arb_pos_inf(x); + arb_sinc_pi(y, x, prec); + if (!arb_is_zero(y)) + { + flint_printf("FAIL: positive infinity\n"); + flint_printf("x = "); arb_printd(x, 30); flint_printf("\n\n"); + flint_printf("y = "); arb_printd(y, 30); flint_printf("\n\n"); + flint_abort(); + } + + arb_zero_pm_inf(x); + arb_sinc_pi(y, x, prec); + if (!(arb_is_finite(y) && + arb_contains_negative(y) && + arb_contains_positive(y) && + arb_contains_zero(y))) + { + flint_printf("FAIL: the whole extended real line\n"); + flint_printf("x = "); arb_printd(x, 30); flint_printf("\n\n"); + flint_printf("y = "); arb_printd(y, 30); flint_printf("\n\n"); + flint_abort(); + } + + arb_clear(x); + arb_clear(y); + } + flint_randclear(state); flint_cleanup(); flint_printf("PASS\n"); diff --git a/arb/test/t-tan_pi.c b/arb/test/t-tan_pi.c new file mode 100644 index 00000000..3fbfad3e --- /dev/null +++ b/arb/test/t-tan_pi.c @@ -0,0 +1,126 @@ +/* + Copyright (C) 2013 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "arb.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("tan_pi...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++) + { + arb_t x, y, a, b, c; + slong prec1, prec2; + + prec1 = 2 + n_randint(state, 1000); + prec2 = prec1 + 30; + + arb_init(x); + arb_init(y); + arb_init(a); + arb_init(b); + arb_init(c); + + arb_randtest(x, state, 1 + n_randint(state, 1000), 2 + n_randint(state, 100)); + + arb_tan_pi(a, x, prec1); + arb_tan_pi(b, x, prec2); + + /* check consistency */ + if (!arb_overlaps(a, b)) + { + flint_printf("FAIL: overlap\n\n"); + flint_printf("x = "); arb_print(x); flint_printf("\n\n"); + flint_printf("a = "); arb_print(a); flint_printf("\n\n"); + flint_printf("b = "); arb_print(b); flint_printf("\n\n"); + flint_abort(); + } + + /* compare with tan */ + arb_const_pi(c, prec1); + arb_mul(y, x, c, prec1); + arb_tan(c, y, prec1); + + if (!arb_overlaps(a, c)) + { + flint_printf("FAIL: functional equation\n\n"); + flint_printf("x = "); arb_print(x); flint_printf("\n\n"); + flint_printf("y = "); arb_print(y); flint_printf("\n\n"); + flint_printf("a = "); arb_print(a); flint_printf("\n\n"); + flint_printf("c = "); arb_print(c); flint_printf("\n\n"); + flint_abort(); + } + + arb_tan_pi(x, x, prec1); + + if (!arb_overlaps(a, x)) + { + flint_printf("FAIL: aliasing\n\n"); + flint_printf("a = "); arb_print(a); flint_printf("\n\n"); + flint_printf("x = "); arb_print(x); flint_printf("\n\n"); + flint_abort(); + } + + arb_clear(x); + arb_clear(y); + arb_clear(a); + arb_clear(b); + arb_clear(c); + } + + /* Check a few special points. */ + { + arb_t x, y; + slong i, prec; + + arb_init(x); + arb_init(y); + + prec = 32; + + for (i = -1; i <= 1; i++) + { + int arbitrary_integer = 7; + + arb_set_d(x, arbitrary_integer*i); + arb_tan_pi(y, x, prec); + if (!arb_is_zero(y)) + { + flint_printf("FAIL: (integer)\n"); + flint_printf("x = "); arb_printd(x, 30); flint_printf("\n\n"); + flint_printf("y = "); arb_printd(y, 30); flint_printf("\n\n"); + flint_abort(); + } + + arb_set_d(x, arbitrary_integer*i + 0.5); + arb_tan_pi(y, x, prec); + if (arb_is_finite(y)) + { + flint_printf("FAIL: (integer + 1/2)\n"); + flint_printf("x = "); arb_printd(x, 30); flint_printf("\n\n"); + flint_printf("y = "); arb_printd(y, 30); flint_printf("\n\n"); + flint_abort(); + } + } + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} +