add arb_contains_interior, acb_contains_interior

acb.h

@@ -215,6 +215,13 @@ acb_contains(const acb_t x, const acb_t y)
             arb_contains(acb_imagref(x), acb_imagref(y));
 }
 
+ACB_INLINE int
+acb_contains_interior(const acb_t x, const acb_t y)
+{
+    return arb_contains_interior(acb_realref(x), acb_realref(y)) &&
+            arb_contains_interior(acb_imagref(x), acb_imagref(y));
+}
+
 ACB_INLINE void
 acb_set_ui(acb_t z, ulong c)
 {

arb.h

@@ -353,6 +353,8 @@ int arb_overlaps(const arb_t x, const arb_t y);
 
 int arb_contains(const arb_t x, const arb_t y);
 
+int arb_contains_interior(const arb_t x, const arb_t y);
+
 int arb_contains_int(const arb_t x);
 
 void arb_get_interval_fmpz_2exp(fmpz_t a, fmpz_t b, fmpz_t exp, const arb_t x);

arb/contains_interior.c

@@ -0,0 +1,71 @@
+/*
+    Copyright (C) 2012, 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 <http://www.gnu.org/licenses/>.
+*/
+
+#include "arb.h"
+
+int
+arb_contains_interior(const arb_t x, const arb_t y)
+{
+    arf_t t;
+    arf_t u;
+    arf_t xr, yr;
+    arf_struct tmp[4];
+    int left_ok, right_ok;
+
+    if (arf_is_nan(arb_midref(x)) || arb_is_exact(x) || !arb_is_finite(y))
+        return 0;
+
+    arf_init(t);
+    arf_init(u);
+
+    arf_init_set_mag_shallow(xr, arb_radref(x));
+    arf_init_set_mag_shallow(yr, arb_radref(y));
+
+    /* fast check */
+    arf_sub(t, arb_midref(x), xr, MAG_BITS, ARF_RND_CEIL);
+    arf_sub(u, arb_midref(y), yr, MAG_BITS, ARF_RND_FLOOR);
+    left_ok = arf_cmp(t, u) < 0;
+
+    /* exact check */
+    if (!left_ok)
+    {
+        arf_init_set_shallow(tmp + 0, arb_midref(x));
+        arf_init_neg_mag_shallow(tmp + 1, arb_radref(x));
+        arf_init_neg_shallow(tmp + 2, arb_midref(y));
+        arf_init_set_mag_shallow(tmp + 3, arb_radref(y));
+
+        arf_sum(t, tmp, 4, MAG_BITS, ARF_RND_DOWN);
+        left_ok = arf_sgn(t) < 0;
+    }
+
+    /* fast check */
+    arf_add(t, arb_midref(x), xr, MAG_BITS, ARF_RND_FLOOR);
+    arf_add(u, arb_midref(y), yr, MAG_BITS, ARF_RND_CEIL);
+    right_ok = (arf_cmp(t, u) > 0);
+
+    /* exact check */
+    if (!right_ok)
+    {
+        arf_init_set_shallow(tmp + 0, arb_midref(x));
+        arf_init_set_mag_shallow(tmp + 1, arb_radref(x));
+        arf_init_neg_shallow(tmp + 2, arb_midref(y));
+        arf_init_neg_mag_shallow(tmp + 3, arb_radref(y));
+
+        arf_sum(t, tmp, 4, MAG_BITS, ARF_RND_DOWN);
+        right_ok = arf_sgn(t) > 0;
+    }
+
+    arf_clear(t);
+    arf_clear(u);
+
+    return left_ok && right_ok;
+}
+

arb/test/t-contains_interior.c

@@ -0,0 +1,88 @@
+/*
+    Copyright (C) 2012 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 <http://www.gnu.org/licenses/>.
+*/
+
+#include "arb.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("contains_interior....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++)
+    {
+        arb_t a, b;
+        fmpq_t am, ar, bm, br, t, u;
+        int c1, c2;
+
+        arb_init(a);
+        arb_init(b);
+
+        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));
+
+        c1 = arb_contains_interior(a, b);
+
+        fmpq_sub(t, am, ar);
+        fmpq_sub(u, bm, br);
+        c2 = fmpq_cmp(t, u) < 0;
+
+        fmpq_add(t, am, ar);
+        fmpq_add(u, bm, br);
+        c2 = c2 && (fmpq_cmp(t, u) > 0);
+
+        if (c1 != c2)
+        {
+            flint_printf("FAIL:\n\n");
+            flint_printf("a = "); arb_print(a); flint_printf("\n\n");
+            flint_printf("b = "); arb_print(b); 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\n\n", c1, c2);
+            flint_abort();
+        }
+
+        arb_clear(a);
+        arb_clear(b);
+
+        fmpq_clear(am);
+        fmpq_clear(ar);
+        fmpq_clear(bm);
+        fmpq_clear(br);
+        fmpq_clear(t);
+        fmpq_clear(u);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

doc/source/acb.rst

@@ -308,6 +308,16 @@ Precision and comparisons
     Returns nonzero iff the complex interval represented by *x* contains
     an integer.
 
+.. function:: int acb_contains_interior(const acb_t x, const acb_t y)
+
+    Tests if *y* is contained in the interior of *x*.
+    This predicate always evaluates to false if *x* and *y* are both
+    real-valued, since an imaginary part of 0 is not considered contained in
+    the interior of the point interval 0. More generally, the same
+    problem occurs for intervals with an exact real or imaginary part.
+    Such intervals must be handled specially by the user where a different
+    interpretation is intended.
+
 .. function:: slong acb_rel_error_bits(const acb_t x)
 
     Returns the effective relative error of *x* measured in bits.

doc/source/arb.rst

@@ -678,6 +678,11 @@ Comparisons
     by *x* satisfying, respectively, `p = 0`, `p < 0`, `p \le 0`, `p > 0`, `p \ge 0`.
     If *x* contains NaN, returns nonzero.
 
+.. function:: int arb_contains_interior(const arb_t x, const arb_t y)
+
+    Tests if *y* is contained in the interior of *x*; that is, contained
+    in *x* and not touching either endpoint.
+
 .. function:: int arb_eq(const arb_t x, const arb_t y)
 
 .. function:: int arb_ne(const arb_t x, const arb_t y)