diff --git a/arb.h b/arb.h index 5cd649a9..5cea6896 100644 --- a/arb.h +++ b/arb.h @@ -535,6 +535,9 @@ void arb_union(arb_t z, const arb_t x, const arb_t y, slong prec); int arb_intersection(arb_t z, const arb_t x, const arb_t y, slong prec); void arb_get_rand_fmpq(fmpq_t q, flint_rand_t state, const arb_t x, slong bits); +void arb_min(arb_t z, const arb_t x, const arb_t y, slong prec); +void arb_max(arb_t z, const arb_t x, const arb_t y, slong prec); + int arb_can_round_arf(const arb_t x, slong prec, arf_rnd_t rnd); int arb_can_round_mpfr(const arb_t x, slong prec, mpfr_rnd_t rnd); diff --git a/arb/max.c b/arb/max.c new file mode 100644 index 00000000..7643ff52 --- /dev/null +++ b/arb/max.c @@ -0,0 +1,59 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2016 Arb authors + +******************************************************************************/ + +#include "arb.h" + +void +arb_max(arb_t z, const arb_t x, const arb_t y, slong prec) +{ + arf_t left, right, t, xr, yr; + + if (arf_is_nan(arb_midref(x)) || arf_is_nan(arb_midref(y))) + { + arb_indeterminate(z); + return; + } + + arf_init(left); + arf_init(right); + arf_init(t); + + arf_init_set_mag_shallow(xr, arb_radref(x)); + arf_init_set_mag_shallow(yr, arb_radref(y)); + + arf_sub(left, arb_midref(x), xr, prec, ARF_RND_FLOOR); + arf_sub(t, arb_midref(y), yr, prec, ARF_RND_FLOOR); + arf_max(left, left, t); + + arf_add(right, arb_midref(x), xr, prec, ARF_RND_CEIL); + arf_add(t, arb_midref(y), yr, prec, ARF_RND_CEIL); + arf_max(right, right, t); + + arb_set_interval_arf(z, left, right, prec); + + arf_clear(left); + arf_clear(right); + arf_clear(t); +} diff --git a/arb/min.c b/arb/min.c new file mode 100644 index 00000000..f752fefc --- /dev/null +++ b/arb/min.c @@ -0,0 +1,59 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2016 Arb authors + +******************************************************************************/ + +#include "arb.h" + +void +arb_min(arb_t z, const arb_t x, const arb_t y, slong prec) +{ + arf_t left, right, t, xr, yr; + + if (arf_is_nan(arb_midref(x)) || arf_is_nan(arb_midref(y))) + { + arb_indeterminate(z); + return; + } + + arf_init(left); + arf_init(right); + arf_init(t); + + arf_init_set_mag_shallow(xr, arb_radref(x)); + arf_init_set_mag_shallow(yr, arb_radref(y)); + + arf_sub(left, arb_midref(x), xr, prec, ARF_RND_FLOOR); + arf_sub(t, arb_midref(y), yr, prec, ARF_RND_FLOOR); + arf_min(left, left, t); + + arf_add(right, arb_midref(x), xr, prec, ARF_RND_CEIL); + arf_add(t, arb_midref(y), yr, prec, ARF_RND_CEIL); + arf_min(right, right, t); + + arb_set_interval_arf(z, left, right, prec); + + arf_clear(left); + arf_clear(right); + arf_clear(t); +} diff --git a/arb/test/t-max.c b/arb/test/t-max.c new file mode 100644 index 00000000..e24da5fc --- /dev/null +++ b/arb/test/t-max.c @@ -0,0 +1,142 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2016 Arb authors + +******************************************************************************/ + +#include "arb.h" + +/* sample (x, y) so that x /in y */ +void _sample_arf_in_arb(arf_t x, arb_t y, flint_rand_t state) +{ + slong bits, prec, expbits; + arf_t a, b; + slong i, n; + + arf_init(a); + arf_init(b); + + bits = 2 + n_randint(state, 1000); + prec = 2 + n_randint(state, 1000); + expbits = n_randint(state, 14); + n = n_randint(state, 3); + + arf_randtest(x, state, bits, expbits); + arf_set(a, x); + arf_set(b, x); + for (i = 0; i < n; i++) + { + arf_randtest(x, state, bits, expbits); + arf_min(a, a, x); + arf_max(b, b, x); + } + + arb_set_interval_arf(y, a, b, prec); + + arf_clear(a); + arf_clear(b); +} + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("max...."); + fflush(stdout); + flint_randinit(state); + + for (iter = 0; iter < 100000; iter++) + { + arf_t a, b, c; + arb_t x, y, z; + slong prec; + + arf_init(a); + arf_init(b); + arf_init(c); + + arb_init(x); + arb_init(y); + arb_init(z); + + _sample_arf_in_arb(a, x, state); + _sample_arf_in_arb(b, y, state); + + prec = 2 + n_randint(state, 200); + + arf_max(c, a, b); + arb_max(z, x, y, prec); + + if (!arb_contains_arf(x, a) || + !arb_contains_arf(y, b) || + !arb_contains_arf(z, c)) + { + flint_printf("FAIL: containment\n\n"); + flint_printf("a = "); arf_print(a); flint_printf("\n\n"); + flint_printf("b = "); arf_print(b); flint_printf("\n\n"); + flint_printf("c = "); arf_print(c); flint_printf("\n\n"); + flint_printf("x = "); arb_print(x); flint_printf("\n\n"); + flint_printf("y = "); arb_print(y); flint_printf("\n\n"); + flint_printf("z = "); arb_print(z); flint_printf("\n\n"); + abort(); + } + + /* aliasing */ + { + int alias; + + if (n_randint(state, 2)) + { + arb_max(x, x, y, prec); + alias = arb_equal(x, z); + } + else + { + arb_max(y, x, y, prec); + alias = arb_equal(y, z); + } + + if (!alias) + { + flint_printf("FAIL: aliasing\n\n"); + flint_printf("x = "); arb_print(x); flint_printf("\n\n"); + flint_printf("y = "); arb_print(y); flint_printf("\n\n"); + flint_printf("z = "); arb_print(z); flint_printf("\n\n"); + abort(); + } + } + + arf_clear(a); + arf_clear(b); + arf_clear(c); + + arb_clear(x); + arb_clear(y); + arb_clear(z); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/arb/test/t-min.c b/arb/test/t-min.c new file mode 100644 index 00000000..4d1d6c73 --- /dev/null +++ b/arb/test/t-min.c @@ -0,0 +1,142 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2016 Arb authors + +******************************************************************************/ + +#include "arb.h" + +/* sample (x, y) so that x /in y */ +void _sample_arf_in_arb(arf_t x, arb_t y, flint_rand_t state) +{ + slong bits, prec, expbits; + arf_t a, b; + slong i, n; + + arf_init(a); + arf_init(b); + + bits = 2 + n_randint(state, 1000); + prec = 2 + n_randint(state, 1000); + expbits = n_randint(state, 14); + n = n_randint(state, 3); + + arf_randtest(x, state, bits, expbits); + arf_set(a, x); + arf_set(b, x); + for (i = 0; i < n; i++) + { + arf_randtest(x, state, bits, expbits); + arf_min(a, a, x); + arf_max(b, b, x); + } + + arb_set_interval_arf(y, a, b, prec); + + arf_clear(a); + arf_clear(b); +} + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("min...."); + fflush(stdout); + flint_randinit(state); + + for (iter = 0; iter < 100000; iter++) + { + arf_t a, b, c; + arb_t x, y, z; + slong prec; + + arf_init(a); + arf_init(b); + arf_init(c); + + arb_init(x); + arb_init(y); + arb_init(z); + + _sample_arf_in_arb(a, x, state); + _sample_arf_in_arb(b, y, state); + + prec = 2 + n_randint(state, 200); + + arf_min(c, a, b); + arb_min(z, x, y, prec); + + if (!arb_contains_arf(x, a) || + !arb_contains_arf(y, b) || + !arb_contains_arf(z, c)) + { + flint_printf("FAIL: containment\n\n"); + flint_printf("a = "); arf_print(a); flint_printf("\n\n"); + flint_printf("b = "); arf_print(b); flint_printf("\n\n"); + flint_printf("c = "); arf_print(c); flint_printf("\n\n"); + flint_printf("x = "); arb_print(x); flint_printf("\n\n"); + flint_printf("y = "); arb_print(y); flint_printf("\n\n"); + flint_printf("z = "); arb_print(z); flint_printf("\n\n"); + abort(); + } + + /* aliasing */ + { + int alias; + + if (n_randint(state, 2)) + { + arb_min(x, x, y, prec); + alias = arb_equal(x, z); + } + else + { + arb_min(y, x, y, prec); + alias = arb_equal(y, z); + } + + if (!alias) + { + flint_printf("FAIL: aliasing\n\n"); + flint_printf("x = "); arb_print(x); flint_printf("\n\n"); + flint_printf("y = "); arb_print(y); flint_printf("\n\n"); + flint_printf("z = "); arb_print(z); flint_printf("\n\n"); + abort(); + } + } + + arf_clear(a); + arf_clear(b); + arf_clear(c); + + arb_clear(x); + arb_clear(y); + arb_clear(z); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/doc/source/arb.rst b/doc/source/arb.rst index 9ff687a3..4a184a29 100644 --- a/doc/source/arb.rst +++ b/doc/source/arb.rst @@ -619,6 +619,12 @@ Arithmetic Sets *y* to the sign function of *x*. The result is `[0 \pm 1]` if *x* contains both zero and nonzero numbers. +.. function:: void arb_min(arb_t z, const arb_t x, const arb_t y, slong prec) + +.. function:: void arb_max(arb_t z, const arb_t x, const arb_t y, slong prec) + + Sets *z* respectively to the minimum and the maximum of *x* and *y*. + .. function:: void arb_add(arb_t z, const arb_t x, const arb_t y, slong prec) .. function:: void arb_add_arf(arb_t z, const arb_t x, const arf_t y, slong prec)