rewrite arb_contains_fmpq to really make it exact, and add test code

arb/contains_fmpq.c

@@ -19,7 +19,7 @@
 =============================================================================*/
 /******************************************************************************
 
-    Copyright (C) 2012 Fredrik Johansson
+    Copyright (C) 2016 Fredrik Johansson
 
 ******************************************************************************/
 
@@ -28,17 +28,49 @@
 int
 arb_contains_fmpq(const arb_t x, const fmpq_t y)
 {
-    fmprb_t t;
+    if (fmpz_is_one(fmpq_denref(y)) || !arb_is_finite(x))
+    {
+        return arb_contains_fmpz(x, fmpq_numref(y));
+    }
+    else
+    {
+        arf_t t, xm, xr, ym;
+        arf_struct tmp[3];
         int result;
 
-    fmprb_init(t);
+        arf_init(t);
+        arf_init(xm);
+        arf_init(xr);
+        arf_init(ym);
 
-    arb_get_fmprb(t, x);
+        /* To compare x with p/q, compare qx with p. */
+        arf_mul_fmpz(xm, arb_midref(x), fmpq_denref(y), ARF_PREC_EXACT, ARF_RND_DOWN);
+        arf_set_mag(xr, arb_radref(x));
+        arf_mul_fmpz(xr, xr, fmpq_denref(y), ARF_PREC_EXACT, ARF_RND_DOWN);
+        arf_set_fmpz(ym, fmpq_numref(y));
 
-    result = fmprb_contains_fmpq(t, y);
+        /* y >= xm - xr  <=>  0 >= xm - xr - y */
+        arf_init_set_shallow(tmp + 0, xm);
+        arf_init_neg_shallow(tmp + 1,  xr);
+        arf_init_neg_shallow(tmp + 2, ym);
 
-    fmprb_clear(t);
+        arf_sum(t, tmp, 3, 30, ARF_RND_DOWN);
+        result = (arf_sgn(t) <= 0);
+
+        if (result)
+        {
+            /* y <= xm + xr  <=>  0 <= xm + xr - y */
+            arf_init_set_shallow(tmp + 1, xr);
+            arf_sum(t, tmp, 3, 30, ARF_RND_DOWN);
+            result = (arf_sgn(t) >= 0);
+        }
+
+        arf_clear(t);
+        arf_clear(xm);
+        arf_clear(xr);
+        arf_clear(ym);
 
         return result;
     }
+}
 

arb/test/t-contains_fmpq.c

@@ -0,0 +1,114 @@
+/*=============================================================================
+
+    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 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arb.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("contains_fmpq....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 10000; iter++)
+    {
+        arb_t x;
+        fmpq_t y;
+        fmpz_t a, b, t;
+        int c1, c2;
+        slong shift;
+
+        arb_init(x);
+        fmpq_init(y);
+        fmpz_init(a);
+        fmpz_init(b);
+        fmpz_init(t);
+
+        arb_randtest_special(x, state, 1 + n_randint(state, 500), 12);
+
+        if (n_randint(state, 2) && arb_is_finite(x))
+        {
+            arb_get_rand_fmpq(y, state, x, 1 + n_randint(state, 500));
+
+            fmpz_add_si(fmpq_numref(y), fmpq_numref(y), n_randint(state, 3) - 1);
+        }
+        else
+            fmpq_randtest(y, state, 1 + n_randint(state, 500));
+
+        c1 = arb_contains_fmpq(x, y);
+
+        if (arb_is_finite(x))
+        {
+            arb_get_interval_fmpz_2exp(a, b, t, x);
+            shift = fmpz_get_si(t);
+
+            fmpz_mul(a, a, fmpq_denref(y));
+            fmpz_mul(b, b, fmpq_denref(y));
+
+            fmpz_set(t, fmpq_numref(y));
+
+            if (shift >= 0)
+            {
+                fmpz_mul_2exp(a, a, shift);
+                fmpz_mul_2exp(b, b, shift);
+            }
+            else
+            {
+                fmpz_mul_2exp(t, t, -shift);
+            }
+
+            c2 = (fmpz_cmp(a, t) <= 0) && (fmpz_cmp(t, b) <= 0);
+        }
+        else
+        {
+            fmpz_tdiv_q(t, fmpq_numref(y), fmpq_denref(y));
+            c2 = arb_contains_fmpz(x, t);
+        }
+
+        if (c1 != c2)
+        {
+            flint_printf("FAIL:\n\n");
+            flint_printf("x = "); arb_print(x); flint_printf("\n\n");
+            flint_printf("y = "); fmpq_print(y); flint_printf("\n\n");
+            flint_printf("c1 = %d, c2 = %d\n\n", c1, c2);
+            abort();
+        }
+
+        arb_clear(x);
+        fmpq_clear(y);
+        fmpz_clear(a);
+        fmpz_clear(b);
+        fmpz_clear(t);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+