diff options
author | Akira Kakuto <kakuto@fuk.kindai.ac.jp> | 2017-12-26 11:50:03 +0000 |
---|---|---|
committer | Akira Kakuto <kakuto@fuk.kindai.ac.jp> | 2017-12-26 11:50:03 +0000 |
commit | 6ddc203e5368cf1d4108e64c18b0bf6c7d53d176 (patch) | |
tree | dda1d9b65325699c60418cfcd02fa5d800cde1f8 /Build/source/libs/mpfr/mpfr-src/src/nrandom.c | |
parent | bbfe0bab9260ac7a06a69f710b053dcc73cecad2 (diff) |
mpfr-4.0.0
git-svn-id: svn://tug.org/texlive/trunk@46140 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/nrandom.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/nrandom.c | 183 |
1 files changed, 183 insertions, 0 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/nrandom.c b/Build/source/libs/mpfr/mpfr-src/src/nrandom.c new file mode 100644 index 00000000000..ed5133bc1fd --- /dev/null +++ b/Build/source/libs/mpfr/mpfr-src/src/nrandom.c @@ -0,0 +1,183 @@ +/* mpfr_nrandom (rop, state, rnd_mode) -- Generate a normal deviate with mean 0 + and variance 1 and round it to the precision of rop according to the given + rounding mode. + +Copyright 2013-2017 Free Software Foundation, Inc. +Contributed by Charles Karney <charles@karney.com>, SRI International. + +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. */ + +/* + * Sampling from the normal distribution with zero mean and unit variance. + * This uses Algorithm N given in: + * Charles F. F. Karney, + * "Sampling exactly from the normal distribution", + * ACM Trans. Math. Software 42(1), 3:1-14 (Jan. 2016). + * https://dx.doi.org/10.1145/2710016 + * http://arxiv.org/abs/1303.6257 + * + * The implementation here closely follows the C++ one given in the paper + * above. However, here, C is simplified by using gmp_urandomm_ui; the initial + * rejection step in H just tests the leading bit of p; and the assignment of + * the sign to the deviate using gmp_urandomb_ui. Finally, the C++ + * implementation benefits from caching temporary random deviates across calls. + * This isn't possible in C without additional machinery which would complicate + * the interface. + * + * There are a few "weasel words" regarding the accuracy of this + * implementation. The algorithm produces exactly rounded normal deviates + * provided that gmp's random number engine delivers truly random bits. If it + * did, the algorithm would be perfect; however, this implementation would have + * problems, e.g., in that the integer part of the normal deviate is + * represented by an unsigned long, whereas in reality the integer part in + * unbounded. In this implementation, asserts catch overflow in the integer + * part and similar (very, very) unlikely events. In reality, of course, gmp's + * random number engine has a finite internal state (19937 bits in the case of + * the MT19937 method). This means that these unlikely events in fact won't + * occur. If the asserts are triggered, then this is an indication that the + * random number engine is defective. (Even if a hardware random number + * generator were used, the most likely explanation for the triggering of the + * asserts would be that the hardware generator was broken.) + */ + +#include "random_deviate.h" + +/* Algorithm H: true with probability exp(-1/2). */ +static int +H (gmp_randstate_t r, mpfr_random_deviate_t p, mpfr_random_deviate_t q) +{ + /* p and q are temporaries */ + mpfr_random_deviate_reset (p); + if (mpfr_random_deviate_tstbit (p, 1, r)) + return 1; + for (;;) + { + mpfr_random_deviate_reset (q); + if (!mpfr_random_deviate_less (q, p, r)) + return 0; + mpfr_random_deviate_reset (p); + if (!mpfr_random_deviate_less (p, q, r)) + return 1; + } +} + +/* Step N1: return n >= 0 with prob. exp(-n/2) * (1 - exp(-1/2)). */ +static unsigned long +G (gmp_randstate_t r, mpfr_random_deviate_t p, mpfr_random_deviate_t q) +{ + /* p and q are temporaries */ + unsigned long n = 0; + + while (H (r, p, q)) + { + ++n; + /* Catch n wrapping around to 0; for a 32-bit unsigned long, the + * probability of this is exp(-2^30)). */ + MPFR_ASSERTN (n != 0UL); + } + return n; +} + +/* Step N2: true with probability exp(-m*n/2). */ +static int +P (unsigned long m, unsigned long n, gmp_randstate_t r, + mpfr_random_deviate_t p, mpfr_random_deviate_t q) +{ + /* p and q are temporaries. m*n is passed as two separate parameters to deal + * with the case where m*n overflows an unsigned long. This may be called + * with m = 0 and n = (unsigned long)(-1) and, because m in handled in to the + * outer loop, this routine will correctly return 1. */ + while (m--) + { + unsigned long k = n; + while (k--) + { + if (!H (r, p, q)) + return 0; + } + } + return 1; +} + +/* Algorithm C: return (-1, 0, 1) with prob (1/m, 1/m, 1-2/m). */ +static int +C (unsigned long m, gmp_randstate_t r) +{ + unsigned long n = gmp_urandomm_ui (r, m); + return n == 0 ? -1 : (n == 1 ? 0 : 1); +} + +/* Algorithm B: true with prob exp(-x * (2*k + x) / (2*k + 2)). */ +static int +B (unsigned long k, mpfr_random_deviate_t x, gmp_randstate_t r, + mpfr_random_deviate_t p, mpfr_random_deviate_t q) +{ + /* p and q are temporaries */ + + unsigned long m = 2 * k + 2; + /* n tracks the parity of the loop; s == 1 on first trip through loop. */ + unsigned n = 0, s = 1; + int f; + + /* Check if 2 * k + 2 would overflow; for a 32-bit unsigned long, the + * probability of this is exp(-2^61)). */ + MPFR_ASSERTN (k < ((unsigned long)(-1) >> 1)); + + for (;; ++n, s = 0) /* overflow of n is innocuous */ + { + if ( ((f = k ? 0 : C (m, r)) < 0) || + (mpfr_random_deviate_reset (q), + !mpfr_random_deviate_less (q, s ? x : p, r)) || + ((f = k ? C (m, r) : f) < 0) || + (f == 0 && + (mpfr_random_deviate_reset (p), + !mpfr_random_deviate_less (p, x, r))) ) + break; + mpfr_random_deviate_swap (p, q); /* an efficient way of doing p = q */ + } + return (n & 1U) == 0; +} + +/* return a normal random deviate with mean 0 and variance 1 as a MPFR */ +int +mpfr_nrandom (mpfr_t z, gmp_randstate_t r, mpfr_rnd_t rnd) +{ + mpfr_random_deviate_t x, p, q; + int inex; + unsigned long k, j; + + mpfr_random_deviate_init (x); + mpfr_random_deviate_init (p); + mpfr_random_deviate_init (q); + for (;;) + { + k = G (r, p, q); /* step 1 */ + if (!P (k, k - 1, r, p, q)) + continue; /* step 2 */ + mpfr_random_deviate_reset (x); /* step 3 */ + for (j = 0; j <= k && B (k, x, r, p, q); ++j); /* step 4 */ + if (j > k) + break; + } + mpfr_random_deviate_clear (q); + mpfr_random_deviate_clear (p); + /* steps 5, 6, 7 */ + inex = mpfr_random_deviate_value (gmp_urandomb_ui (r, 1), k, x, z, r, rnd); + mpfr_random_deviate_clear (x); + return inex; +} |