diff --git a/arb/contains_fmpq.c b/arb/contains_fmpq.c index 955d8b1f..4fff4907 100644 --- a/arb/contains_fmpq.c +++ b/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; - int result; + 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); - return result; + 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; + } } diff --git a/arb/test/t-contains_fmpq.c b/arb/test/t-contains_fmpq.c new file mode 100644 index 00000000..45f94334 --- /dev/null +++ b/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; +} +