GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-44ea597914
lrand48.c
Go to the documentation of this file.
1 /*!
2  * \file lib/gis/lrand48.c
3  *
4  * \brief GIS Library - Pseudo-random number generation
5  *
6  * (C) 2014 by the GRASS Development Team
7  *
8  * This program is free software under the GNU General Public License
9  * (>=v2). Read the file COPYING that comes with GRASS for details.
10  *
11  * \author Glynn Clements
12  */
13 
14 #include <stdio.h>
15 #include <stdlib.h>
16 #include <string.h>
17 #include <errno.h>
18 
19 #include <grass/gis.h>
20 #include <grass/glocale.h>
21 
22 #ifdef HAVE_GETTIMEOFDAY
23 #include <sys/time.h>
24 #else
25 #include <time.h>
26 #endif
27 
28 #include <sys/types.h>
29 #include <unistd.h>
30 
31 typedef unsigned short uint16;
32 typedef unsigned int uint32;
33 typedef signed int int32;
34 
35 static uint16 x0, x1, x2;
36 static const uint32 a0 = 0xE66D;
37 static const uint32 a1 = 0xDEEC;
38 static const uint32 a2 = 0x5;
39 
40 static const uint32 b0 = 0xB;
41 
42 static int seeded;
43 
44 #define LO(x) ((x)&0xFFFFU)
45 #define HI(x) ((x) >> 16)
46 
47 /*!
48  * \brief Seed the pseudo-random number generator
49  *
50  * \param seedval 32-bit integer used to seed the PRNG
51  */
52 
53 void G_srand48(long seedval)
54 {
55  uint32 x = (uint32) * (unsigned long *)&seedval;
56 
57  x2 = (uint16)HI(x);
58  x1 = (uint16)LO(x);
59  x0 = (uint16)0x330E;
60  seeded = 1;
61 }
62 
63 /*!
64  * \brief Seed the pseudo-random number generator from the time and PID
65  *
66  * A weak hash of the current time and PID is generated and used to
67  * seed the PRNG
68  *
69  * \return generated seed value passed to G_srand48()
70  */
71 
72 long G_srand48_auto(void)
73 {
74  unsigned long seed;
75  char *grass_random_seed = getenv("GRASS_RANDOM_SEED");
76 
77  if (!grass_random_seed)
78  grass_random_seed = getenv("SOURCE_DATE_EPOCH");
79  if (grass_random_seed) {
80  seed = strtoull(grass_random_seed, NULL, 10);
81  }
82  else {
83  seed = (unsigned long)getpid();
84 
85 #ifdef HAVE_GETTIMEOFDAY
86  {
87  struct timeval tv;
88 
89  if (gettimeofday(&tv, NULL) < 0)
90  G_fatal_error(_("gettimeofday failed: %s"), strerror(errno));
91  seed += (unsigned long)tv.tv_sec;
92  seed += (unsigned long)tv.tv_usec;
93  }
94 #else
95  {
96  time_t t = time(NULL);
97 
98  seed += (unsigned long)t;
99  }
100 #endif
101  }
102 
103  G_srand48((long)seed);
104  return (long)seed;
105 }
106 
107 static void G__next(void)
108 {
109  uint32 a0x0 = a0 * x0;
110  uint32 a0x1 = a0 * x1;
111  uint32 a0x2 = a0 * x2;
112  uint32 a1x0 = a1 * x0;
113  uint32 a1x1 = a1 * x1;
114  uint32 a2x0 = a2 * x0;
115 
116  uint32 y0 = LO(a0x0) + b0;
117  uint32 y1 = LO(a0x1) + LO(a1x0) + HI(a0x0);
118  uint32 y2 = LO(a0x2) + LO(a1x1) + LO(a2x0) + HI(a0x1) + HI(a1x0);
119 
120  if (!seeded)
121  G_fatal_error(_("Pseudo-random number generator not seeded"));
122 
123  x0 = (uint16)LO(y0);
124  y1 += HI(y0);
125  x1 = (uint16)LO(y1);
126  y2 += HI(y1);
127  x2 = (uint16)LO(y2);
128 }
129 
130 /*!
131  * \brief Generate an integer in the range [0, 2^31)
132  *
133  * \return the generated value
134  */
135 
136 long G_lrand48(void)
137 {
138  uint32 r;
139 
140  G__next();
141  r = ((uint32)x2 << 15) | ((uint32)x1 >> 1);
142  return (long)r;
143 }
144 
145 /*!
146  * \brief Generate an integer in the range [-2^31, 2^31)
147  *
148  * \return the generated value
149  */
150 
151 long G_mrand48(void)
152 {
153  uint32 r;
154 
155  G__next();
156  r = ((uint32)x2 << 16) | ((uint32)x1);
157  return (long)(int32)r;
158 }
159 
160 /*!
161  * \brief Generate a floating-point value in the range [0,1)
162  *
163  * \return the generated value
164  */
165 
166 double G_drand48(void)
167 {
168  double r = 0.0;
169 
170  G__next();
171  r += x2;
172  r *= 0x10000;
173  r += x1;
174  r *= 0x10000;
175  r += x0;
176  r /= 281474976710656.0; /* 2^48 */
177  return r;
178 }
179 
180 /*
181 
182  Test program
183 
184  int main(int argc, char **argv)
185  {
186  long s = (argc > 1) ? atol(argv[1]) : 0;
187  int i;
188 
189  srand48(s);
190  G_srand48(s);
191 
192  for (i = 0; i < 100; i++) {
193  printf("%.50f %.50f\n", drand48(), G_drand48());
194  printf("%lu %lu\n", lrand48(), G_lrand48());
195  printf("%ld %ld\n", mrand48(), G_mrand48());
196  }
197 
198  return 0;
199  }
200 
201  */
#define NULL
Definition: ccmath.h:32
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
#define _(str)
Definition: glocale.h:10
unsigned short uint16
Definition: lrand48.c:31
unsigned int uint32
Definition: lrand48.c:32
long G_mrand48(void)
Generate an integer in the range [-2^31, 2^31)
Definition: lrand48.c:151
long G_lrand48(void)
Generate an integer in the range [0, 2^31)
Definition: lrand48.c:136
long G_srand48_auto(void)
Seed the pseudo-random number generator from the time and PID.
Definition: lrand48.c:72
signed int int32
Definition: lrand48.c:33
#define HI(x)
Definition: lrand48.c:45
void G_srand48(long seedval)
Seed the pseudo-random number generator.
Definition: lrand48.c:53
#define LO(x)
Definition: lrand48.c:44
double G_drand48(void)
Generate a floating-point value in the range [0,1)
Definition: lrand48.c:166
double t
Definition: r_raster.c:39
double r
Definition: r_raster.c:39
#define x