Home | History | Annotate | Line # | Download | only in src 
      1 /* mpfr_erandom (rop, state, rnd_mode) -- Generate an exponential deviate with
      2    mean 1 and round it to the precision of rop according to the given rounding
      3    mode.
      4 
      5 Copyright 2013-2023 Free Software Foundation, Inc.
      6 Contributed by Charles Karney <charles (at) karney.com>, SRI International.
      7 
      8 This file is part of the GNU MPFR Library.
      9 
     10 The GNU MPFR Library is free software; you can redistribute it and/or modify
     11 it under the terms of the GNU Lesser General Public License as published by
     12 the Free Software Foundation; either version 3 of the License, or (at your
     13 option) any later version.
     14 
     15 The GNU MPFR Library is distributed in the hope that it will be useful, but
     16 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
     17 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
     18 License for more details.
     19 
     20 You should have received a copy of the GNU Lesser General Public License
     21 along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
     22 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
     23 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
     24 
     25 /*
     26  * Sampling from the exponential distribution with unit mean using the method
     27  * given in John von Neumann, Various techniques used in connection with random
     28  * digits, in A. S. Householder, G. E. Forsythe, and H. H. Germond, editors,
     29  * "Monte Carlo Method", number 12 in Applied Mathematics Series, pp. 36-38
     30  * (NBS, Washington, DC, 1951), Proceedings of a symposium held June 29-July 1,
     31  * 1949, in Los Angeles.
     32  *
     33  * A modification to this algorithm is given in:
     34  *   Charles F. F. Karney,
     35  *   "Sampling exactly from the normal distribution",
     36  *   ACM Trans. Math. Software 42(1), 3:1-14 (Jan. 2016).
     37  *   https://dx.doi.org/10.1145/2710016
     38  *   https://arxiv.org/abs/1303.6257
     39  * Although this improves the bit efficiency, in practice, it results in
     40  * a slightly slower algorithm for MPFR. So here the original von Neumann
     41  * algorithm is used.
     42  *
     43  * There are a few "weasel words" regarding the accuracy of this
     44  * implementation.  The algorithm produces exactly rounded exponential deviates
     45  * provided that gmp's random number engine delivers truly random bits.  If it
     46  * did, the algorithm would be perfect; however, this implementation would have
     47  * problems, e.g., in that the integer part of the exponential deviate is
     48  * represented by an unsigned long, whereas in reality the integer part in
     49  * unbounded.  In this implementation, asserts catch overflow in the integer
     50  * part and similar (very, very) unlikely events.  In reality, of course, gmp's
     51  * random number engine has a finite internal state (19937 bits in the case of
     52  * the MT19937 method).  This means that these unlikely events in fact won't
     53  * occur.  If the asserts are triggered, then this is an indication that the
     54  * random number engine is defective.  (Even if a hardware random number
     55  * generator were used, the most likely explanation for the triggering of the
     56  * asserts would be that the hardware generator was broken.)
     57  */
     58 
     59 #include "random_deviate.h"
     60 
     61 /* true with prob exp(-x) */
     62 static int
     63 E (mpfr_random_deviate_t x, gmp_randstate_t r,
     64    mpfr_random_deviate_t p, mpfr_random_deviate_t q)
     65 {
     66   /* p and q are temporaries */
     67   mpfr_random_deviate_reset (p);
     68   if (!mpfr_random_deviate_less (p, x, r))
     69     return 1;
     70   for (;;)
     71     {
     72       mpfr_random_deviate_reset (q);
     73       if (!mpfr_random_deviate_less (q, p, r))
     74         return 0;
     75       mpfr_random_deviate_reset (p);
     76       if (!mpfr_random_deviate_less (p, q, r))
     77         return 1;
     78     }
     79 }
     80 
     81 /* return an exponential random deviate with mean 1 as a MPFR  */
     82 int
     83 mpfr_erandom (mpfr_ptr z, gmp_randstate_t r, mpfr_rnd_t rnd)
     84 {
     85   mpfr_random_deviate_t x, p, q;
     86   int inex;
     87   unsigned long k = 0;
     88 
     89   mpfr_random_deviate_init (x);
     90   mpfr_random_deviate_init (p);
     91   mpfr_random_deviate_init (q);
     92   while (!E(x, r, p, q))
     93     {
     94       ++k;
     95       /* Catch k wrapping around to 0; for a 32-bit unsigned long, the
     96        * probability of this is exp(-2^32)). */
     97       MPFR_ASSERTN (k != 0UL);
     98       mpfr_random_deviate_reset (x);
     99     }
    100   mpfr_random_deviate_clear (q);
    101   mpfr_random_deviate_clear (p);
    102   inex = mpfr_random_deviate_value (0, k, x, z, r, rnd);
    103   mpfr_random_deviate_clear (x);
    104   return inex;
    105 }
    106