code for generating random test rationals

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
 

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);
+}

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;
+}