diff options
Diffstat (limited to 'packages')
-rw-r--r-- | packages/base/hmatrix.cabal | 4 | ||||
-rw-r--r-- | packages/base/src/Internal/C/vector-aux.c | 1 | ||||
-rw-r--r-- | packages/base/src/Internal/C/windows_random_r.c | 392 | ||||
-rw-r--r-- | packages/base/src/Internal/C/windows_random_r.h | 35 |
4 files changed, 431 insertions, 1 deletions
diff --git a/packages/base/hmatrix.cabal b/packages/base/hmatrix.cabal index 31bea3e..e248886 100644 --- a/packages/base/hmatrix.cabal +++ b/packages/base/hmatrix.cabal | |||
@@ -25,6 +25,7 @@ build-type: Simple | |||
25 | extra-source-files: THANKS.md CHANGELOG | 25 | extra-source-files: THANKS.md CHANGELOG |
26 | 26 | ||
27 | extra-source-files: src/Internal/C/lapack-aux.h | 27 | extra-source-files: src/Internal/C/lapack-aux.h |
28 | src/Internal/C/windows_random_r.h | ||
28 | 29 | ||
29 | flag openblas | 30 | flag openblas |
30 | description: Link with OpenBLAS (https://github.com/xianyi/OpenBLAS) optimized libraries. | 31 | description: Link with OpenBLAS (https://github.com/xianyi/OpenBLAS) optimized libraries. |
@@ -74,9 +75,10 @@ library | |||
74 | Internal.Util | 75 | Internal.Util |
75 | Internal.Modular | 76 | Internal.Modular |
76 | Internal.Static | 77 | Internal.Static |
77 | 78 | ||
78 | C-sources: src/Internal/C/lapack-aux.c | 79 | C-sources: src/Internal/C/lapack-aux.c |
79 | src/Internal/C/vector-aux.c | 80 | src/Internal/C/vector-aux.c |
81 | src/Internal/C/windows_random_r.c | ||
80 | 82 | ||
81 | 83 | ||
82 | extensions: ForeignFunctionInterface, | 84 | extensions: ForeignFunctionInterface, |
diff --git a/packages/base/src/Internal/C/vector-aux.c b/packages/base/src/Internal/C/vector-aux.c index 5528a9d..a370214 100644 --- a/packages/base/src/Internal/C/vector-aux.c +++ b/packages/base/src/Internal/C/vector-aux.c | |||
@@ -7,6 +7,7 @@ typedef float complex TCF; | |||
7 | #undef complex | 7 | #undef complex |
8 | 8 | ||
9 | #include "lapack-aux.h" | 9 | #include "lapack-aux.h" |
10 | #include "windows_random_r.h" | ||
10 | 11 | ||
11 | #define V(x) x##n,x##p | 12 | #define V(x) x##n,x##p |
12 | 13 | ||
diff --git a/packages/base/src/Internal/C/windows_random_r.c b/packages/base/src/Internal/C/windows_random_r.c new file mode 100644 index 0000000..872cb2e --- /dev/null +++ b/packages/base/src/Internal/C/windows_random_r.c | |||
@@ -0,0 +1,392 @@ | |||
1 | #if defined(_WIN32) || defined(WIN32) | ||
2 | |||
3 | /* | ||
4 | Copyright (C) 1995-2015 Free Software Foundation, Inc. | ||
5 | |||
6 | The GNU C Library is free software; you can redistribute it and/or | ||
7 | modify it under the terms of the GNU Lesser General Public | ||
8 | License as published by the Free Software Foundation; either | ||
9 | version 2.1 of the License, or (at your option) any later version. | ||
10 | |||
11 | The GNU C Library 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 GNU | ||
14 | Lesser General Public License for more details. | ||
15 | |||
16 | You should have received a copy of the GNU Lesser General Public | ||
17 | License along with the GNU C Library; if not, see | ||
18 | <http://www.gnu.org/licenses/>. */ | ||
19 | |||
20 | /* | ||
21 | Copyright (C) 1983 Regents of the University of California. | ||
22 | All rights reserved. | ||
23 | |||
24 | Redistribution and use in source and binary forms, with or without | ||
25 | modification, are permitted provided that the following conditions | ||
26 | are met: | ||
27 | |||
28 | 1. Redistributions of source code must retain the above copyright | ||
29 | notice, this list of conditions and the following disclaimer. | ||
30 | 2. Redistributions in binary form must reproduce the above copyright | ||
31 | notice, this list of conditions and the following disclaimer in the | ||
32 | documentation and/or other materials provided with the distribution. | ||
33 | 4. Neither the name of the University nor the names of its contributors | ||
34 | may be used to endorse or promote products derived from this software | ||
35 | without specific prior written permission. | ||
36 | |||
37 | THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND | ||
38 | ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE | ||
39 | IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE | ||
40 | ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE | ||
41 | FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL | ||
42 | DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS | ||
43 | OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) | ||
44 | HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT | ||
45 | LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY | ||
46 | OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF | ||
47 | SUCH DAMAGE.*/ | ||
48 | |||
49 | /* | ||
50 | * This is derived from the Berkeley source: | ||
51 | * @(#)random.c 5.5 (Berkeley) 7/6/88 | ||
52 | * It was reworked for the GNU C Library by Roland McGrath. | ||
53 | * Rewritten to be reentrant by Ulrich Drepper, 1995 | ||
54 | */ | ||
55 | |||
56 | #include "windows_random_r.h" | ||
57 | |||
58 | /* An improved random number generation package. In addition to the standard | ||
59 | rand()/srand() like interface, this package also has a special state info | ||
60 | interface. The initstate() routine is called with a seed, an array of | ||
61 | bytes, and a count of how many bytes are being passed in; this array is | ||
62 | then initialized to contain information for random number generation with | ||
63 | that much state information. Good sizes for the amount of state | ||
64 | information are 32, 64, 128, and 256 bytes. The state can be switched by | ||
65 | calling the setstate() function with the same array as was initialized | ||
66 | with initstate(). By default, the package runs with 128 bytes of state | ||
67 | information and generates far better random numbers than a linear | ||
68 | congruential generator. If the amount of state information is less than | ||
69 | 32 bytes, a simple linear congruential R.N.G. is used. Internally, the | ||
70 | state information is treated as an array of longs; the zeroth element of | ||
71 | the array is the type of R.N.G. being used (small integer); the remainder | ||
72 | of the array is the state information for the R.N.G. Thus, 32 bytes of | ||
73 | state information will give 7 longs worth of state information, which will | ||
74 | allow a degree seven polynomial. (Note: The zeroth word of state | ||
75 | information also has some other information stored in it; see setstate | ||
76 | for details). The random number generation technique is a linear feedback | ||
77 | shift register approach, employing trinomials (since there are fewer terms | ||
78 | to sum up that way). In this approach, the least significant bit of all | ||
79 | the numbers in the state table will act as a linear feedback shift register, | ||
80 | and will have period 2^deg - 1 (where deg is the degree of the polynomial | ||
81 | being used, assuming that the polynomial is irreducible and primitive). | ||
82 | The higher order bits will have longer periods, since their values are | ||
83 | also influenced by pseudo-random carries out of the lower bits. The | ||
84 | total period of the generator is approximately deg*(2**deg - 1); thus | ||
85 | doubling the amount of state information has a vast influence on the | ||
86 | period of the generator. Note: The deg*(2**deg - 1) is an approximation | ||
87 | only good for large deg, when the period of the shift register is the | ||
88 | dominant factor. With deg equal to seven, the period is actually much | ||
89 | longer than the 7*(2**7 - 1) predicted by this formula. */ | ||
90 | |||
91 | |||
92 | |||
93 | /* For each of the currently supported random number generators, we have a | ||
94 | break value on the amount of state information (you need at least this many | ||
95 | bytes of state info to support this random number generator), a degree for | ||
96 | the polynomial (actually a trinomial) that the R.N.G. is based on, and | ||
97 | separation between the two lower order coefficients of the trinomial. */ | ||
98 | |||
99 | /* Linear congruential. */ | ||
100 | #define TYPE_0 0 | ||
101 | #define BREAK_0 8 | ||
102 | #define DEG_0 0 | ||
103 | #define SEP_0 0 | ||
104 | |||
105 | /* x**7 + x**3 + 1. */ | ||
106 | #define TYPE_1 1 | ||
107 | #define BREAK_1 32 | ||
108 | #define DEG_1 7 | ||
109 | #define SEP_1 3 | ||
110 | |||
111 | /* x**15 + x + 1. */ | ||
112 | #define TYPE_2 2 | ||
113 | #define BREAK_2 64 | ||
114 | #define DEG_2 15 | ||
115 | #define SEP_2 1 | ||
116 | |||
117 | /* x**31 + x**3 + 1. */ | ||
118 | #define TYPE_3 3 | ||
119 | #define BREAK_3 128 | ||
120 | #define DEG_3 31 | ||
121 | #define SEP_3 3 | ||
122 | |||
123 | /* x**63 + x + 1. */ | ||
124 | #define TYPE_4 4 | ||
125 | #define BREAK_4 256 | ||
126 | #define DEG_4 63 | ||
127 | #define SEP_4 1 | ||
128 | |||
129 | |||
130 | /* Array versions of the above information to make code run faster. | ||
131 | Relies on fact that TYPE_i == i. */ | ||
132 | |||
133 | #define MAX_TYPES 5 /* Max number of types above. */ | ||
134 | |||
135 | struct random_poly_info | ||
136 | { | ||
137 | int seps[MAX_TYPES]; | ||
138 | int degrees[MAX_TYPES]; | ||
139 | }; | ||
140 | |||
141 | static const struct random_poly_info random_poly_info = | ||
142 | { | ||
143 | { SEP_0, SEP_1, SEP_2, SEP_3, SEP_4 }, | ||
144 | { DEG_0, DEG_1, DEG_2, DEG_3, DEG_4 } | ||
145 | }; | ||
146 | |||
147 | |||
148 | |||
149 | |||
150 | /* Initialize the random number generator based on the given seed. If the | ||
151 | type is the trivial no-state-information type, just remember the seed. | ||
152 | Otherwise, initializes state[] based on the given "seed" via a linear | ||
153 | congruential generator. Then, the pointers are set to known locations | ||
154 | that are exactly rand_sep places apart. Lastly, it cycles the state | ||
155 | information a given number of times to get rid of any initial dependencies | ||
156 | introduced by the L.C.R.N.G. Note that the initialization of randtbl[] | ||
157 | for default usage relies on values produced by this routine. */ | ||
158 | int srandom_r (unsigned int seed, struct random_data *buf) | ||
159 | { | ||
160 | int type; | ||
161 | int32_t *state; | ||
162 | long int i; | ||
163 | int32_t word; | ||
164 | int32_t *dst; | ||
165 | int kc; | ||
166 | |||
167 | if (buf == NULL) | ||
168 | goto fail; | ||
169 | type = buf->rand_type; | ||
170 | if ((unsigned int) type >= MAX_TYPES) | ||
171 | goto fail; | ||
172 | |||
173 | state = buf->state; | ||
174 | /* We must make sure the seed is not 0. Take arbitrarily 1 in this case. */ | ||
175 | if (seed == 0) | ||
176 | seed = 1; | ||
177 | state[0] = seed; | ||
178 | if (type == TYPE_0) | ||
179 | goto done; | ||
180 | |||
181 | dst = state; | ||
182 | word = seed; | ||
183 | kc = buf->rand_deg; | ||
184 | for (i = 1; i < kc; ++i) | ||
185 | { | ||
186 | /* This does: | ||
187 | state[i] = (16807 * state[i - 1]) % 2147483647; | ||
188 | but avoids overflowing 31 bits. */ | ||
189 | long int hi = word / 127773; | ||
190 | long int lo = word % 127773; | ||
191 | word = 16807 * lo - 2836 * hi; | ||
192 | if (word < 0) | ||
193 | word += 2147483647; | ||
194 | *++dst = word; | ||
195 | } | ||
196 | |||
197 | buf->fptr = &state[buf->rand_sep]; | ||
198 | buf->rptr = &state[0]; | ||
199 | kc *= 10; | ||
200 | while (--kc >= 0) | ||
201 | { | ||
202 | int32_t discard; | ||
203 | (void) random_r (buf, &discard); | ||
204 | } | ||
205 | |||
206 | done: | ||
207 | return 0; | ||
208 | |||
209 | fail: | ||
210 | return -1; | ||
211 | } | ||
212 | |||
213 | /* Initialize the state information in the given array of N bytes for | ||
214 | future random number generation. Based on the number of bytes we | ||
215 | are given, and the break values for the different R.N.G.'s, we choose | ||
216 | the best (largest) one we can and set things up for it. srandom is | ||
217 | then called to initialize the state information. Note that on return | ||
218 | from srandom, we set state[-1] to be the type multiplexed with the current | ||
219 | value of the rear pointer; this is so successive calls to initstate won't | ||
220 | lose this information and will be able to restart with setstate. | ||
221 | Note: The first thing we do is save the current state, if any, just like | ||
222 | setstate so that it doesn't matter when initstate is called. | ||
223 | Returns 0 on success, non-zero on failure. */ | ||
224 | int initstate_r (unsigned int seed, | ||
225 | char *arg_state, | ||
226 | size_t n, | ||
227 | struct random_data *buf) | ||
228 | { | ||
229 | if (buf == NULL) | ||
230 | goto fail; | ||
231 | |||
232 | int32_t *old_state = buf->state; | ||
233 | if (old_state != NULL) | ||
234 | { | ||
235 | int old_type = buf->rand_type; | ||
236 | if (old_type == TYPE_0) | ||
237 | old_state[-1] = TYPE_0; | ||
238 | else | ||
239 | old_state[-1] = (MAX_TYPES * (buf->rptr - old_state)) + old_type; | ||
240 | } | ||
241 | |||
242 | int type; | ||
243 | if (n >= BREAK_3) | ||
244 | type = n < BREAK_4 ? TYPE_3 : TYPE_4; | ||
245 | else if (n < BREAK_1) | ||
246 | { | ||
247 | if (n < BREAK_0) | ||
248 | goto fail; | ||
249 | |||
250 | type = TYPE_0; | ||
251 | } | ||
252 | else | ||
253 | type = n < BREAK_2 ? TYPE_1 : TYPE_2; | ||
254 | |||
255 | int degree = random_poly_info.degrees[type]; | ||
256 | int separation = random_poly_info.seps[type]; | ||
257 | |||
258 | buf->rand_type = type; | ||
259 | buf->rand_sep = separation; | ||
260 | buf->rand_deg = degree; | ||
261 | int32_t *state = &((int32_t *) arg_state)[1]; /* First location. */ | ||
262 | /* Must set END_PTR before srandom. */ | ||
263 | buf->end_ptr = &state[degree]; | ||
264 | |||
265 | buf->state = state; | ||
266 | |||
267 | srandom_r (seed, buf); | ||
268 | |||
269 | state[-1] = TYPE_0; | ||
270 | if (type != TYPE_0) | ||
271 | state[-1] = (buf->rptr - state) * MAX_TYPES + type; | ||
272 | |||
273 | return 0; | ||
274 | |||
275 | fail: | ||
276 | __set_errno (EINVAL); | ||
277 | return -1; | ||
278 | } | ||
279 | |||
280 | /* Restore the state from the given state array. | ||
281 | Note: It is important that we also remember the locations of the pointers | ||
282 | in the current state information, and restore the locations of the pointers | ||
283 | from the old state information. This is done by multiplexing the pointer | ||
284 | location into the zeroth word of the state information. Note that due | ||
285 | to the order in which things are done, it is OK to call setstate with the | ||
286 | same state as the current state | ||
287 | Returns 0 on success, non-zero on failure. */ | ||
288 | int setstate_r (char *arg_state, struct random_data *buf) | ||
289 | { | ||
290 | int32_t *new_state = 1 + (int32_t *) arg_state; | ||
291 | int type; | ||
292 | int old_type; | ||
293 | int32_t *old_state; | ||
294 | int degree; | ||
295 | int separation; | ||
296 | |||
297 | if (arg_state == NULL || buf == NULL) | ||
298 | goto fail; | ||
299 | |||
300 | old_type = buf->rand_type; | ||
301 | old_state = buf->state; | ||
302 | if (old_type == TYPE_0) | ||
303 | old_state[-1] = TYPE_0; | ||
304 | else | ||
305 | old_state[-1] = (MAX_TYPES * (buf->rptr - old_state)) + old_type; | ||
306 | |||
307 | type = new_state[-1] % MAX_TYPES; | ||
308 | if (type < TYPE_0 || type > TYPE_4) | ||
309 | goto fail; | ||
310 | |||
311 | buf->rand_deg = degree = random_poly_info.degrees[type]; | ||
312 | buf->rand_sep = separation = random_poly_info.seps[type]; | ||
313 | buf->rand_type = type; | ||
314 | |||
315 | if (type != TYPE_0) | ||
316 | { | ||
317 | int rear = new_state[-1] / MAX_TYPES; | ||
318 | buf->rptr = &new_state[rear]; | ||
319 | buf->fptr = &new_state[(rear + separation) % degree]; | ||
320 | } | ||
321 | buf->state = new_state; | ||
322 | /* Set end_ptr too. */ | ||
323 | buf->end_ptr = &new_state[degree]; | ||
324 | |||
325 | return 0; | ||
326 | |||
327 | fail: | ||
328 | __set_errno (EINVAL); | ||
329 | return -1; | ||
330 | } | ||
331 | |||
332 | /* If we are using the trivial TYPE_0 R.N.G., just do the old linear | ||
333 | congruential bit. Otherwise, we do our fancy trinomial stuff, which is the | ||
334 | same in all the other cases due to all the global variables that have been | ||
335 | set up. The basic operation is to add the number at the rear pointer into | ||
336 | the one at the front pointer. Then both pointers are advanced to the next | ||
337 | location cyclically in the table. The value returned is the sum generated, | ||
338 | reduced to 31 bits by throwing away the "least random" low bit. | ||
339 | Note: The code takes advantage of the fact that both the front and | ||
340 | rear pointers can't wrap on the same call by not testing the rear | ||
341 | pointer if the front one has wrapped. Returns a 31-bit random number. */ | ||
342 | |||
343 | int random_r (struct random_data *buf, int32_t *result) | ||
344 | { | ||
345 | int32_t *state; | ||
346 | |||
347 | if (buf == NULL || result == NULL) | ||
348 | goto fail; | ||
349 | |||
350 | state = buf->state; | ||
351 | |||
352 | if (buf->rand_type == TYPE_0) | ||
353 | { | ||
354 | int32_t val = state[0]; | ||
355 | val = ((state[0] * 1103515245) + 12345) & 0x7fffffff; | ||
356 | state[0] = val; | ||
357 | *result = val; | ||
358 | } | ||
359 | else | ||
360 | { | ||
361 | int32_t *fptr = buf->fptr; | ||
362 | int32_t *rptr = buf->rptr; | ||
363 | int32_t *end_ptr = buf->end_ptr; | ||
364 | int32_t val; | ||
365 | |||
366 | val = *fptr += *rptr; | ||
367 | /* Chucking least random bit. */ | ||
368 | *result = (val >> 1) & 0x7fffffff; | ||
369 | ++fptr; | ||
370 | if (fptr >= end_ptr) | ||
371 | { | ||
372 | fptr = state; | ||
373 | ++rptr; | ||
374 | } | ||
375 | else | ||
376 | { | ||
377 | ++rptr; | ||
378 | if (rptr >= end_ptr) | ||
379 | rptr = state; | ||
380 | } | ||
381 | buf->fptr = fptr; | ||
382 | buf->rptr = rptr; | ||
383 | } | ||
384 | return 0; | ||
385 | |||
386 | fail: | ||
387 | __set_errno (EINVAL); | ||
388 | return -1; | ||
389 | } | ||
390 | |||
391 | int windows_random_r; | ||
392 | #endif \ No newline at end of file | ||
diff --git a/packages/base/src/Internal/C/windows_random_r.h b/packages/base/src/Internal/C/windows_random_r.h new file mode 100644 index 0000000..4f4d32f --- /dev/null +++ b/packages/base/src/Internal/C/windows_random_r.h | |||
@@ -0,0 +1,35 @@ | |||
1 | #if defined(_WIN32) || defined(WIN32) | ||
2 | |||
3 | #ifndef WINDOWS_RANDOM_R_H | ||
4 | #define WINDOWS_RANDOM_R_H | ||
5 | |||
6 | #include <stddef.h> | ||
7 | |||
8 | #define EINVAL 23 | ||
9 | |||
10 | typedef int int32_t; | ||
11 | |||
12 | struct random_data | ||
13 | { | ||
14 | int *fptr; /* Front pointer. */ | ||
15 | int *rptr; /* Rear pointer. */ | ||
16 | int *state; /* Array of state values. */ | ||
17 | int rand_type; /* Type of random number generator. */ | ||
18 | int rand_deg; /* Degree of random number generator. */ | ||
19 | int rand_sep; /* Distance between front and rear. */ | ||
20 | int *end_ptr; /* Pointer behind state table. */ | ||
21 | }; | ||
22 | |||
23 | extern int windows_random_r_errno; | ||
24 | |||
25 | static inline void __set_errno(int err) { windows_random_r_errno = err; }; | ||
26 | |||
27 | int random_r (struct random_data *buf, int32_t *result); | ||
28 | |||
29 | int initstate_r (unsigned int seed, | ||
30 | char *arg_state, | ||
31 | /*unsigned long*/ size_t n, | ||
32 | struct random_data *buf); | ||
33 | |||
34 | #endif | ||
35 | #endif \ No newline at end of file | ||