1 /*
2 * zrand - subtractive 100 shuffle generator
3 *
4 * Copyright (C) 1999-2007,2021 Landon Curt Noll
5 *
6 * Calc is open software; you can redistribute it and/or modify it under
7 * the terms of the version 2.1 of the GNU Lesser General Public License
8 * as published by the Free Software Foundation.
9 *
10 * Calc is distributed in the hope that it will be useful, but WITHOUT
11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General
13 * Public License for more details.
14 *
15 * A copy of version 2.1 of the GNU Lesser General Public License is
16 * distributed with calc under the filename COPYING-LGPL. You should have
17 * received a copy with calc; if not, write to Free Software Foundation, Inc.
18 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19 *
20 * Under source code control: 1995/01/07 09:45:25
21 * File existed as early as: 1994
22 *
23 * chongo <was here> /\oo/\ http://www.isthe.com/chongo/
24 * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/
25 */
26
27 /*
28 * AN OVERVIEW OF THE FUNCTIONS:
29 *
30 * This module contains an Subtractive 100 shuffle generator wrapped inside
31 * of a shuffle generator.
32 *
33 * We refer to this generator as the s100 generator.
34 *
35 * rand - s100 shuffle generator
36 * srand - seed the s100 shuffle generator
37 *
38 * This generator has two distinct parts, the s100 generator
39 * and the shuffle generator.
40 *
41 * The subtractive 100 generator is described in Knuth's "The Art of
42 * Computer Programming - Seminumerical Algorithms", Vol 2, 3rd edition
43 * (1998), Section 3.6, page 186, formula (2).
44 *
45 * The "use only the first 100 our of every 1009" is described in
46 * Knuth's "The Art of Computer Programming - Seminumerical Algorithms",
47 * Vol 2, 3rd edition (1998), Section 3.6, page 188".
48 *
49 * The period and other properties of this generator make it very
50 * useful to 'seed' other generators.
51 *
52 * The shuffle generator is described in Knuth's "The Art of Computer
53 * Programming - Seminumerical Algorithms", Vol 2, 3rd edition (1998),
54 * Section 3.2.2, page 34, Algorithm B.
55 *
56 * The shuffle generator is fast and serves as a fairly good standard
57 * pseudo-random generator. If you need a fast generator and do not
58 * need a cryptographically strong one, this generator is likely to do
59 * the job.
60 *
61 * The shuffle generator is feed values by the subtractive 100 process.
62 *
63 ******************************************************************************
64 *
65 * GOALS:
66 *
67 * The goals of this package are:
68 *
69 * all magic numbers are explained
70 *
71 * I distrust systems with constants (magic numbers) and tables
72 * that have no justification. I believe that I have
73 * done my best to justify all of the magic numbers used.
74 *
75 * full documentation
76 *
77 * You have this source file, plus background publications,
78 * what more could you ask?
79 *
80 * large selection of seeds
81 *
82 * Seeds are not limited to a small number of bits. A seed
83 * may be of any size.
84 *
85 * the strength of the generators may be tuned to meet the need
86 *
87 * By using the appropriate seed and other arguments, one may
88 * increase the strength of the generator to suit the need of
89 * the application. One does not have just a few levels.
90 *
91 * Even though I have done my best to implement a good system, you still
92 * must use these routines your own risk.
93 *
94 * Share and enjoy! :-)
95 */
96
97 /*
98 * ON THE GENERATORS:
99 *
100 * The subtractive 100 generator has a good period, and is fast. It is
101 * reasonable as generators go, though there are better ones available.
102 * The shuffle generator has a very good period, and is fast. It is
103 * fairly good as generators go, particularly when it is feed reasonably
104 * random numbers. Because of this, we use feed values from the subtractive
105 * 100 process into the shuffle generator.
106 *
107 * The s100 generator uses 2 tables:
108 *
109 * subtractive table - 100 entries of 64 bits used by the subtractive 100
110 * part of the s100 generator
111 *
112 * shuffle table - 256 entries of 64 bits used by the shuffle
113 * part of the s100 generator and feed by the
114 * subtractive table.
115 *
116 * Casual direct use of the shuffle generator may be acceptable. If one
117 * desires cryptographically strong random numbers, or if one is paranoid,
118 * one should use the Blum generator instead (see zrandom.c).
119 *
120 * The s100 generator as the following calc interfaces:
121 *
122 * rand(min,beyond) (where min < beyond)
123 *
124 * Print an s100 generator random value over interval [a,b).
125 *
126 * rand()
127 *
128 * Same as rand(0, 2^64). Print 64 bits.
129 *
130 * rand(lim) (where 0 > lim)
131 *
132 * Same as rand(0, lim).
133 *
134 * randbit(x) (where x > 0)
135 *
136 * Same as rand(0, 2^x). Print x bits.
137 *
138 * randbit(skip) (where skip < 0)
139 *
140 * Skip random bits and return the bit skip count (-skip).
141 */
142
143 /*
144 * INITIALIZATION AND SEEDS:
145 *
146 * All generators come already seeded with precomputed initial constants.
147 * Thus, it is not required to seed a generator before using it.
148 *
149 * The s100 generator may be initialized and seeded via srand().
150 *
151 * Using a seed of '0' will reload generators with their initial states.
152 *
153 * srand(0) restore subtractive 100 generator to the initial state
154 *
155 * The above single arg calls are fairly fast.
156 *
157 * Optimal seed range for the s100 generator:
158 *
159 * There is no limit on the size of a seed. On the other hand,
160 * extremely large seeds require large tables and long seed times.
161 * Using a seed in the range of [2^64, 2^64 * 100!) should be
162 * sufficient for most purposes. An easy way to stay within this
163 * range to to use seeds that are between 21 and 178 digits, or
164 * 64 to 588 bits long.
165 *
166 * To help make the generator produced by seed S, significantly
167 * different from S+1, seeds are scrambled prior to use. The
168 * function randreseed64() maps [0,2^64) into [0,2^64) in a 1-to-1
169 * and onto fashion.
170 *
171 * The purpose of the randreseed64() is not to add security. It
172 * simply helps remove the human perception of the relationship
173 * between the seed and the production of the generator.
174 *
175 * The randreseed64() process does not reduce the security of the
176 * generators. Every seed is converted into a different unique seed.
177 * No seed is ignored or favored.
178 *
179 ******************************************************************************
180 *
181 * srand(seed)
182 *
183 * Seed the s100 generator.
184 *
185 * seed != 0:
186 * ---------
187 * Any buffered random bits are flushed. The subtractive table is loaded
188 * with the default subtractive table. The low order 64 bits of seed is
189 * xor-ed against each table value. The subtractive table is shuffled
190 * according to seed/2^64.
191 *
192 * The following calc code produces the same effect:
193 *
194 * (* reload default subtractive table xor-ed with low 64 seed bits *)
195 * seed_xor = seed & ((1<<64)-1);
196 * for (i=0; i < 100; ++i) {
197 * subtractive[i] = xor(default_subtractive[i], seed_xor);
198 * }
199 *
200 * (* shuffle the subtractive table *)
201 * seed >>= 64;
202 * for (i=100; seed > 0 && i > 0; --i) {
203 * quomod(seed, i+1, seed, j);
204 * swap(subtractive[i], subtractive[j]);
205 * }
206 *
207 * Seed must be >= 0. All seed values < 0 are reserved for future use.
208 *
209 * The subtractive 100 pointers are reset to subtractive[36] and
210 * subtractive[99]. Last the shuffle table is loaded with successive
211 * values from the subtractive 100 generator.
212 *
213 * seed == 0:
214 * ---------
215 * Restore the initial state and modulus of the s100 generator.
216 * After this call, the s100 generator is restored to its initial
217 * state after calc started.
218 *
219 * The subtractive 100 pointers are reset to subtractive[36] and
220 * subtractive[99]. Last the shuffle table is loaded with successive
221 * values from the subtractive 100 generator.
222 *
223 ******************************************************************************
224 *
225 * srand(mat100)
226 *
227 * Seed the s100 generator.
228 *
229 * Any buffered random bits are flushed. The subtractive table with the
230 * first 100 entries of the array mat100, mod 2^64.
231 *
232 * The subtractive 100 pointers are reset to subtractive[36] and
233 * subtractive[99]. Last the shuffle table is loaded with successive
234 * values from the subtractive 100 generator.
235 *
236 ******************************************************************************
237 *
238 * srand()
239 *
240 * Return current s100 generator state. This call does not alter
241 * the generator state.
242 *
243 ******************************************************************************
244 *
245 * srand(state)
246 *
247 * Restore the s100 state and return the previous state. Note that
248 * the argument state is a rand state value (isrand(state) is true).
249 * Any internally buffered random bits are restored.
250 *
251 * The states of the s100 generators can be saved by calling the seed
252 * function with no arguments, and later restored by calling the seed
253 * functions with that same return value.
254 *
255 * rand_state = srand();
256 * ... generate random bits ...
257 * prev_rand_state = srand(rand_state);
258 * ... generate the same random bits ...
259 * srand() == prev_rand_state; (* is true *)
260 *
261 * Saving the state just after seeding a generator and restoring it later
262 * as a very fast way to reseed a generator.
263 */
264
265 /*
266 * TRUTH IN ADVERTISING:
267 *
268 * A "truth in advertising" issue is the use of the term
269 * 'pseudo-random'. All deterministic generators are pseudo random.
270 * This is opposed to true random generators based on some special
271 * physical device.
272 *
273 * A final "truth in advertising" issue deals with how the magic numbers
274 * found in this file were generated. Detains can be found in the
275 * various functions, while a overview can be found in the "SOURCE FOR
276 * MAGIC NUMBERS" section below.
277 */
278
279 /*
280 * SOURCE OF MAGIC NUMBERS:
281 *
282 * Most of the magic constants used on this file ultimately are
283 * based on LavaRnd. LavaRnd produced them via a cryptographic
284 * of the digitization of chaotic system that consisted of a noisy
285 * digital camera and 6 Lava Lite(R) lamps.
286 *
287 * BTW: Lava Lite(R) is a trademark of Haggerty Enterprises, Inc.
288 *
289 * The first 100 groups of 64 bit bits were used to initialize init_s100.slot.
290 *
291 * The function, randreseed64(), uses 2 primes to scramble 64 bits
292 * into 64 bits. These primes were also extracted from the Rand
293 * Book of Random Numbers. See randreseed64() for details.
294 *
295 * The shuffle table size is longer than the 100 entries recommended
296 * by Knuth. We use a power of 2 shuffle table length so that the
297 * shuffle process can select a table entry from a new subtractive 100
298 * value by extracting its low order bits. The value 256 is convenient
299 * in that it is the size of a byte which allows for easy extraction.
300 *
301 * We use the upper byte of the subtractive 100 value to select the shuffle
302 * table entry because it allows all of 64 bits to play a part in the
303 * entry selection. If we were to select a lower 8 bits in the 64 bit
304 * value, carries that propagate above our 8 bits would not impact the
305 * s100 generator output.
306 *
307 ******************************************************************************
308 *
309 * FOR THE PARANOID:
310 *
311 * The truly paranoid might suggest that my claims in the MAGIC NUMBERS
312 * section are a lie intended to entrap people. Well they are not, but
313 * if you that paranoid why would you use a non-cryptographically strong
314 * pseudo-random number generator in the first place? You would be
315 * better off using the random() builtin function.
316 *
317 * The two constants that were picked from the Rand Book of Random Numbers
318 * The random numbers from the Rand Book of Random Numbers can be
319 * verified by anyone who obtains the book. As these numbers were
320 * created before I (Landon Curt Noll) was born (you can look up
321 * my birth record if you want), I claim to have no possible influence
322 * on their generation.
323 *
324 * There is a very slight chance that the electronic copy of the
325 * Rand Book of Random Numbers that I was given access to differs
326 * from the printed text. I am willing to provide access to this
327 * electronic copy should anyone wants to compare it to the printed text.
328 *
329 * When using the s100 generator, one may select your own 100 subtractive
330 * values by calling:
331 *
332 * srand(mat100)
333 *
334 * and avoid using my magic numbers. The randreseed64 process is NOT
335 * applied to the matrix values. Of course, you must pick good subtractive
336 * 100 values yourself!
337 *
338 * One might object to the complexity of the seed scramble/mapping
339 * via the randreseed64() function. The randreseed64() function maps:
340 *
341 * 0 ==> 0
342 * 10239951819489363767 ==> 1363042948800878693
343 *
344 * so that srand(0) does the default action and randreseed64() remains
345 * an 1-to-1 and onto map. Thus calling srand(0) with the randreseed64()
346 * process would be the same as calling srand(4967126403401436567) without
347 * it. No extra security is gained or reduced by using the randreseed64()
348 * process. The meaning of seeds are exchanged, but not lost or favored
349 * (used by more than one input seed).
350 */
351
352 #include <stdio.h>
353
354 #include "alloc.h"
355 #include "zrand.h"
356 #include "have_const.h"
357 #include "have_unused.h"
358
359
360 #include "banned.h" /* include after system header <> includes */
361
362
363 /*
364 * default s100 generator state
365 *
366 * This is the state of the s100 generator after initialization, or srand(0),
367 * or srand(0) is called. The init_s100 value is never changed, only copied.
368 */
369 STATIC CONST RAND init_s100 = {
370 1, /* seeded */
371 0, /* no buffered bits */
372 #if FULL_BITS == SBITS /* buffer */
373 {0},
374 #elif 2*FULL_BITS == SBITS
375 {0, 0},
376 #else
377 /\../\ BASEB is assumed to be 16 or 32 /\../\ !!!
378 #endif
379 INIT_J, /* j */
380 INIT_K, /* k */
381 RAND_CONSEQ_USE, /* use this many before skipping values */
382 #if FULL_BITS == SBITS
383 { /* subtractive 100 table */
384 (FULL)U(0xc8c0370c7db7dc19), (FULL)U(0x738e33b940a06fbb),
385 (FULL)U(0x481abb76a859ed2b), (FULL)U(0x74106bb39ccdccb5),
386 (FULL)U(0x05a8eeb5c3173bfc), (FULL)U(0xefd5100d5a02e577),
387 (FULL)U(0xa69271f74030b24a), (FULL)U(0x641282fc16fe22c5),
388 (FULL)U(0x7aa7267c40438da3), (FULL)U(0x1fdf4abdc2d878d1),
389 (FULL)U(0xd9899e7a95702379), (FULL)U(0x5ea8e217d02d7f08),
390 (FULL)U(0x770587fe4d47a353), (FULL)U(0xde7d1bdd0a33a2b8),
391 (FULL)U(0x4378c3c5900e7c45), (FULL)U(0x77c9447819a514f9),
392 (FULL)U(0xfc5edb22843d1d32), (FULL)U(0x4fc42ce5e8ee5e6e),
393 (FULL)U(0xc938713c8488013e), (FULL)U(0x6a318f0320ab0cac),
394 (FULL)U(0x73e6d1a3ffc8bff3), (FULL)U(0x0cd3232a8ca96aa7),
395 (FULL)U(0x605c8036905f770d), (FULL)U(0x4d037b008b8d04a2),
396 (FULL)U(0x1ed81965cb277294), (FULL)U(0x408d9c477a254ff3),
397 (FULL)U(0x8b68587ae26c7377), (FULL)U(0xcff191a48a48832f),
398 (FULL)U(0x12d3df1d8aeb6fe6), (FULL)U(0xb2bf907e1feda37a),
399 (FULL)U(0x4e5f77193bb5f39f), (FULL)U(0x33ebcf6f8f5d1581),
400 (FULL)U(0x203c8e48d33654eb), (FULL)U(0x68d3656ef19c8a4e),
401 (FULL)U(0x3ec20b04986eb2af), (FULL)U(0x5d73a03b062c3841),
402 (FULL)U(0x836ce7095d4e49eb), (FULL)U(0x2310bc40c3f49221),
403 (FULL)U(0x3868ee48a6d0cbf6), (FULL)U(0x67578aa64a43deb1),
404 (FULL)U(0x6e3426c1150dfc26), (FULL)U(0xc541ccaa3131be30),
405 (FULL)U(0xf7e57432cec7aab2), (FULL)U(0x2b35de998cb3c873),
406 (FULL)U(0x7b9f77648663a5d7), (FULL)U(0x23b00e6aa771e5a6),
407 (FULL)U(0x859c775ca9985d05), (FULL)U(0x99636ea16b692f1f),
408 (FULL)U(0x8700ac703730800d), (FULL)U(0x461425024298a753),
409 (FULL)U(0xea4a411b809e955f), (FULL)U(0x3119ad4033709dfb),
410 (FULL)U(0xb76a6c6e5f01cb7c), (FULL)U(0x6109dc8a15984eaf),
411 (FULL)U(0x5d686db9a5ca9505), (FULL)U(0x8e80d7613b7e6add),
412 (FULL)U(0x79cbd718de6f6fd3), (FULL)U(0x40e9cd151da0f699),
413 (FULL)U(0xe82158bab24f312d), (FULL)U(0x79a4c927f5e5c36b),
414 (FULL)U(0xc25247c9a0039333), (FULL)U(0x936871161766d81d),
415 (FULL)U(0x3c6a03b4a6741327), (FULL)U(0xc8a7b6e8c002f29a),
416 (FULL)U(0x0e2a67c67bbd5ea3), (FULL)U(0x0929042d441eabc1),
417 (FULL)U(0x7dbe232a25e82085), (FULL)U(0x8cfb26e544fbac3d),
418 (FULL)U(0x8e40384d388ab983), (FULL)U(0x48dc1230554632f8),
419 (FULL)U(0xab405048ab492397), (FULL)U(0x21c9e2f5a118e387),
420 (FULL)U(0x484d1a8c343b61b5), (FULL)U(0xd49e3decab256f26),
421 (FULL)U(0xe615c7fd78f2d2e3), (FULL)U(0x8442cc33ce6cc2ed),
422 (FULL)U(0x0a3b93d844d4bbf6), (FULL)U(0x2d7e4efe9301de77),
423 (FULL)U(0x33711b76d8790d8a), (FULL)U(0xc07dc30e44df77e7),
424 (FULL)U(0xb9132ed09ddd508f), (FULL)U(0x45d06cf8c6fb43cc),
425 (FULL)U(0x22bed18ad585dd7b), (FULL)U(0x61c6cced10799ffa),
426 (FULL)U(0xd7f2393be4bd9aa9), (FULL)U(0x706753fbcfd55094),
427 (FULL)U(0xf65a6713ede6e446), (FULL)U(0x8bf6dfae47c0d5c3),
428 (FULL)U(0xfb4dfc179f7927d6), (FULL)U(0x12ebbc16e212c297),
429 (FULL)U(0x43c71283a00a954c), (FULL)U(0x8957087ae7bd40a5),
430 (FULL)U(0xb0859d7108344837), (FULL)U(0xfbf4b9a3aeb313f5),
431 (FULL)U(0x5e66e5bece81823a), (FULL)U(0x09a11c6e58ad6da1),
432 (FULL)U(0xc76f4316c608054f), (FULL)U(0xb582136146084099),
433 (FULL)U(0x4210008f17a725ed), (FULL)U(0xe5ff8912d347c481)
434 },
435 { /* shuffle table */
436 (FULL)U(0x69a2296cec8abd57), (FULL)U(0x867e186999a6df81),
437 (FULL)U(0xc05ab96bd849a48a), (FULL)U(0x7eb3ce0cfa00554b),
438 (FULL)U(0x520d01f65a5a9acd), (FULL)U(0xd4ef1e3336022d81),
439 (FULL)U(0xaf44772bc6f84f70), (FULL)U(0x647e85a6a7c55173),
440 (FULL)U(0x26746cf1959df8d1), (FULL)U(0x98681a904db28abd),
441 (FULL)U(0xb146c969744c5cd2), (FULL)U(0x8ce69d1f706f88c2),
442 (FULL)U(0xfd12eac421b4a748), (FULL)U(0xf12e70fe2710eea5),
443 (FULL)U(0x0b8f78055901f2b5), (FULL)U(0x48860a764f2c115e),
444 (FULL)U(0x0edf6d2a30767e2c), (FULL)U(0x8a6d7dc5fce2713b),
445 (FULL)U(0x46a362ea4e0e2346), (FULL)U(0x6c369a0a359f5aa7),
446 (FULL)U(0xdfca81fe41def54e), (FULL)U(0x4b73381996c2bc4e),
447 (FULL)U(0x659e8b996f3f14f9), (FULL)U(0x8b97b93493d47e6f),
448 (FULL)U(0xa73a8704dfa10a55), (FULL)U(0x8d9eafe9b06503da),
449 (FULL)U(0x2556fb88f32336b0), (FULL)U(0xe71e9f751002a161),
450 (FULL)U(0x27a7be6e200af907), (FULL)U(0x1b9b734ed028e9a3),
451 (FULL)U(0x950cfeed4c0be0d3), (FULL)U(0xf4c416942536d275),
452 (FULL)U(0xf05a58e85687b76e), (FULL)U(0xba53ac0171a62d54),
453 (FULL)U(0x4b14cbcb285adc96), (FULL)U(0xfdf66eddb00a5557),
454 (FULL)U(0xbb43d58d185b6ea1), (FULL)U(0x905db9cdf355c9a6),
455 (FULL)U(0xfc3a07fc04429c8a), (FULL)U(0x65d7e365aa3a4f7e),
456 (FULL)U(0x2d284c18b243ac65), (FULL)U(0x72fba65d44e417fd),
457 (FULL)U(0x422d50b45c934805), (FULL)U(0xb62a6053d1587441),
458 (FULL)U(0xa5e71ce96f7ae035), (FULL)U(0x93abca2e595c8dd8),
459 (FULL)U(0x534231afe39afad5), (FULL)U(0x08d26cac12eaad56),
460 (FULL)U(0xec18bf8d7fb1b1c2), (FULL)U(0x3d28ea16faf6f09b),
461 (FULL)U(0xea357a7816697dd6), (FULL)U(0x51471ea1420f3f51),
462 (FULL)U(0x5e051aeb7f8946b4), (FULL)U(0x881be0970cf0524c),
463 (FULL)U(0xd558b25b1b31489e), (FULL)U(0x707d1a943a8b065c),
464 (FULL)U(0x37017e66568ff836), (FULL)U(0xb9cd627c24c2f747),
465 (FULL)U(0x1485549ffb1d9ff6), (FULL)U(0x308d32d9bdf2dc6f),
466 (FULL)U(0x4d4142cad543818a), (FULL)U(0x5d9c7aee87ebba43),
467 (FULL)U(0x81c5bdd8e17adb2f), (FULL)U(0x3dc9752ec8d8677a),
468 (FULL)U(0x66b086e6c34e4212), (FULL)U(0x3af7a90dc62b25e3),
469 (FULL)U(0xf8349f7935539315), (FULL)U(0x6bcfd9d5a22917f0),
470 (FULL)U(0x8639bb765f5ee517), (FULL)U(0xd3c5e3698095b092),
471 (FULL)U(0x8a33851e7eb44748), (FULL)U(0x5e29d443ea54bbcf),
472 (FULL)U(0x0f84651f4d59a834), (FULL)U(0x85040beaf1a5f951),
473 (FULL)U(0x3dba1c7498002078), (FULL)U(0x5d70712bf0b2cc15),
474 (FULL)U(0xfa3af8ebcce8e5a7), (FULL)U(0xfb3e223704bba57d),
475 (FULL)U(0x5d3b87858a950434), (FULL)U(0xce3112bdba3f8dcf),
476 (FULL)U(0x44904f55860d3051), (FULL)U(0xcec8fed44ed3e98b),
477 (FULL)U(0x4581698d25d01ea4), (FULL)U(0x11eb68289a9548e0),
478 (FULL)U(0x796cb4c6e911fac8), (FULL)U(0x2164cf26b5fd813e),
479 (FULL)U(0x4ac8e0f5d5de640f), (FULL)U(0xe9e757d78802ab4e),
480 (FULL)U(0x3c97de26f49dfcbd), (FULL)U(0xc604881b6ee6dbe6),
481 (FULL)U(0xa7c22a6e57d6154e), (FULL)U(0x234e2370877b3cc7),
482 (FULL)U(0xc0bdb72bdf1f8358), (FULL)U(0x6522e0fca95b7b55),
483 (FULL)U(0xba174c9022344162), (FULL)U(0x712c9b2d75d48867),
484 (FULL)U(0x240f7e92e59f3700), (FULL)U(0xe83cc2d4ad95d763),
485 (FULL)U(0x8509445a4336d717), (FULL)U(0xf1e572c5dfff1804),
486 (FULL)U(0xed10eb5d623232dd), (FULL)U(0x9205ea1bd4f957e8),
487 (FULL)U(0x4973a54f2ff062f5), (FULL)U(0x26b018f1e8c48cd5),
488 (FULL)U(0x56908401d1c7ed9f), (FULL)U(0x2e48937bdf89a247),
489 (FULL)U(0x9d53069b2be47129), (FULL)U(0x98069e3bc048a2b0),
490 (FULL)U(0xf25b7d651cd83f93), (FULL)U(0x2b004e6ce6f886c8),
491 (FULL)U(0xf618442a5c635935), (FULL)U(0xa502ab5c7198e052),
492 (FULL)U(0xc14241a4a6c41b0b), (FULL)U(0x720e845a7db9b18e),
493 (FULL)U(0x2abb13e94b713918), (FULL)U(0x90fc0c207f52467d),
494 (FULL)U(0x799c8ccd7868d348), (FULL)U(0xf4817ced912a0ea4),
495 (FULL)U(0xd68c0f4cc4903a57), (FULL)U(0xa3171f29e2b7934c),
496 (FULL)U(0xb1158baa0b4ccc22), (FULL)U(0xf5d8555349a29eda),
497 (FULL)U(0x59d1a078959442ef), (FULL)U(0xdb9b4a96a67fd518),
498 (FULL)U(0xcc7ca9eed2870636), (FULL)U(0x548f021cecf59920),
499 (FULL)U(0x25b7f4b6571bc8c5), (FULL)U(0x4fa527473a44f536),
500 (FULL)U(0xb246845fdf0ebdc2), (FULL)U(0xdd8d68ae42058793),
501 (FULL)U(0x3ba133289f6c39fb), (FULL)U(0x8bfdfbf37b6b42af),
502 (FULL)U(0xfb34c5ca7fb2b3b0), (FULL)U(0x2345dcecd428e32a),
503 (FULL)U(0x6891e850ad42b63e), (FULL)U(0x930642c8362c1381),
504 (FULL)U(0x13871e9b1886aff5), (FULL)U(0xd0cf2407482bda55),
505 (FULL)U(0x125b5fc95069bc31), (FULL)U(0x9b71d0a9f07dfa5d),
506 (FULL)U(0x55c044cc6712e524), (FULL)U(0xf0377358bb601978),
507 (FULL)U(0x152ad5f87fa51e8b), (FULL)U(0xe5ebf4789fcdd9af),
508 (FULL)U(0x3d78e18c66ebce7e), (FULL)U(0x8246db72f36aa83f),
509 (FULL)U(0xcc6ddc6d2c64c0a3), (FULL)U(0xa758d6870d91851e),
510 (FULL)U(0x24b20a6f9488ee36), (FULL)U(0xbe11ccdf09798197),
511 (FULL)U(0x11aca01599c1f4e3), (FULL)U(0x40e89e366437ac05),
512 (FULL)U(0xc8bfc7625af675f8), (FULL)U(0x6367c578b577e759),
513 (FULL)U(0x00380346615f0b74), (FULL)U(0xee964cc48de07d81),
514 (FULL)U(0x17f6ac16859d9261), (FULL)U(0x092f4a173a6e2f6c),
515 (FULL)U(0x79981a3db9024b95), (FULL)U(0x36db166004f7f540),
516 (FULL)U(0xc36252cf65a2f1c8), (FULL)U(0x705b6fde124c9bd2),
517 (FULL)U(0x31e58dda85db40ce), (FULL)U(0x6342b1a59f5e8d6d),
518 (FULL)U(0x5c2c67d0bd6d1d4d), (FULL)U(0x1fe5b46fba7e069d),
519 (FULL)U(0x21c46c6cac72e13c), (FULL)U(0xb80c5fd59eb8f52a),
520 (FULL)U(0x56c3aebfa74c92bc), (FULL)U(0xc1aff1fcbf8c4196),
521 (FULL)U(0x2b1df645754ad208), (FULL)U(0x5c734600d46eeb50),
522 (FULL)U(0xe0ff1b126a70a765), (FULL)U(0xd54164977a94547c),
523 (FULL)U(0x67b59d7c4ea35206), (FULL)U(0x53be7146779203b4),
524 (FULL)U(0x6b589fe5414026b8), (FULL)U(0x9e81016c3083bfee),
525 (FULL)U(0xb23526b93b4b7671), (FULL)U(0x4fa9ffb17ee300ba),
526 (FULL)U(0x6217e212ad05fb21), (FULL)U(0xf5b3fcd3b294e6c2),
527 (FULL)U(0xac040bbe216beb2a), (FULL)U(0x1f8d8a5471d0e78c),
528 (FULL)U(0xb6d15b419cfec96b), (FULL)U(0xc5477845d0508c78),
529 (FULL)U(0x5b486e81b4bba621), (FULL)U(0x90c35c94ef4c4121),
530 (FULL)U(0xefce7346f6a6bc55), (FULL)U(0xa27828d925bdb9bb),
531 (FULL)U(0xe3a53095a1f0b205), (FULL)U(0x1bfa6093d9f208ab),
532 (FULL)U(0xfb078f6a6842cdf4), (FULL)U(0x07806d7297133a38),
533 (FULL)U(0x2c6c901ba3ce9592), (FULL)U(0x1f0ab2cfebc1b789),
534 (FULL)U(0x2ce81415e2d03d5e), (FULL)U(0x7da45d5baa9f2417),
535 (FULL)U(0x3be4f76ddd800682), (FULL)U(0xdbf4e4a3364d72d3),
536 (FULL)U(0xb538cccf4fc59da5), (FULL)U(0xb0aa39d5487f66ec),
537 (FULL)U(0x2fd28dfd87927d3d), (FULL)U(0xd14e77f05900c6b1),
538 (FULL)U(0x2523fad25330c7b4), (FULL)U(0x991b5938d82368a4),
539 (FULL)U(0xb7c114432b9c1302), (FULL)U(0xdb842db61394b116),
540 (FULL)U(0x3641548d78ed26d8), (FULL)U(0x274fa8ef0a61dacf),
541 (FULL)U(0xa554ba63112df6f1), (FULL)U(0x7b7fe9856b50438d),
542 (FULL)U(0xc9fa0042bb63bbad), (FULL)U(0x3abf45d0e27f00da),
543 (FULL)U(0xd95faa159f87aabb), (FULL)U(0x4a95012e3488e7ae),
544 (FULL)U(0x1be2bdb90c642d04), (FULL)U(0x145c88818b4abf3e),
545 (FULL)U(0x7f9fb635544cf17f), (FULL)U(0xb8ab2f62cc78db70),
546 (FULL)U(0x8ee64bcdb4242f9a), (FULL)U(0xabd5285895dad129),
547 (FULL)U(0xbe722c2fccf31141), (FULL)U(0x7c330703575e26a9),
548 (FULL)U(0x45d3e3b3361b79e4), (FULL)U(0x241163a754b2e6a6),
549 (FULL)U(0x8f678d7df7cacb77), (FULL)U(0x988a68a483211d19),
550 (FULL)U(0x79599598ba7836f6), (FULL)U(0x4850c887eeda68bf),
551 (FULL)U(0xafa69a718052ce25), (FULL)U(0x8b21efc6bdd73573),
552 (FULL)U(0x89dbae18d0972493), (FULL)U(0x560776bf537d9454),
553 (FULL)U(0x3c009f78165310f2), (FULL)U(0xa36800210160c3af),
554 (FULL)U(0x3353ec3ca643bd40), (FULL)U(0x7e593f99911dab02),
555 (FULL)U(0x72d1ddd94f676e89), (FULL)U(0xfd18b8bd6b43c0ea),
556 (FULL)U(0x43cacef2ddbd697d), (FULL)U(0x2868a4d0acefe884),
557 (FULL)U(0x5f377b63a506f013), (FULL)U(0xeaa0975e05ca662b),
558 (FULL)U(0x3740e6b8eb433931), (FULL)U(0xce85df0008557948),
559 (FULL)U(0x784745fb547e33f9), (FULL)U(0x4a1fc5d4e5c6f598),
560 (FULL)U(0x85fa6fec768430a7), (FULL)U(0x990d0c24d2332a51),
561 (FULL)U(0x55245c2c33b676d5), (FULL)U(0xb1091519e2bcfa71),
562 (FULL)U(0x38521478d23a28d8), (FULL)U(0x9b794f899a573010),
563 (FULL)U(0x61d225e8699bb486), (FULL)U(0x21476d241c2158b0)
564 }
565 #elif 2*FULL_BITS == SBITS
566 { /* subtractive 100 table */
567 (FULL)0x7db7dc19,(FULL)0xc8c0370c,(FULL)0x40a06fbb,(FULL)0x738e33b9,
568 (FULL)0xa859ed2b,(FULL)0x481abb76,(FULL)0x9ccdccb5,(FULL)0x74106bb3,
569 (FULL)0xc3173bfc,(FULL)0x05a8eeb5,(FULL)0x5a02e577,(FULL)0xefd5100d,
570 (FULL)0x4030b24a,(FULL)0xa69271f7,(FULL)0x16fe22c5,(FULL)0x641282fc,
571 (FULL)0x40438da3,(FULL)0x7aa7267c,(FULL)0xc2d878d1,(FULL)0x1fdf4abd,
572 (FULL)0x95702379,(FULL)0xd9899e7a,(FULL)0xd02d7f08,(FULL)0x5ea8e217,
573 (FULL)0x4d47a353,(FULL)0x770587fe,(FULL)0x0a33a2b8,(FULL)0xde7d1bdd,
574 (FULL)0x900e7c45,(FULL)0x4378c3c5,(FULL)0x19a514f9,(FULL)0x77c94478,
575 (FULL)0x843d1d32,(FULL)0xfc5edb22,(FULL)0xe8ee5e6e,(FULL)0x4fc42ce5,
576 (FULL)0x8488013e,(FULL)0xc938713c,(FULL)0x20ab0cac,(FULL)0x6a318f03,
577 (FULL)0xffc8bff3,(FULL)0x73e6d1a3,(FULL)0x8ca96aa7,(FULL)0x0cd3232a,
578 (FULL)0x905f770d,(FULL)0x605c8036,(FULL)0x8b8d04a2,(FULL)0x4d037b00,
579 (FULL)0xcb277294,(FULL)0x1ed81965,(FULL)0x7a254ff3,(FULL)0x408d9c47,
580 (FULL)0xe26c7377,(FULL)0x8b68587a,(FULL)0x8a48832f,(FULL)0xcff191a4,
581 (FULL)0x8aeb6fe6,(FULL)0x12d3df1d,(FULL)0x1feda37a,(FULL)0xb2bf907e,
582 (FULL)0x3bb5f39f,(FULL)0x4e5f7719,(FULL)0x8f5d1581,(FULL)0x33ebcf6f,
583 (FULL)0xd33654eb,(FULL)0x203c8e48,(FULL)0xf19c8a4e,(FULL)0x68d3656e,
584 (FULL)0x986eb2af,(FULL)0x3ec20b04,(FULL)0x062c3841,(FULL)0x5d73a03b,
585 (FULL)0x5d4e49eb,(FULL)0x836ce709,(FULL)0xc3f49221,(FULL)0x2310bc40,
586 (FULL)0xa6d0cbf6,(FULL)0x3868ee48,(FULL)0x4a43deb1,(FULL)0x67578aa6,
587 (FULL)0x150dfc26,(FULL)0x6e3426c1,(FULL)0x3131be30,(FULL)0xc541ccaa,
588 (FULL)0xcec7aab2,(FULL)0xf7e57432,(FULL)0x8cb3c873,(FULL)0x2b35de99,
589 (FULL)0x8663a5d7,(FULL)0x7b9f7764,(FULL)0xa771e5a6,(FULL)0x23b00e6a,
590 (FULL)0xa9985d05,(FULL)0x859c775c,(FULL)0x6b692f1f,(FULL)0x99636ea1,
591 (FULL)0x3730800d,(FULL)0x8700ac70,(FULL)0x4298a753,(FULL)0x46142502,
592 (FULL)0x809e955f,(FULL)0xea4a411b,(FULL)0x33709dfb,(FULL)0x3119ad40,
593 (FULL)0x5f01cb7c,(FULL)0xb76a6c6e,(FULL)0x15984eaf,(FULL)0x6109dc8a,
594 (FULL)0xa5ca9505,(FULL)0x5d686db9,(FULL)0x3b7e6add,(FULL)0x8e80d761,
595 (FULL)0xde6f6fd3,(FULL)0x79cbd718,(FULL)0x1da0f699,(FULL)0x40e9cd15,
596 (FULL)0xb24f312d,(FULL)0xe82158ba,(FULL)0xf5e5c36b,(FULL)0x79a4c927,
597 (FULL)0xa0039333,(FULL)0xc25247c9,(FULL)0x1766d81d,(FULL)0x93687116,
598 (FULL)0xa6741327,(FULL)0x3c6a03b4,(FULL)0xc002f29a,(FULL)0xc8a7b6e8,
599 (FULL)0x7bbd5ea3,(FULL)0x0e2a67c6,(FULL)0x441eabc1,(FULL)0x0929042d,
600 (FULL)0x25e82085,(FULL)0x7dbe232a,(FULL)0x44fbac3d,(FULL)0x8cfb26e5,
601 (FULL)0x388ab983,(FULL)0x8e40384d,(FULL)0x554632f8,(FULL)0x48dc1230,
602 (FULL)0xab492397,(FULL)0xab405048,(FULL)0xa118e387,(FULL)0x21c9e2f5,
603 (FULL)0x343b61b5,(FULL)0x484d1a8c,(FULL)0xab256f26,(FULL)0xd49e3dec,
604 (FULL)0x78f2d2e3,(FULL)0xe615c7fd,(FULL)0xce6cc2ed,(FULL)0x8442cc33,
605 (FULL)0x44d4bbf6,(FULL)0x0a3b93d8,(FULL)0x9301de77,(FULL)0x2d7e4efe,
606 (FULL)0xd8790d8a,(FULL)0x33711b76,(FULL)0x44df77e7,(FULL)0xc07dc30e,
607 (FULL)0x9ddd508f,(FULL)0xb9132ed0,(FULL)0xc6fb43cc,(FULL)0x45d06cf8,
608 (FULL)0xd585dd7b,(FULL)0x22bed18a,(FULL)0x10799ffa,(FULL)0x61c6cced,
609 (FULL)0xe4bd9aa9,(FULL)0xd7f2393b,(FULL)0xcfd55094,(FULL)0x706753fb,
610 (FULL)0xede6e446,(FULL)0xf65a6713,(FULL)0x47c0d5c3,(FULL)0x8bf6dfae,
611 (FULL)0x9f7927d6,(FULL)0xfb4dfc17,(FULL)0xe212c297,(FULL)0x12ebbc16,
612 (FULL)0xa00a954c,(FULL)0x43c71283,(FULL)0xe7bd40a5,(FULL)0x8957087a,
613 (FULL)0x08344837,(FULL)0xb0859d71,(FULL)0xaeb313f5,(FULL)0xfbf4b9a3,
614 (FULL)0xce81823a,(FULL)0x5e66e5be,(FULL)0x58ad6da1,(FULL)0x09a11c6e,
615 (FULL)0xc608054f,(FULL)0xc76f4316,(FULL)0x46084099,(FULL)0xb5821361,
616 (FULL)0x17a725ed,(FULL)0x4210008f,(FULL)0xd347c481,(FULL)0xe5ff8912
617 },
618 { /* shuffle table */
619 (FULL)0xec8abd57,(FULL)0x69a2296c,(FULL)0x99a6df81,(FULL)0x867e1869,
620 (FULL)0xd849a48a,(FULL)0xc05ab96b,(FULL)0xfa00554b,(FULL)0x7eb3ce0c,
621 (FULL)0x5a5a9acd,(FULL)0x520d01f6,(FULL)0x36022d81,(FULL)0xd4ef1e33,
622 (FULL)0xc6f84f70,(FULL)0xaf44772b,(FULL)0xa7c55173,(FULL)0x647e85a6,
623 (FULL)0x959df8d1,(FULL)0x26746cf1,(FULL)0x4db28abd,(FULL)0x98681a90,
624 (FULL)0x744c5cd2,(FULL)0xb146c969,(FULL)0x706f88c2,(FULL)0x8ce69d1f,
625 (FULL)0x21b4a748,(FULL)0xfd12eac4,(FULL)0x2710eea5,(FULL)0xf12e70fe,
626 (FULL)0x5901f2b5,(FULL)0x0b8f7805,(FULL)0x4f2c115e,(FULL)0x48860a76,
627 (FULL)0x30767e2c,(FULL)0x0edf6d2a,(FULL)0xfce2713b,(FULL)0x8a6d7dc5,
628 (FULL)0x4e0e2346,(FULL)0x46a362ea,(FULL)0x359f5aa7,(FULL)0x6c369a0a,
629 (FULL)0x41def54e,(FULL)0xdfca81fe,(FULL)0x96c2bc4e,(FULL)0x4b733819,
630 (FULL)0x6f3f14f9,(FULL)0x659e8b99,(FULL)0x93d47e6f,(FULL)0x8b97b934,
631 (FULL)0xdfa10a55,(FULL)0xa73a8704,(FULL)0xb06503da,(FULL)0x8d9eafe9,
632 (FULL)0xf32336b0,(FULL)0x2556fb88,(FULL)0x1002a161,(FULL)0xe71e9f75,
633 (FULL)0x200af907,(FULL)0x27a7be6e,(FULL)0xd028e9a3,(FULL)0x1b9b734e,
634 (FULL)0x4c0be0d3,(FULL)0x950cfeed,(FULL)0x2536d275,(FULL)0xf4c41694,
635 (FULL)0x5687b76e,(FULL)0xf05a58e8,(FULL)0x71a62d54,(FULL)0xba53ac01,
636 (FULL)0x285adc96,(FULL)0x4b14cbcb,(FULL)0xb00a5557,(FULL)0xfdf66edd,
637 (FULL)0x185b6ea1,(FULL)0xbb43d58d,(FULL)0xf355c9a6,(FULL)0x905db9cd,
638 (FULL)0x04429c8a,(FULL)0xfc3a07fc,(FULL)0xaa3a4f7e,(FULL)0x65d7e365,
639 (FULL)0xb243ac65,(FULL)0x2d284c18,(FULL)0x44e417fd,(FULL)0x72fba65d,
640 (FULL)0x5c934805,(FULL)0x422d50b4,(FULL)0xd1587441,(FULL)0xb62a6053,
641 (FULL)0x6f7ae035,(FULL)0xa5e71ce9,(FULL)0x595c8dd8,(FULL)0x93abca2e,
642 (FULL)0xe39afad5,(FULL)0x534231af,(FULL)0x12eaad56,(FULL)0x08d26cac,
643 (FULL)0x7fb1b1c2,(FULL)0xec18bf8d,(FULL)0xfaf6f09b,(FULL)0x3d28ea16,
644 (FULL)0x16697dd6,(FULL)0xea357a78,(FULL)0x420f3f51,(FULL)0x51471ea1,
645 (FULL)0x7f8946b4,(FULL)0x5e051aeb,(FULL)0x0cf0524c,(FULL)0x881be097,
646 (FULL)0x1b31489e,(FULL)0xd558b25b,(FULL)0x3a8b065c,(FULL)0x707d1a94,
647 (FULL)0x568ff836,(FULL)0x37017e66,(FULL)0x24c2f747,(FULL)0xb9cd627c,
648 (FULL)0xfb1d9ff6,(FULL)0x1485549f,(FULL)0xbdf2dc6f,(FULL)0x308d32d9,
649 (FULL)0xd543818a,(FULL)0x4d4142ca,(FULL)0x87ebba43,(FULL)0x5d9c7aee,
650 (FULL)0xe17adb2f,(FULL)0x81c5bdd8,(FULL)0xc8d8677a,(FULL)0x3dc9752e,
651 (FULL)0xc34e4212,(FULL)0x66b086e6,(FULL)0xc62b25e3,(FULL)0x3af7a90d,
652 (FULL)0x35539315,(FULL)0xf8349f79,(FULL)0xa22917f0,(FULL)0x6bcfd9d5,
653 (FULL)0x5f5ee517,(FULL)0x8639bb76,(FULL)0x8095b092,(FULL)0xd3c5e369,
654 (FULL)0x7eb44748,(FULL)0x8a33851e,(FULL)0xea54bbcf,(FULL)0x5e29d443,
655 (FULL)0x4d59a834,(FULL)0x0f84651f,(FULL)0xf1a5f951,(FULL)0x85040bea,
656 (FULL)0x98002078,(FULL)0x3dba1c74,(FULL)0xf0b2cc15,(FULL)0x5d70712b,
657 (FULL)0xcce8e5a7,(FULL)0xfa3af8eb,(FULL)0x04bba57d,(FULL)0xfb3e2237,
658 (FULL)0x8a950434,(FULL)0x5d3b8785,(FULL)0xba3f8dcf,(FULL)0xce3112bd,
659 (FULL)0x860d3051,(FULL)0x44904f55,(FULL)0x4ed3e98b,(FULL)0xcec8fed4,
660 (FULL)0x25d01ea4,(FULL)0x4581698d,(FULL)0x9a9548e0,(FULL)0x11eb6828,
661 (FULL)0xe911fac8,(FULL)0x796cb4c6,(FULL)0xb5fd813e,(FULL)0x2164cf26,
662 (FULL)0xd5de640f,(FULL)0x4ac8e0f5,(FULL)0x8802ab4e,(FULL)0xe9e757d7,
663 (FULL)0xf49dfcbd,(FULL)0x3c97de26,(FULL)0x6ee6dbe6,(FULL)0xc604881b,
664 (FULL)0x57d6154e,(FULL)0xa7c22a6e,(FULL)0x877b3cc7,(FULL)0x234e2370,
665 (FULL)0xdf1f8358,(FULL)0xc0bdb72b,(FULL)0xa95b7b55,(FULL)0x6522e0fc,
666 (FULL)0x22344162,(FULL)0xba174c90,(FULL)0x75d48867,(FULL)0x712c9b2d,
667 (FULL)0xe59f3700,(FULL)0x240f7e92,(FULL)0xad95d763,(FULL)0xe83cc2d4,
668 (FULL)0x4336d717,(FULL)0x8509445a,(FULL)0xdfff1804,(FULL)0xf1e572c5,
669 (FULL)0x623232dd,(FULL)0xed10eb5d,(FULL)0xd4f957e8,(FULL)0x9205ea1b,
670 (FULL)0x2ff062f5,(FULL)0x4973a54f,(FULL)0xe8c48cd5,(FULL)0x26b018f1,
671 (FULL)0xd1c7ed9f,(FULL)0x56908401,(FULL)0xdf89a247,(FULL)0x2e48937b,
672 (FULL)0x2be47129,(FULL)0x9d53069b,(FULL)0xc048a2b0,(FULL)0x98069e3b,
673 (FULL)0x1cd83f93,(FULL)0xf25b7d65,(FULL)0xe6f886c8,(FULL)0x2b004e6c,
674 (FULL)0x5c635935,(FULL)0xf618442a,(FULL)0x7198e052,(FULL)0xa502ab5c,
675 (FULL)0xa6c41b0b,(FULL)0xc14241a4,(FULL)0x7db9b18e,(FULL)0x720e845a,
676 (FULL)0x4b713918,(FULL)0x2abb13e9,(FULL)0x7f52467d,(FULL)0x90fc0c20,
677 (FULL)0x7868d348,(FULL)0x799c8ccd,(FULL)0x912a0ea4,(FULL)0xf4817ced,
678 (FULL)0xc4903a57,(FULL)0xd68c0f4c,(FULL)0xe2b7934c,(FULL)0xa3171f29,
679 (FULL)0x0b4ccc22,(FULL)0xb1158baa,(FULL)0x49a29eda,(FULL)0xf5d85553,
680 (FULL)0x959442ef,(FULL)0x59d1a078,(FULL)0xa67fd518,(FULL)0xdb9b4a96,
681 (FULL)0xd2870636,(FULL)0xcc7ca9ee,(FULL)0xecf59920,(FULL)0x548f021c,
682 (FULL)0x571bc8c5,(FULL)0x25b7f4b6,(FULL)0x3a44f536,(FULL)0x4fa52747,
683 (FULL)0xdf0ebdc2,(FULL)0xb246845f,(FULL)0x42058793,(FULL)0xdd8d68ae,
684 (FULL)0x9f6c39fb,(FULL)0x3ba13328,(FULL)0x7b6b42af,(FULL)0x8bfdfbf3,
685 (FULL)0x7fb2b3b0,(FULL)0xfb34c5ca,(FULL)0xd428e32a,(FULL)0x2345dcec,
686 (FULL)0xad42b63e,(FULL)0x6891e850,(FULL)0x362c1381,(FULL)0x930642c8,
687 (FULL)0x1886aff5,(FULL)0x13871e9b,(FULL)0x482bda55,(FULL)0xd0cf2407,
688 (FULL)0x5069bc31,(FULL)0x125b5fc9,(FULL)0xf07dfa5d,(FULL)0x9b71d0a9,
689 (FULL)0x6712e524,(FULL)0x55c044cc,(FULL)0xbb601978,(FULL)0xf0377358,
690 (FULL)0x7fa51e8b,(FULL)0x152ad5f8,(FULL)0x9fcdd9af,(FULL)0xe5ebf478,
691 (FULL)0x66ebce7e,(FULL)0x3d78e18c,(FULL)0xf36aa83f,(FULL)0x8246db72,
692 (FULL)0x2c64c0a3,(FULL)0xcc6ddc6d,(FULL)0x0d91851e,(FULL)0xa758d687,
693 (FULL)0x9488ee36,(FULL)0x24b20a6f,(FULL)0x09798197,(FULL)0xbe11ccdf,
694 (FULL)0x99c1f4e3,(FULL)0x11aca015,(FULL)0x6437ac05,(FULL)0x40e89e36,
695 (FULL)0x5af675f8,(FULL)0xc8bfc762,(FULL)0xb577e759,(FULL)0x6367c578,
696 (FULL)0x615f0b74,(FULL)0x00380346,(FULL)0x8de07d81,(FULL)0xee964cc4,
697 (FULL)0x859d9261,(FULL)0x17f6ac16,(FULL)0x3a6e2f6c,(FULL)0x092f4a17,
698 (FULL)0xb9024b95,(FULL)0x79981a3d,(FULL)0x04f7f540,(FULL)0x36db1660,
699 (FULL)0x65a2f1c8,(FULL)0xc36252cf,(FULL)0x124c9bd2,(FULL)0x705b6fde,
700 (FULL)0x85db40ce,(FULL)0x31e58dda,(FULL)0x9f5e8d6d,(FULL)0x6342b1a5,
701 (FULL)0xbd6d1d4d,(FULL)0x5c2c67d0,(FULL)0xba7e069d,(FULL)0x1fe5b46f,
702 (FULL)0xac72e13c,(FULL)0x21c46c6c,(FULL)0x9eb8f52a,(FULL)0xb80c5fd5,
703 (FULL)0xa74c92bc,(FULL)0x56c3aebf,(FULL)0xbf8c4196,(FULL)0xc1aff1fc,
704 (FULL)0x754ad208,(FULL)0x2b1df645,(FULL)0xd46eeb50,(FULL)0x5c734600,
705 (FULL)0x6a70a765,(FULL)0xe0ff1b12,(FULL)0x7a94547c,(FULL)0xd5416497,
706 (FULL)0x4ea35206,(FULL)0x67b59d7c,(FULL)0x779203b4,(FULL)0x53be7146,
707 (FULL)0x414026b8,(FULL)0x6b589fe5,(FULL)0x3083bfee,(FULL)0x9e81016c,
708 (FULL)0x3b4b7671,(FULL)0xb23526b9,(FULL)0x7ee300ba,(FULL)0x4fa9ffb1,
709 (FULL)0xad05fb21,(FULL)0x6217e212,(FULL)0xb294e6c2,(FULL)0xf5b3fcd3,
710 (FULL)0x216beb2a,(FULL)0xac040bbe,(FULL)0x71d0e78c,(FULL)0x1f8d8a54,
711 (FULL)0x9cfec96b,(FULL)0xb6d15b41,(FULL)0xd0508c78,(FULL)0xc5477845,
712 (FULL)0xb4bba621,(FULL)0x5b486e81,(FULL)0xef4c4121,(FULL)0x90c35c94,
713 (FULL)0xf6a6bc55,(FULL)0xefce7346,(FULL)0x25bdb9bb,(FULL)0xa27828d9,
714 (FULL)0xa1f0b205,(FULL)0xe3a53095,(FULL)0xd9f208ab,(FULL)0x1bfa6093,
715 (FULL)0x6842cdf4,(FULL)0xfb078f6a,(FULL)0x97133a38,(FULL)0x07806d72,
716 (FULL)0xa3ce9592,(FULL)0x2c6c901b,(FULL)0xebc1b789,(FULL)0x1f0ab2cf,
717 (FULL)0xe2d03d5e,(FULL)0x2ce81415,(FULL)0xaa9f2417,(FULL)0x7da45d5b,
718 (FULL)0xdd800682,(FULL)0x3be4f76d,(FULL)0x364d72d3,(FULL)0xdbf4e4a3,
719 (FULL)0x4fc59da5,(FULL)0xb538cccf,(FULL)0x487f66ec,(FULL)0xb0aa39d5,
720 (FULL)0x87927d3d,(FULL)0x2fd28dfd,(FULL)0x5900c6b1,(FULL)0xd14e77f0,
721 (FULL)0x5330c7b4,(FULL)0x2523fad2,(FULL)0xd82368a4,(FULL)0x991b5938,
722 (FULL)0x2b9c1302,(FULL)0xb7c11443,(FULL)0x1394b116,(FULL)0xdb842db6,
723 (FULL)0x78ed26d8,(FULL)0x3641548d,(FULL)0x0a61dacf,(FULL)0x274fa8ef,
724 (FULL)0x112df6f1,(FULL)0xa554ba63,(FULL)0x6b50438d,(FULL)0x7b7fe985,
725 (FULL)0xbb63bbad,(FULL)0xc9fa0042,(FULL)0xe27f00da,(FULL)0x3abf45d0,
726 (FULL)0x9f87aabb,(FULL)0xd95faa15,(FULL)0x3488e7ae,(FULL)0x4a95012e,
727 (FULL)0x0c642d04,(FULL)0x1be2bdb9,(FULL)0x8b4abf3e,(FULL)0x145c8881,
728 (FULL)0x544cf17f,(FULL)0x7f9fb635,(FULL)0xcc78db70,(FULL)0xb8ab2f62,
729 (FULL)0xb4242f9a,(FULL)0x8ee64bcd,(FULL)0x95dad129,(FULL)0xabd52858,
730 (FULL)0xccf31141,(FULL)0xbe722c2f,(FULL)0x575e26a9,(FULL)0x7c330703,
731 (FULL)0x361b79e4,(FULL)0x45d3e3b3,(FULL)0x54b2e6a6,(FULL)0x241163a7,
732 (FULL)0xf7cacb77,(FULL)0x8f678d7d,(FULL)0x83211d19,(FULL)0x988a68a4,
733 (FULL)0xba7836f6,(FULL)0x79599598,(FULL)0xeeda68bf,(FULL)0x4850c887,
734 (FULL)0x8052ce25,(FULL)0xafa69a71,(FULL)0xbdd73573,(FULL)0x8b21efc6,
735 (FULL)0xd0972493,(FULL)0x89dbae18,(FULL)0x537d9454,(FULL)0x560776bf,
736 (FULL)0x165310f2,(FULL)0x3c009f78,(FULL)0x0160c3af,(FULL)0xa3680021,
737 (FULL)0xa643bd40,(FULL)0x3353ec3c,(FULL)0x911dab02,(FULL)0x7e593f99,
738 (FULL)0x4f676e89,(FULL)0x72d1ddd9,(FULL)0x6b43c0ea,(FULL)0xfd18b8bd,
739 (FULL)0xddbd697d,(FULL)0x43cacef2,(FULL)0xacefe884,(FULL)0x2868a4d0,
740 (FULL)0xa506f013,(FULL)0x5f377b63,(FULL)0x05ca662b,(FULL)0xeaa0975e,
741 (FULL)0xeb433931,(FULL)0x3740e6b8,(FULL)0x08557948,(FULL)0xce85df00,
742 (FULL)0x547e33f9,(FULL)0x784745fb,(FULL)0xe5c6f598,(FULL)0x4a1fc5d4,
743 (FULL)0x768430a7,(FULL)0x85fa6fec,(FULL)0xd2332a51,(FULL)0x990d0c24,
744 (FULL)0x33b676d5,(FULL)0x55245c2c,(FULL)0xe2bcfa71,(FULL)0xb1091519,
745 (FULL)0xd23a28d8,(FULL)0x38521478,(FULL)0x9a573010,(FULL)0x9b794f89,
746 (FULL)0x699bb486,(FULL)0x61d225e8,(FULL)0x1c2158b0,(FULL)0x21476d24
747 }
748 #else
749 /\../\ FULL_BITS must be 32 or 64 /\../\ !!!
750 #endif
751 };
752
753 /*
754 * default subtractive 100 table
755 *
756 * The subtractive 100 table in init_s100 has been processed 256 times in order
757 * to preload the shuffle table. The array below is the table before
758 * this processing. These values have came from LavaRnd.
759 *
760 * This array is never changed, only copied.
761 */
762 STATIC CONST FULL def_subtract[SCNT] = {
763 #if FULL_BITS == SBITS
764 (FULL)U(0xc8c0370c7db7dc19), (FULL)U(0x738e33b940a06fbb),
765 (FULL)U(0x481abb76a859ed2b), (FULL)U(0x74106bb39ccdccb5),
766 (FULL)U(0x05a8eeb5c3173bfc), (FULL)U(0xefd5100d5a02e577),
767 (FULL)U(0xa69271f74030b24a), (FULL)U(0x641282fc16fe22c5),
768 (FULL)U(0x7aa7267c40438da3), (FULL)U(0x1fdf4abdc2d878d1),
769 (FULL)U(0xd9899e7a95702379), (FULL)U(0x5ea8e217d02d7f08),
770 (FULL)U(0x770587fe4d47a353), (FULL)U(0xde7d1bdd0a33a2b8),
771 (FULL)U(0x4378c3c5900e7c45), (FULL)U(0x77c9447819a514f9),
772 (FULL)U(0xfc5edb22843d1d32), (FULL)U(0x4fc42ce5e8ee5e6e),
773 (FULL)U(0xc938713c8488013e), (FULL)U(0x6a318f0320ab0cac),
774 (FULL)U(0x73e6d1a3ffc8bff3), (FULL)U(0x0cd3232a8ca96aa7),
775 (FULL)U(0x605c8036905f770d), (FULL)U(0x4d037b008b8d04a2),
776 (FULL)U(0x1ed81965cb277294), (FULL)U(0x408d9c477a254ff3),
777 (FULL)U(0x8b68587ae26c7377), (FULL)U(0xcff191a48a48832f),
778 (FULL)U(0x12d3df1d8aeb6fe6), (FULL)U(0xb2bf907e1feda37a),
779 (FULL)U(0x4e5f77193bb5f39f), (FULL)U(0x33ebcf6f8f5d1581),
780 (FULL)U(0x203c8e48d33654eb), (FULL)U(0x68d3656ef19c8a4e),
781 (FULL)U(0x3ec20b04986eb2af), (FULL)U(0x5d73a03b062c3841),
782 (FULL)U(0x836ce7095d4e49eb), (FULL)U(0x2310bc40c3f49221),
783 (FULL)U(0x3868ee48a6d0cbf6), (FULL)U(0x67578aa64a43deb1),
784 (FULL)U(0x6e3426c1150dfc26), (FULL)U(0xc541ccaa3131be30),
785 (FULL)U(0xf7e57432cec7aab2), (FULL)U(0x2b35de998cb3c873),
786 (FULL)U(0x7b9f77648663a5d7), (FULL)U(0x23b00e6aa771e5a6),
787 (FULL)U(0x859c775ca9985d05), (FULL)U(0x99636ea16b692f1f),
788 (FULL)U(0x8700ac703730800d), (FULL)U(0x461425024298a753),
789 (FULL)U(0xea4a411b809e955f), (FULL)U(0x3119ad4033709dfb),
790 (FULL)U(0xb76a6c6e5f01cb7c), (FULL)U(0x6109dc8a15984eaf),
791 (FULL)U(0x5d686db9a5ca9505), (FULL)U(0x8e80d7613b7e6add),
792 (FULL)U(0x79cbd718de6f6fd3), (FULL)U(0x40e9cd151da0f699),
793 (FULL)U(0xe82158bab24f312d), (FULL)U(0x79a4c927f5e5c36b),
794 (FULL)U(0xc25247c9a0039333), (FULL)U(0x936871161766d81d),
795 (FULL)U(0x3c6a03b4a6741327), (FULL)U(0xc8a7b6e8c002f29a),
796 (FULL)U(0x0e2a67c67bbd5ea3), (FULL)U(0x0929042d441eabc1),
797 (FULL)U(0x7dbe232a25e82085), (FULL)U(0x8cfb26e544fbac3d),
798 (FULL)U(0x8e40384d388ab983), (FULL)U(0x48dc1230554632f8),
799 (FULL)U(0xab405048ab492397), (FULL)U(0x21c9e2f5a118e387),
800 (FULL)U(0x484d1a8c343b61b5), (FULL)U(0xd49e3decab256f26),
801 (FULL)U(0xe615c7fd78f2d2e3), (FULL)U(0x8442cc33ce6cc2ed),
802 (FULL)U(0x0a3b93d844d4bbf6), (FULL)U(0x2d7e4efe9301de77),
803 (FULL)U(0x33711b76d8790d8a), (FULL)U(0xc07dc30e44df77e7),
804 (FULL)U(0xb9132ed09ddd508f), (FULL)U(0x45d06cf8c6fb43cc),
805 (FULL)U(0x22bed18ad585dd7b), (FULL)U(0x61c6cced10799ffa),
806 (FULL)U(0xd7f2393be4bd9aa9), (FULL)U(0x706753fbcfd55094),
807 (FULL)U(0xf65a6713ede6e446), (FULL)U(0x8bf6dfae47c0d5c3),
808 (FULL)U(0xfb4dfc179f7927d6), (FULL)U(0x12ebbc16e212c297),
809 (FULL)U(0x43c71283a00a954c), (FULL)U(0x8957087ae7bd40a5),
810 (FULL)U(0xb0859d7108344837), (FULL)U(0xfbf4b9a3aeb313f5),
811 (FULL)U(0x5e66e5bece81823a), (FULL)U(0x09a11c6e58ad6da1),
812 (FULL)U(0xc76f4316c608054f), (FULL)U(0xb582136146084099),
813 (FULL)U(0x4210008f17a725ed), (FULL)U(0xe5ff8912d347c481)
814 #elif 2*FULL_BITS == SBITS
815 (FULL)0x7db7dc19,(FULL)0xc8c0370c,(FULL)0x40a06fbb,(FULL)0x738e33b9,
816 (FULL)0xa859ed2b,(FULL)0x481abb76,(FULL)0x9ccdccb5,(FULL)0x74106bb3,
817 (FULL)0xc3173bfc,(FULL)0x05a8eeb5,(FULL)0x5a02e577,(FULL)0xefd5100d,
818 (FULL)0x4030b24a,(FULL)0xa69271f7,(FULL)0x16fe22c5,(FULL)0x641282fc,
819 (FULL)0x40438da3,(FULL)0x7aa7267c,(FULL)0xc2d878d1,(FULL)0x1fdf4abd,
820 (FULL)0x95702379,(FULL)0xd9899e7a,(FULL)0xd02d7f08,(FULL)0x5ea8e217,
821 (FULL)0x4d47a353,(FULL)0x770587fe,(FULL)0x0a33a2b8,(FULL)0xde7d1bdd,
822 (FULL)0x900e7c45,(FULL)0x4378c3c5,(FULL)0x19a514f9,(FULL)0x77c94478,
823 (FULL)0x843d1d32,(FULL)0xfc5edb22,(FULL)0xe8ee5e6e,(FULL)0x4fc42ce5,
824 (FULL)0x8488013e,(FULL)0xc938713c,(FULL)0x20ab0cac,(FULL)0x6a318f03,
825 (FULL)0xffc8bff3,(FULL)0x73e6d1a3,(FULL)0x8ca96aa7,(FULL)0x0cd3232a,
826 (FULL)0x905f770d,(FULL)0x605c8036,(FULL)0x8b8d04a2,(FULL)0x4d037b00,
827 (FULL)0xcb277294,(FULL)0x1ed81965,(FULL)0x7a254ff3,(FULL)0x408d9c47,
828 (FULL)0xe26c7377,(FULL)0x8b68587a,(FULL)0x8a48832f,(FULL)0xcff191a4,
829 (FULL)0x8aeb6fe6,(FULL)0x12d3df1d,(FULL)0x1feda37a,(FULL)0xb2bf907e,
830 (FULL)0x3bb5f39f,(FULL)0x4e5f7719,(FULL)0x8f5d1581,(FULL)0x33ebcf6f,
831 (FULL)0xd33654eb,(FULL)0x203c8e48,(FULL)0xf19c8a4e,(FULL)0x68d3656e,
832 (FULL)0x986eb2af,(FULL)0x3ec20b04,(FULL)0x062c3841,(FULL)0x5d73a03b,
833 (FULL)0x5d4e49eb,(FULL)0x836ce709,(FULL)0xc3f49221,(FULL)0x2310bc40,
834 (FULL)0xa6d0cbf6,(FULL)0x3868ee48,(FULL)0x4a43deb1,(FULL)0x67578aa6,
835 (FULL)0x150dfc26,(FULL)0x6e3426c1,(FULL)0x3131be30,(FULL)0xc541ccaa,
836 (FULL)0xcec7aab2,(FULL)0xf7e57432,(FULL)0x8cb3c873,(FULL)0x2b35de99,
837 (FULL)0x8663a5d7,(FULL)0x7b9f7764,(FULL)0xa771e5a6,(FULL)0x23b00e6a,
838 (FULL)0xa9985d05,(FULL)0x859c775c,(FULL)0x6b692f1f,(FULL)0x99636ea1,
839 (FULL)0x3730800d,(FULL)0x8700ac70,(FULL)0x4298a753,(FULL)0x46142502,
840 (FULL)0x809e955f,(FULL)0xea4a411b,(FULL)0x33709dfb,(FULL)0x3119ad40,
841 (FULL)0x5f01cb7c,(FULL)0xb76a6c6e,(FULL)0x15984eaf,(FULL)0x6109dc8a,
842 (FULL)0xa5ca9505,(FULL)0x5d686db9,(FULL)0x3b7e6add,(FULL)0x8e80d761,
843 (FULL)0xde6f6fd3,(FULL)0x79cbd718,(FULL)0x1da0f699,(FULL)0x40e9cd15,
844 (FULL)0xb24f312d,(FULL)0xe82158ba,(FULL)0xf5e5c36b,(FULL)0x79a4c927,
845 (FULL)0xa0039333,(FULL)0xc25247c9,(FULL)0x1766d81d,(FULL)0x93687116,
846 (FULL)0xa6741327,(FULL)0x3c6a03b4,(FULL)0xc002f29a,(FULL)0xc8a7b6e8,
847 (FULL)0x7bbd5ea3,(FULL)0x0e2a67c6,(FULL)0x441eabc1,(FULL)0x0929042d,
848 (FULL)0x25e82085,(FULL)0x7dbe232a,(FULL)0x44fbac3d,(FULL)0x8cfb26e5,
849 (FULL)0x388ab983,(FULL)0x8e40384d,(FULL)0x554632f8,(FULL)0x48dc1230,
850 (FULL)0xab492397,(FULL)0xab405048,(FULL)0xa118e387,(FULL)0x21c9e2f5,
851 (FULL)0x343b61b5,(FULL)0x484d1a8c,(FULL)0xab256f26,(FULL)0xd49e3dec,
852 (FULL)0x78f2d2e3,(FULL)0xe615c7fd,(FULL)0xce6cc2ed,(FULL)0x8442cc33,
853 (FULL)0x44d4bbf6,(FULL)0x0a3b93d8,(FULL)0x9301de77,(FULL)0x2d7e4efe,
854 (FULL)0xd8790d8a,(FULL)0x33711b76,(FULL)0x44df77e7,(FULL)0xc07dc30e,
855 (FULL)0x9ddd508f,(FULL)0xb9132ed0,(FULL)0xc6fb43cc,(FULL)0x45d06cf8,
856 (FULL)0xd585dd7b,(FULL)0x22bed18a,(FULL)0x10799ffa,(FULL)0x61c6cced,
857 (FULL)0xe4bd9aa9,(FULL)0xd7f2393b,(FULL)0xcfd55094,(FULL)0x706753fb,
858 (FULL)0xede6e446,(FULL)0xf65a6713,(FULL)0x47c0d5c3,(FULL)0x8bf6dfae,
859 (FULL)0x9f7927d6,(FULL)0xfb4dfc17,(FULL)0xe212c297,(FULL)0x12ebbc16,
860 (FULL)0xa00a954c,(FULL)0x43c71283,(FULL)0xe7bd40a5,(FULL)0x8957087a,
861 (FULL)0x08344837,(FULL)0xb0859d71,(FULL)0xaeb313f5,(FULL)0xfbf4b9a3,
862 (FULL)0xce81823a,(FULL)0x5e66e5be,(FULL)0x58ad6da1,(FULL)0x09a11c6e,
863 (FULL)0xc608054f,(FULL)0xc76f4316,(FULL)0x46084099,(FULL)0xb5821361,
864 (FULL)0x17a725ed,(FULL)0x4210008f,(FULL)0xd347c481,(FULL)0xe5ff8912
865 #else
866 /\../\ FULL_BITS must be 32 or 64 /\../\ !!!
867 #endif
868 };
869
870
871 /*
872 * Linear Congruential Constants
873 *
874 * a = 6316878969928993981 = 0x57aa0ff473c0ccbd
875 * c = 1363042948800878693 = 0x12ea805718e09865
876 *
877 * These constants are used in the randreseed64(). See below.
878 */
879 #if FULL_BITS == SBITS
880 STATIC CONST HALF a_vec[SHALFS] = { (HALF)0x73c0ccbd, (HALF)0x57aa0ff4 };
881 STATIC CONST HALF c_vec[SHALFS] = { (HALF)0x18e09865, (HALF)0x12ea8057 };
882 #elif 2*FULL_BITS == SBITS
883 STATIC CONST HALF a_vec[SHALFS] = { (HALF)0xccbd, (HALF)0x73c0,
884 (HALF)0x0ff4, (HALF)0x57aa };
885 STATIC CONST HALF c_vec[SHALFS] = { (HALF)0x9865, (HALF)0x18e0,
886 (HALF)0x8057, (HALF)0x12ea };
887 #else
888 /\../\ FULL_BITS must be 32 or 64 /\../\ !!!
889 #endif
890 STATIC CONST ZVALUE a_val = {(HALF *)a_vec, SHALFS, 0};
891 STATIC CONST ZVALUE c_val = {(HALF *)c_vec, SHALFS, 0};
892
893
894 /*
895 * current s100 generator state
896 */
897 STATIC RAND s100;
898
899
900 /*
901 * declare static functions
902 */
903 S_FUNC void randreseed64(ZVALUE seed, ZVALUE *res);
904 S_FUNC int slotcp(BITSTR *bitstr, FULL *src, int count);
905 S_FUNC void slotcp64(BITSTR *bitstr, FULL *src);
906
907
908 /*
909 * randreseed64 - scramble seed in 64 bit chunks
910 *
911 * given:
912 * a seed
913 *
914 * returns:
915 * a scrambled seed, or 0 if seed was 0
916 *
917 * It is 'nice' when a seed of "n" produces a 'significantly different'
918 * sequence than a seed of "n+1". Generators, by convention, assign
919 * special significance to the seed of '0'. It is an unfortunate that
920 * people often pick small seed values, particularly when large seed
921 * are of significance to the generators found in this file.
922 *
923 * We will process seed 64 bits at a time until the entire seed has been
924 * exhausted. If a 64 bit chunk is 0, then 0 is returned. If the 64 bit
925 * chunk is non-zero, we will produce a different and unique new scrambled
926 * chunk. In particular, if the seed is 0 we will return 0. If the seed
927 * is non-zero, we will return a different value (though chunks of 64
928 * zeros will remain zero). This scrambling will effectively eliminate
929 * the human perceptions that are noted above.
930 *
931 * It should be noted that the purpose of this process is to scramble a seed
932 * ONLY. We do not care if these generators produce good random numbers.
933 * We only want to help eliminate the human factors and perceptions
934 * noted above.
935 *
936 * This function scrambles all 64 bit chunks of a seed, by mapping [0,2^64)
937 * into [0,2^64). This map is one-to-one and onto. Mapping is performed
938 * using a linear congruence generator of the form:
939 *
940 * X1 <-- (a*X0 + c) % m
941 *
942 * with the exception that:
943 *
944 * 0 ==> 0 (so that srand(0) acts as default)
945 * randreseed64() is an 1-to-1 and onto map
946 *
947 * The generator are based on the linear congruential generators found in
948 * Knuth's "The Art of Computer Programming - Seminumerical Algorithms",
949 * vol 2, 2nd edition (1981), Section 3.6, pages 170-171.
950 *
951 * Because we process 64 bits we will take:
952 *
953 * m = 2^64 (based on note ii)
954 *
955 * We will scan the Rand Book of Random Numbers to look for an 'a' such that:
956 *
957 * a mod 8 == 5 (based on note iii)
958 * 0.01*m < a < 0.99*m (based on note iv)
959 * 0.01*2^64 < a < 0.99*2^64
960 *
961 * To help keep the generators independent, we want:
962 *
963 * a is prime
964 *
965 * The choice of an adder 'c' is considered immaterial according (based
966 * in note v). Knuth suggests 'c==1' or 'c==a'. We elect to select 'c'
967 * using the same process as we used to select 'a'. The choice is
968 * 'immaterial' after all, and as long as:
969 *
970 * gcd(c, m) == 1 (based on note v)
971 * gcd(c, 2^64) == 1
972 *
973 * the concerns are met. It can be shown that if we have:
974 *
975 * gcd(a, c) == 1
976 *
977 * then the adders and multipliers will be more independent.
978 *
979 * We will obtain the values 'a' and 'c for our generator from the
980 * Rand Book of Random Numbers. Because m=2^64 is 20 decimal digits long,
981 * we will search the Rand Book of Random Numbers 20 at a time. We will
982 * skip any of the 100 values that were used to initialize the subtractive 100
983 * generators. The values obtained from the Rand Book of Random Numbers are:
984 *
985 * a = 6316878969928993981
986 * c = 1363042948800878693
987 *
988 * As we stated before, we must map 0 ==> 0. The linear congruence
989 * generator would normally map as follows:
990 *
991 * 0 ==> 1363042948800878693 (0 ==> c)
992 *
993 * We can determine which 0<=y<m will produce a given value z by inverting the
994 * linear congruence generator:
995 *
996 * z = (a*y + c) % m
997 *
998 * z = a*y % m + c % m
999 *
1000 * z-c % m = a*y % m
1001 *
1002 * (z-c % m) * minv(a,m) = (a*y % m) * minv(a,m)
1003 * [minv(a,m) is the multiplicative inverse of a % m]
1004 *
1005 * ((z-c % m) * minv(a,m)) % m = ((a*y % m) * minv(a,m)) % m
1006 *
1007 * ((z-c % m) * minv(a,m)) % m = y % m
1008 * [a % m * minv(a,m) = 1 % m by definition]
1009 *
1010 * ((z-c) * minv(a,m)) % m = y % m
1011 *
1012 * ((z-c) * minv(a,m)) % m = y
1013 * [we defined y to be 0<=y<m]
1014 *
1015 * To determine which value maps back into 0, we let z = 0 and compute:
1016 *
1017 * ((0-c) * minv(a,m)) % m ==> 10239951819489363767
1018 *
1019 * and thus we find that the congruence generator would also normally map:
1020 *
1021 * 10239951819489363767 ==> 0
1022 *
1023 * To overcome this, and preserve the 1-to-1 and onto map, we force:
1024 *
1025 * 0 ==> 0
1026 * 10239951819489363767 ==> 1363042948800878693
1027 *
1028 * To repeat, this function converts a values into a seed value. With the
1029 * except of 'seed == 0', every value is mapped into a unique seed value.
1030 * This mapping need not be complex, random or secure. All we attempt
1031 * to do here is to allow humans who pick small or successive seed values
1032 * to obtain reasonably different sequences from the generators below.
1033 *
1034 * NOTE: This is NOT a pseudo random number generator. This function is
1035 * intended to be used internally by ss100rand() and sshufrand().
1036 */
1037 S_FUNC void
randreseed64(ZVALUE seed,ZVALUE * res)1038 randreseed64(ZVALUE seed, ZVALUE *res)
1039 {
1040 ZVALUE t; /* temp value */
1041 ZVALUE chunk; /* processed 64 bit chunk value */
1042 ZVALUE seed64; /* seed mod 2^64 */
1043 HALF *v64; /* 64 bit array of HALFs */
1044 long chunknum; /* 64 bit chunk number */
1045
1046 /*
1047 * quickly return 0 if seed is 0
1048 */
1049 if (ziszero(seed) || seed.len <= 0) {
1050 itoz(0, res);
1051 return;
1052 }
1053
1054 /*
1055 * allocate result
1056 */
1057 seed.sign = 0; /* use the value of seed */
1058 res->len = (int)(((seed.len+SHALFS-1) / SHALFS) * SHALFS);
1059 res->v = alloc(res->len);
1060 res->sign = 0;
1061 memset(res->v, 0, res->len*sizeof(HALF)); /* default value is 0 */
1062
1063 /*
1064 * process 64 bit chunks until done
1065 */
1066 chunknum = 0;
1067 while (!zislezero(seed)) {
1068
1069 /*
1070 * grab the lower 64 bits of seed
1071 */
1072 if (zge64b(seed)) {
1073 v64 = alloc(SHALFS);
1074 memcpy(v64, seed.v, SHALFS*sizeof(HALF));
1075 seed64.v = v64;
1076 seed64.len = SHALFS;
1077 seed64.sign = 0;
1078 } else {
1079 zcopy(seed, &seed64);
1080 }
1081 zshiftr(seed, SBITS);
1082 ztrim(&seed);
1083 ztrim(&seed64);
1084
1085 /*
1086 * do nothing if chunk is zero
1087 */
1088 if (ziszero(seed64)) {
1089 ++chunknum;
1090 zfree(seed64);
1091 continue;
1092 }
1093
1094 /*
1095 * Compute the linear congruence generator map:
1096 *
1097 * X1 <-- (a*X0 + c) mod m
1098 *
1099 * in other words:
1100 *
1101 * chunk == (a_val*seed + c_val) mod 2^64
1102 */
1103 zmul(seed64, a_val, &t);
1104 zfree(seed64);
1105 zadd(t, c_val, &chunk);
1106 zfree(t);
1107
1108 /*
1109 * form chunk mod 2^64
1110 */
1111 if (chunk.len > (SB32)SHALFS) {
1112 /* result is too large, reduce to 64 bits */
1113 v64 = alloc(SHALFS);
1114 memcpy(v64, chunk.v, SHALFS*sizeof(HALF));
1115 free(chunk.v);
1116 chunk.v = v64;
1117 chunk.len = SHALFS;
1118 ztrim(&chunk);
1119 }
1120
1121 /*
1122 * Normally, the above equation would map:
1123 *
1124 * f(0) == 1363042948800878693
1125 * f(10239951819489363767) == 0
1126 *
1127 * However, we have already forced f(0) == 0. To preserve the
1128 * 1-to-1 and onto map property, we force:
1129 *
1130 * f(10239951819489363767) ==> 1363042948800878693
1131 */
1132 if (ziszero(chunk)) {
1133 /* load 1363042948800878693 instead of 0 */
1134 zcopy(c_val, &chunk);
1135 memcpy(res->v+(chunknum*SHALFS), c_val.v,
1136 c_val.len*sizeof(HALF));
1137
1138 /*
1139 * load the 64 bit chunk into the result
1140 */
1141 } else {
1142 memcpy(res->v+(chunknum*SHALFS), chunk.v,
1143 chunk.len*sizeof(HALF));
1144 }
1145 ++chunknum;
1146 zfree(chunk);
1147 }
1148 ztrim(res);
1149 }
1150
1151
1152 /*
1153 * zsrand - seed the s100 generator
1154 *
1155 * given:
1156 * pseed - ptr to seed of the generator or NULL
1157 * pmat100 - subtractive 100 state table or NULL
1158 *
1159 * returns:
1160 * previous s100 state
1161 */
1162 RAND *
zsrand(CONST ZVALUE * pseed,CONST MATRIX * pmat100)1163 zsrand(CONST ZVALUE *pseed, CONST MATRIX *pmat100)
1164 {
1165 RAND *ret; /* previous s100 state */
1166 CONST VALUE *v; /* value from a passed matrix */
1167 ZVALUE zscram; /* scrambled 64 bit seed */
1168 ZVALUE ztmp; /* temp holding value for zscram */
1169 ZVALUE seed; /* to hold *pseed */
1170 FULL shufxor[SLEN]; /* zshufxor as an 64 bit array of FULLs */
1171 long indx; /* index to shuffle slots for seeding */
1172 int i;
1173
1174 /*
1175 * firewall
1176 */
1177 if (pseed != NULL && zisneg(*pseed)) {
1178 math_error("neg seeds for srand reserved for future use");
1179 /*NOTREACHED*/
1180 }
1181
1182 /*
1183 * initialize state if first call
1184 */
1185 if (!s100.seeded) {
1186 s100 = init_s100;
1187 }
1188
1189 /*
1190 * save the current state to return later
1191 */
1192 ret = (RAND *)malloc(sizeof(RAND));
1193 if (ret == NULL) {
1194 math_error("cannot allocate RAND state");
1195 /*NOTREACHED*/
1196 }
1197 *ret = s100;
1198
1199 /*
1200 * if call was srand(), just return current state
1201 */
1202 if (pseed == NULL && pmat100 == NULL) {
1203 return ret;
1204 }
1205
1206 /*
1207 * if call is srand(0), initialize and return quickly
1208 */
1209 if (pmat100 == NULL && ziszero(*pseed)) {
1210 s100 = init_s100;
1211 return ret;
1212 }
1213
1214 /*
1215 * clear buffered bits, initialize pointers
1216 */
1217 s100.seeded = 0; /* not seeded now */
1218 s100.j = INIT_J;
1219 s100.k = INIT_K;
1220 s100.bits = 0;
1221 memset(s100.buffer, 0, sizeof(s100.buffer));
1222
1223 /*
1224 * load subtractive table
1225 *
1226 * We will load the default subtractive table unless we are passed a
1227 * matrix. If we are passed a matrix, we will load the first 100
1228 * values mod 2^64 instead.
1229 */
1230 if (pmat100 == NULL) {
1231 memcpy(s100.slot, def_subtract, sizeof(def_subtract));
1232 } else {
1233
1234 /*
1235 * use the first 100 entries from the matrix arg
1236 */
1237 if (pmat100->m_size < S100) {
1238 math_error("matrix for srand has < 100 elements");
1239 /*NOTREACHED*/
1240 }
1241 for (v=pmat100->m_table, i=0; i < S100; ++i, ++v) {
1242
1243 /* reject if not integer */
1244 if (v->v_type != V_NUM || qisfrac(v->v_num)) {
1245 math_error("matrix for srand must contain ints");
1246 /*NOTREACHED*/
1247 }
1248
1249 /* load table element from matrix element mod 2^64 */
1250 SLOAD(s100, i, v->v_num->num);
1251 }
1252 }
1253
1254 /*
1255 * scramble the seed in 64 bit chunks
1256 */
1257 if (pseed != NULL) {
1258 seed.sign = pseed->sign;
1259 seed.len = pseed->len;
1260 seed.v = alloc(seed.len);
1261 zcopyval(*pseed, seed);
1262 randreseed64(seed, &zscram);
1263 zfree(seed);
1264 }
1265
1266 /*
1267 * xor subtractive table with the rehashed lower 64 bits of seed
1268 */
1269 if (pseed != NULL && !ziszero(zscram)) {
1270
1271 /* xor subtractive table with lower 64 bits of seed */
1272 SMOD64(shufxor, zscram);
1273 for (i=0; i < S100; ++i) {
1274 SXOR(s100, i, shufxor);
1275 }
1276 }
1277
1278 /*
1279 * shuffle subtractive 100 table according to seed, if passed
1280 */
1281 if (pseed != NULL && zge64b(zscram)) {
1282
1283 /* prepare the seed for subtractive slot shuffling */
1284 zshiftr(zscram, 64);
1285 ztrim(&zscram);
1286
1287 /* shuffle subtractive table */
1288 for (i=S100-1; i > 0 && !zislezero(zscram); --i) {
1289
1290 /* determine what we will swap with */
1291 indx = zdivi(zscram, i+1, &ztmp);
1292 zfree(zscram);
1293 zscram = ztmp;
1294
1295 /* do nothing if swap with itself */
1296 if (indx == i) {
1297 continue;
1298 }
1299
1300 /* swap slot[i] with slot[indx] */
1301 SSWAP(s100, i, indx);
1302 }
1303 zfree(zscram);
1304 } else if (pseed != NULL) {
1305 zfree(zscram);
1306 }
1307
1308 /*
1309 * load the shuffle table
1310 *
1311 * We will generate SHUFCNT entries from the subtractive 100 slots
1312 * and fill the shuffle table in consecutive order.
1313 */
1314 for (i=0; i < SHUFCNT; ++i) {
1315
1316 /*
1317 * skip values if required
1318 *
1319 * See:
1320 * Knuth's "The Art of Computer Programming -
1321 * Seminumerical Algorithms", Vol 2, 3rd edition (1998),
1322 * Section 3.6, page 188".
1323 */
1324 if (s100.need_to_skip <= 0) {
1325 int sk;
1326
1327 /* skip the require number of values */
1328 for (sk=0; sk < RAND_SKIP; ++sk) {
1329
1330 /* bump j and k */
1331 if (++s100.j >= S100) {
1332 s100.j = 0;
1333 }
1334 if (++s100.k >= S100) {
1335 s100.k = 0;
1336 }
1337
1338 /* slot[k] -= slot[j] */
1339 SSUB(s100, s100.k, s100.j);
1340
1341 /* NOTE: don't shuffle, no shuffle table yet */
1342 }
1343
1344 /* reset the skip count */
1345 s100.need_to_skip = RAND_CONSEQ_USE;
1346 if (conf->calc_debug & CALCDBG_RAND) {
1347 printf("rand: skipped %d states\n", RAND_SKIP);
1348 }
1349
1350 /* no skip, just decrement use counter */
1351 } else {
1352 --s100.need_to_skip;
1353 }
1354
1355 /* bump j and k */
1356 if (++s100.j >= S100) {
1357 s100.j = 0;
1358 }
1359 if (++s100.k >= S100) {
1360 s100.k = 0;
1361 }
1362
1363 /* slot[k] -= slot[j] */
1364 SSUB(s100, s100.k, s100.j);
1365
1366 /* shuf[i] = slot[k] */
1367 SSHUF(s100, i, s100.k);
1368 }
1369
1370 /*
1371 * note that we are seeded
1372 */
1373 s100.seeded = 1;
1374
1375 /*
1376 * return the previous state
1377 */
1378 return ret;
1379 }
1380
1381
1382 /*
1383 * zsetrand - set the s100 generator state
1384 *
1385 * given:
1386 * state - the state to copy
1387 *
1388 * returns:
1389 * previous s100 state
1390 */
1391 RAND *
zsetrand(CONST RAND * state)1392 zsetrand(CONST RAND *state)
1393 {
1394 RAND *ret; /* previous s100 state */
1395
1396 /*
1397 * initialize state if first call
1398 */
1399 if (!s100.seeded) {
1400 s100 = init_s100;
1401 }
1402
1403 /*
1404 * save the current state to return later
1405 */
1406 ret = randcopy(&s100);
1407
1408 /*
1409 * load the new state
1410 */
1411 s100 = *state;
1412
1413 /*
1414 * return the previous state
1415 */
1416 return ret;
1417 }
1418
1419
1420 /*
1421 * slotcp - copy up to 64 bits from a 64 bit array of FULLs to some HALFs
1422 *
1423 * We will copy data from an array of FULLs into an array of HALFs.
1424 * The destination within the HALFs is some bit location found in bitstr.
1425 * We will attempt to copy 64 bits, but if there is not enough room
1426 * (bits not yet loaded) in the destination bit string we will transfer
1427 * what we can.
1428 *
1429 * The src slot is 64 bits long and is stored as an array of FULLs.
1430 * When FULL_BITS is 64 the element is 1 FULL, otherwise FULL_BITS
1431 * is 32 bits and the element is 2 FULLs. The most significant bit
1432 * in the array (highest bit in the last FULL of the array) is to
1433 * be transferred to the most significant bit in the destination.
1434 *
1435 * given:
1436 * bitstr - most significant destination bit in a bit string
1437 * src - low order FULL in a 64 bit slot
1438 * count - number of bits to transfer (must be 0 < count <= 64)
1439 *
1440 * returns:
1441 * number of bits transferred
1442 */
1443 S_FUNC int
slotcp(BITSTR * bitstr,FULL * src,int count)1444 slotcp(BITSTR *bitstr, FULL *src, int count)
1445 {
1446 HALF *dh; /* most significant HALF in dest */
1447 int dnxtbit; /* next bit beyond most significant in dh */
1448 int need; /* number of bits we need to transfer */
1449 int ret; /* bits transferred */
1450
1451 /*
1452 * determine how many bits we actually need to transfer
1453 */
1454 dh = bitstr->loc;
1455 dnxtbit = bitstr->bit+1;
1456 count &= (SBITS-1);
1457 need = (bitstr->len < count) ? bitstr->len : count;
1458
1459 /*
1460 * prepare for the return
1461 *
1462 * Note the final bitstr location after we have moved the
1463 * position down 'need' bits.
1464 */
1465 bitstr->len -= need;
1466 bitstr->loc -= need / BASEB;
1467 bitstr->bit -= need % BASEB;
1468 if (bitstr->bit < 0) {
1469 --bitstr->loc;
1470 bitstr->bit += BASEB;
1471 }
1472 ret = need;
1473
1474 /*
1475 * deal with aligned copies quickly
1476 */
1477 if (dnxtbit == BASEB) {
1478 if (need == SBITS) {
1479 #if 2*FULL_BITS == SBITS
1480 *dh-- = (HALF)(src[SLEN-1] >> BASEB);
1481 *dh-- = (HALF)(src[SLEN-1]);
1482 #endif
1483 *dh-- = (HALF)(src[0] >> BASEB);
1484 *dh = (HALF)(src[0]);
1485 #if 2*FULL_BITS == SBITS
1486 } else if (need > FULL_BITS+BASEB) {
1487 *dh-- = (HALF)(src[SLEN-1] >> BASEB);
1488 *dh-- = (HALF)(src[SLEN-1]);
1489 *dh-- = (HALF)(src[0] >> BASEB);
1490 *dh = ((HALF)src[0] &
1491 highhalf[need-FULL_BITS-BASEB]);
1492 } else if (need > FULL_BITS) {
1493 *dh-- = (HALF)(src[SLEN-1] >> BASEB);
1494 *dh-- = (HALF)(src[SLEN-1]);
1495 *dh = ((HALF)(src[0] >> BASEB) &
1496 highhalf[need-FULL_BITS]);
1497 #endif
1498 } else if (need > BASEB) {
1499 *dh-- = (HALF)(src[SLEN-1] >> BASEB);
1500 *dh = ((HALF)(src[SLEN-1]) & highhalf[need-BASEB]);
1501 } else {
1502 *dh = ((HALF)(src[SLEN-1] >> BASEB) & highhalf[need]);
1503 }
1504 return ret;
1505 }
1506
1507 /*
1508 * load the most significant HALF
1509 */
1510 if (need >= dnxtbit) {
1511 /* fill up the most significant HALF */
1512 *dh-- |= (HALF)(src[SLEN-1] >> (FULL_BITS-dnxtbit));
1513 need -= dnxtbit;
1514 } else if (need > 0) {
1515 /* we exhaust our need before 1st half is filled */
1516 *dh |= (HALF)((src[SLEN-1] >> (FULL_BITS-need)) <<
1517 (dnxtbit-need));
1518 return ret; /* our need has been filled */
1519 } else {
1520 return ret; /* our need has been filled */
1521 }
1522
1523 /*
1524 * load the 2nd most significant HALF
1525 */
1526 if (need > BASEB) {
1527 /* fill up the 2nd most significant HALF */
1528 *dh-- = (HALF)(src[SLEN-1] >> (BASEB-dnxtbit));
1529 need -= BASEB;
1530 } else if (need > 0) {
1531 /* we exhaust our need before 2nd half is filled */
1532 *dh |= ((HALF)(src[SLEN-1] >> (BASEB-dnxtbit)) &
1533 highhalf[need]);
1534 return ret; /* our need has been filled */
1535 } else {
1536 return ret; /* our need has been filled */
1537 }
1538
1539 /*
1540 * load the 3rd most significant HALF
1541 *
1542 * At this point we know that our 3rd HALF will force us
1543 * to cross into a second FULL for systems with 32 bit FULLs.
1544 * We know this because the aligned case has been previously
1545 * taken care of above.
1546 *
1547 * For systems that have 64 bit FULLs (and 32 bit HALFs) this
1548 * is will be our least significant HALF. We also know that
1549 * the need must be < BASEB.
1550 */
1551 #if FULL_BITS == SBITS
1552 *dh |= (((HALF)src[0] & highhalf[dnxtbit+need]) << dnxtbit);
1553 #else
1554 if (need > BASEB) {
1555 /* load the remaining bits from the most significant FULL */
1556 *dh-- = ((((HALF)src[SLEN-1] & lowhalf[BASEB-dnxtbit])
1557 << dnxtbit) | (HALF)(src[0] >> (FULL_BITS-dnxtbit)));
1558 need -= BASEB;
1559 } else if (need > 0) {
1560 /* load the remaining bits from the most significant FULL */
1561 *dh-- |= (((((HALF)src[SLEN-1] & lowhalf[BASEB-dnxtbit])
1562 << dnxtbit) | (HALF)(src[0] >> (FULL_BITS-dnxtbit))) &
1563 highhalf[need]);
1564 return ret; /* our need has been filled */
1565 } else {
1566 return ret; /* our need has been filled */
1567 }
1568
1569 /*
1570 * load the 4th most significant HALF
1571 *
1572 * At this point, only 32 bit FULLs are operating.
1573 */
1574 if (need > BASEB) {
1575 /* fill up the 2nd most significant HALF */
1576 *dh-- = (HALF)(src[0] >> (BASEB-dnxtbit));
1577 /* no need TODO: need -= BASEB, because we are nearly done */
1578 } else if (need > 0) {
1579 /* we exhaust our need before 2nd half is filled */
1580 *dh |= ((HALF)(src[0] >> (BASEB-dnxtbit)) &
1581 highhalf[need]);
1582 return ret; /* our need has been filled */
1583 } else {
1584 return ret; /* our need has been filled */
1585 }
1586
1587 /*
1588 * load the 5th and least significant HALF
1589 *
1590 * At this point we know that the need will be satisfied.
1591 */
1592 *dh |= (((HALF)src[0] & lowhalf[BASEB-dnxtbit]) << dnxtbit);
1593 #endif
1594 return ret; /* our need has been filled */
1595 }
1596
1597
1598 /*
1599 * slotcp64 - copy 64 bits from a 64 bit array of FULLs to some HALFs
1600 *
1601 * We will copy data from an array of FULLs into an array of HALFs.
1602 * The destination within the HALFs is some bit location found in bitstr.
1603 * Unlike slotcp(), this function will always copy 64 bits.
1604 *
1605 * The src slot is 64 bits long and is stored as an array of FULLs.
1606 * When FULL_BITS is 64 this array is 1 FULL, otherwise FULL_BITS
1607 * is 32 bits and the array is 2 FULLs. The most significant bit
1608 * in the array (highest bit in the last FULL of the array) is to
1609 * be transferred to the most significant bit in the destination.
1610 *
1611 * given:
1612 * bitstr - most significant destination bit in a bit string
1613 * src - low order FULL in a 64 bit slot
1614 *
1615 * returns:
1616 * number of bits transferred
1617 */
1618 S_FUNC void
slotcp64(BITSTR * bitstr,FULL * src)1619 slotcp64(BITSTR *bitstr, FULL *src)
1620 {
1621 HALF *dh = bitstr->loc; /* most significant HALF in dest */
1622 int dnxtbit = bitstr->bit+1; /* next dh bit beyond most significant */
1623
1624 /*
1625 * prepare for the return
1626 *
1627 * Since we are moving the point 64 bits down, we know that
1628 * the bit location (bitstr->bit) will remain the same.
1629 */
1630 bitstr->len -= SBITS;
1631 bitstr->loc -= SBITS/BASEB;
1632
1633 /*
1634 * deal with aligned copies quickly
1635 */
1636 if (dnxtbit == BASEB) {
1637 #if 2*FULL_BITS == SBITS
1638 *dh-- = (HALF)(src[SLEN-1] >> BASEB);
1639 *dh-- = (HALF)(src[SLEN-1]);
1640 #endif
1641 *dh-- = (HALF)(src[0] >> BASEB);
1642 *dh = (HALF)(src[0]);
1643 return;
1644 }
1645
1646 /*
1647 * load the most significant HALF
1648 */
1649 *dh-- |= (HALF)(src[SLEN-1] >> (FULL_BITS-dnxtbit));
1650
1651 /*
1652 * load the 2nd most significant HALF
1653 */
1654 *dh-- = (HALF)(src[SLEN-1] >> (BASEB-dnxtbit));
1655
1656 /*
1657 * load the 3rd most significant HALF
1658 *
1659 * At this point we know that our 3rd HALF will force us
1660 * to cross into a second FULL for systems with 32 bit FULLs.
1661 * We know this because the aligned case has been previously
1662 * taken care of above.
1663 *
1664 * For systems that have 64 bit FULLs (and 32 bit HALFs) this
1665 * is will be our least significant HALF.
1666 */
1667 #if FULL_BITS == SBITS
1668 *dh |= (((HALF)src[0] & lowhalf[BASEB-dnxtbit]) << dnxtbit);
1669 #else
1670 /* load the remaining bits from the most significant FULL */
1671 *dh-- = ((((HALF)src[SLEN-1] & lowhalf[BASEB-dnxtbit])
1672 << dnxtbit) | (HALF)(src[0] >> (FULL_BITS-dnxtbit)));
1673
1674 /*
1675 * load the 4th most significant HALF
1676 *
1677 * At this point, only 32 bit FULLs are operating.
1678 */
1679 *dh-- = (HALF)(src[0] >> (BASEB-dnxtbit));
1680
1681 /*
1682 * load the 5th and least significant HALF
1683 *
1684 * At this point we know that the need will be satisfied.
1685 */
1686 *dh |= (((HALF)src[0] & lowhalf[BASEB-dnxtbit]) << dnxtbit);
1687 #endif
1688 }
1689
1690
1691 /*
1692 * zrandskip - skip s s100 bits
1693 *
1694 * given:
1695 * count - number of bits to be skipped
1696 */
1697 void
zrandskip(long cnt)1698 zrandskip(long cnt)
1699 {
1700 int indx; /* shuffle entry index */
1701
1702 /*
1703 * initialize state if first call
1704 */
1705 if (!s100.seeded) {
1706 s100 = init_s100;
1707 }
1708
1709 /*
1710 * skip required bits in the buffer
1711 */
1712 if (s100.bits > 0 && s100.bits <= cnt) {
1713
1714 /* just toss the buffer bits */
1715 cnt -= s100.bits;
1716 s100.bits = 0;
1717 memset(s100.buffer, 0, sizeof(s100.buffer));
1718
1719 } else if (s100.bits > 0 && s100.bits > cnt) {
1720
1721 /* buffer contains more bits than we need to toss */
1722 #if FULL_BITS == SBITS
1723 s100.buffer[0] <<= cnt;
1724 #else
1725 if (cnt >= FULL_BITS) {
1726 s100.buffer[SLEN-1] = (s100.buffer[0] << cnt);
1727 s100.buffer[0] = 0;
1728 } else {
1729 s100.buffer[SLEN-1] =
1730 ((s100.buffer[SLEN-1] << cnt) |
1731 (s100.buffer[0] >> (FULL_BITS-cnt)));
1732 s100.buffer[0] <<= cnt;
1733 }
1734 #endif
1735 s100.bits -= cnt;
1736 return; /* skip need satisfied */
1737 }
1738
1739 /*
1740 * skip 64 bits at a time
1741 */
1742 while (cnt >= SBITS) {
1743
1744 /*
1745 * skip values if required
1746 *
1747 * NOTE: This skip loop is part of the algorithm, not part
1748 * of the builtin function request.
1749 *
1750 * See:
1751 * Knuth's "The Art of Computer Programming -
1752 * Seminumerical Algorithms", Vol 2, 3rd edition (1998),
1753 * Section 3.6, page 188".
1754 */
1755 if (s100.need_to_skip <= 0) {
1756 int sk;
1757
1758 /* skip the require number of values */
1759 for (sk=0; sk < RAND_SKIP; ++sk) {
1760 /* bump j and k */
1761 if (++s100.j >= S100) {
1762 s100.j = 0;
1763 }
1764 if (++s100.k >= S100) {
1765 s100.k = 0;
1766 }
1767
1768 /* slot[k] -= slot[j] */
1769 SSUB(s100, s100.k, s100.j);
1770
1771 /* store s100.k into s100.slot[indx] */
1772 indx = SINDX(s100, s100.k);
1773 SSHUF(s100, indx, s100.k);
1774 }
1775
1776 /* reset the skip count */
1777 s100.need_to_skip = RAND_CONSEQ_USE;
1778 if (conf->calc_debug & CALCDBG_RAND) {
1779 printf("rand: skipped %d states\n", RAND_SKIP);
1780 }
1781
1782 /* no skip, just decrement use counter */
1783 } else {
1784 --s100.need_to_skip;
1785 }
1786
1787 /* bump j and k */
1788 if (++s100.j >= S100) {
1789 s100.j = 0;
1790 }
1791 if (++s100.k >= S100) {
1792 s100.k = 0;
1793 }
1794
1795 /* slot[k] -= slot[j] */
1796 SSUB(s100, s100.k, s100.j);
1797
1798 /* we will ignore the output value of s100.slot[indx] */
1799 indx = SINDX(s100, s100.k);
1800 cnt -= SBITS;
1801
1802 /* store s100.k into s100.slot[indx] */
1803 SSHUF(s100, indx, s100.k);
1804 }
1805
1806 /*
1807 * skip the final bits
1808 */
1809 if (cnt > 0) {
1810
1811 /*
1812 * skip values if required
1813 *
1814 * NOTE: This skip loop is part of the algorithm, not part
1815 * of the builtin function request.
1816 *
1817 * See:
1818 * Knuth's "The Art of Computer Programming -
1819 * Seminumerical Algorithms", Vol 2, 3rd edition (1998),
1820 * Section 3.6, page 188".
1821 */
1822 if (s100.need_to_skip <= 0) {
1823 int sk;
1824
1825 /* skip the require number of values */
1826 for (sk=0; sk < RAND_SKIP; ++sk) {
1827
1828 /* bump j and k */
1829 if (++s100.j >= S100) {
1830 s100.j = 0;
1831 }
1832 if (++s100.k >= S100) {
1833 s100.k = 0;
1834 }
1835
1836 /* slot[k] -= slot[j] */
1837 SSUB(s100, s100.k, s100.j);
1838
1839 /* store s100.k into s100.slot[indx] */
1840 indx = SINDX(s100, s100.k);
1841 SSHUF(s100, indx, s100.k);
1842 }
1843
1844 /* reset the skip count */
1845 s100.need_to_skip = RAND_CONSEQ_USE;
1846 if (conf->calc_debug & CALCDBG_RAND) {
1847 printf("rand: skipped %d states\n", RAND_SKIP);
1848 }
1849
1850 /* no skip, just decrement use counter */
1851 } else {
1852 --s100.need_to_skip;
1853 }
1854
1855 /* bump j and k */
1856 if (++s100.j >= S100) {
1857 s100.j = 0;
1858 }
1859 if (++s100.k >= S100) {
1860 s100.k = 0;
1861 }
1862
1863 /* slot[k] -= slot[j] */
1864 SSUB(s100, s100.k, s100.j);
1865
1866 /* we will ignore the output value of s100.slot[indx] */
1867 indx = SINDX(s100, s100.k);
1868
1869 /*
1870 * We know the buffer is empty, so fill it
1871 * with any unused bits. Copy SBITS-trans bits
1872 * from slot[indx] into buffer.
1873 */
1874 s100.bits = (int)(SBITS-cnt);
1875 memcpy(s100.buffer, &s100.shuf[indx*SLEN],
1876 sizeof(s100.buffer));
1877
1878 /*
1879 * shift the buffer bits all the way up to
1880 * the most significant bit
1881 */
1882 #if FULL_BITS == SBITS
1883 s100.buffer[0] <<= cnt;
1884 #else
1885 if (cnt >= FULL_BITS) {
1886 s100.buffer[SLEN-1] = (s100.buffer[0] << cnt);
1887 s100.buffer[0] = 0;
1888 } else {
1889 s100.buffer[SLEN-1] =
1890 ((s100.buffer[SLEN-1] << cnt) |
1891 (s100.buffer[0] >> (FULL_BITS-cnt)));
1892 s100.buffer[0] <<= cnt;
1893 }
1894 #endif
1895
1896 /* store s100.k into s100.slot[indx] */
1897 SSHUF(s100, indx, s100.k);
1898 }
1899 }
1900
1901
1902 /*
1903 * zrand - crank the s100 generator for some bits
1904 *
1905 * We will load the ZVALUE with random bits starting from the
1906 * most significant and ending with the lowest bit in the
1907 * least significant HALF.
1908 *
1909 * given:
1910 * count - number of bits required
1911 * res - where to place the random bits as ZVALUE
1912 */
1913 void
zrand(long cnt,ZVALUE * res)1914 zrand(long cnt, ZVALUE *res)
1915 {
1916 BITSTR dest; /* destination bit string */
1917 int trans; /* bits transferred */
1918 int indx; /* shuffle entry index */
1919
1920 /*
1921 * firewall
1922 */
1923 if (cnt <= 0) {
1924 if (cnt == 0) {
1925 /* zero length random number is always 0 */
1926 itoz(0, res);
1927 return;
1928 } else {
1929 math_error("negative zrand bit count");
1930 /*NOTREACHED*/
1931 }
1932 #if LONG_BITS > 32
1933 } else if (cnt > (1L<<31)) {
1934 math_error("huge rand bit count in internal zrand function");
1935 /*NOTREACHED*/
1936 #endif
1937 }
1938
1939 /*
1940 * initialize state if first call
1941 */
1942 if (!s100.seeded) {
1943 s100 = init_s100;
1944 }
1945
1946 /*
1947 * allocate storage
1948 */
1949 res->len = (LEN)((cnt+BASEB-1)/BASEB);
1950 res->v = alloc((LEN)((cnt+BASEB-1)/BASEB));
1951
1952 /*
1953 * dest bit string
1954 */
1955 dest.len = (int)cnt;
1956 dest.loc = res->v + (((cnt+BASEB-1)/BASEB)-1);
1957 dest.bit = (int)((cnt-1) % BASEB);
1958 memset(res->v, 0, (LEN)((cnt+BASEB-1)/BASEB)*sizeof(HALF));
1959
1960 /*
1961 * load from buffer first
1962 */
1963 if (s100.bits > 0) {
1964
1965 /*
1966 * We know there are only s100.bits in the buffer, so
1967 * transfer as much as we can (treating it as a slot)
1968 * and return the bit transfer count.
1969 */
1970 trans = slotcp(&dest, s100.buffer, s100.bits);
1971
1972 /*
1973 * If we need to keep bits in the buffer,
1974 * shift the buffer bits all the way up to
1975 * the most significant unused bit.
1976 */
1977 if (trans < s100.bits) {
1978 #if FULL_BITS == SBITS
1979 s100.buffer[0] <<= trans;
1980 #else
1981 if (trans >= FULL_BITS) {
1982 s100.buffer[SLEN-1] =
1983 (s100.buffer[0] << (trans-FULL_BITS));
1984 s100.buffer[0] = 0;
1985 } else {
1986 s100.buffer[SLEN-1] =
1987 ((s100.buffer[SLEN-1] << trans) |
1988 (s100.buffer[0] >> (FULL_BITS-trans)));
1989 s100.buffer[0] <<= trans;
1990 }
1991 #endif
1992 }
1993 /* note that we have fewer bits in the buffer */
1994 s100.bits -= trans;
1995 }
1996
1997 /*
1998 * spin the generator until we need less than 64 bits
1999 *
2000 * The buffer did not contain enough bits, so we crank the
2001 * s100 generator and load then 64 bits at a time.
2002 */
2003 while (dest.len >= SBITS) {
2004
2005 /*
2006 * skip values if required
2007 *
2008 * See:
2009 * Knuth's "The Art of Computer Programming -
2010 * Seminumerical Algorithms", Vol 2, 3rd edition (1998),
2011 * Section 3.6, page 188".
2012 */
2013 if (s100.need_to_skip <= 0) {
2014 int sk;
2015
2016 /* skip the require number of values */
2017 for (sk=0; sk < RAND_SKIP; ++sk) {
2018
2019 /* bump j and k */
2020 if (++s100.j >= S100) {
2021 s100.j = 0;
2022 }
2023 if (++s100.k >= S100) {
2024 s100.k = 0;
2025 }
2026
2027 /* slot[k] -= slot[j] */
2028 SSUB(s100, s100.k, s100.j);
2029
2030 /* select slot index to output */
2031 indx = SINDX(s100, s100.k);
2032
2033 /* store s100.k into s100.slot[indx] */
2034 SSHUF(s100, indx, s100.k);
2035 }
2036
2037 /* reset the skip count */
2038 s100.need_to_skip = RAND_CONSEQ_USE;
2039 if (conf->calc_debug & CALCDBG_RAND) {
2040 printf("rand: skipped %d states\n", RAND_SKIP);
2041 }
2042
2043 /* no skip, just decrement use counter */
2044 } else {
2045 --s100.need_to_skip;
2046 }
2047
2048 /* bump j and k */
2049 if (++s100.j >= S100) {
2050 s100.j = 0;
2051 }
2052 if (++s100.k >= S100) {
2053 s100.k = 0;
2054 }
2055
2056 /* slot[k] -= slot[j] */
2057 SSUB(s100, s100.k, s100.j);
2058
2059 /* select slot index to output */
2060 indx = SINDX(s100, s100.k);
2061
2062 /* move up to 64 bits from slot[indx] to dest */
2063 slotcp64(&dest, &s100.shuf[indx*SLEN]);
2064
2065 /* store s100.k into s100.slot[indx] */
2066 SSHUF(s100, indx, s100.k);
2067 }
2068
2069 /*
2070 * spin the generator one last time to fill out the remaining need
2071 */
2072 if (dest.len > 0) {
2073
2074 /*
2075 * skip values if required
2076 *
2077 * See:
2078 * Knuth's "The Art of Computer Programming -
2079 * Seminumerical Algorithms", Vol 2, 3rd edition (1998),
2080 * Section 3.6, page 188".
2081 */
2082 if (s100.need_to_skip <= 0) {
2083 int sk;
2084
2085 /* skip the require number of values */
2086 for (sk=0; sk < RAND_SKIP; ++sk) {
2087 /* bump j and k */
2088 if (++s100.j >= S100) {
2089 s100.j = 0;
2090 }
2091 if (++s100.k >= S100) {
2092 s100.k = 0;
2093 }
2094
2095 /* slot[k] -= slot[j] */
2096 SSUB(s100, s100.k, s100.j);
2097
2098 /* select slot index to output */
2099 indx = SINDX(s100, s100.k);
2100
2101 /* store s100.k into s100.slot[indx] */
2102 SSHUF(s100, indx, s100.k);
2103 }
2104
2105 /* reset the skip count */
2106 s100.need_to_skip = RAND_CONSEQ_USE;
2107 if (conf->calc_debug & CALCDBG_RAND) {
2108 printf("rand: skipped %d states\n", RAND_SKIP);
2109 }
2110
2111 /* no skip, just decrement use counter */
2112 } else {
2113 --s100.need_to_skip;
2114 }
2115
2116 /* bump j and k */
2117 if (++s100.j >= S100) {
2118 s100.j = 0;
2119 }
2120 if (++s100.k >= S100) {
2121 s100.k = 0;
2122 }
2123
2124 /* slot[k] -= slot[j] */
2125 SSUB(s100, s100.k, s100.j);
2126
2127 /* select slot index to output */
2128 indx = SINDX(s100, s100.k);
2129
2130 /* move up to 64 bits from slot[indx] to dest */
2131 trans = slotcp(&dest, &s100.shuf[indx*SLEN], dest.len);
2132
2133 /* buffer up unused bits if we are done */
2134 if (trans != SBITS) {
2135
2136 /*
2137 * We know the buffer is empty, so fill it
2138 * with any unused bits. Copy SBITS-trans bits
2139 * from slot[indx] into buffer.
2140 */
2141 s100.bits = SBITS-trans;
2142 memcpy(s100.buffer, &s100.shuf[indx*SLEN],
2143 sizeof(s100.buffer));
2144
2145 /*
2146 * shift the buffer bits all the way up to
2147 * the most significant bit
2148 */
2149 #if FULL_BITS == SBITS
2150 s100.buffer[0] <<= trans;
2151 #else
2152 if (trans >= FULL_BITS) {
2153 s100.buffer[SLEN-1] =
2154 (s100.buffer[0] << (trans-FULL_BITS));
2155 s100.buffer[0] = 0;
2156 } else {
2157 s100.buffer[SLEN-1] =
2158 ((s100.buffer[SLEN-1] << trans) |
2159 (s100.buffer[0] >> (FULL_BITS-trans)));
2160 s100.buffer[0] <<= trans;
2161 }
2162 #endif
2163 }
2164
2165 /* store s100.k into s100.slot[indx] */
2166 SSHUF(s100, indx, s100.k);
2167 }
2168 res->sign = 0;
2169 ztrim(res);
2170 }
2171
2172
2173 /*
2174 * zrandrange - generate an s100 random value in the range [low, beyond)
2175 *
2176 * given:
2177 * low - low value of range
2178 * beyond - beyond end of range
2179 * res - where to place the random bits as ZVALUE
2180 */
2181 void
zrandrange(CONST ZVALUE low,CONST ZVALUE beyond,ZVALUE * res)2182 zrandrange(CONST ZVALUE low, CONST ZVALUE beyond, ZVALUE *res)
2183 {
2184 ZVALUE range; /* beyond-low */
2185 ZVALUE rval; /* random value [0, 2^bitlen) */
2186 ZVALUE rangem1; /* range - 1 */
2187 long bitlen; /* smallest power of 2 >= diff */
2188
2189 /*
2190 * firewall
2191 */
2192 if (zrel(low, beyond) >= 0) {
2193 math_error("srand low range >= beyond range");
2194 /*NOTREACHED*/
2195 }
2196
2197 /*
2198 * determine the size of the random number needed
2199 */
2200 zsub(beyond, low, &range);
2201 if (zisone(range)) {
2202 zfree(range);
2203 *res = low;
2204 return;
2205 }
2206 zsub(range, _one_, &rangem1);
2207 bitlen = 1+zhighbit(rangem1);
2208 zfree(rangem1);
2209
2210 /*
2211 * generate a random value between [0, diff)
2212 *
2213 * We will not fall into the trap of thinking that we can simply take
2214 * a value mod 'range'. Consider the case where 'range' is '80'
2215 * and we are given pseudo-random numbers [0,100). If we took them
2216 * mod 80, then the numbers [0,20) would be produced more frequently
2217 * because the numbers [81,100) mod 80 wrap back into [0,20).
2218 */
2219 rval.v = NULL;
2220 do {
2221 if (rval.v != NULL) {
2222 zfree(rval);
2223 }
2224 zrand(bitlen, &rval);
2225 } while (zrel(rval, range) >= 0);
2226
2227 /*
2228 * add in low value to produce the range [0+low, diff+low)
2229 * which is the range [low, beyond)
2230 */
2231 zadd(rval, low, res);
2232 zfree(rval);
2233 zfree(range);
2234 }
2235
2236
2237 /*
2238 * irand - generate an s100 random long in the range [0, s)
2239 *
2240 * given:
2241 * s - limit of the range
2242 *
2243 * returns:
2244 * random long in the range [0, s)
2245 */
2246 long
irand(long s)2247 irand(long s)
2248 {
2249 ZVALUE z1, z2;
2250 long res;
2251
2252 if (s <= 0) {
2253 math_error("Non-positive argument for irand()");
2254 /*NOTREACHED*/
2255 }
2256 if (s == 1)
2257 return 0;
2258 itoz(s, &z1);
2259 zrandrange(_zero_, z1, &z2);
2260 res = ztoi(z2);
2261 zfree(z1);
2262 zfree(z2);
2263 return res;
2264 }
2265
2266
2267 /*
2268 * randcopy - make a copy of an s100 state
2269 *
2270 * given:
2271 * state - the state to copy
2272 *
2273 * returns:
2274 * a malloced copy of the state
2275 */
2276 RAND *
randcopy(CONST RAND * state)2277 randcopy(CONST RAND *state)
2278 {
2279 RAND *ret; /* return copy of state */
2280
2281 /*
2282 * malloc state
2283 */
2284 ret = (RAND *)malloc(sizeof(RAND));
2285 if (ret == NULL) {
2286 math_error("can't allocate RAND state");
2287 /*NOTREACHED*/
2288 }
2289 memcpy(ret, state, sizeof(RAND));
2290
2291 /*
2292 * return copy
2293 */
2294 return ret;
2295 }
2296
2297
2298 /*
2299 * randfree - free an s100 state
2300 *
2301 * given:
2302 * state - the state to free
2303 */
2304 void
randfree(RAND * state)2305 randfree(RAND *state)
2306 {
2307 /* free it */
2308 free(state);
2309 }
2310
2311
2312 /*
2313 * randcmp - compare two s100 states
2314 *
2315 * given:
2316 * s1 - first state to compare
2317 * s2 - second state to compare
2318 *
2319 * return:
2320 * TRUE if states differ
2321 */
2322 BOOL
randcmp(CONST RAND * s1,CONST RAND * s2)2323 randcmp(CONST RAND *s1, CONST RAND *s2)
2324 {
2325 /*
2326 * assume uninitialized state == the default seeded state
2327 */
2328 if (!s1->seeded) {
2329 if (!s2->seeded) {
2330 /* uninitialized == uninitialized */
2331 return FALSE;
2332 } else {
2333 /* uninitialized only equals default state */
2334 return randcmp(s2, &init_s100);
2335 }
2336 } else if (!s2->seeded) {
2337 /* uninitialized only equals default state */
2338 return randcmp(s1, &init_s100);
2339 }
2340
2341 /* compare states */
2342 return (BOOL)(memcmp(s1, s2, sizeof(RAND)) != 0);
2343 }
2344
2345
2346 /*
2347 * randprint - print an s100 state
2348 *
2349 * given:
2350 * state - state to print
2351 * flags - print flags passed from printvalue() in value.c
2352 */
2353 /*ARGSUSED*/
2354 void
randprint(CONST RAND * UNUSED (state),int UNUSED (flags))2355 randprint(CONST RAND *UNUSED(state), int UNUSED(flags))
2356 {
2357 math_str("RAND state");
2358 }
2359