diff --git a/arb.h b/arb.h index e40f3bda..14c113cc 100644 --- a/arb.h +++ b/arb.h @@ -532,6 +532,7 @@ void arb_get_interval_arf(arf_t a, arf_t b, const arb_t x, slong prec); void arb_get_interval_mpfr(mpfr_t a, mpfr_t b, const arb_t x); 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); int arb_can_round_arf(const arb_t x, slong prec, arf_rnd_t rnd); diff --git a/arb/intersection.c b/arb/intersection.c new file mode 100644 index 00000000..63e9b6dc --- /dev/null +++ b/arb/intersection.c @@ -0,0 +1,73 @@ +/*============================================================================= + + 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" + +int +arb_intersection(arb_t z, const arb_t x, const arb_t y, slong prec) +{ + arf_t left, right, t, xr, yr; + int result; + + if (arf_is_nan(arb_midref(x)) || arf_is_nan(arb_midref(y))) + { + arb_indeterminate(z); + return 1; + } + + if (mag_is_inf(arb_radref(x)) && mag_is_inf(arb_radref(y))) + { + arb_zero_pm_inf(z); + return 1; + } + + result = arb_overlaps(x, y); + + if (result) + { + 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_min(right, right, t); + + arb_set_interval_arf(z, left, right, prec); + + arf_clear(left); + arf_clear(right); + arf_clear(t); + } + + return result; +} diff --git a/arb/test/t-intersection.c b/arb/test/t-intersection.c new file mode 100644 index 00000000..3d717b8a --- /dev/null +++ b/arb/test/t-intersection.c @@ -0,0 +1,213 @@ +/*============================================================================= + + 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" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("intersection...."); + fflush(stdout); + flint_randinit(state); + + /* check a containment requirement */ + for (iter = 0; iter < 100000; iter++) + { + arb_t x, y, z, w; + arb_t xy, yz; + slong pa, pb, pc; + int alias; + + arb_init(x); + arb_init(y); + arb_init(z); + arb_init(w); + arb_init(xy); + arb_init(yz); + + arb_randtest_special(x, state, 200, 10); + arb_randtest_special(y, state, 200, 10); + arb_randtest_special(z, state, 200, 10); + arb_randtest_special(w, state, 200, 10); + arb_randtest_special(xy, state, 200, 10); + arb_randtest_special(yz, state, 200, 10); + + pa = 2 + n_randint(state, 200); + pb = 2 + n_randint(state, 200); + pc = 2 + n_randint(state, 200); + + arb_union(xy, x, y, pa); + arb_union(yz, y, z, pb); + arb_intersection(w, xy, yz, pc); + + if (!arb_contains(w, y)) + { + flint_printf("FAIL (containment):\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"); + flint_printf("w = "); arb_print(w); flint_printf("\n\n"); + abort(); + } + + if (n_randint(state, 2)) + { + arb_intersection(xy, xy, yz, pc); + alias = arb_equal(xy, w); + } + else + { + arb_intersection(yz, xy, yz, pc); + alias = arb_equal(yz, w); + } + + 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"); + flint_printf("w = "); arb_print(w); flint_printf("\n\n"); + abort(); + } + + arb_clear(x); + arb_clear(y); + arb_clear(z); + arb_clear(w); + arb_clear(xy); + arb_clear(yz); + } + + /* require that the return value is the same as for arb_overlaps */ + for (iter = 0; iter < 10000; iter++) + { + arb_t a, b, y; + fmpq_t am, ar, bm, br, t, u; + int c1, c2, c3; + slong prec; + + prec = 2 + n_randint(state, 200); + + arb_init(a); + arb_init(b); + arb_init(y); + + fmpq_init(am); + fmpq_init(ar); + fmpq_init(bm); + fmpq_init(br); + fmpq_init(t); + fmpq_init(u); + + arb_randtest(a, state, 1 + n_randint(state, 500), 14); + arb_randtest(b, state, 1 + n_randint(state, 500), 14); + + arf_get_fmpq(am, arb_midref(a)); + mag_get_fmpq(ar, arb_radref(a)); + arf_get_fmpq(bm, arb_midref(b)); + mag_get_fmpq(br, arb_radref(b)); + + fmpq_sub(t, am, bm); + fmpz_abs(fmpq_numref(t), fmpq_numref(t)); + fmpq_add(u, ar, br); + + c1 = arb_overlaps(a, b); + + c2 = (fmpq_cmp(t, u) <= 0); + + c3 = arb_intersection(y, a, b, prec); + + if (c1 != c2 || c1 != c3) + { + flint_printf("FAIL (compatibility with arb_overlaps):\n\n"); + flint_printf("a = "); arb_print(a); flint_printf("\n\n"); + flint_printf("b = "); arb_print(b); flint_printf("\n\n"); + flint_printf("y = "); arb_print(y); flint_printf("\n\n"); + flint_printf("am = "); fmpq_print(am); flint_printf("\n\n"); + flint_printf("ar = "); fmpq_print(ar); flint_printf("\n\n"); + flint_printf("bm = "); fmpq_print(bm); flint_printf("\n\n"); + flint_printf("br = "); fmpq_print(br); flint_printf("\n\n"); + flint_printf("t = "); fmpq_print(t); flint_printf("\n\n"); + flint_printf("u = "); fmpq_print(u); flint_printf("\n\n"); + flint_printf("c1 = %d, c2 = %d, c3 = %d\n\n", c1, c2, c3); + abort(); + } + + arb_clear(a); + arb_clear(b); + arb_clear(y); + + fmpq_clear(am); + fmpq_clear(ar); + fmpq_clear(bm); + fmpq_clear(br); + fmpq_clear(t); + fmpq_clear(u); + } + + /* check a simple hardcoded example */ + { + slong prec; + arb_t xy, yz, y, v, w; + + prec = 32; + + arb_init(xy); + arb_init(yz); + arb_init(y); + arb_init(v); + arb_init(w); + + arb_set_str(xy, "1 +/- 1", prec); + arb_set_str(yz, "2 +/- 1", prec); + arb_set_str(y, "1.5 +/- 0.6", prec); + arb_set_str(v, "1.5 +/- 0.4", prec); + + arb_intersection(w, xy, yz, prec); + + if (!arb_contains(y, w) || !arb_contains(w, v)) + { + flint_printf("FAIL (hardcoded example)\n\n"); + flint_printf("xy = "); arb_print(xy); flint_printf("\n\n"); + flint_printf("yx = "); arb_print(yz); flint_printf("\n\n"); + flint_printf("y = "); arb_print(y); flint_printf("\n\n"); + flint_printf("w = "); arb_print(w); flint_printf("\n\n"); + abort(); + } + + arb_clear(xy); + arb_clear(yz); + arb_clear(y); + arb_clear(w); + } + + 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 30ac5d23..da8be1e4 100644 --- a/doc/source/arb.rst +++ b/doc/source/arb.rst @@ -313,6 +313,14 @@ Radius and interval operations Sets *z* to a ball containing both *x* and *y*. +.. function:: int arb_intersection(arb_t z, const arb_t x, const arb_t y, slong prec) + + If *x* and *y* overlap according to :func:`arb_overlaps`, + then *z* is set to a ball containing the intersection of *x* and *y* + and a nonzero value is returned. + Otherwise zero is returned and the value of *z* is undefined. + If *x* or *y* contains NaN, the result is NaN. + .. function:: void arb_get_abs_ubound_arf(arf_t u, const arb_t x, slong prec) Sets *u* to the upper bound for the absolute value of *x*,