1 /*============================================================================
2 * Random number generation.
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *----------------------------------------------------------------------------*/
32
33 #include <assert.h>
34 #include <math.h>
35
36 /*----------------------------------------------------------------------------
37 * Local headers
38 *----------------------------------------------------------------------------*/
39
40 #include "bft_error.h"
41
42 /*----------------------------------------------------------------------------
43 * Header for the current file
44 *----------------------------------------------------------------------------*/
45
46 #include "cs_random.h"
47
48 /*----------------------------------------------------------------------------*/
49
50 BEGIN_C_DECLS
51
52 /*=============================================================================
53 * Additional doxygen documentation
54 *============================================================================*/
55
56 /*!
57 \file cs_random.c
58 Random number generation.
59
60 Based on the uniform, gaussian, and poisson random number generation code
61 from netlib.org: lagged (-273,-607) Fibonacci; Box-Muller;
62 by W.P. Petersen, IPS, ETH Zuerich.
63 */
64
65 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
66
67 /*=============================================================================
68 * Macro definitions
69 *============================================================================*/
70
71 /*============================================================================
72 * Type definitions
73 *============================================================================*/
74
75 /*============================================================================
76 * Static global variables
77 *============================================================================*/
78
79 static struct {
80 int fill_1[1214];
81 double *buff;
82 int ptr;
83 int e_2;
84 } klotz0_1 = {{0}, NULL, 0, 0};
85
86 static struct {
87 double xbuff[1024];
88 int first;
89 int xptr;
90 double e_3;
91 } klotz1_1 = {{0}, 0, 0, 0.};
92
93 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
94
95 /*============================================================================
96 * Global variables
97 *============================================================================*/
98
99 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
100
101 /*============================================================================
102 * Private function definitions
103 *============================================================================*/
104
105 /*----------------------------------------------------------------------------*/
106 /*!
107 * \brief Save static variables used by uniform number generator.
108 *
109 * \param[out] save_block saved state values
110 */
111 /*----------------------------------------------------------------------------*/
112
113 static void
_random_uniform_save(cs_real_t save_block[608])114 _random_uniform_save(cs_real_t save_block[608])
115 {
116 /* Saves blocks klotz0, containing seeds and
117 pointer to position in seed block. The entire contents
118 of klotz0 (pointer in buff, and buff) must be saved. */
119
120 save_block[0] = (double) klotz0_1.ptr;
121 # pragma omp simd
122 for (int i = 0; i < 607; ++i) {
123 save_block[i + 1] = klotz0_1.buff[i];
124 }
125 }
126
127 /*----------------------------------------------------------------------------*/
128 /*!
129 * \brief Restore static variables used by uniform number generator.
130 *
131 * \param[in] save_block saved state values
132 */
133 /*----------------------------------------------------------------------------*/
134
135 static void
_random_uniform_restore(cs_real_t save_block[608])136 _random_uniform_restore(cs_real_t save_block[608])
137 {
138 /* Restores block klotz0, containing seeds and pointer
139 to position in seed block. The entire contents
140 of klotz0 must be restored. */
141
142 klotz0_1.ptr = (int) save_block[0];
143 # pragma omp simd
144 for (int i = 0; i < 607; ++i) {
145 klotz0_1.buff[i] = save_block[i + 1];
146 }
147 }
148
149 /*----------------------------------------------------------------------------*/
150 /*!
151 * \brief Internal function for Box_Muller normal distribution.
152 */
153 /*----------------------------------------------------------------------------*/
154
155 static void
_normal00(void)156 _normal00(void)
157 {
158 /* Local variables */
159 double twopi, r1, r2, t1, t2;
160
161 twopi = 6.2831853071795862;
162 cs_random_uniform(1024, klotz1_1.xbuff);
163 #pragma omp simd
164 for (int i = 0; i < 1023; i += 2) {
165 r1 = twopi * klotz1_1.xbuff[i];
166 t1 = cos(r1);
167 t2 = sin(r1);
168 r2 = sqrt(-2.*(log(1. - klotz1_1.xbuff[i+1])));
169 klotz1_1.xbuff[i] = t1 * r2;
170 klotz1_1.xbuff[i+1] = t2 * r2;
171 }
172 }
173
174 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
175
176 /*=============================================================================
177 * Public function definitions
178 *============================================================================*/
179
180 /*----------------------------------------------------------------------------*/
181 /*!
182 * \brief Initialize random number generator.
183 *
184 * Generates initial seed buffer by linear congruential method.
185 * Taken from Marsaglia, FSU report FSU-SCRI-87-50.
186 *
187 * \param[in] seed variable seed, with 0 < seed < 31328
188 */
189 /*----------------------------------------------------------------------------*/
190
191 void
cs_random_seed(int seed)192 cs_random_seed(int seed)
193 {
194 int kl = 9373;
195 int ij = 1802;
196
197 /* Set pointer */
198
199 klotz0_1.buff = (double *)(klotz0_1.fill_1);
200
201 /* Seed should be > 0 and < 31328 */
202 if (seed > 0)
203 ij = seed % 31328;
204
205 int i = ij / 177 % 177 + 2;
206 int j = ij % 177 + 2;
207 int k = kl / 169 % 178 + 1;
208 int l = kl % 169;
209
210 for (int ii = 0; ii < 607; ++ii) {
211 double s = 0.;
212 double t = .5;
213 for (int jj = 1; jj <= 24; ++jj) {
214 int m = i * j % 179 * k % 179;
215 i = j;
216 j = k;
217 k = m;
218 l = (l * 53 + 1) % 169;
219 if (l * m % 64 >= 32) {
220 s += t;
221 }
222 t *= (double).5;
223 }
224 klotz0_1.buff[ii] = s;
225 }
226 }
227
228 /*----------------------------------------------------------------------------*/
229 /*!
230 * \brief Uniform distribution random number generator.
231 *
232 * Portable lagged Fibonacci series uniform random number generator
233 * with "lags" -273 und -607:
234 * W.P. Petersen, IPS, ETH Zuerich, 19 Mar. 92
235 *
236 * \param[in] n number of values to compute
237 * \param[out] a pseudo-random numbers following uniform distribution
238 */
239 /*----------------------------------------------------------------------------*/
240
241 void
cs_random_uniform(cs_lnum_t n,cs_real_t a[])242 cs_random_uniform(cs_lnum_t n,
243 cs_real_t a[])
244 {
245 int buffsz = 607;
246
247 int bptr, aptr0, i, k, q, kptr;
248 double t;
249 int left, vl, qq, k273, k607;
250
251 int aptr = 0;
252 int nn = n;
253
254 L1:
255
256 if (nn <= 0)
257 return;
258
259 /* factor nn = q*607 + r */
260
261 q = (nn - 1) / 607;
262 left = buffsz - klotz0_1.ptr;
263
264 if (q <= 1) {
265
266 /* only one or fewer full segments */
267
268 if (nn < left) {
269 kptr = klotz0_1.ptr;
270 for (i = 0; i < nn; ++i) {
271 a[i + aptr] = klotz0_1.buff[kptr + i];
272 }
273 klotz0_1.ptr += nn;
274 return;
275 }
276 else {
277 kptr = klotz0_1.ptr;
278 # pragma omp simd
279 for (i = 0; i < left; ++i) {
280 a[i + aptr] = klotz0_1.buff[kptr + i];
281 }
282 klotz0_1.ptr = 0;
283 aptr += left;
284 nn -= left;
285 /* buff -> buff case */
286 vl = 273;
287 k273 = 334;
288 k607 = 0;
289 for (k = 0; k < 3; ++k) {
290 # pragma omp simd
291 for (i = 0; i < vl; ++i) {
292 t = klotz0_1.buff[k273+i]+klotz0_1.buff[k607+i];
293 klotz0_1.buff[k607+i] = t - (double) ((int) t);
294 }
295 k607 += vl;
296 k273 += vl;
297 vl = 167;
298 if (k == 0) {
299 k273 = 0;
300 }
301 }
302 goto L1;
303 }
304 } else {
305
306 /* more than 1 full segment */
307
308 kptr = klotz0_1.ptr;
309 # pragma omp simd
310 for (i = 0; i < left; ++i) {
311 a[i + aptr] = klotz0_1.buff[kptr + i];
312 }
313 nn -= left;
314 klotz0_1.ptr = 0;
315 aptr += left;
316
317 /* buff -> a(aptr0) */
318
319 vl = 273;
320 k273 = 334;
321 k607 = 0;
322 for (k = 0; k < 3; ++k) {
323 if (k == 0) {
324 # pragma omp simd
325 for (i = 0; i < vl; ++i) {
326 t = klotz0_1.buff[k273+i]+klotz0_1.buff[k607+i];
327 a[aptr + i] = t - (double) ((int) t);
328 }
329 k273 = aptr;
330 k607 += vl;
331 aptr += vl;
332 vl = 167;
333 } else {
334 # pragma omp simd
335 for (i = 0; i < vl; ++i) {
336 t = a[k273 + i] + klotz0_1.buff[k607 + i];
337 a[aptr + i] = t - (double) ((int) t);
338 }
339 k607 += vl;
340 k273 += vl;
341 aptr += vl;
342 }
343 }
344 nn += -607;
345
346 /* a(aptr-607) -> a(aptr) for last of the q-1 segments */
347
348 aptr0 = aptr - 607;
349 vl = 607;
350
351 for (qq = 0; qq < q-2; ++qq) {
352 k273 = aptr0 + 334;
353 # pragma omp simd
354 for (i = 0; i < vl; ++i) {
355 t = a[k273 + i] + a[aptr0 + i];
356 a[aptr + i] = t - (double) ((int) t);
357 }
358 nn += -607;
359 aptr += vl;
360 aptr0 += vl;
361 }
362
363 /* a(aptr0) -> buff, last segment before residual */
364
365 vl = 273;
366 k273 = aptr0 + 334;
367 k607 = aptr0;
368 bptr = 0;
369 for (k = 0; k < 3; ++k) {
370 if (k == 0) {
371 # pragma omp simd
372 for (i = 0; i < vl; ++i) {
373 t = a[k273 + i] + a[k607 + i];
374 klotz0_1.buff[bptr + i] = t - (double) ((int) t);
375 }
376 k273 = 0;
377 k607 += vl;
378 bptr += vl;
379 vl = 167;
380 } else {
381 # pragma omp simd
382 for (i = 0; i < vl; ++i) {
383 t = klotz0_1.buff[k273 + i] + a[k607 + i];
384 klotz0_1.buff[bptr + i] = t - (double) ((int) t);
385 }
386 k607 += vl;
387 k273 += vl;
388 bptr += vl;
389 }
390 }
391 goto L1;
392 }
393 }
394
395 /*----------------------------------------------------------------------------*/
396 /*!
397 * \brief Normal distribution random number generator.
398 *
399 * Box-Muller method for Gaussian random numbers.
400 *
401 * \param[in] n number of values to compute
402 * \param[out] x pseudo-random numbers following normal distribution
403 */
404 /*----------------------------------------------------------------------------*/
405
406 void
cs_random_normal(cs_lnum_t n,cs_real_t x[])407 cs_random_normal(cs_lnum_t n,
408 cs_real_t x[])
409 {
410 int buffsz = 1024;
411
412 /* Local variables */
413 int left, i, nn, ptr, kptr;
414
415 /* Box-Muller method for Gaussian random numbers */
416
417 nn = n;
418 if (nn <= 0)
419 return;
420 if (klotz1_1.first == 0) {
421 _normal00();
422 klotz1_1.first = 1;
423 }
424 ptr = 0;
425
426 L1:
427 left = buffsz - klotz1_1.xptr;
428 if (nn < left) {
429 kptr = klotz1_1.xptr;
430 # pragma omp simd
431 for (i = 0; i < nn; ++i) {
432 x[i + ptr] = klotz1_1.xbuff[kptr + i];
433 }
434 klotz1_1.xptr += nn;
435 return;
436 } else {
437 kptr = klotz1_1.xptr;
438 # pragma omp simd
439 for (i = 0; i < left; ++i) {
440 x[i + ptr] = klotz1_1.xbuff[kptr + i];
441 }
442 klotz1_1.xptr = 0;
443 ptr += left;
444 nn -= left;
445 _normal00();
446 goto L1;
447 }
448 }
449
450 /*----------------------------------------------------------------------------*/
451 /*!
452 * \brief Poisson distribution random number generator.
453 *
454 * q(mu,p) = exp(-mu) mu**p/p!
455 *
456 * \param[in] n number of values to compute
457 * \param[in] mu Poisson distribution parameter
458 * \param[out] p pseudo-random numbers following Poisson distribution
459 */
460 /*----------------------------------------------------------------------------*/
461
462 void
cs_random_poisson(cs_lnum_t n,cs_real_t mu,int p[])463 cs_random_poisson(cs_lnum_t n,
464 cs_real_t mu,
465 int p[])
466 {
467 /* Local variables */
468 int left, indx[1024], i, k;
469 double q[1024], u[1024];
470 int nsegs, p0;
471 double q0;
472 int ii, jj;
473 int nl0;
474 double pmu;
475
476 if (n <= 0)
477 return;
478
479 pmu = exp(-mu);
480 p0 = 0;
481
482 nsegs = (n - 1) / 1024;
483 left = n - (nsegs << 10);
484 ++nsegs;
485 nl0 = left;
486
487 for (k = 0; k < nsegs; ++k) {
488
489 for (i = 0; i < left; ++i) {
490 indx[i] = i;
491 p[p0 + i] = 0;
492 q[i] = 1.;
493 }
494
495 /* Begin iterative loop on segment of p's */
496
497 L1:
498
499 /* Get the needed uniforms */
500
501 cs_random_uniform(left, u);
502
503 jj = 0;
504
505 for (i = 0; i < left; ++i) {
506 ii = indx[i];
507 q0 = q[ii] * u[i];
508 q[ii] = q0;
509 if (q0 > pmu) {
510 indx[jj++] = ii;
511 ++p[p0 + ii];
512 }
513 }
514
515 /* any left in this segment? */
516
517 left = jj;
518 if (left > 0) {
519 goto L1;
520 }
521
522 p0 += nl0;
523 nl0 = 1024;
524 left = 1024;
525 }
526 }
527
528 /*----------------------------------------------------------------------------*/
529 /*!
530 * \brief Save static variables used by random number generator.
531 *
532 * \param[out] save_block saved state values
533 */
534 /*----------------------------------------------------------------------------*/
535
536 void
cs_random_save(cs_real_t save_block[1634])537 cs_random_save(cs_real_t save_block[1634])
538 {
539 int i, k;
540
541 /* The entire contents of blocks klotz0 and klotz1 must be saved. */
542
543 if (klotz1_1.first == 0) {
544 _normal00();
545 klotz1_1.first = 1;
546 }
547
548 _random_uniform_save(save_block);
549
550 save_block[608] = (double) klotz1_1.first;
551 save_block[609] = (double) klotz1_1.xptr;
552 k = 610;
553 # pragma omp simd
554 for (i = 0; i < 1024; ++i) {
555 save_block[i + k] = klotz1_1.xbuff[i];
556 }
557
558 }
559
560 /*----------------------------------------------------------------------------*/
561 /*!
562 * \brief Restore static variables used by random number generator.
563 *
564 * \param[out] save_block saved state values
565 */
566 /*----------------------------------------------------------------------------*/
567
568 void
cs_random_restore(cs_real_t save_block[1634])569 cs_random_restore(cs_real_t save_block[1634])
570 {
571 int i, k;
572
573 /* The entire contents of klotz0 and klotz1 must be restored. */
574
575 _random_uniform_restore(save_block);
576 klotz1_1.first = (int) save_block[608];
577 if (klotz1_1.first == 0)
578 bft_error(__FILE__, __LINE__, 0,
579 "In %s, restore of uninitialized block.", __func__);
580 klotz1_1.xptr = (int) save_block[609];
581 k = 610;
582 # pragma omp simd
583 for (i = 0; i < 1024; ++i) {
584 klotz1_1.xbuff[i] = save_block[i + k];
585 }
586 }
587
588 /*----------------------------------------------------------------------------*/
589
590 END_C_DECLS
591