efros_freeman  0.1
random_number.c
Go to the documentation of this file.
1 /*
2  Copyright (c) 2016 rafael grompone von gioi <grompone@gmail.com>
3 
4  Quilting is free software: you can redistribute it and/or modify
5  it under the terms of the GNU Affero General Public License as
6  published by the Free Software Foundation, either version 3 of the
7  License, or (at your option) any later version.
8 
9  This program is distributed in the hope that it will be useful,
10  but WITHOUT ANY WARRANTY; without even the implied warranty of
11  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  GNU Affero General Public License for more details.
13 
14  You should have received a copy of the GNU Affero General Public License
15  along with this program. If not, see <http://www.gnu.org/licenses/>.
16 */
17 
18 /**
19 * @file random_number.c
20 * @brief this header file contains all definitions.
21 *
22 * @author rafael grompone von gioi
23 *
24 * @date 05/02/2011
25 */
26 
27 /*----------------------------------------------------------------------------*/
28 /* random_number generator limites. */
29 #define RANDOM_MIN 1
30 #define RANDOM_MAX 2147483562
31 
32 #include <stdlib.h>
33 
34 /*----------------------------------------------------------------------------*/
35 /* Uniform random number generator in [1,2147483562].
36 
37  This algorithm is described and analyzed on:
38 
39  "Efficient and Portable Combined Random Number Generators"
40  by Pierre L'Ecuyer, Communications of the ACM, vol. 31, num. 6,
41  pp. 742-749 and 774, June 1988.
42 
43  This code is a modification from the code available at
44  http://cg.scs.carleton.ca/~luc/lecuyer.c
45  on February 5, 2011.
46 
47  The algorithm is based on two combined Multiplicative Linear
48  Congruential Generators (MLCGs). The generator produces double
49  numbers in the range (0,1), (it will never return 0.0 or 1.0).
50  The period of the generator is about 2.30584 x 10^18,
51  in the order of 2^61.
52 
53  A normal call is in the form 'uniform_rand(NULL)'.
54  But it can also be called giving a pointer to long with a seed:
55  long seed = 263287632;
56  uniform_rand(&seed);
57  This is usually done only once. The function 'rand_time_seed'
58  performs this initilization with the current time.
59 
60  The state of the generator is determined by the static variables
61  s1 and s2, that must take values in the following intervals:
62 
63  s1 in [1,2147483562]
64  s2 in [1,2147483398]
65  */
66 long random_number(long * seed)
67 {
68  static long s1 = 55555;
69  static long s2 = 99999;
70  #pragma omp threadprivate(s1)
71  #pragma omp threadprivate(s2)
72  long k,z;
73 
74  /* Initialization of the seed, if demanded. */
75  if( seed != NULL ) s1 = 1 + (*seed % 2147483561); /* s1 in [1,2147483562] */
76 
77  /* First MLCG: s1 = (40014 * s1) mod 2147483563, efficient implementation */
78  k = s1 / 53668;
79  s1 = 40014 * ( s1 % 53668 ) - k * 12211;
80  if( s1 < 0 ) s1 += 2147483563;
81 
82  /* Second MLCG: s2 = (40692 * s2) mod 2147483399, efficient implementation */
83  k = s2 / 52774;
84  s2 = 40692 * ( s2 % 52774 ) - k * 3791;
85  if( s2 < 0 ) s2 += 2147483399;
86 
87  /* Combination of both MLCGs */
88  z = ( s1 - 2147483563 ) + s2;
89  if(z < 1) z += 2147483562;
90 
91  return z;
92 }
93