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@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  }