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
6Copyright 2011, 2012, 2013 Free Software Foundation, Inc.
7Contributed by the AriC and Caramel projects, INRIA.
8
9This file is part of the GNU MPFR Library.
10
11The GNU MPFR Library is free software; you can redistribute it and/or modify
12it under the terms of the GNU Lesser General Public License as published by
13the Free Software Foundation; either version 3 of the License, or (at your
14option) any later version.
15
16The GNU MPFR Library is distributed in the hope that it will be useful, but
17WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
19License for more details.
20
21You should have received a copy of the GNU Lesser General Public License
22along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
23http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2451 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
25
26
27/* #define MPFR_NEED_LONGLONG_H */
28#include "mpfr-impl.h"
29
30
31int
32mpfr_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}
199