diff --git a/fmprb.h b/fmprb.h index e39002ca..4f4f14d2 100644 --- a/fmprb.h +++ b/fmprb.h @@ -320,5 +320,7 @@ fmprb_rel_accuracy_bits(const fmprb_t x) void fmprb_randtest(fmprb_t x, flint_rand_t state, long prec, long mag_bits); +void fmprb_get_rand_fmpq(fmpq_t q, flint_rand_t state, const fmprb_t x, long bits); + #endif diff --git a/fmprb/get_rand_fmpq.c b/fmprb/get_rand_fmpq.c new file mode 100644 index 00000000..6e165b17 --- /dev/null +++ b/fmprb/get_rand_fmpq.c @@ -0,0 +1,125 @@ +/*============================================================================= + + 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) 2012 Fredrik Johansson + +******************************************************************************/ + +#include "fmprb.h" + +/* +Returns a random rational in x. The binary denominator of x is multiplied +by den_mult to select a base denominator to choose random integers +between. + +The outcome is undefined if the midpoint or radius of x is non-finite, +or if the exponent of the midpoint or radius is so large or small +that representing the endpoints as exact rational numbers would +cause overflows. +*/ + +void +_fmprb_get_rand_fmpq(fmpz_t num, fmpz_t den, flint_rand_t state, + const fmpz_t den_mult, const fmprb_t x) +{ + fmpz_t a, b, exp; + + fmpz_init(a); + fmpz_init(b); + fmpz_init(exp); + + fmprb_get_interval_fmpz_2exp(a, b, exp, x); + + if (COEFF_IS_MPZ(*exp)) + { + printf("exception: fmprb_get_rand_fmpq: too large exponent\n"); + abort(); + } + + if (*exp >= 0) + { + fmpz_mul_2exp(a, a, *exp); + fmpz_mul_2exp(b, b, *exp); + } + + /* generate random integer in [a*den, b*den] */ + fmpz_mul(a, a, den_mult); + fmpz_mul(b, b, den_mult); + fmpz_add_ui(b, b, 1UL); + fmpz_sub(b, b, a); + + /* return one endpoint with high probability (used for stress + testing rounding) */ + if (n_randint(state, 6) == 0) + { + if (n_randint(state, 2)) + fmpz_zero(num); + else + fmpz_sub_ui(num, b, 1UL); + } + else + { + fmpz_randtest_mod(num, state, b); + } + + fmpz_add(num, num, a); + + fmpz_set(den, den_mult); + + if (*exp < 0) + fmpz_mul_2exp(den, den, -(*exp)); + + fmpz_clear(a); + fmpz_clear(b); + fmpz_clear(exp); +} + + +/* +Chooses a random rational number from the interval represented by x. +A denominator is chosen by multiplying the binary denominator of x +by a random integer up to size bits. + +The outcome is undefined if the midpoint or radius of x is non-finite, +or if the exponent of the midpoint or radius is so large or small +that representing the endpoints as exact rational numbers would +cause overflows. +*/ + +void +fmprb_get_rand_fmpq(fmpq_t q, flint_rand_t state, const fmprb_t x, long bits) +{ + /* there is only one rational */ + if (fmprb_is_exact(x)) + { + fmpr_get_fmpq(q, fmprb_midref(x)); + return; + } + + /* pick a denominator */ + fmpz_randbits(fmpq_denref(q), state, n_randint(state, bits + 1)); + fmpz_abs(fmpq_denref(q), fmpq_denref(q)); + if (fmpz_is_zero(fmpq_denref(q))) + fmpz_one(fmpq_denref(q)); + + _fmprb_get_rand_fmpq(fmpq_numref(q), fmpq_denref(q), state, fmpq_denref(q), x); + fmpq_canonicalise(q); +} diff --git a/fmprb/test/t-get_rand_fmpq.c b/fmprb/test/t-get_rand_fmpq.c new file mode 100644 index 00000000..e0de0604 --- /dev/null +++ b/fmprb/test/t-get_rand_fmpq.c @@ -0,0 +1,64 @@ +/*============================================================================= + + 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) 2012 Fredrik Johansson + +******************************************************************************/ + +#include "fmprb.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("get_rand_fmpq...."); + fflush(stdout); + flint_randinit(state); + + for (iter = 0; iter < 100000; iter++) + { + fmprb_t x; + fmpq_t q; + + fmprb_init(x); + fmpq_init(q); + + fmprb_randtest(x, state, 1 + n_randint(state, 200), 10); + fmprb_get_rand_fmpq(q, state, x, 1 + n_randint(state, 200)); + + if (!fmprb_contains_fmpq(x, q) || !fmpq_is_canonical(q)) + { + printf("FAIL:\n\n"); + printf("x = "); fmprb_print(x); printf("\n\n"); + printf("q = "); fmpq_print(q); printf("\n\n"); + abort(); + } + + fmprb_clear(x); + fmpq_clear(q); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +}