diff options
| -rw-r--r-- | Makefile.am | 41 | ||||
| -rw-r--r-- | NEWS | 1 | ||||
| -rw-r--r-- | mpfr.h | 2 | ||||
| -rw-r--r-- | mpfr.texi | 14 | ||||
| -rw-r--r-- | tests/Makefile.am | 44 | ||||
| -rw-r--r-- | tests/turandom.c | 187 | ||||
| -rw-r--r-- | urandom.c | 121 |
7 files changed, 368 insertions, 42 deletions
diff --git a/Makefile.am b/Makefile.am index 939978d80..f8002bf60 100644 --- a/Makefile.am +++ b/Makefile.am @@ -35,26 +35,27 @@ const_pi.c pow.c pow_si.c pow_ui.c print_raw.c print_rnd_mode.c \ reldiff.c round_prec.c set.c setmax.c setmin.c set_d.c set_dfl_prec.c \ set_exp.c set_rnd.c set_f.c set_prc_raw.c set_prec.c set_q.c set_si.c \ set_str.c set_str_raw.c set_ui.c set_z.c sqrt.c sqrt_ui.c sub.c sub1.c \ -sub_ui.c rint.c ui_div.c ui_sub.c urandomb.c get_z_exp.c swap.c \ -factorial.c cosh.c sinh.c tanh.c sinh_cosh.c acosh.c asinh.c atanh.c \ -atan.c cmp2.c exp_2.c asin.c const_euler.c cos.c sin.c tan.c fma.c \ -fms.c hypot.c log1p.c expm1.c log2.c log10.c ui_pow.c ui_pow_ui.c \ -minmax.c dim.c signbit.c copysign.c setsign.c gmp_op.c init2.c acos.c \ -sin_cos.c set_nan.c set_inf.c set_zero.c powerof2.c gamma.c set_ld.c \ -get_ld.c cbrt.c volatile.c fits_s.h fits_sshort.c fits_sint.c \ -fits_slong.c fits_u.h fits_ushort.c fits_uint.c fits_ulong.c \ -fits_uintmax.c fits_intmax.c get_si.c get_ui.c zeta.c cmp_d.c erf.c \ -inits.c inits2.c clears.c sgn.c check.c sub1sp.c version.c mpn_exp.c \ -mpfr-gmp.c mp_clz_tab.c sum.c add1sp.c free_cache.c si_op.c cmp_ld.c \ -set_ui_2exp.c set_si_2exp.c set_uj.c set_sj.c get_sj.c get_uj.c \ -get_z.c iszero.c cache.c sqr.c int_ceil_log2.c isqrt.c strtofr.c \ -pow_z.c logging.c mulders.c get_f.c round_p.c erfc.c atan2.c \ -subnormal.c const_catalan.c root.c gen_inverse.h sec.c csc.c cot.c \ -eint.c sech.c csch.c coth.c round_near_x.c constant.c abort_prec_max.c \ -stack_interface.c lngamma.c zeta_ui.c set_d64.c get_d64.c jn.c yn.c \ -rem1.c get_patches.c add_d.c sub_d.c d_sub.c mul_d.c div_d.c d_div.c \ -li2.c rec_sqrt.c min_prec.c buildopt.c digamma.c bernoulli.c \ -isregular.c set_flt.c get_flt.c scale2.c +sub_ui.c rint.c ui_div.c ui_sub.c urandom.c urandomb.c get_z_exp.c \ +swap.c factorial.c cosh.c sinh.c tanh.c sinh_cosh.c acosh.c asinh.c \ +atanh.c atan.c cmp2.c exp_2.c asin.c const_euler.c cos.c sin.c tan.c \ +fma.c fms.c hypot.c log1p.c expm1.c log2.c log10.c ui_pow.c \ +ui_pow_ui.c minmax.c dim.c signbit.c copysign.c setsign.c gmp_op.c \ +init2.c acos.c sin_cos.c set_nan.c set_inf.c set_zero.c powerof2.c \ +gamma.c set_ld.c get_ld.c cbrt.c volatile.c fits_s.h fits_sshort.c \ +fits_sint.c fits_slong.c fits_u.h fits_ushort.c fits_uint.c \ +fits_ulong.c fits_uintmax.c fits_intmax.c get_si.c get_ui.c zeta.c \ +cmp_d.c erf.c inits.c inits2.c clears.c sgn.c check.c sub1sp.c \ +version.c mpn_exp.c mpfr-gmp.c mp_clz_tab.c sum.c add1sp.c \ +free_cache.c si_op.c cmp_ld.c set_ui_2exp.c set_si_2exp.c set_uj.c \ +set_sj.c get_sj.c get_uj.c get_z.c iszero.c cache.c sqr.c \ +int_ceil_log2.c isqrt.c strtofr.c pow_z.c logging.c mulders.c get_f.c \ +round_p.c erfc.c atan2.c subnormal.c const_catalan.c root.c \ +gen_inverse.h sec.c csc.c cot.c eint.c sech.c csch.c coth.c \ +round_near_x.c constant.c abort_prec_max.c stack_interface.c lngamma.c \ +zeta_ui.c set_d64.c get_d64.c jn.c yn.c rem1.c get_patches.c add_d.c \ +sub_d.c d_sub.c mul_d.c div_d.c d_div.c li2.c rec_sqrt.c min_prec.c \ +buildopt.c digamma.c bernoulli.c isregular.c set_flt.c get_flt.c \ +scale2.c libmpfr_la_LIBADD = @LIBOBJS@ @@ -34,6 +34,7 @@ _ The rounding mode type mp_rnd_t is now mpfr_rnd_t (mp_rnd_t kept for - New function mpfr_set_zero. - New function mpfr_digamma. - New functions mpfr_set_flt and mpfr_get_flt to convert from/to the float type +- New function mpfr_urandom. - Speed improvement for large operands in the trigonometric functions (mpfr_sin, mpfr_cos, mpfr_tan, mpfr_sin_cos): speedup of about 2.5 for 10^5 digits, of about 5 for 10^6 digits. @@ -367,6 +367,8 @@ __MPFR_DECLSPEC int mpfr_get_z _MPFR_PROTO ((mpz_ptr z, mpfr_srcptr f, __MPFR_DECLSPEC void mpfr_free_str _MPFR_PROTO ((char *)); +__MPFR_DECLSPEC int mpfr_urandom _MPFR_PROTO ((mpfr_ptr, gmp_randstate_t, + mpfr_rnd_t)); __MPFR_DECLSPEC int mpfr_urandomb _MPFR_PROTO ((mpfr_ptr, gmp_randstate_t)); __MPFR_DECLSPEC void mpfr_nextabove _MPFR_PROTO ((mpfr_ptr)); @@ -2475,6 +2475,7 @@ Generate a uniformly distributed random float in the interval float with a random non-normalized significand and exponent 0, which is then normalized (thus if @var{e} denotes the exponent after normalization, then the least @math{-@var{e}} significant bits of the significand are always 0). + Return 0, unless the exponent is not in the current exponent range, in which case @var{rop} is set to NaN and a non-zero value is returned (this should never happen in practice, except in very specific cases). The @@ -2482,6 +2483,19 @@ second argument is a @code{gmp_randstate_t} structure which should be created using the GMP @code{gmp_randinit} function, see the GMP manual. @end deftypefun +@deftypefun int mpfr_urandom (mpfr_t @var{rop}, gmp_randstate_t @var{state}, mpfr_rnd_t @var{rnd}) +Generate a uniformly distributed random float. +The floating-point number @var{rop} can be seen as if a random real number is +generated according to the continuous uniform distribution on the interval +[0, 1] and then rounded in the direction @var{rnd}. + +Return 0, unless the exponent is not in the current exponent range, in +which case @var{rop} is set to NaN and a non-zero value is returned (this +should never happen in practice, except in very specific cases). +The second argument is a @code{gmp_randstate_t} structure which should be +created using the GMP @code{gmp_randinit} function, see the GMP manual. +@end deftypefun + @deftypefun mp_exp_t mpfr_get_exp (mpfr_t @var{x}) Get the exponent of @var{x}, assuming that @var{x} is a non-zero ordinary number and the significand is chosen in [1/2,1). The behavior for NaN, diff --git a/tests/Makefile.am b/tests/Makefile.am index 0e679490e..33251e098 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -11,28 +11,28 @@ AUTOMAKE_OPTIONS = 1.6 gnu $(top_builddir)/ansi2knr -check_PROGRAMS = tversion tinternals tinits tisqrt tsgn \ - tcheck tisnan texceptions tset_exp tset mpf_compat \ - mpfr_compat reuse tabs tacos tacosh tadd tadd1sp tadd_d \ - tadd_ui tagm tasin tasinh tatan tatanh taway tbuildopt \ - tcan_round tcbrt tcmp tcmp2 tcmp_d tcmp_ld tcmp_ui \ - tcmpabs tcomparisons tconst_catalan tconst_euler \ - tconst_log2 tconst_pi tcopysign tcos tcosh tcot tcoth \ - tcsc tcsch td_div td_sub tdigamma tdim tdiv tdiv_d tdiv_ui \ - teint teq terf texp texp10 texp2 texpm1 tfactorial tfits \ - tfma tfmod tfms tfprintf tfrac tgamma tget_flt tget_d \ - tget_d_2exp tget_f tget_ld_2exp tget_set_d64 tget_sj \ - tget_str tget_z tgmpop thyperbolic thypot tinp_str tj0 tj1 \ - tjn tl2b tlgamma tli2 tlngamma tlog tlog10 tlog1p tlog2 \ - tmin_prec tminmax tmodf tmul tmul_2exp tmul_d tmul_ui \ - tnext tout_str toutimpl tpow tpow3 tpow_all tpow_z \ - tprintf trandom trec_sqrt tremquo trint troot \ - tround_prec tsec tsech tset_d tset_f tset_ld tset_q \ - tset_si tset_sj tset_str tset_z tsi_op tsin tsin_cos \ - tsinh tsinh_cosh tsprintf tsqr tsqrt tsqrt_ui tstckintc \ - tstrtofr tsub tsub1sp tsub_d tsub_ui tsubnormal tsum \ - tswap ttan ttanh ttrunc tui_div tui_pow tui_sub ty0 ty1 \ - tyn tzeta tzeta_ui +check_PROGRAMS = tversion tinternals tinits tisqrt tsgn tcheck \ + tisnan texceptions tset_exp tset mpf_compat mpfr_compat \ + reuse tabs tacos tacosh tadd tadd1sp tadd_d tadd_ui tagm \ + tasin tasinh tatan tatanh taway tbuildopt tcan_round \ + tcbrt tcmp tcmp2 tcmp_d tcmp_ld tcmp_ui tcmpabs \ + tcomparisons tconst_catalan tconst_euler tconst_log2 \ + tconst_pi tcopysign tcos tcosh tcot tcoth tcsc tcsch \ + td_div td_sub tdigamma tdim tdiv tdiv_d tdiv_ui teint teq \ + terf texp texp10 texp2 texpm1 tfactorial tfits tfma tfmod \ + tfms tfprintf tfrac tgamma tget_flt tget_d tget_d_2exp \ + tget_f tget_ld_2exp tget_set_d64 tget_sj tget_str tget_z \ + tgmpop thyperbolic thypot tinp_str tj0 tj1 tjn tl2b \ + tlgamma tli2 tlngamma tlog tlog10 tlog1p tlog2 tmin_prec \ + tminmax tmodf tmul tmul_2exp tmul_d tmul_ui tnext \ + tout_str toutimpl tpow tpow3 tpow_all tpow_z tprintf \ + trandom trec_sqrt tremquo trint troot tround_prec tsec \ + tsech tset_d tset_f tset_ld tset_q tset_si tset_sj \ + tset_str tset_z tsi_op tsin tsin_cos tsinh tsinh_cosh \ + tsprintf tsqr tsqrt tsqrt_ui tstckintc tstrtofr tsub \ + tsub1sp tsub_d tsub_ui tsubnormal tsum tswap turandom \ + ttan ttanh ttrunc tui_div tui_pow tui_sub ty0 ty1 tyn \ + tzeta tzeta_ui EXTRA_DIST = tgeneric.c tgeneric_ui.c mpf_compat.h inp_str.data tmul.dat diff --git a/tests/turandom.c b/tests/turandom.c new file mode 100644 index 000000000..7c3af4595 --- /dev/null +++ b/tests/turandom.c @@ -0,0 +1,187 @@ +/* Test file for mpfr_urandom + +Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2009 Free Software Foundation, Inc. +Contributed by the Arenaire and Cacao projects, INRIA. + +This file is part of the GNU MPFR Library. + +The GNU MPFR Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +The GNU MPFR Library 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 Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see +http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., +51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ + +#include <stdio.h> +#include <stdlib.h> + +#include "mpfr-test.h" + +static void +test_urandom (long nbtests, mp_prec_t prec, mpfr_rnd_t rnd, long bit_index, + int verbose) +{ + mpfr_t x; + int *tab, size_tab, k, sh, xn; + double d, av = 0, var = 0, chi2 = 0, th; + mp_exp_t emin; + mp_size_t limb_index = 0; + mp_limb_t limb_mask = 0; + long count = 0; + int i; + + size_tab = (nbtests >= 1000 ? nbtests / 50 : 20); + tab = (int *) calloc (size_tab, sizeof(int)); + if (tab == NULL) + { + fprintf (stderr, "trandom: can't allocate memory in test_urandom\n"); + exit (1); + } + + mpfr_init2 (x, prec); + xn = 1 + (prec - 1) / mp_bits_per_limb; + sh = xn * mp_bits_per_limb - prec; + if (bit_index >= 0 && bit_index < prec) + { + /* compute the limb index and limb mask to fetch the bit #bit_index */ + limb_index = (prec - bit_index) / mp_bits_per_limb; + i = 1 + bit_index - (bit_index / mp_bits_per_limb) * mp_bits_per_limb; + limb_mask = MPFR_LIMB_ONE << (mp_bits_per_limb - i); + } + + for (k = 0; k < nbtests; k++) + { + mpfr_urandom (x, RANDS, rnd); + /* check that lower bits are zero */ + if (MPFR_MANT(x)[0] & MPFR_LIMB_MASK(sh) && !MPFR_IS_ZERO (x)) + { + printf ("Error: mpfr_urandom() returns invalid numbers:\n"); + mpfr_print_binary (x); puts (""); + exit (1); + } + d = mpfr_get_d1 (x); av += d; var += d*d; + i = (int)(size_tab * d); + if (d == 1.0) i --; + tab[i]++; + + if (limb_mask && (MPFR_MANT (x)[limb_index] & limb_mask)) + count ++; + } + + /* coverage test */ + emin = mpfr_get_emin (); + set_emin (1); + /* the generated number in [0,1] is not in the exponent range, except: + - if it is zero and rnd is toward zero or to nearest + - if it is one and rnd is not toward zero and not to nearest */ + k = mpfr_urandom (x, RANDS, rnd); + if ((k == 0 || mpfr_nan_p (x) == 0) + && (((rnd != MPFR_RNDZ && rnd != MPFR_RNDN) && MPFR_IS_ZERO (x)) + || ((rnd==MPFR_RNDZ || rnd==MPFR_RNDN) && mpfr_cmp_ui (x, 1)==0))) + { + printf ("Error in mpfr_urandom, expected NaN, got "); + mpfr_dump (x); + exit (1); + } + set_emin (emin); + + mpfr_clear (x); + if (!verbose) + { + free(tab); + return; + } + + av /= nbtests; + var = (var / nbtests) - av * av; + + th = (double)nbtests / size_tab; + printf("Average = %.5f\nVariance = %.5f\n", av, var); + printf("Repartition for urandom with rounding mode %s. " + "Each integer should be close to %d.\n", + mpfr_print_rnd_mode (rnd), (int)th); + + for (k = 0; k < size_tab; k++) + { + chi2 += (tab[k] - th) * (tab[k] - th) / th; + printf("%d ", tab[k]); + if (((k+1) & 7) == 0) + printf("\n"); + } + + printf("\nChi2 statistics value (with %d degrees of freedom) : %.5f\n", + size_tab - 1, chi2); + + if (limb_mask) + printf ("Bit #%ld is set %ld/%ld = %.1f %% of time\n", + bit_index, count, nbtests, count * 100.0 / nbtests); + + puts (""); + + free(tab); + return; +} + + +int +main (int argc, char *argv[]) +{ + long nbtests; + mp_prec_t prec; + int verbose = 0; + mpfr_rnd_t rnd; + long bit_index; + + tests_start_mpfr (); + + if (argc > 1) + verbose = 1; + + nbtests = 10000; + if (argc > 1) + { + long a = atol(argv[1]); + if (a != 0) + nbtests = a; + } + + if (argc <= 2) + prec = 1000; + else + prec = atol(argv[2]); + + if (argc <= 3) + bit_index = -1; + else + { + bit_index = atol(argv[3]); + if (bit_index >= prec) + { + printf ("Warning. Cannot compute the bit frequency: the given bit " + "index (= %ld) is not less than the precision (= %ld).\n", + bit_index, prec); + bit_index = -1; + } + } + + RND_LOOP(rnd) + { + test_urandom (nbtests, prec, rnd, bit_index, verbose); + + if (argc == 1) /* check also small precision */ + { + test_urandom (nbtests, 2, rnd, -1, 0); + } + } + + tests_end_mpfr (); + return 0; +} diff --git a/urandom.c b/urandom.c new file mode 100644 index 000000000..376857fbe --- /dev/null +++ b/urandom.c @@ -0,0 +1,121 @@ +/* mpfr_urandom (rop, state, rnd_mode) -- Generate a uniform pseudorandom + real number between 0 and 1 (exclusive) and round it to the precision of rop + according to the given rounding mode. + +Copyright 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2009 Free Software Foundation, Inc. +Contributed by the Arenaire and Cacao projects, INRIA. + +This file is part of the GNU MPFR Library. + +The GNU MPFR Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +The GNU MPFR Library 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 Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see +http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., +51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ + + +#define MPFR_NEED_LONGLONG_H +#include "mpfr-impl.h" + +int +mpfr_urandom (mpfr_ptr rop, gmp_randstate_t rstate, mpfr_rnd_t rnd_mode) +{ + mp_ptr rp; + mp_prec_t nbits; + mp_size_t nlimbs; + mp_size_t n; + mp_size_t k; /* number of high zero limbs */ + mp_exp_t exp; + int rndbit = 0; + int cnt; + + rp = MPFR_MANT (rop); + nbits = MPFR_PREC (rop); + nlimbs = MPFR_LIMB_SIZE (rop); + MPFR_SET_POS (rop); + exp = 0; + + + /* Rounding bit */ + if (rnd_mode == MPFR_RNDN) + { + mpfr_rand_raw (rp, rstate, GMP_NUMB_BITS); + rndbit = rp[0] & 1; + } + + + /* Exponent */ + k = nlimbs; + while (k == nlimbs) + { + mpfr_rand_raw (rp, rstate, nlimbs * GMP_NUMB_BITS); + + /* Count the null significant limbs k and remaining limbs n */ + k = 0; + n = nlimbs; + while (k != nlimbs && rp[n - 1] == 0) + { + k ++; + n --; + } + + if (exp < MPFR_EXP_MIN + k * GMP_NUMB_BITS) + goto tiny; + + exp -= k * GMP_NUMB_BITS; + } + count_leading_zeros (cnt, rp[n - 1]); + if (mpfr_set_exp (rop, exp - cnt)) + goto outofrange; + + + /* Significand */ + mpfr_rand_raw (rp, rstate, nlimbs * GMP_NUMB_BITS); + + /* Set the msb to 1 since it was fixed by the exponent choice */ + rp[nlimbs - 1] |= MPFR_LIMB_HIGHBIT; + + /* If nbits isn't a multiple of GMP_NUMB_BITS, mask the low bits */ + n = nlimbs * GMP_NUMB_BITS - nbits; + if (MPFR_LIKELY (n != 0)) + rp[0] &= ~MPFR_LIMB_MASK (n); + + + /* Rounding */ + if (rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDA + || (rnd_mode == MPFR_RNDN && rndbit)) + mpfr_nextabove (rop); + + normalexit: + return 0; + + tiny: + /* To get here, we have been drawing more than 2^31 zeros in a raw + (very unlucky). */ + MPFR_SET_ZERO (rop); + if (rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDA + || (rnd_mode == MPFR_RNDN && rndbit)) + mpfr_nextabove (rop); + + return 0; + + outofrange: + /* If the exponent is not in the current exponent range, we choose + to return a NaN as this is probably a user error. + Indeed this can happen only if the exponent range has been + reduced to a very small interval and/or the precision is huge + (very unlikely). */ + MPFR_SET_NAN (rop); + __gmpfr_flags |= MPFR_FLAGS_NAN; /* Can't use MPFR_RET_NAN */ + + return 1; +} |