/src/postgis/liblwgeom/lwrandom.c
Line | Count | Source |
1 | | /********************************************************************** |
2 | | * |
3 | | * PostGIS - Spatial Types for PostgreSQL |
4 | | * http://postgis.net |
5 | | * |
6 | | * PostGIS is free software: you can redistribute it and/or modify |
7 | | * it under the terms of the GNU General Public License as published by |
8 | | * the Free Software Foundation, either version 2 of the License, or |
9 | | * (at your option) any later version. |
10 | | * |
11 | | * PostGIS is distributed in the hope that it will be useful, |
12 | | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
13 | | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
14 | | * GNU General Public License for more details. |
15 | | * |
16 | | * You should have received a copy of the GNU General Public License |
17 | | * along with PostGIS. If not, see <http://www.gnu.org/licenses/>. |
18 | | * |
19 | | ********************************************************************** |
20 | | * |
21 | | * Copyright 2019 Mike Taves |
22 | | * |
23 | | **********************************************************************/ |
24 | | |
25 | | #include "lwrandom.h" |
26 | | |
27 | | #include <stdlib.h> |
28 | | #include <stdint.h> |
29 | | #include <stddef.h> |
30 | | #include <time.h> |
31 | | |
32 | | #ifdef _WIN32 |
33 | | #include <process.h> |
34 | | #define getpid _getpid |
35 | | #else |
36 | | #include <unistd.h> |
37 | | #endif |
38 | | |
39 | | static unsigned char _lwrandom_seed_set = 0; |
40 | | static int32_t _lwrandom_seed[3] = {0x330e, 0xabcd, 0x1234}; |
41 | | |
42 | | /* |
43 | | * Set seed for a random number generator. |
44 | | * Repeatable numbers are generated with seed values >= 1. |
45 | | * When seed is zero and has not previously been set, it is based on |
46 | | * Unix time (seconds) and process ID. */ |
47 | | void |
48 | | lwrandom_set_seed(int32_t seed) |
49 | 0 | { |
50 | 0 | if (seed == 0) |
51 | 0 | { |
52 | 0 | if (_lwrandom_seed_set == 0) |
53 | 0 | seed = (unsigned int)time(NULL) + (unsigned int)getpid() - 0xbadd; |
54 | 0 | else |
55 | 0 | return; |
56 | 0 | } |
57 | | /* s1 value between 1 and 2147483562 */ |
58 | 0 | _lwrandom_seed[1] = (((int64_t)seed + 0xfeed) % 2147483562) + 1; |
59 | | /* s2 value between 1 and 2147483398 */ |
60 | 0 | _lwrandom_seed[2] = ((((int64_t)seed + 0xdefeb) << 5) % 2147483398) + 1; |
61 | 0 | _lwrandom_seed_set = 1; |
62 | 0 | } |
63 | | |
64 | | /* for low-level external debugging */ |
65 | | void |
66 | | _lwrandom_set_seeds(int32_t s1, int32_t s2) |
67 | 0 | { |
68 | | /* _lwrandom_seed[0] not used */ |
69 | 0 | _lwrandom_seed[1] = s1; |
70 | 0 | _lwrandom_seed[2] = s2; |
71 | 0 | _lwrandom_seed_set = 1; |
72 | 0 | } |
73 | | int32_t |
74 | | _lwrandom_get_seed(size_t idx) |
75 | 0 | { |
76 | 0 | return _lwrandom_seed[idx]; |
77 | 0 | } |
78 | | |
79 | | /* |
80 | | * Generate a random floating-point value. |
81 | | * Values are uniformly distributed between 0 and 1. |
82 | | * |
83 | | * Authors: |
84 | | * Pierre L'Ecuyer (1988), see source code in Figure 3. |
85 | | * C version by John Burkardt, modified by Mike Taves. |
86 | | * |
87 | | * Reference: |
88 | | * Pierre L'Ecuyer, |
89 | | * Efficient and Portable Combined Random Number Generators, |
90 | | * Communications of the ACM, Volume 31, Number 6, June 1988, |
91 | | * pages 742-751. doi:10.1145/62959.62969 |
92 | | */ |
93 | | double |
94 | | lwrandom_uniform(void) |
95 | 0 | { |
96 | 0 | double value; |
97 | 0 | int32_t k; |
98 | 0 | int32_t z; |
99 | 0 | int32_t *s1 = &_lwrandom_seed[1]; |
100 | 0 | int32_t *s2 = &_lwrandom_seed[2]; |
101 | |
|
102 | 0 | k = *s1 / 53668; |
103 | 0 | *s1 = 40014 * (*s1 - k * 53668) - k * 12211; |
104 | 0 | if (*s1 < 0) |
105 | 0 | *s1 += 2147483563; |
106 | |
|
107 | 0 | k = *s2 / 52774; |
108 | 0 | *s2 = 40692 * (*s2 - k * 52774) - k * 3791; |
109 | 0 | if (*s2 < 0) |
110 | 0 | *s2 += 2147483399; |
111 | |
|
112 | 0 | z = *s1 - *s2; |
113 | 0 | if (z < 1) |
114 | 0 | z += 2147483562; |
115 | |
|
116 | 0 | value = (double)(z) / 2147483563.0; |
117 | |
|
118 | 0 | return value; |
119 | 0 | } |