contrib/mpfr/src/grandom.c

   1 /* mpfr_grandom (rop1, rop2, state, rnd_mode) -- Generate up to two
   2    pseudorandom real numbers according to a standard normal gaussian
   3    distribution and round it to the precision of rop1, rop2 according
   4    to the given rounding mode.
   5
   6 Copyright 2011 Free Software Foundation, Inc.
   7 Contributed by the Arenaire and Caramel projects, INRIA.
   8
   9 This file is part of the GNU MPFR Library.
  10
  11 The GNU MPFR Library is free software; you can redistribute it and/or modify
  12 it under the terms of the GNU Lesser General Public License as published by
  13 the Free Software Foundation; either version 3 of the License, or (at your
  14 option) any later version.
  15
  16 The GNU MPFR Library is distributed in the hope that it will be useful, but
  17 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  18 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  19 License for more details.
  20
  21 You should have received a copy of the GNU Lesser General Public License
  22 along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
  23 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
  24 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
  25
  26
  27 /* #define MPFR_NEED_LONGLONG_H */
  28 #include "mpfr-impl.h"
  29
  30
  31 int
  32 mpfr_grandom (mpfr_ptr rop1, mpfr_ptr rop2, gmp_randstate_t rstate,
  33               mpfr_rnd_t rnd)
  34 {
  35   int inex1, inex2, s1, s2;
  36   mpz_t x, y, xp, yp, t, a, b, s;
  37   mpfr_t sfr, l, r1, r2;
  38   mpfr_prec_t tprec, tprec0;
  39
  40   inex2 = inex1 = 0;
  41
  42   if (rop2 == NULL) /* only one output requested. */
  43     {
  44       tprec0 = MPFR_PREC (rop1);
  45     }
  46   else
  47     {
  48       tprec0 = MAX (MPFR_PREC (rop1), MPFR_PREC (rop2));
  49     }
  50
  51   tprec0 += 11;
  52
  53   /* We use "Marsaglia polar method" here (cf.
  54      George Marsaglia, Normal (Gaussian) random variables for supercomputers
  55      The Journal of Supercomputing, Volume 5, Number 1, 49–55
  56      DOI: 10.1007/BF00155857).
  57
  58      First we draw uniform x and y in [0,1] using mpz_urandomb (in
  59      fixed precision), and scale them to [-1, 1].
  60   */
  61
  62   mpz_init (xp);
  63   mpz_init (yp);
  64   mpz_init (x);
  65   mpz_init (y);
  66   mpz_init (t);
  67   mpz_init (s);
  68   mpz_init (a);
  69   mpz_init (b);
  70   mpfr_init2 (sfr, MPFR_PREC_MIN);
  71   mpfr_init2 (l, MPFR_PREC_MIN);
  72   mpfr_init2 (r1, MPFR_PREC_MIN);
  73   if (rop2 != NULL)
  74     mpfr_init2 (r2, MPFR_PREC_MIN);
  75
  76   mpz_set_ui (xp, 0);
  77   mpz_set_ui (yp, 0);
  78
  79   for (;;)
  80     {
  81       tprec = tprec0;
  82       do
  83         {
  84           mpz_urandomb (xp, rstate, tprec);
  85           mpz_urandomb (yp, rstate, tprec);
  86           mpz_mul (a, xp, xp);
  87           mpz_mul (b, yp, yp);
  88           mpz_add (s, a, b);
  89         }
  90       while (mpz_sizeinbase (s, 2) > tprec * 2); /* x^2 + y^2 <= 2^{2tprec} */
  91
  92       for (;;)
  93         {
  94           /* FIXME: compute s as s += 2x + 2y + 2 */
  95           mpz_add_ui (a, xp, 1);
  96           mpz_add_ui (b, yp, 1);
  97           mpz_mul (a, a, a);
  98           mpz_mul (b, b, b);
  99           mpz_add (s, a, b);
 100           if ((mpz_sizeinbase (s, 2) <= 2 * tprec) ||
 101               ((mpz_sizeinbase (s, 2) == 2 * tprec + 1) &&
 102                (mpz_scan1 (s, 0) == 2 * tprec)))
 103             goto yeepee;
 104           /* Extend by 32 bits */
 105           mpz_mul_2exp (xp, xp, 32);
 106           mpz_mul_2exp (yp, yp, 32);
 107           mpz_urandomb (x, rstate, 32);
 108           mpz_urandomb (y, rstate, 32);
 109           mpz_add (xp, xp, x);
 110           mpz_add (yp, yp, y);
 111           tprec += 32;
 112
 113           mpz_mul (a, xp, xp);
 114           mpz_mul (b, yp, yp);
 115           mpz_add (s, a, b);
 116           if (mpz_sizeinbase (s, 2) > tprec * 2)
 117             break;
 118         }
 119     }
 120  yeepee:
 121
 122   /* FIXME: compute s with s -= 2x + 2y + 2 */
 123   mpz_mul (a, xp, xp);
 124   mpz_mul (b, yp, yp);
 125   mpz_add (s, a, b);
 126   /* Compute the signs of the output */
 127   mpz_urandomb (x, rstate, 2);
 128   s1 = mpz_tstbit (x, 0);
 129   s2 = mpz_tstbit (x, 1);
 130   for (;;)
 131     {
 132       /* s = xp^2 + yp^2 (loop invariant) */
 133       mpfr_set_prec (sfr, 2 * tprec);
 134       mpfr_set_prec (l, tprec);
 135       mpfr_set_z (sfr, s, MPFR_RNDN); /* exact */
 136       mpfr_mul_2si (sfr, sfr, -2 * tprec, MPFR_RNDN); /* exact */
 137       mpfr_log (l, sfr, MPFR_RNDN);
 138       mpfr_neg (l, l, MPFR_RNDN);
 139       mpfr_mul_2si (l, l, 1, MPFR_RNDN);
 140       mpfr_div (l, l, sfr, MPFR_RNDN);
 141       mpfr_sqrt (l, l, MPFR_RNDN);
 142
 143       mpfr_set_prec (r1, tprec);
 144       mpfr_mul_z (r1, l, xp, MPFR_RNDN);
 145       mpfr_div_2ui (r1, r1, tprec, MPFR_RNDN); /* exact */
 146       if (s1)
 147         mpfr_neg (r1, r1, MPFR_RNDN);
 148       if (MPFR_CAN_ROUND (r1, tprec - 2, MPFR_PREC (rop1), rnd))
 149         {
 150           if (rop2 != NULL)
 151             {
 152               mpfr_set_prec (r2, tprec);
 153               mpfr_mul_z (r2, l, yp, MPFR_RNDN);
 154               mpfr_div_2ui (r2, r2, tprec, MPFR_RNDN); /* exact */
 155               if (s2)
 156                 mpfr_neg (r2, r2, MPFR_RNDN);
 157               if (MPFR_CAN_ROUND (r2, tprec - 2, MPFR_PREC (rop2), rnd))
 158                 break;
 159             }
 160           else
 161             break;
 162         }
 163       /* Extend by 32 bits */
 164       mpz_mul_2exp (xp, xp, 32);
 165       mpz_mul_2exp (yp, yp, 32);
 166       mpz_urandomb (x, rstate, 32);
 167       mpz_urandomb (y, rstate, 32);
 168       mpz_add (xp, xp, x);
 169       mpz_add (yp, yp, y);
 170       tprec += 32;
 171       mpz_mul (a, xp, xp);
 172       mpz_mul (b, yp, yp);
 173       mpz_add (s, a, b);
 174     }
 175   inex1 = mpfr_set (rop1, r1, rnd);
 176   if (rop2 != NULL)
 177     {
 178       inex2 = mpfr_set (rop2, r2, rnd);
 179       inex2 = mpfr_check_range (rop2, inex2, rnd);
 180     }
 181   inex1 = mpfr_check_range (rop1, inex1, rnd);
 182
 183   if (rop2 != NULL)
 184     mpfr_clear (r2);
 185   mpfr_clear (r1);
 186   mpfr_clear (l);
 187   mpfr_clear (sfr);
 188   mpz_clear (b);
 189   mpz_clear (a);
 190   mpz_clear (s);
 191   mpz_clear (t);
 192   mpz_clear (y);
 193   mpz_clear (x);
 194   mpz_clear (yp);
 195   mpz_clear (xp);
 196
 197   return INEX (inex1, inex2);
 198 }