1 /*
2 * zrandom - Blum-Blum-Shub pseudo-random 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: 1997/02/15 04:01:56
21 * File existed as early as: 1997
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 a Blum-Blum-Shub generator:
31 *
32 * To shorten the name, We refer to this generator as the Blum generator.
33 *
34 * This generator is described in the papers:
35 *
36 * Blum, Blum, and Shub, "Comparison of Two Pseudorandom Number
37 * Generators", in Chaum, D. et. al., "Advances in Cryptology:
38 * Proceedings Crypto 82", pp. 61-79, Plenum Press, 1983.
39 *
40 * Blum, Blum, and Shub, "A Simple Unpredictable Pseudo-Random
41 * Number Generator", SIAM Journal of Computing, v. 15, n. 2,
42 * 1986, pp. 364-383.
43 *
44 * U. V. Vazirani and V. V. Vazirani, "Trapdoor Pseudo-Random
45 * Number Generators with Applications to Protocol Design",
46 * Proceedings of the 24th IEEE Symposium on the Foundations
47 * of Computer Science, 1983, pp. 23-30.
48 *
49 * U. V. Vazirani and V. V. Vazirani, "Efficient and Secure
50 * Pseudo-Random Number Generation", Proceedings of the 24th
51 * IEEE Symposium on the Foundations of Computer Science,
52 * 1984, pp. 458-463.
53 *
54 * U. V. Vazirani and V. V. Vazirani, "Efficient and Secure
55 * Pseudo-Random Number Generation", Advances in Cryptology -
56 * Proceedings of CRYPTO '84, Berlin: Springer-Verlag, 1985,
57 * pp. 193-202.
58 *
59 * Sciences 28, pp. 270-299.
60 *
61 * Bruce Schneier, "Applied Cryptography", John Wiley & Sons,
62 * 1st edition (1994), pp 365-366.
63 *
64 * This generator is considered 'strong' in that it passes all
65 * polynomial-time statistical tests. The sequences produced
66 * are random in an absolutely precise way. There is absolutely
67 * no better way to predict the sequence than by tossing a coin
68 * (as with TRULY random numbers) EVEN IF YOU KNOW THE MODULUS!
69 * Furthermore, having a large chunk of output from the sequence
70 * does not help. The BITS THAT FOLLOW OR PRECEDE A SEQUENCE
71 * ARE UNPREDICTABLE!
72 *
73 * Of course the Blum modulus should have a long period. The default
74 * Blum modulus as well as the compiled in Blum moduli have very long
75 * periods. When using your own Blum modulus, a little care is needed
76 * to avoid generators with very short periods. (see below)
77 *
78 * To compromise the generator, an adversary must either factor the
79 * modulus or perform an exhaustive search just to determine the next
80 * (or previous) bit. If we make the modulus hard to factor
81 * (such as the product of two large well chosen primes) breaking
82 * the sequence could be intractable for todays computers and methods.
83 *
84 ******************************************************************************
85 *
86 * GOALS:
87 *
88 * The goals of this package are:
89 *
90 * All magic numbers are explained:
91 *
92 * I distrust systems with constants (magic numbers) and tables
93 * that have no justification. I believe that I have
94 * done my best to justify all of the magic numbers used.
95 *
96 * Full documentation:
97 *
98 * You have this source file, plus background publications,
99 * what more could you ask?
100 *
101 * Large selection of seeds:
102 *
103 * Seeds are not limited to a small number of bits. A seed
104 * may be of any size.
105 *
106 * The strength of the generators may be tuned to meet the need:
107 *
108 * By using the appropriate seed and other arguments, one may
109 * increase the strength of the generator to suit the need of
110 * the application. One does not have just a few levels.
111 *
112 * Even though I have done my best to implement a good system, you still
113 * must use these routines your own risk.
114 *
115 * Share and enjoy! :-)
116 */
117
118 /*
119 * ON THE GENERATOR:
120 *
121 * The Blum generator is the best generator in this package. It
122 * produces a cryptographically strong pseudo-random bit sequence.
123 * Internally, a fixed number of bits are generated after each
124 * generator iteration. Any unused bits are saved for the next call
125 * to the generator. The Blum generator is not too slow, though
126 * seeding the generator via srandom(seed,plen,qlen) can be slow.
127 * Shortcuts and pre-defined generators have been provided for this reason.
128 * Use of Blum should be more than acceptable for many applications.
129 *
130 * The Blum generator as the following calc interfaces:
131 *
132 * random(min, beyond) (where min < beyond)
133 *
134 * Print a Blum generator random value over interval [min,beyond).
135 *
136 * This form returns integers of the form:
137 *
138 * min <= value < beyond
139 *
140 * random()
141 *
142 * Same as random(0, 2^64). Print 64 bits.
143 *
144 * This form returns integers of the form:
145 *
146 * 0 <= value < 2^64
147 *
148 * random(lim) (where lim > 0)
149 *
150 * Same as random(0, lim).
151 *
152 * This form returns integers of the form:
153 *
154 * 0 <= value < lim
155 *
156 * randombit(x) (where x > 0)
157 *
158 * Same as random(0, 2^x). Print x bits.
159 *
160 * This form returns integers of the form:
161 *
162 * 0 <= value < 2^x
163 *
164 * randombit(neg) (where neg < 0)
165 *
166 * Skip neg random bits and return the bit skip count.
167 */
168
169 /*
170 * INITIALIZATION AND SEEDS:
171 *
172 * All generators come already seeded with precomputed initial constants.
173 * Thus, it is not required to seed a generator before using it.
174 *
175 * The Blum generator may be initialized and seeded via srandom().
176 *
177 * Using a seed of '0' will reload generators with their initial states.
178 *
179 * srandom(0) restore Blum generator to the initial state
180 *
181 * The above single arg calls are fairly fast.
182 *
183 * The call:
184 *
185 * srandom(seed, newn)
186 *
187 * is fast when the config value "srandom" is 0, 1 or 2.
188 *
189 * Optimal seed range for the Blum generator:
190 *
191 * There is no limit on the size of a seed. On the other hand,
192 * in most cases the seed is taken modulo the Blum modulus.
193 * Using a seed that is too small (except for 0) results in
194 * an internal generator be used to increase its size.
195 *
196 * It is faster to use seeds that are in the half open internal
197 * [sqrt(n), n) where n is the Blum modulus.
198 *
199 * The default Blum modulus is 260 bits long, so when using a the
200 * single arg call, a seed of between 128 and 256 bits is reasonable.
201 *
202 ******************************************************************************
203 *
204 * srandom(seed)
205 *
206 * We attempt to set the quadratic residue and possibly the Blum modulus.
207 *
208 * Any internally buffered random bits are flushed.
209 *
210 * The Blum modulus is only set if seed == 0.
211 *
212 * The initial quadratic residue is set according to the seed
213 * arg as defined below.
214 *
215 * seed >= 2^32:
216 * -------------
217 * Use seed to compute a new quadratic residue for use with
218 * the current Blum modulus. We will successively square mod Blum
219 * modulus until we get a smaller value (modulus wrap).
220 *
221 * The follow calc resource file produces an equivalent effect:
222 *
223 * n = default_modulus; (* n is the new Blum modulus *)
224 * r = seed;
225 * do {
226 * last_r = r;
227 * r = pmod(last_r, 2, n);
228 * } while (r > last_r); (* r is the new quadratic residue *)
229 *
230 * NOTE: The Blum modulus is not set by this call.
231 *
232 * 0 < seed < 2^32:
233 * ----------------
234 * Reserved for future use.
235 *
236 * seed == 0:
237 * ----------
238 * Restore the initial state and modulus of the Blum generator.
239 * After this call, the Blum generator is restored to its initial
240 * state after calc started.
241 *
242 * The Blum prime factors of the modulus have been disclosed (see
243 * "SOURCE OF MAGIC NUMBERS" below). If you want to use moduli that
244 * have not been disclosed, use srandom(seed, newn) with the
245 * appropriate args as noted below.
246 *
247 * The follow calc resource file produces an equivalent effect:
248 *
249 * n = default_modulus; (* as used by the initial state *)
250 * r = default_residue; (* as used by the initial state *)
251 *
252 * NOTE: The Blum modulus is changed by this call.
253 *
254 * seed < 0:
255 * ---------
256 * Reserved for future use.
257 *
258 ******************************************************************************
259 *
260 * srandom(seed, newn)
261 *
262 * We attempt to set the Blum modulus and quadratic residue.
263 * Any internally buffered random bits are flushed.
264 *
265 * If newn == 1 mod 4, then we will assume that it is the
266 * product of two Blum primes (primes == 3 mod 4) and use it
267 * as the Blum modulus.
268 *
269 * The new quadratic residue is set according to the seed
270 * arg as defined below.
271 *
272 * seed >= 2^32, newn >= 2^32:
273 * ---------------------------
274 * Assuming that 'newn' == 3 mod 4, then we will use it as
275 * the Blum modulus.
276 *
277 * We will use the seed arg to compute a new quadratic residue.
278 * We will successively square it mod Blum modulus until we get
279 * a smaller value (modulus wrap).
280 *
281 * The follow calc resource file produces an equivalent effect:
282 *
283 * if (newn % 4 == 1) {
284 * n = newn; (* n is the new Blum modulus *)
285 * r = seed;
286 * do {
287 * last_r = r;
288 * r = pmod(last_r, 2, n);
289 * } while (r > last_r); (* r is the new quadratic residue *)
290 * } else {
291 * quit "newn (2nd arg) must be 3 mod 4";
292 * }
293 *
294 * 0 < seed < 2^32, newn >= 2^32:
295 * ------------------------------
296 * Reserved for future use.
297 *
298 * seed == 0, newn >= 2^32:
299 * ------------------------
300 * Assuming that 'newn' == 3 mod 4, then we will use it as
301 * the Blum modulus.
302 *
303 * The initial quadratic residue will be as if the default initial
304 * quadratic residue arg was given.
305 *
306 * The follow calc resource file produces an equivalent effect:
307 *
308 * srandom(default_residue, newn)
309 *
310 * or in other words:
311 *
312 * if (newn % 4 == 1) {
313 * n = newn; (* n is the new Blum modulus *)
314 * r = default_residue; (* as used by the initial state *)
315 * do {
316 * last_r = r;
317 * r = pmod(last_r, 2, n);
318 * } while (r > last_r); (* r is the new quadratic residue *)
319 * } else {
320 * quit "newn (2nd arg) must be 3 mod 4";
321 * }
322 *
323 * seed < 0, newn >= 2^32:
324 * -----------------------
325 * Reserved for future use.
326 *
327 * any seed, 20 < newn < 2^32:
328 * ---------------------------
329 * Reserved for future use.
330 *
331 * seed >= 2^32, 0 < newn <= 20:
332 * -----------------------------
333 * Set the Blum modulus to one of the pre-defined Blum moduli.
334 * See below for the values of these pre-defined Blum moduli and how
335 * they were computed.
336 *
337 * We will use the seed arg to compute a new quadratic residue.
338 * We will successively square it mod Blum modulus until we get
339 * a smaller value (modulus wrap).
340 *
341 * The follow calc resource file produces an equivalent effect:
342 *
343 * n = n[newn]; (* n is new Blum modulus, see below *)
344 * r = seed;
345 * do {
346 * last_r = r;
347 * r = pmod(last_r, 2, n);
348 * } while (r > last_r); (* r is the new quadratic residue *)
349 *
350 * 0 < seed < 2^32, 0 < newn <= 20:
351 * --------------------------------
352 * Reserved for future use.
353 *
354 * seed == 0, 0 < newn <= 20:
355 * --------------------------
356 * Set the Blum modulus to one of the pre-defined Blum moduli.
357 * The new quadratic residue will also be set to one of
358 * the pre-defined quadratic residues.
359 *
360 * The follow calc resource file produces an equivalent effect:
361 *
362 * srandom(r[newn], n[newn])
363 *
364 * or in other words:
365 *
366 * n = n[newn]; (* n is the new Blum modulus, see below *)
367 * r = r[newn]; (* r is the new quadratic residue *)
368 *
369 * The pre-defined Blum moduli was computed by searching for Blum
370 * primes (primes == 3 mod 4) starting from new values that
371 * were selected by LavaRnd, a hardware random number generator.
372 * See the URL:
373 *
374 * http://www.LavaRnd.org/
375 *
376 * for an explanation of how the LavaRnd random number generator works.
377 *
378 * For a given newn, we select a given bit length. For 0 < newn <= 20,
379 * the bit length selected was by:
380 *
381 * bitlen = 2^(int((newn-1)/4)+7) + small_random_value;
382 *
383 * where small_random_value is also generated by LavaRnd. For
384 * 1 <= newn <= 16, small_random_value is a random value in [0,40).
385 * For 17 < newn <= 20, small_random_value is a random value in [0,120).
386 * Given two random integers generated by LavaRnd, we used the following
387 * to compute Blum primes:
388 *
389 * (* find the first Blum prime *)
390 * fp = int((ip-1)/2); (* ip was generated by LavaRnd *)
391 * do {
392 * fp = nextcand(fp+2, 25, 0, 3, 4);
393 * p = 2*fp+1;
394 * } while (ptest(p, 25) == 0);
395 *
396 * (* find the 2nd Blum prime *)
397 * fq = int((iq-1)/2); (* iq was generated by LavaRnd *)
398 * do {
399 * fq = nextcand(fq+2, 25, 0, 3, 4);
400 * q = 2*fq+1;
401 * } while (ptest(q, 25) == 0);
402 *
403 * (* compute the Blum modulus *)
404 * n[newn] = p * q;
405 *
406 * The pre-defined quadratic residues was also generated by LavaRnd.
407 * The value produced by LavaRnd was squared mod the Blum moduli
408 * that was previously computed.
409 *
410 * The purpose of these pre-defined Blum moduli is to provide users with
411 * an easy way to use a generator where the individual Blum primes used
412 * are not well known. True, these values are in some way "MAGIC", on
413 * the other hand that is their purpose! If this bothers you, don't
414 * use them. See the section "FOR THE PARANOID" below for details.
415 *
416 * The value 'newn' determines which pre-defined generator is used.
417 * For a given 'newn' the Blum modulus 'n[newn]' (product of 2 Blum
418 * (primes) and new quadratic residue 'r[newn]' is set as follows:
419 *
420 * newn == 1: (Blum modulus bit length 130)
421 * n[ 1] = 0x5049440736fe328caf0db722d83de9361
422 * r[ 1] = 0xb226980f11d952e74e5dbb01a4cc42ec
423 *
424 * newn == 2: (Blum modulus bit length 137)
425 * n[ 2] = 0x2c5348a2555dd374a18eb286ea9353443f1
426 * r[ 2] = 0x40f3d643446cd710e3e893616b21e3a218
427 *
428 * newn == 3: (Blum modulus bit length 147)
429 * n[ 3] = 0x9cfd959d6ce4e3a81f1e0f2ca661f11d001f1
430 * r[ 3] = 0xfae5b44d9b64ff5cea4f3e142de2a0d7d76a
431 *
432 * newn == 4: (Blum modulus bit length 157)
433 * n[ 4] = 0x3070f9245c894ed75df12a1a2decc680dfcc0751
434 * r[ 4] = 0x20c2d8131b2bdca2c0af8aa220ddba4b984570
435 *
436 * newn == 5: (Blum modulus bit length 257)
437 * n[ 5] = 0x2109b1822db81a85b38f75aac680bc2fa5d3fe1118769a0108b99e5e799
438 * 166ef1
439 * r[ 5] = 0x5e9b890eae33b792e821a9605f5df6db234f7b7d1e70aeed0e6c77c859e
440 * 2efa9
441 *
442 * newn == 6: (Blum modulus bit length 259)
443 * n[ 6] = 0xa7bfd9d7d9ada2c79f2dbf2185c6440263a38db775ee732dad85557f1e1
444 * ddf431
445 * r[ 6] = 0x5e94a02f88667154e097aedece1c925ce1f3495d2c98eccfc5dc2e80c94
446 * 04daf
447 *
448 * newn == 7: (Blum modulus bit length 286)
449 * n[ 7] = 0x43d87de8f2399ef237801cd5628643fcff569d6b0dcf53ce52882e7f602
450 * f9125cf9ec751
451 * r[ 7] = 0x13522d1ee014c7bfbe90767acced049d876aefcf18d4dd64f0b58c3992d
452 * 2e5098d25e6
453 *
454 * newn == 8: (Blum modulus bit length 294)
455 * n[ 8] = 0x5847126ca7eb4699b7f13c9ce7bdc91fed5bdbd2f99ad4a6c2b59cd9f0b
456 * c42e66a26742f11
457 * r[ 8] = 0x853016dca3269116b7e661fa3d344f9a28e9c9475597b4b8a35da929aae
458 * 95f3a489dc674
459 *
460 * newn == 9: (Blum modulus bit length 533)
461 * n[ 9] = 0x39e8be52322fd3218d923814e81b003d267bb0562157a3c1797b4f4a867
462 * 52a84d895c3e08eb61c36a6ff096061c6fd0fdece0d62b16b66b980f95112
463 * 745db4ab27e3d1
464 * r[ 9] = 0xb458f8ad1e6bbab915bfc01508864b787343bc42a8aa82d9d2880107e3f
465 * d8357c0bd02de3222796b2545e5ab7d81309a89baedaa5d9e8e59f959601e
466 * f2b87d4ed20d
467 *
468 * newn == 10: (Blum modulus bit length 537)
469 * n[10] = 0x25f2435c9055666c23ef596882d7f98bd1448bf23b50e88250d3cc952c8
470 * 1b3ba524a02fd38582de74511c4008d4957302abe36c6092ce222ef9c73cc
471 * 3cdc363b7e64b89
472 * r[10] = 0x66bb7e47b20e0c18401468787e2b707ca81ec9250df8cfc24b5ffbaaf2c
473 * f3008ed8b408d075d56f62c669fadc4f1751baf950d145f40ce23442aee59
474 * 4f5ad494cfc482
475 *
476 * newn == 11: (Blum modulus bit length 542)
477 * n[11] = 0x497864de82bdb3094217d56b874ecd7769a791ea5ec5446757f3f9b6286
478 * e58704499daa2dd37a74925873cfa68f27533920ee1a9a729cf522014dab2
479 * 2e1a530c546ee069
480 * r[11] = 0x8684881cb5e630264a4465ae3af8b69ce3163f806549a7732339eea2c54
481 * d5c590f47fbcedfa07c1ef5628134d918fee5333fed9c094d65461d88b13a
482 * 0aded356e38b04
483 *
484 * newn == 12: (Blum modulus bit length 549)
485 * n[12] = 0x3457582ab3c0ccb15f08b8911665b18ca92bb7c2a12b4a1a66ee4251da1
486 * 90b15934c94e315a1bf41e048c7c7ce812fdd25d653416557d3f09887efad
487 * 2b7f66d151f14c7b99
488 * r[12] = 0xdf719bd1f648ed935870babd55490137758ca3b20add520da4c5e8cdcbf
489 * c4333a13f72a10b604eb7eeb07c573dd2c0208e736fe56ed081aa9488fbc4
490 * 5227dd68e207b4a0
491 *
492 * newn == 13: (Blum modulus bit length 1048)
493 * n[13] = 0x1517c19166b7dd21b5af734ed03d833daf66d82959a553563f4345bd439
494 * 510a7bda8ee0cb6bf6a94286bfd66e49e25678c1ee99ceec891da8b18e843
495 * 7575113aaf83c638c07137fdd3a76c3a49322a11b5a1a84c32d99cbb2b056
496 * 671589917ed14cc7f1b5915f6495dd1892b4ed7417d79a63cc8aaa503a208
497 * e3420cca200323314fc49
498 * r[13] = 0xd42e8e9a560d1263fa648b04f6a69b706d2bc4918c3317ddd162cb4be7a
499 * 5e3bbdd1564a4aadae9fd9f00548f730d5a68dc146f05216fe509f0b8f404
500 * 902692de080bbeda0a11f445ff063935ce78a67445eae5c9cea5a8f6b9883
501 * faeda1bbe5f1ad3ef6409600e2f67b92ed007aba432b567cc26cf3e965e20
502 * 722407bfe46b7736f5
503 *
504 * newn == 14: (Blum modulus bit length 1054)
505 * n[14] = 0x5e56a00e93c6f4e87479ac07b9d983d01f564618b314b4bfec7931eee85
506 * eb909179161e23e78d32110560b22956b22f3bc7e4a034b0586e463fd40c6
507 * f01a33e30ede912acb86a0c1e03483c45f289a271d14bd52792d0a076fdfe
508 * fe32159054b217092237f0767434b3db112fee83005b33f925bacb3185cc4
509 * 409a1abdef8c0fc116af01
510 * r[14] = 0xf7aa7cb67335096ef0c5d09b18f15415b9a564b609913f75f627fc6b0c5
511 * b686c86563fe86134c5a0ea19d243350dfc6b9936ba1512abafb81a0a6856
512 * c9ae7816bf2073c0fb58d8138352b261a704b3ce64d69dee6339010186b98
513 * 3677c84167d4973444194649ad6d71f8fa8f1f1c313edfbbbb6b1b220913c
514 * c8ea47a4db680ff9f190
515 *
516 * newn == 15: (Blum modulus bit length 1055)
517 * n[15] = 0x97dd840b9edfbcdb02c46c175ba81ca845352ebe470be6075326a26770c
518 * ab84bfc0f2e82aa95aac14f40de42a0590445b902c2b8ebb916753e72ab86
519 * c3278cccc1a783b3e962d81b80df03e4380a8fa08b0d86ed0caa515c196a5
520 * 30e49c558ddb53082310b1d0c7aee6f92b619798624ffe6c337299bc51ff5
521 * d2c721061e7597c8d97079
522 * r[15] = 0xb8220703b8c75869ab99f9b50025daa8d77ca6df8cef423ede521f55b1c
523 * 25d74fbf6d6cc31f5ef45e3b29660ef43797f226860a4aa1023dbe522b1fe
524 * 6224d01eb77dee9ad97e8970e4a9e28e7391a6a70557fa0e46eca78866241
525 * ba3c126fc0c5469f8a2f65c33db95d1749d3f0381f401b9201e6abd43d98d
526 * b92e808f0aaa6c3e2110
527 *
528 * newn == 16: (Blum modulus bit length 1062)
529 * n[16] = 0x456e348549b82fbb12b56f84c39f544cb89e43536ae8b2b497d426512c7
530 * f3c9cc2311e0503928284391959e379587bc173e6bc51ba51c856ba557fee
531 * 8dd69cee4bd40845bd34691046534d967e40fe15b6d7cf61e30e283c05be9
532 * 93c44b6a2ea8ade0f5578bd3f618336d9731fed1f1c5996a5828d4ca857ac
533 * 2dc9bd36184183f6d84346e1
534 * r[16] = 0xb0d7dcb19fb27a07973e921a4a4b6dcd7895ae8fced828de8a81a3dbf25
535 * 24def719225404bfd4977a1508c4bac0f3bc356e9d83b9404b5bf86f6d19f
536 * f75645dffc9c5cc153a41772670a5e1ae87a9521416e117a0c0d415fb15d2
537 * 454809bad45d6972f1ab367137e55ad0560d29ada9a2bcda8f4a70fbe04a1
538 * abe4a570605db87b4e8830
539 *
540 * newn == 17: (Blum modulus bit length 2062)
541 * n[17] = 0x6177813aeac0ffa3040b33be3c0f96e0faf97ca54266bfedd7be68494f7
542 * 6a7a91144598bf28b3a5a9dc35a6c9f58d0e5fb19839814bc9d456bff7f29
543 * 953bdac7cafd66e2fc30531b8d544d2720b97025e22b1c71fa0b2eb9a499d
544 * 49484615d07af7a3c23b568531e9b8507543362027ec5ebe0209b4647b7ff
545 * 54be530e9ef50aa819c8ff11f6d7d0a00b25e88f2e6e9de4a7747022b949a
546 * b2c2e1ab0876e2f1177105718c60196f6c3ac0bde26e6cd4e5b8a20e9f0f6
547 * 0974f0b3868ff772ab2ceaf77f328d7244c9ad30e11a2700a120a314aff74
548 * c7f14396e2a39cc14a9fa6922ca0fce40304166b249b574ffd9cbb927f766
549 * c9b150e970a8d1edc24ebf72b72051
550 * r[17] = 0x53720b6eaf3bc3b8adf1dd665324c2d2fc5b2a62f32920c4e167537284d
551 * a802fc106be4b0399caf97519486f31e0fa45a3a677c6cb265c5551ba4a51
552 * 68a7ce3c29731a4e9345eac052ee1b84b7b3a82f906a67aaf7b35949fd7fc
553 * 2f9f4fbc8c18689694c8d30810fff31ebee99b1cf029a33bd736750e7fe0a
554 * 56f7e1d2a9b5321b5117fe9a10e46bf43c896e4a33faebd584f7431e7edbe
555 * bd1703ccee5771b44f0c149888af1a4264cb9cf2e0294ea7719ed6fda1b09
556 * fa6e016c039aeb6d02a03281bcea8c278dd2a807eacae6e52ade048f58f2e
557 * b5193f4ffb9dd68467bc6f8e9d14286bfef09b0aec414c9dadfbf5c46d945
558 * d147b52aa1e0cbd625800522b41dac
559 *
560 * newn == 18: (Blum modulus bit length 2074)
561 * n[18] = 0x68f2a38fb61b42af07cb724fec0c7c65378efcbafb3514e268d7ee38e21
562 * a5680de03f4e63e1e52bde1218f689900be4e5407950539b9d28e9730e8e6
563 * ad6438008aa956b259cd965f3a9d02e1711e6b344b033de6425625b6346d2
564 * ca62e41605e8eae0a7e2f45c25119ef9eece4d3b18369e753419d94118d51
565 * 803842f4de5956b8349e6a0a330145aa4cd1a72afd4ef9db5d8233068e691
566 * 18ff4b93bcc67859f211886bb660033f8170640c6e3d61471c3b7dd62c595
567 * b156d77f317dc272d6b7e7f4fdc20ed82f172fe29776f3bddf697fb673c70
568 * defd6476198a408642ed62081447886a625812ac6576310f23036a7cd3c93
569 * 1c96f7df128ad4ed841351b18c8b78629
570 * r[18] = 0x4735e921f1ac6c3f0d5cda84cd835d75358be8966b99ff5e5d36bdb4be1
571 * 2c5e1df70ac249c0540a99113a8962778dc75dac65af9f3ab4672b4c575c4
572 * 9926f7f3f306fd122ac033961d042c416c3aa43b13ef51b764d505bb1f369
573 * ac7340f8913ddd812e9e75e8fde8c98700e1d3353da18f255e7303db3bcbb
574 * eda4bc5b8d472fbc9697f952cfc243c6f32f3f1bb4541e73ca03f5109df80
575 * 37219a06430e88a6e94be870f8d36dbcc381a1c449c357753a535aa5666db
576 * 92af2aaf1f50a3ddde95024d9161548c263973665a909bd325441a3c18fc7
577 * 0502f2c9a1c944adda164e84a8f3f0230ff2aef8304b5af333077e04920db
578 * a179158f6a2b3afb78df2ef9735ea3c63
579 *
580 * newn == 19: (Blum modulus bit length 2133)
581 * n[19] = 0x230d7ab23bb9e8d6788b252ad6534bdde276540721c3152e410ad4244de
582 * b0df28f4a6de063ba1e51d7cd1736c3d8410e2516b4eb903b8d9206b92026
583 * 64cacbd0425c516833770d118bd5011f3de57e8f607684088255bf7da7530
584 * 56bf373715ed9a7ab85f698b965593fe2b674225fa0a02ebd87402ffb3d97
585 * 172acadaa841664c361f7c11b2af47a472512ee815c970af831f95b737c34
586 * 2508e4c23f3148f3cdf622744c1dcfb69a43fd535e55eebcdc992ee62f2b5
587 * 2c94ac02e0921884fe275b3a528bdb14167b7dec3f3f390cd5a82d80c6c30
588 * 6624cc7a7814fb567cd4d687eede573358f43adfcf1e32f4ee7a2dc4af029
589 * 6435ade8099bf0001d4ae0c7d204df490239c12d6b659a79
590 * r[19] = 0x8f1725f21e245e4fc17982196605b999518b4e21f65126fa6fa759332c8
591 * e27d80158b7537da39d001cc62b83bbef0713b1e82f8293dad522993f86d1
592 * 761015414b2900e74fa23f3eaaa55b31cffd2e801fefb0ac73fd99b5d0cf9
593 * a635c3f4c73d8892d36ad053fc17a423cdcbcf07967a8608c7735e287d784
594 * ae089b3ddea9f2d2bb5d43d2ee25be346832e8dd186fc7a88d82847c03d1c
595 * 05ee52c1f2a51a85f733338547fdbab657cb64b43d44d41148eb32ea68c7e
596 * 66a8d47806f460cd6573b6ca1dd3eeaf1ce8db9621f1e121d2bb4a1878621
597 * dd2dbdd7b5390ab06a5dcd9307d6662eb4248dff2ee263ef2ab778e77724a
598 * 14c62406967daa0d9ad4445064483193d53a5b7698ef473
599 *
600 * newn == 20: (Blum modulus bit length 2166)
601 * n[20] = 0x4fd2b820e0d8b13322e890dddc63a0267e5b3a648b03276066a3f356d79
602 * 660c67704c1be6803b8e7590ee8a962c8331a05778d010e9ba10804d661f3
603 * 354be1932f90babb741bd4302a07a92c42253fd4921864729fb0f0b1e0a42
604 * d66b6777893195abd2ee2141925624bf71ad7328360135c565064ee502773
605 * 6f42a78b988f47407ba4f7996892ffdc5cf9e7ab78ac95734dbf4e3a3def1
606 * 615b5b4341cfbf6c3d0a61b75f4974080bbac03ee9de55221302b40da0c50
607 * ded31d28a2f04921a532b3a486ae36e0bb5273e811d119adf90299a74e623
608 * 3ccce7069676db00a3e8ce255a82fd9748b26546b98c8f4430a8db2a4b230
609 * fa365c51e0985801abba4bbcf3727f7c8765cc914d262fcec3c1d081
610 * r[20] = 0x46ef0184445feaa3099293ee960da14b0f8b046fa9f608241bc08ddeef1
611 * 7ee49194fd9bb2c302840e8da88c4e88df810ce387cc544209ec67656bd1d
612 * a1e9920c7b1aad69448bb58455c9ae4e9cd926911b30d6b5843ff3d306d56
613 * 54a41dc20e2de4eb174ec5ac3e6e70849de5d5f9166961207e2d8b31014cf
614 * 35f801de8372881ae1ba79e58942e5bef0a7e40f46387bf775c54b1d15a14
615 * 40e84beb39cd9e931f5638234ea730ed81d6fca1d7cea9e8ffb171f6ca228
616 * 56264a36a2a783fd7ac39361a6598ed3a565d58acf1f5759bd294e5f53131
617 * bc8e4ee3750794df727b29b1f5788ae14e6a1d1a5b26c2947ed46f49e8377
618 * 3292d7dd5650580faebf85fd126ac98d98f47cf895abdc7ba048bd1a
619 *
620 * NOTE: The Blum moduli associated with 1 <= newn < 9 are subject
621 * to having their Blum moduli factored, depending in their size,
622 * by small PCs in a reasonable to large supercomputers/highly
623 * parallel processors over a long time. Their value lies in their
624 * speed relative the default Blum generator. As of Jan 1997,
625 * the Blum moduli associated with 13 <= newn < 20 appear to
626 * be well beyond the scope of hardware and algorithms,
627 * and 9 <= newn < 12 might be factorable with extreme difficulty.
628 *
629 * The following table may be useful as a guide for how easy it
630 * is to factor the modulus:
631 *
632 * 1 <= newn <= 4 PC using ECM in a short amount of time
633 * 5 <= newn <= 8 Workstation using MPQS in a short amount of time
634 * 8 <= newn <= 12 High end supercomputer or high parallel processor
635 * using state of the art factoring over a long time
636 * 12 <= newn <= 16 Beyond Feb 1997 systems and factoring methods
637 * 17 <= newn <= 20 Well beyond Feb 1997 systems and factoring methods
638 *
639 * See the section titled 'FOR THE PARANOID' for more details.
640 *
641 * seed < 0, 0 < newn <= 20:
642 * -------------------------
643 * Reserved for future use.
644 *
645 ******************************************************************************
646 *
647 * srandom(seed, ip, iq, trials)
648 *
649 * We attempt to set the Blum modulus and quadratic residue.
650 * Any internally buffered random bits are flushed.
651 *
652 * Use the ip and iq args as starting points for Blum primes.
653 * The trials arg determines how many ptest cycles are performed
654 * on each candidate.
655 *
656 * The new quadratic residue is set according to the seed
657 * arg as defined below.
658 *
659 * seed >= 2^32, ip >= 2^16, iq >= 2^16:
660 * -------------------------------------
661 * Set the Blum modulus by searching from the ip and iq search points.
662 *
663 * We will use the seed arg to compute a new quadratic residue.
664 * We will successively square it mod Blum modulus until we get
665 * a smaller value (modulus wrap).
666 *
667 * The follow calc resource file produces an equivalent effect:
668 *
669 * p = nextcand(ip-2, trials, 0, 3, 4); (* find the 1st Blum prime *)
670 * q = nextcand(iq-2, trials, 0, 3, 4); (* find the 2nd Blum prime *)
671 * n = p * q; (* n is the new Blum modulus *)
672 * r = seed;
673 * do {
674 * last_r = r;
675 * r = pmod(last_r, 2, n);
676 * } while (r > last_r); (* r is the new quadratic residue *)
677 * srandom(r, n);
678 *
679 * any seed, ip <= 2^16 or iq <= 2^16:
680 * -----------------------------------
681 * Reserved for future use.
682 *
683 * 0 < seed < 2^32, any ip, any iq:
684 * --------------------------------
685 * Reserved for future use.
686 *
687 * seed == 0, ip > 2^16, iq > 2^16:
688 * --------------------------------
689 * Set the Blum modulus by searching from the ip and iq search points.
690 * If trials is omitted, 1 is assumed.
691 *
692 * The initial quadratic residue will be as if the default initial
693 * quadratic residue arg was given.
694 *
695 * The follow calc resource file produces an equivalent effect:
696 *
697 * srandom(default_residue, ip, iq, trials)
698 *
699 * or in other words:
700 *
701 * (* trials, if omitted, is assumed to be 1 *)
702 * p = nextcand(ip-2, trials, 0, 3, 4); (* find the 1st Blum prime *)
703 * q = nextcand(iq-2, trials, 0, 3, 4); (* find the 2nd Blum prime *)
704 * n = p * q; (* n is the new Blum modulus *)
705 * r = default_residue; (* as used by the initial state *)
706 * do {
707 * last_r = r;
708 * r = pmod(last_r, 2, n);
709 * } while (r > last_r); (* r is the new quadratic residue *)
710 *
711 * seed < 0, any ip, any iq:
712 * -------------------------
713 * Reserved for future use.
714 *
715 ******************************************************************************
716 *
717 * srandom()
718 *
719 * Return current Blum generator state. This call does not alter
720 * the generator state.
721 *
722 ******************************************************************************
723 *
724 * srandom(state)
725 *
726 * Restore the Blum state and return the previous state. Note that
727 * the argument state is a random state value (israndom(state) is true).
728 * Any internally buffered random bits are restored.
729 *
730 * The states of the Blum generators can be saved by calling the seed
731 * function with no arguments, and later restored by calling the seed
732 * functions with that same return value.
733 *
734 * random_state = srandom();
735 * ... generate random bits ...
736 * prev_random_state = srandom(random_state);
737 * ... generate the same random bits ...
738 * srandom() == prev_random_state; (* is true *)
739 *
740 * Saving the state just after seeding a generator and restoring it later
741 * as a very fast way to reseed a generator.
742 */
743
744 /*
745 * TRUTH IN ADVERTISING:
746 *
747 * When the call:
748 *
749 * srandom(seed, nextcand(ip,25,0,3,4)*nextcand(iq,25,0,3,4))
750 *
751 * probable primes from nextcand are used. We use the word
752 * 'probable' because of an extremely extremely small chance that a
753 * composite (a non-prime) may be returned.
754 *
755 * We use the builtin function nextcand in its 5 arg form:
756 *
757 * nextcand(value, 25, 0, 3, 4)
758 *
759 * The odds that a number returned by the above call is not prime is
760 * less than 1 in 4^25. For our purposes, this is sufficient as the
761 * chance of returning a composite is much smaller than the chance that
762 * a hardware glitch will cause nextcand() to return a bogus result.
763 *
764 * Another "truth in advertising" issue is the use of the term
765 * 'pseudo-random'. All deterministic generators are pseudo random.
766 * This is opposed to true random generators based on some special
767 * physical device.
768 *
769 * Even though Blum generator is 'pseudo-random', there is no statistical
770 * test, which runs in polynomial time, that can distinguish the Blum
771 * generator from a truly random source. See the comment under
772 * the "Blum-Blum-Shub generator" section above.
773 *
774 * A final "truth in advertising" issue deals with how the magic numbers
775 * found in this file were generated. Detains can be found in the
776 * various functions, while a overview can be found in the "SOURCE FOR
777 * MAGIC NUMBERS" section below.
778 */
779
780 /*
781 * SOURCE OF MAGIC NUMBERS:
782 *
783 * When seeding the Blum generator, we disallow seeds < 2^32 in an
784 * effort to avoid trivial seed values such as 0, 1 or other small values.
785 * The 2^32 lower bound limit was also selected because it provides a
786 * large reserved value space for special seeds. Currently the
787 * [1,2^32) seed range is reserved for future use.
788 *
789 ***
790 *
791 * When using the 2 arg srandom(seed, newn), we require newn > 2^32
792 * to avoid trivial Blum moduli. The 2^32 lower bound limit was also
793 * selected because it provides a large reserved value space for special
794 * moduli. Currently the [21,2^32) newn range is reserved for future use.
795 *
796 * When using the 3 or 4 arg srandom(seed, ip, iq [, trials]) form,
797 * we require ip>2^16 and ip>2^16. This ensures that the resulting
798 * Blum modulus is > 2^32.
799 *
800 ***
801 *
802 * Taking some care to select a good initial residue helps eliminate cheap
803 * search attacks. It is true that a subsequent residue could be one of the
804 * residues that we would first avoid. However such an occurrence will
805 * happen after the generator is well underway and any such seed information
806 * has been lost.
807 *
808 * In an effort to avoid trivial seed values, we force the seed arg
809 * to srandom() to be > 2^32. We then square this value mod the
810 * Blum modulus until it is less than the previous value. This ensures
811 * that the previous seed value is large enough that its square is > Blum
812 * modulus, and this the square mod Blum modulus is non-trivial.
813 *
814 ***
815 *
816 * The size of default Blum modulus 'n=p*q' was taken to be > 2^259, or
817 * 260 bits (79 digits) long. A modulus > 2^256 will generate 8 bits
818 * per crank of the generator. The period of this generator is long
819 * enough to be reasonable, and the modulus is small enough to be fast.
820 *
821 * The default Blum modulus is not a secure modulus because it can
822 * be factored with ease. As if Feb 1997, the upper reach of the
823 * state of the art for factoring general numbers was around 2^512.
824 * Clearly factoring a 260 bit number if well within the reach of even
825 * a low life Pentium.
826 *
827 * The fact that the default modulus can be factored with ease is
828 * not a drawback. Afterall, if we are going to keep to the goal
829 * of disclosing the source magic numbers, we need to disclose how
830 * the Blum Modulus was produced ... including its factors. Knowing
831 * the facotrs of the Blum modulus does not reduce its quality,
832 * only the ability for someone to determine where you are in the
833 * sequence. But heck, the default seed is well known anyway so
834 * there is no real loss if the factors are also known.
835 *
836 ***
837 *
838 * The default Blum modulus is the product of two Blum probable primes
839 * that were selected by the Rand Book of Random Numbers. Using the '6% rule',
840 * a default Blum modulus n=p*q > 2^256 would be satisfied if p were
841 * 38 decimal digits and q were 42 decimal digits in length. We restate
842 * the sizes in decimal digits because the Rand Book of Random Numbers is a
843 * book of decimal digits. Using the first 38 rand digits as a
844 * starting search point for 'p', and the next 42 digits for a starting
845 * search point for 'q'.
846 *
847 * (*
848 * * setup the search points (lines split for readability)
849 * *)
850 * ip = 10097325337652013586346735487680959091;
851 * iq = 173929274945375420480564894742962480524037;
852 *
853 * (*
854 * * find the first Blum prime
855 * *)
856 * fp = int((ip-1)/2);
857 * do {
858 * fp = nextcand(fp+2, 25, 0, 3, 4);
859 * p = 2*fp+1;
860 * } while (ptest(p, 25) == 0);
861 *
862 * (*
863 * * find the 2nd Blum prime
864 * *)
865 * fq = int((iq-1)/2);
866 * do {
867 * fq = nextcand(fq+2, 25, 0, 3, 4);
868 * q = 2*fq+1;
869 * } while (ptest(q, 25) == 0);
870 *
871 * The above resource file produces the Blum probable primes and initial
872 * quadratic residue (line wrapped for readability):
873 *
874 * p= 0x798ac934c7a3318ad446190f3474e57
875 *
876 * q= 0x1ff21d7e1dd7d5965e224d485d84c3ef44f
877 *
878 * These Blum primes were found after 1.81s of CPU time on a 195 Mhz IP28
879 * R10000 version 2.5 processor. The first Blum prime 'p' was 31716 higher
880 * than the initial search value 'ip'. The second Blum prime 'q' was 18762
881 * higher than the initial starting 'iq'.
882 *
883 * The product of the two Blum primes results in a 260 bit Blum modulus of:
884 *
885 * n = 0xf2ac1903156af9e373d78613ed0e8d30284f34b644a9027d9ba55a689d6be18d9
886 *
887 * The selection if the initial quadratic residue comes from the next
888 * unused digits of the Rand Book of Random Numbers. Now the two initial
889 * search values 'ip' and 'iq' used above needed the first 38 digits and
890 * the next 42 digits. Thus we will skip the first 38+42=80 digits
891 * and begin to build in initial search value for a quadratic residue (most
892 * significant digit first) from the Rand Book of Numbers digits until we
893 * have a value whose square mod n > 4th power mod n. In other words, we
894 * need to build ir up from new Rand Book of Random Numbers digits until
895 * we find a value in which srandom(ir), for the Blum Modulus 'n' produces
896 * an initial quadratic residue on the first loop.
897 *
898 * Clearly we need to find an ir that is > sqrt(n). The first ir:
899 *
900 * ir = 2063610402008229166508422689531964509303
901 *
902 * fails the single loop criteria. So we add the next digit:
903 *
904 * ir = 20636104020082291665084226895319645093032
905 *
906 * Here we find that:
907 *
908 * pmod(ir,2,n) > pmod(pmod(ir,2,n),2,n)
909 *
910 * Thus, for the Blum modulus 'n', the method outlined for srandom(ir) yields
911 * the initial quadratic residue of:
912 *
913 * r = 0x748b6d882ff4b074e2f1e93a8627d626506c73ca5a62546c90f23fd7ed3e7b11e
914 *
915 ***
916 *
917 * In the above process of finding the Blum primes used in the default
918 * Blum modulus, we selected primes of the form:
919 *
920 * 2*x + 1 x is also prime
921 *
922 * because Blum generators with modulus 'n=p*q' have a period:
923 *
924 * lambda(n) = lcm(factors of p-1 & q-1)
925 *
926 * since 'p' and 'q' are both odd, 'p-1' and 'q-1' have 2 as
927 * a factor. The calc resource file above ensures that '(p-1)/2' and
928 * '(q-1)/2' are probable prime, thus maximizing the period
929 * of the default generator to:
930 *
931 * lambda(n) = lcm(2,2,fp,fq) = 2*fp*fq = ~2*(p/2)*(q/2) = ~n/2
932 *
933 * The process above resulted in a default generator Blum modulus n > 2^259
934 * with period of at least 2^258 bits. To be exact, the period of the
935 * default Blum generator is:
936 *
937 * 0x79560c818ab57cf1b9ebc309f68746881adc15e79c05e476f741e5f904b9beb1a
938 *
939 * which is approximately:
940 *
941 * ~8.781 * 10^77
942 *
943 * This period is more than long enough for computationally tractable tasks.
944 *
945 ***
946 *
947 * The 20 builtin generators, on the other hand, were selected
948 * with more care. While the lower order 20 generators have
949 * factorable generators, the higher order are likely to be
950 * be beyond the reach for a while.
951 *
952 * The lengths of the two Blum probable primes 'p' and 'q' used to make up
953 * the 20 Blum moduli 'n=p*q' differ slightly to avoid certain
954 * factorization attacks that work on numbers that are a perfect square,
955 * or where the two primes are nearly the same. I elected to have the
956 * sizes differ by up to 6% of the product size to avoid such attacks.
957 * Clearly one does not want the size of the two factors to differ
958 * by a large percentage: p=3 and q large would result in a easy
959 * to factor Blum modulus. Thus we select sizes that differ by
960 * up to 6% but not (significantly) greater than 6%.
961 *
962 * Using the '6% rule' above, a Blum modulus n=p*q > 2^1024 would have two
963 * Blum factors p > 2^482 and q > 2^542. This is because 482+542 = 1024.
964 * The difference 542-482 is ~5.86% of 1024, and is the largest difference
965 * that is < 6%.
966 *
967 ******************************************************************************
968 *
969 * FOR THE PARANOID:
970 *
971 * The truly paranoid might suggest that my claims in the MAGIC NUMBERS
972 * section are a lie intended to entrap people. Well they are not, but
973 * you need not take my word for it.
974 *
975 ***
976 *
977 * One could take issue with the above resource file that produced a 260 bit
978 * Blum modulus. So if that bothers you, then seed your generator
979 * with your own Blum modulus and initial quadratic residue. And
980 * if you are truly paranoid, why would you want to use the default seed,
981 * which is well known?
982 *
983 ***
984 *
985 * If all the above fails to pacify the truly paranoid, then one may
986 * select your own modulus and initial quadratic residue by calling:
987 *
988 * srandom(r, n);
989 *
990 * Of course, you will need to select a correct Blum modulus 'n' as the
991 * product of two Blum primes (both 3 mod 4) and with a long period (where
992 * lcm(factors of one less than the two Blum primes) is large) and an
993 * initial quadratic residue 'r' that is hard to guess selected at random.
994 *
995 * A simple way to seed the generator would be to:
996 *
997 * srandom(ir, ip, iq, 25)
998 *
999 * where 'ip', 'iq' and 'ir' are large integers that are unlikely to be
1000 * 'guessed' and where numbers around the size of iq*ir are beyond
1001 * the current reach of the best factoring methods on the fastest
1002 * SGI/Cray supercomputers.
1003 *
1004 * Of course you can increase the '25' value if 1 of 4^25 odds of a
1005 * non-prime are too probable for you.
1006 *
1007 * The problem with the above call is that the period of the Blum generator
1008 * could be small if 'p' and 'q' have lots of small prime factors in common.
1009 *
1010 * A better way to do seed the Blum generator yourself is to use the
1011 * seedrandom(seed1, seed2, size [,tests]) function from "seedrandom.cal"
1012 * with the args:
1013 *
1014 * seed1 - seed rand() to search for 'p', select from [2^64, 2^308)
1015 * seed2 - seed rand() to search for 'q', select from [2^64, 2^308)
1016 * size - minimum bit size of the Blum modulus 'n=p*q'
1017 * tests - optional arg for number of pseudo prime tests (default is 25)
1018 *
1019 * Last, one could use some external random source to select starting
1020 * search points for 'p', 'q' and the quadratic residue. One way to
1021 * do this is:
1022 *
1023 * fp = int((ip-1)/2);
1024 * do {
1025 * fp = nextcand(fp+2, tests, 0, 3, 4);
1026 * p = 2*fp+1;
1027 * } while (ptest(p, tests) == 0);
1028 * fq = int((iq-1)/2);
1029 * do {
1030 * fq = nextcand(fq+2, tests, 0, 3, 4);
1031 * q = 2*fq+1;
1032 * } while (ptest(q, tests) == 0);
1033 * srandom(pmod(ir,2,p*q), p*q);
1034 *
1035 * where 'tests' is the number of pseudo prime tests that a candidate must
1036 * pass before being considered a probable prime (must be >0, perhaps 25), and
1037 * where 'ip' is the initial search location for the Blum prime 'p', and
1038 * where 'iq' is the initial search location for the Blum prime 'q', and
1039 * where 'ir' is the initial Blum quadratic residue generator. The 'ir'
1040 * value should be a random value in the range [2^(binsize*4/5), 2^(binsize-2))
1041 * where 2^(binsize-1) < n=p*q <= 2^binsize.
1042 *
1043 * Your external generator would need to generate 'ip', 'iq' and 'ir'.
1044 * While any value for 'ip' and 'iq will do (provided that their product
1045 * is large enough to meet your modulus needs), 'ir' should be selected
1046 * to avoid values near 0 or near 'n' (or ip*iq).
1047 *
1048 ***
1049 *
1050 * The Blum moduli used with the pre-defined generators (via the call
1051 * srandom(seed, 0<x<=16)) were generated using the above process. The
1052 * 'ip', 'iq' and 'ir' values were produced by special purpose hardware
1053 * that produced cryptographically strong random numbers. The Blum
1054 * primes used in these special pre-defined generators are unknown.
1055 *
1056 * Not being able to factor 'n=p*q' into 'p' and 'q' does not directly
1057 * improve the quality Blum generator. On the other hand, it does
1058 * improve the security of it.
1059 *
1060 * I (Landon Curt Noll) did not keep the search values of these 20 special
1061 * pre-defined generators. While some of the smaller Blum moduli is
1062 * within the range of some factoring methods, others are not. As of
1063 * Feb 1997, the following is the estimate of what can factor the
1064 * pre-defined moduli:
1065 *
1066 * 1 <= newn < 4 PC using ECM in a short amount of time
1067 * 5 <= newn < 8 Workstation using MPQS in a short amount of time
1068 * 8 <= newn < 12 High end supercomputer or high parallel processor
1069 * using state of the art factoring over a long time
1070 * 12 <= newn < 16 Beyond Feb 1997 systems and factoring methods
1071 * 17 <= newn < 20 Well beyond Feb 1997 systems and factoring methods
1072 *
1073 * This is not to say that in the future things will not change. One
1074 * can say that faster hardware will not help in the factoring of 1024+
1075 * bit Blum moduli found in 12 <= newn as well as in the default
1076 * Blum generator. A combination of better algorithms, helped by faster
1077 * hardware will be needed.
1078 *
1079 * It is true that the pre-defined moduli are 'magic'. On the other
1080 * hand, there purpose was in part to produce users with pre-seeded
1081 * generators where the individual Blum primes are not well known. If
1082 * this bothers you, don't use them and seed your own!
1083 *
1084 ***
1085 *
1086 * The output of the Blum generator has been proven to be cryptographically
1087 * strong. I.e., there is absolutely no better way to predict the next
1088 * bit in the sequence than by tossing a coin (as with TRULY random numbers)
1089 * EVEN IF YOU HAVE A LARGE NUMBER OF PREVIOUSLY GENERATED BITS AND KNOW
1090 * A LARGE NUMBER OF BITS THAT FOLLOW THE NEXT BIT! An adversary would
1091 * be far better advised to try to factor the modulus or exhaustively search
1092 * for the quadratic residue in use.
1093 *
1094 */
1095
1096
1097 #include "zrandom.h"
1098 #include "have_const.h"
1099 #include "have_unused.h"
1100
1101
1102 #include "banned.h" /* include after system header <> includes */
1103
1104
1105 /*
1106 * current Blum generator state
1107 */
1108 STATIC RANDOM blum;
1109
1110
1111 /*
1112 * static constants 3 and 4 used by zsrandom4
1113 */
1114 STATIC HALF _threeval_[] = { 3 };
1115 STATIC ZVALUE _three_ = { _threeval_, 1, 0 };
1116 STATIC HALF _fourval_[] = { 4 };
1117 STATIC ZVALUE _four_ = { _fourval_, 1, 0 };
1118
1119
1120 /*
1121 * Default Blum generator
1122 *
1123 * The init_blum generator is established via the srandom(0) call.
1124 *
1125 * The z_rdefault ZVALUE is the 'r' (quadratic residue) of init_blum.
1126 */
1127 #if FULL_BITS == 64
1128 STATIC CONST HALF h_ndefvec[] = {
1129 (HALF)0xd6be18d9, (HALF)0xba55a689, (HALF)0x4a9027d9, (HALF)0x84f34b64,
1130 (HALF)0xd0e8d302, (HALF)0x3d78613e, (HALF)0x56af9e37, (HALF)0x2ac19031,
1131 (HALF)0xf
1132 };
1133 STATIC CONST HALF h_rdefvec[] = {
1134 (HALF)0xd3e7b11e, (HALF)0x0f23fd7e, (HALF)0xa62546c9, (HALF)0x06c73ca5,
1135 (HALF)0x627d6265, (HALF)0x2f1e93a8, (HALF)0xff4b074e, (HALF)0x48b6d882,
1136 (HALF)0x7
1137 };
1138 #elif 2*FULL_BITS == 64
1139 STATIC CONST HALF h_ndefvec[] = {
1140 (HALF)0x18d9, (HALF)0xd6be, (HALF)0xa689, (HALF)0xba55,
1141 (HALF)0x27d9, (HALF)0x4a90, (HALF)0x4b64, (HALF)0x84f3,
1142 (HALF)0xd302, (HALF)0xd0e8, (HALF)0x613e, (HALF)0x3d78,
1143 (HALF)0x9e37, (HALF)0x56af, (HALF)0x9031, (HALF)0x2ac1,
1144 (HALF)0xf
1145 };
1146 STATIC CONST HALF h_rdefvec[] = {
1147 (HALF)0xb11e, (HALF)0xd3e7, (HALF)0xfd7e, (HALF)0x0f23,
1148 (HALF)0x46c9, (HALF)0xa625, (HALF)0x3ca5, (HALF)0x06c7,
1149 (HALF)0x6265, (HALF)0x627d, (HALF)0x93a8, (HALF)0x2f1e,
1150 (HALF)0x074e, (HALF)0xff4b, (HALF)0xd882, (HALF)0x48b6,
1151 (HALF)0x7
1152 };
1153 #else
1154 /\../\ FULL_BITS must be 32 or 64 /\../\ !!!
1155 #endif
1156 STATIC CONST RANDOM init_blum = {1, 0, 8, (HALF)0, (HALF)0xff,
1157 {(HALF *)h_ndefvec,
1158 sizeof(h_ndefvec)/sizeof(HALF), 0},
1159 {(HALF *)h_rdefvec,
1160 sizeof(h_rdefvec)/sizeof(HALF), 0}};
1161 #if FULL_BITS == 64
1162 STATIC CONST HALF h_rdefvec_2[] = {
1163 (HALF)0xd3e7b11e, (HALF)0x0f23fd7e, (HALF)0xa62546c9, (HALF)0x06c73ca5,
1164 (HALF)0x627d6265, (HALF)0x2f1e93a8, (HALF)0xff4b074e, (HALF)0x48b6d882,
1165 (HALF)0x7
1166 };
1167 #elif 2*FULL_BITS == 64
1168 STATIC CONST HALF h_rdefvec_2[] = {
1169 (HALF)0xb11e, (HALF)0xd3e7, (HALF)0xfd7e, (HALF)0x0f23,
1170 (HALF)0x46c9, (HALF)0xa625, (HALF)0x3ca5, (HALF)0x06c7,
1171 (HALF)0x6265, (HALF)0x627d, (HALF)0x93a8, (HALF)0x2f1e,
1172 (HALF)0x074e, (HALF)0xff4b, (HALF)0xd882, (HALF)0x48b6,
1173 (HALF)0x7
1174 };
1175 #else
1176 /\../\ FULL_BITS must be 32 or 64 /\../\ !!!
1177 #endif
1178 STATIC CONST ZVALUE z_rdefault = {
1179 (HALF *)h_rdefvec_2, sizeof(h_rdefvec_2)/sizeof(HALF), 0
1180 };
1181
1182
1183 /*
1184 * Pre-defined Blum generators
1185 *
1186 * These generators are seeded via the srandom(0, newn) where
1187 * 1 <= newn < BLUM_PREGEN.
1188 */
1189 #if FULL_BITS == 64
1190 STATIC CONST HALF h_nvec01[] = {
1191 (HALF)0x83de9361, (HALF)0xf0db722d, (HALF)0x6fe328ca, (HALF)0x04944073,
1192 (HALF)0x5
1193 };
1194 STATIC CONST HALF h_rvec01[] = {
1195 (HALF)0xa4cc42ec, (HALF)0x4e5dbb01, (HALF)0x11d952e7, (HALF)0xb226980f
1196 };
1197 STATIC CONST HALF h_nvec02[] = {
1198 (HALF)0x353443f1, (HALF)0xeb286ea9, (HALF)0xdd374a18, (HALF)0x348a2555,
1199 (HALF)0x2c5
1200 };
1201 STATIC CONST HALF h_rvec02[] = {
1202 (HALF)0x21e3a218, (HALF)0xe893616b, (HALF)0x6cd710e3, (HALF)0xf3d64344,
1203 (HALF)0x40
1204 };
1205 STATIC CONST HALF h_nvec03[] = {
1206 (HALF)0x11d001f1, (HALF)0xf2ca661f, (HALF)0x3a81f1e0, (HALF)0x59d6ce4e,
1207 (HALF)0x0009cfd9
1208 };
1209 STATIC CONST HALF h_rvec03[] = {
1210 (HALF)0xa0d7d76a, (HALF)0x3e142de2, (HALF)0xff5cea4f, (HALF)0xb44d9b64,
1211 (HALF)0xfae5
1212 };
1213 STATIC CONST HALF h_nvec04[] = {
1214 (HALF)0xdfcc0751, (HALF)0x2decc680, (HALF)0x5df12a1a, (HALF)0x5c894ed7,
1215 (HALF)0x3070f924
1216 };
1217 STATIC CONST HALF h_rvec04[] = {
1218 (HALF)0x4b984570, (HALF)0xa220ddba, (HALF)0xa2c0af8a, (HALF)0x131b2bdc,
1219 (HALF)0x0020c2d8
1220 };
1221 STATIC CONST HALF h_nvec05[] = {
1222 (HALF)0x99166ef1, (HALF)0x8b99e5e7, (HALF)0x8769a010, (HALF)0x5d3fe111,
1223 (HALF)0x680bc2fa, (HALF)0x38f75aac, (HALF)0xdb81a85b, (HALF)0x109b1822,
1224 (HALF)0x2
1225 };
1226 STATIC CONST HALF h_rvec05[] = {
1227 (HALF)0x59e2efa9, (HALF)0x0e6c77c8, (HALF)0x1e70aeed, (HALF)0x234f7b7d,
1228 (HALF)0x5f5df6db, (HALF)0xe821a960, (HALF)0xae33b792, (HALF)0x5e9b890e
1229 };
1230 STATIC CONST HALF h_nvec06[] = {
1231 (HALF)0xe1ddf431, (HALF)0xd85557f1, (HALF)0x5ee732da, (HALF)0x3a38db77,
1232 (HALF)0x5c644026, (HALF)0xf2dbf218, (HALF)0x9ada2c79, (HALF)0x7bfd9d7d,
1233 (HALF)0xa
1234 };
1235 STATIC CONST HALF h_rvec06[] = {
1236 (HALF)0xc9404daf, (HALF)0xc5dc2e80, (HALF)0x2c98eccf, (HALF)0xe1f3495d,
1237 (HALF)0xce1c925c, (HALF)0xe097aede, (HALF)0x88667154, (HALF)0x5e94a02f
1238 };
1239 STATIC CONST HALF h_nvec07[] = {
1240 (HALF)0xcf9ec751, (HALF)0x602f9125, (HALF)0x52882e7f, (HALF)0x0dcf53ce,
1241 (HALF)0xff569d6b, (HALF)0x628643fc, (HALF)0x37801cd5, (HALF)0xf2399ef2,
1242 (HALF)0x43d87de8
1243 };
1244 STATIC CONST HALF h_rvec07[] = {
1245 (HALF)0x098d25e6, (HALF)0x3992d2e5, (HALF)0x64f0b58c, (HALF)0xcf18d4dd,
1246 (HALF)0x9d876aef, (HALF)0x7acced04, (HALF)0xbfbe9076, (HALF)0x1ee014c7,
1247 (HALF)0x0013522d
1248 };
1249 STATIC CONST HALF h_nvec08[] = {
1250 (HALF)0x26742f11, (HALF)0xbc42e66a, (HALF)0xb59cd9f0, (HALF)0x9ad4a6c2,
1251 (HALF)0x5bdbd2f9, (HALF)0xbdc91fed, (HALF)0xf13c9ce7, (HALF)0xeb4699b7,
1252 (HALF)0x47126ca7, (HALF)0x58
1253 };
1254 STATIC CONST HALF h_rvec08[] = {
1255 (HALF)0x489dc674, (HALF)0xaae95f3a, (HALF)0xa35da929, (HALF)0x5597b4b8,
1256 (HALF)0x28e9c947, (HALF)0x3d344f9a, (HALF)0xb7e661fa, (HALF)0xa3269116,
1257 (HALF)0x853016dc
1258 };
1259 STATIC CONST HALF h_nvec09[] = {
1260 (HALF)0xab27e3d1, (HALF)0x12745db4, (HALF)0xb980f951, (HALF)0x62b16b66,
1261 (HALF)0x0fdece0d, (HALF)0x6061c6fd, (HALF)0x36a6ff09, (HALF)0xe08eb61c,
1262 (HALF)0x84d895c3, (HALF)0x4a86752a, (HALF)0xc1797b4f, (HALF)0x562157a3,
1263 (HALF)0x3d267bb0, (HALF)0x14e81b00, (HALF)0x218d9238, (HALF)0x52322fd3,
1264 (HALF)0x0039e8be
1265 };
1266 STATIC CONST HALF h_rvec09[] = {
1267 (HALF)0x7d4ed20d, (HALF)0x601ef2b8, (HALF)0x8e59f959, (HALF)0xedaa5d9e,
1268 (HALF)0x309a89ba, (HALF)0xe5ab7d81, (HALF)0x796b2545, (HALF)0x02de3222,
1269 (HALF)0x8357c0bd, (HALF)0x0107e3fd, (HALF)0x82d9d288, (HALF)0xbc42a8aa,
1270 (HALF)0x4b787343, (HALF)0xc0150886, (HALF)0xbab915bf, (HALF)0xf8ad1e6b,
1271 (HALF)0xb458
1272 };
1273 STATIC CONST HALF h_nvec10[] = {
1274 (HALF)0xb7e64b89, (HALF)0xc3cdc363, (HALF)0x2ef9c73c, (HALF)0x6092ce22,
1275 (HALF)0x02abe36c, (HALF)0x08d49573, (HALF)0x74511c40, (HALF)0xd38582de,
1276 (HALF)0xa524a02f, (HALF)0x52c81b3b, (HALF)0x250d3cc9, (HALF)0x23b50e88,
1277 (HALF)0xbd1448bf, (HALF)0x882d7f98, (HALF)0xc23ef596, (HALF)0xc9055666,
1278 (HALF)0x025f2435
1279 };
1280 STATIC CONST HALF h_rvec10[] = {
1281 (HALF)0x94cfc482, (HALF)0x594f5ad4, (HALF)0x23442aee, (HALF)0x145f40ce,
1282 (HALF)0x1baf950d, (HALF)0xadc4f175, (HALF)0xf62c669f, (HALF)0x8d075d56,
1283 (HALF)0x08ed8b40, (HALF)0xaaf2cf30, (HALF)0xc24b5ffb, (HALF)0x250df8cf,
1284 (HALF)0x7ca81ec9, (HALF)0x787e2b70, (HALF)0x18401468, (HALF)0x47b20e0c,
1285 (HALF)0x0066bb7e
1286 };
1287 STATIC CONST HALF h_nvec11[] = {
1288 (HALF)0x546ee069, (HALF)0x2e1a530c, (HALF)0x2014dab2, (HALF)0xa729cf52,
1289 (HALF)0x920ee1a9, (HALF)0x68f27533, (HALF)0x25873cfa, (HALF)0xdd37a749,
1290 (HALF)0x4499daa2, (HALF)0x286e5870, (HALF)0x57f3f9b6, (HALF)0x5ec54467,
1291 (HALF)0x69a791ea, (HALF)0x874ecd77, (HALF)0x4217d56b, (HALF)0x82bdb309,
1292 (HALF)0x497864de
1293 };
1294 STATIC CONST HALF h_rvec11[] = {
1295 (HALF)0x56e38b04, (HALF)0x3a0aded3, (HALF)0x461d88b1, (HALF)0x9c094d65,
1296 (HALF)0xe5333fed, (HALF)0x34d918fe, (HALF)0x1ef56281, (HALF)0xcedfa07c,
1297 (HALF)0x590f47fb, (HALF)0xa2c54d5c, (HALF)0x732339ee, (HALF)0x806549a7,
1298 (HALF)0x9ce3163f, (HALF)0xae3af8b6, (HALF)0x264a4465, (HALF)0x1cb5e630,
1299 (HALF)0x00868488
1300 };
1301 STATIC CONST HALF h_nvec12[] = {
1302 (HALF)0xf14c7b99, (HALF)0x7f66d151, (HALF)0x87efad2b, (HALF)0x57d3f098,
1303 (HALF)0xd6534165, (HALF)0x812fdd25, (HALF)0x48c7c7ce, (HALF)0xa1bf41e0,
1304 (HALF)0x4c94e315, (HALF)0x190b1593, (HALF)0xee4251da, (HALF)0x2b4a1a66,
1305 (HALF)0x2bb7c2a1, (HALF)0x65b18ca9, (HALF)0x08b89116, (HALF)0xc0ccb15f,
1306 (HALF)0x57582ab3, (HALF)0x34
1307 };
1308 STATIC CONST HALF h_rvec12[] = {
1309 (HALF)0xe207b4a0, (HALF)0x5227dd68, (HALF)0x9488fbc4, (HALF)0x6ed081aa,
1310 (HALF)0x8e736fe5, (HALF)0x3dd2c020, (HALF)0xeeb07c57, (HALF)0x0b604eb7,
1311 (HALF)0xa13f72a1, (HALF)0xcbfc4333, (HALF)0xa4c5e8cd, (HALF)0x0add520d,
1312 (HALF)0x758ca3b2, (HALF)0x55490137, (HALF)0x5870babd, (HALF)0xf648ed93,
1313 (HALF)0xdf719bd1
1314 };
1315 STATIC CONST HALF h_nvec13[] = {
1316 (HALF)0x3314fc49, (HALF)0xcca20032, (HALF)0x208e3420, (HALF)0x8aaa503a,
1317 (HALF)0xd79a63cc, (HALF)0xb4ed7417, (HALF)0x95dd1892, (HALF)0xb5915f64,
1318 (HALF)0xd14cc7f1, (HALF)0x1589917e, (HALF)0xb2b05667, (HALF)0xc32d99cb,
1319 (HALF)0x1b5a1a84, (HALF)0xa49322a1, (HALF)0xdd3a76c3, (HALF)0x8c07137f,
1320 (HALF)0xaaf83c63, (HALF)0x37575113, (HALF)0xa8b18e84, (HALF)0xceec891d,
1321 (HALF)0x78c1ee99, (HALF)0x6e49e256, (HALF)0x4286bfd6, (HALF)0xcb6bf6a9,
1322 (HALF)0x7bda8ee0, (HALF)0xd439510a, (HALF)0x63f4345b, (HALF)0x959a5535,
1323 (HALF)0xdaf66d82, (HALF)0xed03d833, (HALF)0x1b5af734, (HALF)0x166b7dd2,
1324 (HALF)0x01517c19
1325 };
1326 STATIC CONST HALF h_rvec13[] = {
1327 (HALF)0x6b7736f5, (HALF)0x2407bfe4, (HALF)0x965e2072, (HALF)0xcc26cf3e,
1328 (HALF)0xa432b567, (HALF)0x2ed007ab, (HALF)0x0e2f67b9, (HALF)0xef640960,
1329 (HALF)0xbe5f1ad3, (HALF)0x3faeda1b, (HALF)0xa8f6b988, (HALF)0xe5c9cea5,
1330 (HALF)0xa67445ea, (HALF)0x3935ce78, (HALF)0xf445ff06, (HALF)0xbeda0a11,
1331 (HALF)0x92de080b, (HALF)0xf4049026, (HALF)0xe509f0b8, (HALF)0x6f05216f,
1332 (HALF)0x5a68dc14, (HALF)0x548f730d, (HALF)0xe9fd9f00, (HALF)0x64a4aada,
1333 (HALF)0xe3bbdd15, (HALF)0xcb4be7a5, (HALF)0x17ddd162, (HALF)0xc4918c33,
1334 (HALF)0x9b706d2b, (HALF)0x8b04f6a6, (HALF)0x1263fa64, (HALF)0x8e9a560d,
1335 (HALF)0xd42e
1336 };
1337 STATIC CONST HALF h_nvec14[] = {
1338 (HALF)0xc116af01, (HALF)0xbdef8c0f, (HALF)0xc4409a1a, (HALF)0xacb3185c,
1339 (HALF)0xb33f925b, (HALF)0xfee83005, (HALF)0x4b3db112, (HALF)0x7f076743,
1340 (HALF)0x21709223, (HALF)0x2159054b, (HALF)0x6fdfefe3, (HALF)0x792d0a07,
1341 (HALF)0x1d14bd52, (HALF)0x5f289a27, (HALF)0xe03483c4, (HALF)0xcb86a0c1,
1342 (HALF)0x0ede912a, (HALF)0xf01a33e3, (HALF)0x63fd40c6, (HALF)0x4b0586e4,
1343 (HALF)0xbc7e4a03, (HALF)0x956b22f3, (HALF)0x10560b22, (HALF)0x3e78d321,
1344 (HALF)0x179161e2, (HALF)0xe85eb909, (HALF)0xec7931ee, (HALF)0xb314b4bf,
1345 (HALF)0x1f564618, (HALF)0xb9d983d0, (HALF)0x7479ac07, (HALF)0x93c6f4e8,
1346 (HALF)0x5e56a00e
1347 };
1348 STATIC CONST HALF h_rvec14[] = {
1349 (HALF)0x0ff9f190, (HALF)0x47a4db68, (HALF)0x913cc8ea, (HALF)0xb6b1b220,
1350 (HALF)0x13edfbbb, (HALF)0xa8f1f1c3, (HALF)0xd6d71f8f, (HALF)0x4194649a,
1351 (HALF)0x7d497344, (HALF)0x677c8416, (HALF)0x0186b983, (HALF)0xee633901,
1352 (HALF)0xce64d69d, (HALF)0x61a704b3, (HALF)0x138352b2, (HALF)0xc0fb58d8,
1353 (HALF)0x16bf2073, (HALF)0x56c9ae78, (HALF)0xb81a0a68, (HALF)0x1512abaf,
1354 (HALF)0x6b9936ba, (HALF)0x43350dfc, (HALF)0xa0ea19d2, (HALF)0xe86134c5,
1355 (HALF)0x6c86563f, (HALF)0x6b0c5b68, (HALF)0x75f627fc, (HALF)0xb609913f,
1356 (HALF)0x15b9a564, (HALF)0x9b18f154, (HALF)0x6ef0c5d0, (HALF)0xb6733509,
1357 (HALF)0x00f7aa7c
1358 };
1359 STATIC CONST HALF h_nvec15[] = {
1360 (HALF)0xc8d97079, (HALF)0x061e7597, (HALF)0xf5d2c721, (HALF)0x299bc51f,
1361 (HALF)0xffe6c337, (HALF)0x19798624, (HALF)0xee6f92b6, (HALF)0x0b1d0c7a,
1362 (HALF)0xb5308231, (HALF)0x49c558dd, (HALF)0x196a530e, (HALF)0x0caa515c,
1363 (HALF)0x8b0d86ed, (HALF)0x380a8fa0, (HALF)0x80df03e4, (HALF)0xe962d81b,
1364 (HALF)0xc1a783b3, (HALF)0xc3278ccc, (HALF)0x3e72ab86, (HALF)0xebb91675,
1365 (HALF)0xb902c2b8, (HALF)0xa0590445, (HALF)0x4f40de42, (HALF)0xaa95aac1,
1366 (HALF)0xfc0f2e82, (HALF)0x70cab84b, (HALF)0x5326a267, (HALF)0x470be607,
1367 (HALF)0x45352ebe, (HALF)0x5ba81ca8, (HALF)0x02c46c17, (HALF)0x9edfbcdb,
1368 (HALF)0x97dd840b
1369 };
1370 STATIC CONST HALF h_rvec15[] = {
1371 (HALF)0x6c3e2110, (HALF)0x808f0aaa, (HALF)0xd98db92e, (HALF)0x1e6abd43,
1372 (HALF)0xf401b920, (HALF)0x9d3f0381, (HALF)0xdb95d174, (HALF)0xa2f65c33,
1373 (HALF)0x0c5469f8, (HALF)0xa3c126fc, (HALF)0x8866241b, (HALF)0x0e46eca7,
1374 (HALF)0xa70557fa, (HALF)0x8e7391a6, (HALF)0x70e4a9e2, (HALF)0x9ad97e89,
1375 (HALF)0x1eb77dee, (HALF)0xfe6224d0, (HALF)0xdbe522b1, (HALF)0xa4aa1023,
1376 (HALF)0x7f226860, (HALF)0x60ef4379, (HALF)0x45e3b296, (HALF)0xcc31f5ef,
1377 (HALF)0x74fbf6d6, (HALF)0x55b1c25d, (HALF)0x3ede521f, (HALF)0xdf8cef42,
1378 (HALF)0xa8d77ca6, (HALF)0xb50025da, (HALF)0x69ab99f9, (HALF)0x03b8c758,
1379 (HALF)0x00b82207
1380 };
1381 STATIC CONST HALF h_nvec16[] = {
1382 (HALF)0xd84346e1, (HALF)0x184183f6, (HALF)0x2dc9bd36, (HALF)0x4ca857ac,
1383 (HALF)0x96a5828d, (HALF)0xed1f1c59, (HALF)0x36d9731f, (HALF)0xbd3f6183,
1384 (HALF)0xde0f5578, (HALF)0xb6a2ea8a, (HALF)0xbe993c44, (HALF)0x0e283c05,
1385 (HALF)0xd7cf61e3, (HALF)0x40fe15b6, (HALF)0x534d967e, (HALF)0x34691046,
1386 (HALF)0xd40845bd, (HALF)0xd69cee4b, (HALF)0x557fee8d, (HALF)0x51c856ba,
1387 (HALF)0xe6bc51ba, (HALF)0x587bc173, (HALF)0x1959e379, (HALF)0x92828439,
1388 (HALF)0x311e0503, (HALF)0x7f3c9cc2, (HALF)0xd426512c, (HALF)0xe8b2b497,
1389 (HALF)0x9e43536a, (HALF)0x9f544cb8, (HALF)0xb56f84c3, (HALF)0xb82fbb12,
1390 (HALF)0x6e348549, (HALF)0x45
1391 };
1392 STATIC CONST HALF h_rvec16[] = {
1393 (HALF)0x7b4e8830, (HALF)0x70605db8, (HALF)0xa1abe4a5, (HALF)0xa70fbe04,
1394 (HALF)0x2bcda8f4, (HALF)0xd29ada9a, (HALF)0x55ad0560, (HALF)0xb367137e,
1395 (HALF)0xd6972f1a, (HALF)0x809bad45, (HALF)0xb15d2454, (HALF)0x0c0d415f,
1396 (HALF)0x416e117a, (HALF)0xe87a9521, (HALF)0x670a5e1a, (HALF)0x53a41772,
1397 (HALF)0xfc9c5cc1, (HALF)0xf75645df, (HALF)0x86f6d19f, (HALF)0x9404b5bf,
1398 (HALF)0x56e9d83b, (HALF)0xac0f3bc3, (HALF)0xa1508c4b, (HALF)0x4bfd4977,
1399 (HALF)0x71922540, (HALF)0xf2524def, (HALF)0x8a81a3db, (HALF)0xced828de,
1400 (HALF)0x7895ae8f, (HALF)0x4a4b6dcd, (HALF)0x973e921a, (HALF)0x9fb27a07,
1401 (HALF)0xb0d7dcb1
1402 };
1403 STATIC CONST HALF h_nvec17[] = {
1404 (HALF)0x72b72051, (HALF)0xedc24ebf, (HALF)0xe970a8d1, (HALF)0x66c9b150,
1405 (HALF)0xcbb927f7, (HALF)0xb574ffd9, (HALF)0x4166b249, (HALF)0x0fce4030,
1406 (HALF)0xfa6922ca, (HALF)0x39cc14a9, (HALF)0x14396e2a, (HALF)0xaff74c7f,
1407 (HALF)0xa120a314, (HALF)0xe11a2700, (HALF)0x44c9ad30, (HALF)0x7f328d72,
1408 (HALF)0xab2ceaf7, (HALF)0x868ff772, (HALF)0x0974f0b3, (HALF)0x20e9f0f6,
1409 (HALF)0xcd4e5b8a, (HALF)0x0bde26e6, (HALF)0x96f6c3ac, (HALF)0x5718c601,
1410 (HALF)0x2f117710, (HALF)0x1ab0876e, (HALF)0x49ab2c2e, (HALF)0x747022b9,
1411 (HALF)0x6e9de4a7, (HALF)0x25e88f2e, (HALF)0xd7d0a00b, (HALF)0xc8ff11f6,
1412 (HALF)0xf50aa819, (HALF)0xbe530e9e, (HALF)0x47b7ff54, (HALF)0xe0209b46,
1413 (HALF)0x027ec5eb, (HALF)0x07543362, (HALF)0x531e9b85, (HALF)0x3c23b568,
1414 (HALF)0x5d07af7a, (HALF)0xd4948461, (HALF)0x2eb9a499, (HALF)0x1c71fa0b,
1415 (HALF)0x7025e22b, (HALF)0x4d2720b9, (HALF)0x531b8d54, (HALF)0x66e2fc30,
1416 (HALF)0xdac7cafd, (HALF)0x7f29953b, (HALF)0x9d456bff, (HALF)0x839814bc,
1417 (HALF)0xd0e5fb19, (HALF)0x5a6c9f58, (HALF)0x3a5a9dc3, (HALF)0x598bf28b,
1418 (HALF)0xa7a91144, (HALF)0x68494f76, (HALF)0xbfedd7be, (HALF)0x7ca54266,
1419 (HALF)0x96e0faf9, (HALF)0x33be3c0f, (HALF)0xffa3040b, (HALF)0x813aeac0,
1420 (HALF)0x6177
1421 };
1422 STATIC CONST HALF h_rvec17[] = {
1423 (HALF)0x22b41dac, (HALF)0xd6258005, (HALF)0x2aa1e0cb, (HALF)0x45d147b5,
1424 (HALF)0xbf5c46d9, (HALF)0x14c9dadf, (HALF)0x09b0aec4, (HALF)0x4286bfef,
1425 (HALF)0xc6f8e9d1, (HALF)0xdd68467b, (HALF)0x93f4ffb9, (HALF)0x58f2eb51,
1426 (HALF)0x2ade048f, (HALF)0xeacae6e5, (HALF)0x8dd2a807, (HALF)0xbcea8c27,
1427 (HALF)0x02a03281, (HALF)0x039aeb6d, (HALF)0xfa6e016c, (HALF)0x6fda1b09,
1428 (HALF)0xea7719ed, (HALF)0xcf2e0294, (HALF)0xa4264cb9, (HALF)0x49888af1,
1429 (HALF)0x1b44f0c1, (HALF)0x3ccee577, (HALF)0xdbebd170, (HALF)0xf7431e7e,
1430 (HALF)0xfaebd584, (HALF)0x896e4a33, (HALF)0xe46bf43c, (HALF)0x17fe9a10,
1431 (HALF)0xb5321b51, (HALF)0xf7e1d2a9, (HALF)0xe7fe0a56, (HALF)0xbd736750,
1432 (HALF)0xcf029a33, (HALF)0xebee99b1, (HALF)0x810fff31, (HALF)0x694c8d30,
1433 (HALF)0xc8c18689, (HALF)0xc2f9f4fb, (HALF)0x5949fd7f, (HALF)0x67aaf7b3,
1434 (HALF)0xa82f906a, (HALF)0x1b84b7b3, (HALF)0xeac052ee, (HALF)0x1a4e9345,
1435 (HALF)0xce3c2973, (HALF)0x4a5168a7, (HALF)0x5c5551ba, (HALF)0x77c6cb26,
1436 (HALF)0xfa45a3a6, (HALF)0x486f31e0, (HALF)0xcaf97519, (HALF)0xbe4b0399,
1437 (HALF)0x802fc106, (HALF)0x537284da, (HALF)0x20c4e167, (HALF)0x2a62f329,
1438 (HALF)0xc2d2fc5b, (HALF)0xdd665324, (HALF)0xc3b8adf1, (HALF)0x0b6eaf3b,
1439 (HALF)0x5372
1440 };
1441 STATIC CONST HALF h_nvec18[] = {
1442 (HALF)0xc8b78629, (HALF)0x41351b18, (HALF)0x28ad4ed8, (HALF)0xc96f7df1,
1443 (HALF)0x7cd3c931, (HALF)0x0f23036a, (HALF)0xac657631, (HALF)0x6a625812,
1444 (HALF)0x08144788, (HALF)0x8642ed62, (HALF)0x76198a40, (HALF)0x70defd64,
1445 (HALF)0x97fb673c, (HALF)0x6f3bddf6, (HALF)0x72fe2977, (HALF)0x20ed82f1,
1446 (HALF)0x7e7f4fdc, (HALF)0xdc272d6b, (HALF)0x6d77f317, (HALF)0x2c595b15,
1447 (HALF)0x1c3b7dd6, (HALF)0x6e3d6147, (HALF)0x8170640c, (HALF)0xb660033f,
1448 (HALF)0xf211886b, (HALF)0xbcc67859, (HALF)0x18ff4b93, (HALF)0x3068e691,
1449 (HALF)0x9db5d823, (HALF)0x72afd4ef, (HALF)0x5aa4cd1a, (HALF)0xa0a33014,
1450 (HALF)0x6b8349e6, (HALF)0x2f4de595, (HALF)0xd5180384, (HALF)0x19d94118,
1451 (HALF)0x369e7534, (HALF)0xce4d3b18, (HALF)0x119ef9ee, (HALF)0xe2f45c25,
1452 (HALF)0xe8eae0a7, (HALF)0x62e41605, (HALF)0x6346d2ca, (HALF)0x6425625b,
1453 (HALF)0x44b033de, (HALF)0x1711e6b3, (HALF)0xf3a9d02e, (HALF)0x259cd965,
1454 (HALF)0x08aa956b, (HALF)0x6ad64380, (HALF)0xe9730e8e, (HALF)0x539b9d28,
1455 (HALF)0xe5407950, (HALF)0x89900be4, (HALF)0xde1218f6, (HALF)0x63e1e52b,
1456 (HALF)0x0de03f4e, (HALF)0x8e21a568, (HALF)0x268d7ee3, (HALF)0xafb3514e,
1457 (HALF)0x5378efcb, (HALF)0xfec0c7c6, (HALF)0xf07cb724, (HALF)0xfb61b42a,
1458 (HALF)0x068f2a38
1459 };
1460 STATIC CONST HALF h_rvec18[] = {
1461 (HALF)0x35ea3c63, (HALF)0x8df2ef97, (HALF)0xa2b3afb7, (HALF)0x179158f6,
1462 (HALF)0x04920dba, (HALF)0xf333077e, (HALF)0xf8304b5a, (HALF)0x230ff2ae,
1463 (HALF)0x84a8f3f0, (HALF)0xadda164e, (HALF)0xc9a1c944, (HALF)0xc70502f2,
1464 (HALF)0x41a3c18f, (HALF)0x09bd3254, (HALF)0x973665a9, (HALF)0x1548c263,
1465 (HALF)0x5024d916, (HALF)0x0a3ddde9, (HALF)0xf2aaf1f5, (HALF)0x666db92a,
1466 (HALF)0x3a535aa5, (HALF)0x49c35775, (HALF)0xc381a1c4, (HALF)0xf8d36dbc,
1467 (HALF)0xe94be870, (HALF)0x430e88a6, (HALF)0x37219a06, (HALF)0x5109df80,
1468 (HALF)0xe73ca03f, (HALF)0xf1bb4541, (HALF)0x3c6f32f3, (HALF)0x952cfc24,
1469 (HALF)0xfbc9697f, (HALF)0xc5b8d472, (HALF)0xcbbeda4b, (HALF)0x7303db3b,
1470 (HALF)0xa18f255e, (HALF)0xe1d3353d, (HALF)0xe8c98700, (HALF)0x9e75e8fd,
1471 (HALF)0x3ddd812e, (HALF)0x7340f891, (HALF)0xb1f369ac, (HALF)0x764d505b,
1472 (HALF)0xb13ef51b, (HALF)0x16c3aa43, (HALF)0x61d042c4, (HALF)0x22ac0339,
1473 (HALF)0x3f306fd1, (HALF)0x49926f7f, (HALF)0x2b4c575c, (HALF)0x9f3ab467,
1474 (HALF)0x5dac65af, (HALF)0x62778dc7, (HALF)0x99113a89, (HALF)0x49c0540a,
1475 (HALF)0x1df70ac2, (HALF)0x4be12c5e, (HALF)0xe5d36bdb, (HALF)0x66b99ff5,
1476 (HALF)0x5358be89, (HALF)0x4cd835d7, (HALF)0xf0d5cda8, (HALF)0x1f1ac6c3,
1477 (HALF)0x04735e92
1478 };
1479 STATIC CONST HALF h_nvec19[] = {
1480 (HALF)0x6b659a79, (HALF)0x0239c12d, (HALF)0xd204df49, (HALF)0x1d4ae0c7,
1481 (HALF)0x099bf000, (HALF)0x6435ade8, (HALF)0xdc4af029, (HALF)0x2f4ee7a2,
1482 (HALF)0xadfcf1e3, (HALF)0x73358f43, (HALF)0x687eede5, (HALF)0xb567cd4d,
1483 (HALF)0xc7a7814f, (HALF)0xc306624c, (HALF)0xa82d80c6, (HALF)0x3f390cd5,
1484 (HALF)0x7b7dec3f, (HALF)0x8bdb1416, (HALF)0x275b3a52, (HALF)0x921884fe,
1485 (HALF)0x94ac02e0, (HALF)0x62f2b52c, (HALF)0xcdc992ee, (HALF)0x35e55eeb,
1486 (HALF)0x69a43fd5, (HALF)0x44c1dcfb, (HALF)0x3cdf6227, (HALF)0x23f3148f,
1487 (HALF)0x42508e4c, (HALF)0x95b737c3, (HALF)0x70af831f, (HALF)0x2ee815c9,
1488 (HALF)0x47a47251, (HALF)0x7c11b2af, (HALF)0x664c361f, (HALF)0xcadaa841,
1489 (HALF)0x3d97172a, (HALF)0x87402ffb, (HALF)0xa0a02ebd, (HALF)0xb674225f,
1490 (HALF)0x65593fe2, (HALF)0x85f698b9, (HALF)0x5ed9a7ab, (HALF)0x6bf37371,
1491 (HALF)0x7da75305, (HALF)0x088255bf, (HALF)0x8f607684, (HALF)0x1f3de57e,
1492 (HALF)0x118bd501, (HALF)0x6833770d, (HALF)0xd0425c51, (HALF)0x2664cacb,
1493 (HALF)0x9206b920, (HALF)0xeb903b8d, (HALF)0x0e2516b4, (HALF)0x36c3d841,
1494 (HALF)0x51d7cd17, (HALF)0xe063ba1e, (HALF)0xf28f4a6d, (HALF)0x244deb0d,
1495 (HALF)0x2e410ad4, (HALF)0x0721c315, (HALF)0xdde27654, (HALF)0x2ad6534b,
1496 (HALF)0xd6788b25, (HALF)0xb23bb9e8, (HALF)0x00230d7a
1497 };
1498 STATIC CONST HALF h_rvec19[] = {
1499 (HALF)0x698ef473, (HALF)0x3d53a5b7, (HALF)0x06448319, (HALF)0xd9ad4445,
1500 (HALF)0x6967daa0, (HALF)0xa14c6240, (HALF)0x78e77724, (HALF)0x63ef2ab7,
1501 (HALF)0x8dff2ee2, (HALF)0x662eb424, (HALF)0xcd9307d6, (HALF)0x0ab06a5d,
1502 (HALF)0xbdd7b539, (HALF)0x8621dd2d, (HALF)0x2bb4a187, (HALF)0x1f1e121d,
1503 (HALF)0xce8db962, (HALF)0xdd3eeaf1, (HALF)0x573b6ca1, (HALF)0x6f460cd6,
1504 (HALF)0x6a8d4780, (HALF)0xea68c7e6, (HALF)0x1148eb32, (HALF)0xb43d44d4,
1505 (HALF)0xb657cb64, (HALF)0x8547fdba, (HALF)0x85f73333, (HALF)0xc1f2a51a,
1506 (HALF)0x1c05ee52, (HALF)0x2847c03d, (HALF)0xfc7a88d8, (HALF)0x2e8dd186,
1507 (HALF)0x5be34683, (HALF)0xd43d2ee2, (HALF)0x9f2d2bb5, (HALF)0x89b3ddea,
1508 (HALF)0x7d784ae0, (HALF)0xc7735e28, (HALF)0x967a8608, (HALF)0xcdcbcf07,
1509 (HALF)0xfc17a423, (HALF)0xd36ad053, (HALF)0xc73d8892, (HALF)0xa635c3f4,
1510 (HALF)0x9b5d0cf9, (HALF)0x0ac73fd9, (HALF)0xe801fefb, (HALF)0xb31cffd2,
1511 (HALF)0xf3eaaa55, (HALF)0x0e74fa23, (HALF)0x5414b290, (HALF)0x6d176101,
1512 (HALF)0x522993f8, (HALF)0xf8293dad, (HALF)0x713b1e82, (HALF)0xb83bbef0,
1513 (HALF)0xd001cc62, (HALF)0x7537da39, (HALF)0x7d80158b, (HALF)0x9332c8e2,
1514 (HALF)0x6fa6fa75, (HALF)0xe21f6512, (HALF)0x999518b4, (HALF)0x2196605b,
1515 (HALF)0xe4fc1798, (HALF)0x5f21e245, (HALF)0x0008f172
1516 };
1517 STATIC CONST HALF h_nvec20[] = {
1518 (HALF)0xc3c1d081, (HALF)0x4d262fce, (HALF)0x8765cc91, (HALF)0xf3727f7c,
1519 (HALF)0xabba4bbc, (HALF)0xe0985801, (HALF)0xfa365c51, (HALF)0xb2a4b230,
1520 (HALF)0xf4430a8d, (HALF)0x546b98c8, (HALF)0xd9748b26, (HALF)0xe255a82f,
1521 (HALF)0xb00a3e8c, (HALF)0x7069676d, (HALF)0x6233ccce, (HALF)0x0299a74e,
1522 (HALF)0xd119adf9, (HALF)0x5273e811, (HALF)0xae36e0bb, (HALF)0x32b3a486,
1523 (HALF)0xf04921a5, (HALF)0xd31d28a2, (HALF)0xda0c50de, (HALF)0x21302b40,
1524 (HALF)0xee9de552, (HALF)0x80bbac03, (HALF)0x75f49740, (HALF)0xc3d0a61b,
1525 (HALF)0x341cfbf6, (HALF)0x1615b5b4, (HALF)0x4e3a3def, (HALF)0x95734dbf,
1526 (HALF)0xe7ab78ac, (HALF)0xffdc5cf9, (HALF)0xf7996892, (HALF)0x47407ba4,
1527 (HALF)0xa78b988f, (HALF)0x27736f42, (HALF)0x5064ee50, (HALF)0x60135c56,
1528 (HALF)0x1ad73283, (HALF)0x25624bf7, (HALF)0x2ee21419, (HALF)0x93195abd,
1529 (HALF)0x66b67778, (HALF)0xb1e0a42d, (HALF)0x729fb0f0, (HALF)0xd4921864,
1530 (HALF)0x2c42253f, (HALF)0x302a07a9, (HALF)0xbb741bd4, (HALF)0x932f90ba,
1531 (HALF)0xf3354be1, (HALF)0x0804d661, (HALF)0x010e9ba1, (HALF)0x1a05778d,
1532 (HALF)0xa962c833, (HALF)0xe7590ee8, (HALF)0xbe6803b8, (HALF)0xc67704c1,
1533 (HALF)0x56d79660, (HALF)0x6066a3f3, (HALF)0x648b0327, (HALF)0x267e5b3a,
1534 (HALF)0xdddc63a0, (HALF)0x3322e890, (HALF)0x20e0d8b1, (HALF)0x004fd2b8
1535 };
1536 STATIC CONST HALF h_rvec20[] = {
1537 (HALF)0xa048bd1a, (HALF)0x95abdc7b, (HALF)0x98f47cf8, (HALF)0x126ac98d,
1538 (HALF)0xaebf85fd, (HALF)0x5650580f, (HALF)0x3292d7dd, (HALF)0xf49e8377,
1539 (HALF)0x2947ed46, (HALF)0xd1a5b26c, (HALF)0xae14e6a1, (HALF)0x9b1f5788,
1540 (HALF)0x4df727b2, (HALF)0xee375079, (HALF)0x131bc8e4, (HALF)0x294e5f53,
1541 (HALF)0x1f5759bd, (HALF)0x65d58acf, (HALF)0x598ed3a5, (HALF)0xc39361a6,
1542 (HALF)0xa783fd7a, (HALF)0x264a36a2, (HALF)0x6ca22856, (HALF)0x8ffb171f,
1543 (HALF)0x1d7cea9e, (HALF)0xd81d6fca, (HALF)0x34ea730e, (HALF)0x31f56382,
1544 (HALF)0xb39cd9e9, (HALF)0x440e84be, (HALF)0x4b1d15a1, (HALF)0x7bf775c5,
1545 (HALF)0xe40f4638, (HALF)0xe5bef0a7, (HALF)0x79e58942, (HALF)0x881ae1ba,
1546 (HALF)0x01de8372, (HALF)0x14cf35f8, (HALF)0xe2d8b310, (HALF)0x66961207,
1547 (HALF)0xde5d5f91, (HALF)0xe6e70849, (HALF)0x74ec5ac3, (HALF)0xe2de4eb1,
1548 (HALF)0x4a41dc20, (HALF)0xd306d565, (HALF)0xb5843ff3, (HALF)0x911b30d6,
1549 (HALF)0x4e9cd926, (HALF)0x8455c9ae, (HALF)0x69448bb5, (HALF)0x0c7b1aad,
1550 (HALF)0x1da1e992, (HALF)0xc67656bd, (HALF)0xc544209e, (HALF)0x10ce387c,
1551 (HALF)0xc4e88df8, (HALF)0x40e8da88, (HALF)0xbb2c3028, (HALF)0x49194fd9,
1552 (HALF)0xdeef17ee, (HALF)0x241bc08d, (HALF)0x6fa9f608, (HALF)0x4b0f8b04,
1553 (HALF)0xee960da1, (HALF)0xa3099293, (HALF)0x84445fea, (HALF)0x0046ef01
1554 };
1555 #elif 2*FULL_BITS == 64
1556 STATIC CONST HALF h_nvec01[] = {
1557 (HALF)0x9361, (HALF)0x83de, (HALF)0x722d, (HALF)0xf0db,
1558 (HALF)0x28ca, (HALF)0x6fe3, (HALF)0x4073, (HALF)0x0494,
1559 (HALF)0x5
1560 };
1561 STATIC CONST HALF h_rvec01[] = {
1562 (HALF)0x42ec, (HALF)0xa4cc, (HALF)0xbb01, (HALF)0x4e5d,
1563 (HALF)0x52e7, (HALF)0x11d9, (HALF)0x980f, (HALF)0xb226
1564 };
1565 STATIC CONST HALF h_nvec02[] = {
1566 (HALF)0x43f1, (HALF)0x3534, (HALF)0x6ea9, (HALF)0xeb28,
1567 (HALF)0x4a18, (HALF)0xdd37, (HALF)0x2555, (HALF)0x348a,
1568 (HALF)0x2c5
1569 };
1570 STATIC CONST HALF h_rvec02[] = {
1571 (HALF)0xa218, (HALF)0x21e3, (HALF)0x616b, (HALF)0xe893,
1572 (HALF)0x10e3, (HALF)0x6cd7, (HALF)0x4344, (HALF)0xf3d6,
1573 (HALF)0x40
1574 };
1575 STATIC CONST HALF h_nvec03[] = {
1576 (HALF)0x01f1, (HALF)0x11d0, (HALF)0x661f, (HALF)0xf2ca,
1577 (HALF)0xf1e0, (HALF)0x3a81, (HALF)0xce4e, (HALF)0x59d6,
1578 (HALF)0xcfd9, (HALF)0x0009
1579 };
1580 STATIC CONST HALF h_rvec03[] = {
1581 (HALF)0xd76a, (HALF)0xa0d7, (HALF)0x2de2, (HALF)0x3e14,
1582 (HALF)0xea4f, (HALF)0xff5c, (HALF)0x9b64, (HALF)0xb44d,
1583 (HALF)0xfae5
1584 };
1585 STATIC CONST HALF h_nvec04[] = {
1586 (HALF)0x0751, (HALF)0xdfcc, (HALF)0xc680, (HALF)0x2dec,
1587 (HALF)0x2a1a, (HALF)0x5df1, (HALF)0x4ed7, (HALF)0x5c89,
1588 (HALF)0xf924, (HALF)0x3070
1589 };
1590 STATIC CONST HALF h_rvec04[] = {
1591 (HALF)0x4570, (HALF)0x4b98, (HALF)0xddba, (HALF)0xa220,
1592 (HALF)0xaf8a, (HALF)0xa2c0, (HALF)0x2bdc, (HALF)0x131b,
1593 (HALF)0xc2d8, (HALF)0x0020
1594 };
1595 STATIC CONST HALF h_nvec05[] = {
1596 (HALF)0x6ef1, (HALF)0x9916, (HALF)0xe5e7, (HALF)0x8b99,
1597 (HALF)0xa010, (HALF)0x8769, (HALF)0xe111, (HALF)0x5d3f,
1598 (HALF)0xc2fa, (HALF)0x680b, (HALF)0x5aac, (HALF)0x38f7,
1599 (HALF)0xa85b, (HALF)0xdb81, (HALF)0x1822, (HALF)0x109b,
1600 (HALF)0x2
1601 };
1602 STATIC CONST HALF h_rvec05[] = {
1603 (HALF)0xefa9, (HALF)0x59e2, (HALF)0x77c8, (HALF)0x0e6c,
1604 (HALF)0xaeed, (HALF)0x1e70, (HALF)0x7b7d, (HALF)0x234f,
1605 (HALF)0xf6db, (HALF)0x5f5d, (HALF)0xa960, (HALF)0xe821,
1606 (HALF)0xb792, (HALF)0xae33, (HALF)0x890e, (HALF)0x5e9b
1607 };
1608 STATIC CONST HALF h_nvec06[] = {
1609 (HALF)0xf431, (HALF)0xe1dd, (HALF)0x57f1, (HALF)0xd855,
1610 (HALF)0x32da, (HALF)0x5ee7, (HALF)0xdb77, (HALF)0x3a38,
1611 (HALF)0x4026, (HALF)0x5c64, (HALF)0xf218, (HALF)0xf2db,
1612 (HALF)0x2c79, (HALF)0x9ada, (HALF)0x9d7d, (HALF)0x7bfd,
1613 (HALF)0xa
1614 };
1615 STATIC CONST HALF h_rvec06[] = {
1616 (HALF)0x4daf, (HALF)0xc940, (HALF)0x2e80, (HALF)0xc5dc,
1617 (HALF)0xeccf, (HALF)0x2c98, (HALF)0x495d, (HALF)0xe1f3,
1618 (HALF)0x925c, (HALF)0xce1c, (HALF)0xaede, (HALF)0xe097,
1619 (HALF)0x7154, (HALF)0x8866, (HALF)0xa02f, (HALF)0x5e94
1620 };
1621 STATIC CONST HALF h_nvec07[] = {
1622 (HALF)0xc751, (HALF)0xcf9e, (HALF)0x9125, (HALF)0x602f,
1623 (HALF)0x2e7f, (HALF)0x5288, (HALF)0x53ce, (HALF)0x0dcf,
1624 (HALF)0x9d6b, (HALF)0xff56, (HALF)0x43fc, (HALF)0x6286,
1625 (HALF)0x1cd5, (HALF)0x3780, (HALF)0x9ef2, (HALF)0xf239,
1626 (HALF)0x7de8, (HALF)0x43d8
1627 };
1628 STATIC CONST HALF h_rvec07[] = {
1629 (HALF)0x25e6, (HALF)0x098d, (HALF)0xd2e5, (HALF)0x3992,
1630 (HALF)0xb58c, (HALF)0x64f0, (HALF)0xd4dd, (HALF)0xcf18,
1631 (HALF)0x6aef, (HALF)0x9d87, (HALF)0xed04, (HALF)0x7acc,
1632 (HALF)0x9076, (HALF)0xbfbe, (HALF)0x14c7, (HALF)0x1ee0,
1633 (HALF)0x522d, (HALF)0x0013
1634 };
1635 STATIC CONST HALF h_nvec08[] = {
1636 (HALF)0x2f11, (HALF)0x2674, (HALF)0xe66a, (HALF)0xbc42,
1637 (HALF)0xd9f0, (HALF)0xb59c, (HALF)0xa6c2, (HALF)0x9ad4,
1638 (HALF)0xd2f9, (HALF)0x5bdb, (HALF)0x1fed, (HALF)0xbdc9,
1639 (HALF)0x9ce7, (HALF)0xf13c, (HALF)0x99b7, (HALF)0xeb46,
1640 (HALF)0x6ca7, (HALF)0x4712, (HALF)0x58
1641 };
1642 STATIC CONST HALF h_rvec08[] = {
1643 (HALF)0xc674, (HALF)0x489d, (HALF)0x5f3a, (HALF)0xaae9,
1644 (HALF)0xa929, (HALF)0xa35d, (HALF)0xb4b8, (HALF)0x5597,
1645 (HALF)0xc947, (HALF)0x28e9, (HALF)0x4f9a, (HALF)0x3d34,
1646 (HALF)0x61fa, (HALF)0xb7e6, (HALF)0x9116, (HALF)0xa326,
1647 (HALF)0x16dc, (HALF)0x8530
1648 };
1649 STATIC CONST HALF h_nvec09[] = {
1650 (HALF)0xe3d1, (HALF)0xab27, (HALF)0x5db4, (HALF)0x1274,
1651 (HALF)0xf951, (HALF)0xb980, (HALF)0x6b66, (HALF)0x62b1,
1652 (HALF)0xce0d, (HALF)0x0fde, (HALF)0xc6fd, (HALF)0x6061,
1653 (HALF)0xff09, (HALF)0x36a6, (HALF)0xb61c, (HALF)0xe08e,
1654 (HALF)0x95c3, (HALF)0x84d8, (HALF)0x752a, (HALF)0x4a86,
1655 (HALF)0x7b4f, (HALF)0xc179, (HALF)0x57a3, (HALF)0x5621,
1656 (HALF)0x7bb0, (HALF)0x3d26, (HALF)0x1b00, (HALF)0x14e8,
1657 (HALF)0x9238, (HALF)0x218d, (HALF)0x2fd3, (HALF)0x5232,
1658 (HALF)0xe8be, (HALF)0x0039
1659 };
1660 STATIC CONST HALF h_rvec09[] = {
1661 (HALF)0xd20d, (HALF)0x7d4e, (HALF)0xf2b8, (HALF)0x601e,
1662 (HALF)0xf959, (HALF)0x8e59, (HALF)0x5d9e, (HALF)0xedaa,
1663 (HALF)0x89ba, (HALF)0x309a, (HALF)0x7d81, (HALF)0xe5ab,
1664 (HALF)0x2545, (HALF)0x796b, (HALF)0x3222, (HALF)0x02de,
1665 (HALF)0xc0bd, (HALF)0x8357, (HALF)0xe3fd, (HALF)0x0107,
1666 (HALF)0xd288, (HALF)0x82d9, (HALF)0xa8aa, (HALF)0xbc42,
1667 (HALF)0x7343, (HALF)0x4b78, (HALF)0x0886, (HALF)0xc015,
1668 (HALF)0x15bf, (HALF)0xbab9, (HALF)0x1e6b, (HALF)0xf8ad,
1669 (HALF)0xb458
1670 };
1671 STATIC CONST HALF h_nvec10[] = {
1672 (HALF)0x4b89, (HALF)0xb7e6, (HALF)0xc363, (HALF)0xc3cd,
1673 (HALF)0xc73c, (HALF)0x2ef9, (HALF)0xce22, (HALF)0x6092,
1674 (HALF)0xe36c, (HALF)0x02ab, (HALF)0x9573, (HALF)0x08d4,
1675 (HALF)0x1c40, (HALF)0x7451, (HALF)0x82de, (HALF)0xd385,
1676 (HALF)0xa02f, (HALF)0xa524, (HALF)0x1b3b, (HALF)0x52c8,
1677 (HALF)0x3cc9, (HALF)0x250d, (HALF)0x0e88, (HALF)0x23b5,
1678 (HALF)0x48bf, (HALF)0xbd14, (HALF)0x7f98, (HALF)0x882d,
1679 (HALF)0xf596, (HALF)0xc23e, (HALF)0x5666, (HALF)0xc905,
1680 (HALF)0x2435, (HALF)0x025f
1681 };
1682 STATIC CONST HALF h_rvec10[] = {
1683 (HALF)0xc482, (HALF)0x94cf, (HALF)0x5ad4, (HALF)0x594f,
1684 (HALF)0x2aee, (HALF)0x2344, (HALF)0x40ce, (HALF)0x145f,
1685 (HALF)0x950d, (HALF)0x1baf, (HALF)0xf175, (HALF)0xadc4,
1686 (HALF)0x669f, (HALF)0xf62c, (HALF)0x5d56, (HALF)0x8d07,
1687 (HALF)0x8b40, (HALF)0x08ed, (HALF)0xcf30, (HALF)0xaaf2,
1688 (HALF)0x5ffb, (HALF)0xc24b, (HALF)0xf8cf, (HALF)0x250d,
1689 (HALF)0x1ec9, (HALF)0x7ca8, (HALF)0x2b70, (HALF)0x787e,
1690 (HALF)0x1468, (HALF)0x1840, (HALF)0x0e0c, (HALF)0x47b2,
1691 (HALF)0xbb7e, (HALF)0x0066
1692 };
1693 STATIC CONST HALF h_nvec11[] = {
1694 (HALF)0xe069, (HALF)0x546e, (HALF)0x530c, (HALF)0x2e1a,
1695 (HALF)0xdab2, (HALF)0x2014, (HALF)0xcf52, (HALF)0xa729,
1696 (HALF)0xe1a9, (HALF)0x920e, (HALF)0x7533, (HALF)0x68f2,
1697 (HALF)0x3cfa, (HALF)0x2587, (HALF)0xa749, (HALF)0xdd37,
1698 (HALF)0xdaa2, (HALF)0x4499, (HALF)0x5870, (HALF)0x286e,
1699 (HALF)0xf9b6, (HALF)0x57f3, (HALF)0x4467, (HALF)0x5ec5,
1700 (HALF)0x91ea, (HALF)0x69a7, (HALF)0xcd77, (HALF)0x874e,
1701 (HALF)0xd56b, (HALF)0x4217, (HALF)0xb309, (HALF)0x82bd,
1702 (HALF)0x64de, (HALF)0x4978
1703 };
1704 STATIC CONST HALF h_rvec11[] = {
1705 (HALF)0x8b04, (HALF)0x56e3, (HALF)0xded3, (HALF)0x3a0a,
1706 (HALF)0x88b1, (HALF)0x461d, (HALF)0x4d65, (HALF)0x9c09,
1707 (HALF)0x3fed, (HALF)0xe533, (HALF)0x18fe, (HALF)0x34d9,
1708 (HALF)0x6281, (HALF)0x1ef5, (HALF)0xa07c, (HALF)0xcedf,
1709 (HALF)0x47fb, (HALF)0x590f, (HALF)0x4d5c, (HALF)0xa2c5,
1710 (HALF)0x39ee, (HALF)0x7323, (HALF)0x49a7, (HALF)0x8065,
1711 (HALF)0x163f, (HALF)0x9ce3, (HALF)0xf8b6, (HALF)0xae3a,
1712 (HALF)0x4465, (HALF)0x264a, (HALF)0xe630, (HALF)0x1cb5,
1713 (HALF)0x8488, (HALF)0x0086
1714 };
1715 STATIC CONST HALF h_nvec12[] = {
1716 (HALF)0x7b99, (HALF)0xf14c, (HALF)0xd151, (HALF)0x7f66,
1717 (HALF)0xad2b, (HALF)0x87ef, (HALF)0xf098, (HALF)0x57d3,
1718 (HALF)0x4165, (HALF)0xd653, (HALF)0xdd25, (HALF)0x812f,
1719 (HALF)0xc7ce, (HALF)0x48c7, (HALF)0x41e0, (HALF)0xa1bf,
1720 (HALF)0xe315, (HALF)0x4c94, (HALF)0x1593, (HALF)0x190b,
1721 (HALF)0x51da, (HALF)0xee42, (HALF)0x1a66, (HALF)0x2b4a,
1722 (HALF)0xc2a1, (HALF)0x2bb7, (HALF)0x8ca9, (HALF)0x65b1,
1723 (HALF)0x9116, (HALF)0x08b8, (HALF)0xb15f, (HALF)0xc0cc,
1724 (HALF)0x2ab3, (HALF)0x5758, (HALF)0x34
1725 };
1726 STATIC CONST HALF h_rvec12[] = {
1727 (HALF)0xb4a0, (HALF)0xe207, (HALF)0xdd68, (HALF)0x5227,
1728 (HALF)0xfbc4, (HALF)0x9488, (HALF)0x81aa, (HALF)0x6ed0,
1729 (HALF)0x6fe5, (HALF)0x8e73, (HALF)0xc020, (HALF)0x3dd2,
1730 (HALF)0x7c57, (HALF)0xeeb0, (HALF)0x4eb7, (HALF)0x0b60,
1731 (HALF)0x72a1, (HALF)0xa13f, (HALF)0x4333, (HALF)0xcbfc,
1732 (HALF)0xe8cd, (HALF)0xa4c5, (HALF)0x520d, (HALF)0x0add,
1733 (HALF)0xa3b2, (HALF)0x758c, (HALF)0x0137, (HALF)0x5549,
1734 (HALF)0xbabd, (HALF)0x5870, (HALF)0xed93, (HALF)0xf648,
1735 (HALF)0x9bd1, (HALF)0xdf71
1736 };
1737 STATIC CONST HALF h_nvec13[] = {
1738 (HALF)0xfc49, (HALF)0x3314, (HALF)0x0032, (HALF)0xcca2,
1739 (HALF)0x3420, (HALF)0x208e, (HALF)0x503a, (HALF)0x8aaa,
1740 (HALF)0x63cc, (HALF)0xd79a, (HALF)0x7417, (HALF)0xb4ed,
1741 (HALF)0x1892, (HALF)0x95dd, (HALF)0x5f64, (HALF)0xb591,
1742 (HALF)0xc7f1, (HALF)0xd14c, (HALF)0x917e, (HALF)0x1589,
1743 (HALF)0x5667, (HALF)0xb2b0, (HALF)0x99cb, (HALF)0xc32d,
1744 (HALF)0x1a84, (HALF)0x1b5a, (HALF)0x22a1, (HALF)0xa493,
1745 (HALF)0x76c3, (HALF)0xdd3a, (HALF)0x137f, (HALF)0x8c07,
1746 (HALF)0x3c63, (HALF)0xaaf8, (HALF)0x5113, (HALF)0x3757,
1747 (HALF)0x8e84, (HALF)0xa8b1, (HALF)0x891d, (HALF)0xceec,
1748 (HALF)0xee99, (HALF)0x78c1, (HALF)0xe256, (HALF)0x6e49,
1749 (HALF)0xbfd6, (HALF)0x4286, (HALF)0xf6a9, (HALF)0xcb6b,
1750 (HALF)0x8ee0, (HALF)0x7bda, (HALF)0x510a, (HALF)0xd439,
1751 (HALF)0x345b, (HALF)0x63f4, (HALF)0x5535, (HALF)0x959a,
1752 (HALF)0x6d82, (HALF)0xdaf6, (HALF)0xd833, (HALF)0xed03,
1753 (HALF)0xf734, (HALF)0x1b5a, (HALF)0x7dd2, (HALF)0x166b,
1754 (HALF)0x7c19, (HALF)0x0151
1755 };
1756 STATIC CONST HALF h_rvec13[] = {
1757 (HALF)0x36f5, (HALF)0x6b77, (HALF)0xbfe4, (HALF)0x2407,
1758 (HALF)0x2072, (HALF)0x965e, (HALF)0xcf3e, (HALF)0xcc26,
1759 (HALF)0xb567, (HALF)0xa432, (HALF)0x07ab, (HALF)0x2ed0,
1760 (HALF)0x67b9, (HALF)0x0e2f, (HALF)0x0960, (HALF)0xef64,
1761 (HALF)0x1ad3, (HALF)0xbe5f, (HALF)0xda1b, (HALF)0x3fae,
1762 (HALF)0xb988, (HALF)0xa8f6, (HALF)0xcea5, (HALF)0xe5c9,
1763 (HALF)0x45ea, (HALF)0xa674, (HALF)0xce78, (HALF)0x3935,
1764 (HALF)0xff06, (HALF)0xf445, (HALF)0x0a11, (HALF)0xbeda,
1765 (HALF)0x080b, (HALF)0x92de, (HALF)0x9026, (HALF)0xf404,
1766 (HALF)0xf0b8, (HALF)0xe509, (HALF)0x216f, (HALF)0x6f05,
1767 (HALF)0xdc14, (HALF)0x5a68, (HALF)0x730d, (HALF)0x548f,
1768 (HALF)0x9f00, (HALF)0xe9fd, (HALF)0xaada, (HALF)0x64a4,
1769 (HALF)0xdd15, (HALF)0xe3bb, (HALF)0xe7a5, (HALF)0xcb4b,
1770 (HALF)0xd162, (HALF)0x17dd, (HALF)0x8c33, (HALF)0xc491,
1771 (HALF)0x6d2b, (HALF)0x9b70, (HALF)0xf6a6, (HALF)0x8b04,
1772 (HALF)0xfa64, (HALF)0x1263, (HALF)0x560d, (HALF)0x8e9a,
1773 (HALF)0xd42e
1774 };
1775 STATIC CONST HALF h_nvec14[] = {
1776 (HALF)0xaf01, (HALF)0xc116, (HALF)0x8c0f, (HALF)0xbdef,
1777 (HALF)0x9a1a, (HALF)0xc440, (HALF)0x185c, (HALF)0xacb3,
1778 (HALF)0x925b, (HALF)0xb33f, (HALF)0x3005, (HALF)0xfee8,
1779 (HALF)0xb112, (HALF)0x4b3d, (HALF)0x6743, (HALF)0x7f07,
1780 (HALF)0x9223, (HALF)0x2170, (HALF)0x054b, (HALF)0x2159,
1781 (HALF)0xefe3, (HALF)0x6fdf, (HALF)0x0a07, (HALF)0x792d,
1782 (HALF)0xbd52, (HALF)0x1d14, (HALF)0x9a27, (HALF)0x5f28,
1783 (HALF)0x83c4, (HALF)0xe034, (HALF)0xa0c1, (HALF)0xcb86,
1784 (HALF)0x912a, (HALF)0x0ede, (HALF)0x33e3, (HALF)0xf01a,
1785 (HALF)0x40c6, (HALF)0x63fd, (HALF)0x86e4, (HALF)0x4b05,
1786 (HALF)0x4a03, (HALF)0xbc7e, (HALF)0x22f3, (HALF)0x956b,
1787 (HALF)0x0b22, (HALF)0x1056, (HALF)0xd321, (HALF)0x3e78,
1788 (HALF)0x61e2, (HALF)0x1791, (HALF)0xb909, (HALF)0xe85e,
1789 (HALF)0x31ee, (HALF)0xec79, (HALF)0xb4bf, (HALF)0xb314,
1790 (HALF)0x4618, (HALF)0x1f56, (HALF)0x83d0, (HALF)0xb9d9,
1791 (HALF)0xac07, (HALF)0x7479, (HALF)0xf4e8, (HALF)0x93c6,
1792 (HALF)0xa00e, (HALF)0x5e56
1793 };
1794 STATIC CONST HALF h_rvec14[] = {
1795 (HALF)0xf190, (HALF)0x0ff9, (HALF)0xdb68, (HALF)0x47a4,
1796 (HALF)0xc8ea, (HALF)0x913c, (HALF)0xb220, (HALF)0xb6b1,
1797 (HALF)0xfbbb, (HALF)0x13ed, (HALF)0xf1c3, (HALF)0xa8f1,
1798 (HALF)0x1f8f, (HALF)0xd6d7, (HALF)0x649a, (HALF)0x4194,
1799 (HALF)0x7344, (HALF)0x7d49, (HALF)0x8416, (HALF)0x677c,
1800 (HALF)0xb983, (HALF)0x0186, (HALF)0x3901, (HALF)0xee63,
1801 (HALF)0xd69d, (HALF)0xce64, (HALF)0x04b3, (HALF)0x61a7,
1802 (HALF)0x52b2, (HALF)0x1383, (HALF)0x58d8, (HALF)0xc0fb,
1803 (HALF)0x2073, (HALF)0x16bf, (HALF)0xae78, (HALF)0x56c9,
1804 (HALF)0x0a68, (HALF)0xb81a, (HALF)0xabaf, (HALF)0x1512,
1805 (HALF)0x36ba, (HALF)0x6b99, (HALF)0x0dfc, (HALF)0x4335,
1806 (HALF)0x19d2, (HALF)0xa0ea, (HALF)0x34c5, (HALF)0xe861,
1807 (HALF)0x563f, (HALF)0x6c86, (HALF)0x5b68, (HALF)0x6b0c,
1808 (HALF)0x27fc, (HALF)0x75f6, (HALF)0x913f, (HALF)0xb609,
1809 (HALF)0xa564, (HALF)0x15b9, (HALF)0xf154, (HALF)0x9b18,
1810 (HALF)0xc5d0, (HALF)0x6ef0, (HALF)0x3509, (HALF)0xb673,
1811 (HALF)0xaa7c, (HALF)0x00f7
1812 };
1813 STATIC CONST HALF h_nvec15[] = {
1814 (HALF)0x7079, (HALF)0xc8d9, (HALF)0x7597, (HALF)0x061e,
1815 (HALF)0xc721, (HALF)0xf5d2, (HALF)0xc51f, (HALF)0x299b,
1816 (HALF)0xc337, (HALF)0xffe6, (HALF)0x8624, (HALF)0x1979,
1817 (HALF)0x92b6, (HALF)0xee6f, (HALF)0x0c7a, (HALF)0x0b1d,
1818 (HALF)0x8231, (HALF)0xb530, (HALF)0x58dd, (HALF)0x49c5,
1819 (HALF)0x530e, (HALF)0x196a, (HALF)0x515c, (HALF)0x0caa,
1820 (HALF)0x86ed, (HALF)0x8b0d, (HALF)0x8fa0, (HALF)0x380a,
1821 (HALF)0x03e4, (HALF)0x80df, (HALF)0xd81b, (HALF)0xe962,
1822 (HALF)0x83b3, (HALF)0xc1a7, (HALF)0x8ccc, (HALF)0xc327,
1823 (HALF)0xab86, (HALF)0x3e72, (HALF)0x1675, (HALF)0xebb9,
1824 (HALF)0xc2b8, (HALF)0xb902, (HALF)0x0445, (HALF)0xa059,
1825 (HALF)0xde42, (HALF)0x4f40, (HALF)0xaac1, (HALF)0xaa95,
1826 (HALF)0x2e82, (HALF)0xfc0f, (HALF)0xb84b, (HALF)0x70ca,
1827 (HALF)0xa267, (HALF)0x5326, (HALF)0xe607, (HALF)0x470b,
1828 (HALF)0x2ebe, (HALF)0x4535, (HALF)0x1ca8, (HALF)0x5ba8,
1829 (HALF)0x6c17, (HALF)0x02c4, (HALF)0xbcdb, (HALF)0x9edf,
1830 (HALF)0x840b, (HALF)0x97dd
1831 };
1832 STATIC CONST HALF h_rvec15[] = {
1833 (HALF)0x2110, (HALF)0x6c3e, (HALF)0x0aaa, (HALF)0x808f,
1834 (HALF)0xb92e, (HALF)0xd98d, (HALF)0xbd43, (HALF)0x1e6a,
1835 (HALF)0xb920, (HALF)0xf401, (HALF)0x0381, (HALF)0x9d3f,
1836 (HALF)0xd174, (HALF)0xdb95, (HALF)0x5c33, (HALF)0xa2f6,
1837 (HALF)0x69f8, (HALF)0x0c54, (HALF)0x26fc, (HALF)0xa3c1,
1838 (HALF)0x241b, (HALF)0x8866, (HALF)0xeca7, (HALF)0x0e46,
1839 (HALF)0x57fa, (HALF)0xa705, (HALF)0x91a6, (HALF)0x8e73,
1840 (HALF)0xa9e2, (HALF)0x70e4, (HALF)0x7e89, (HALF)0x9ad9,
1841 (HALF)0x7dee, (HALF)0x1eb7, (HALF)0x24d0, (HALF)0xfe62,
1842 (HALF)0x22b1, (HALF)0xdbe5, (HALF)0x1023, (HALF)0xa4aa,
1843 (HALF)0x6860, (HALF)0x7f22, (HALF)0x4379, (HALF)0x60ef,
1844 (HALF)0xb296, (HALF)0x45e3, (HALF)0xf5ef, (HALF)0xcc31,
1845 (HALF)0xf6d6, (HALF)0x74fb, (HALF)0xc25d, (HALF)0x55b1,
1846 (HALF)0x521f, (HALF)0x3ede, (HALF)0xef42, (HALF)0xdf8c,
1847 (HALF)0x7ca6, (HALF)0xa8d7, (HALF)0x25da, (HALF)0xb500,
1848 (HALF)0x99f9, (HALF)0x69ab, (HALF)0xc758, (HALF)0x03b8,
1849 (HALF)0x2207, (HALF)0x00b8
1850 };
1851 STATIC CONST HALF h_nvec16[] = {
1852 (HALF)0x46e1, (HALF)0xd843, (HALF)0x83f6, (HALF)0x1841,
1853 (HALF)0xbd36, (HALF)0x2dc9, (HALF)0x57ac, (HALF)0x4ca8,
1854 (HALF)0x828d, (HALF)0x96a5, (HALF)0x1c59, (HALF)0xed1f,
1855 (HALF)0x731f, (HALF)0x36d9, (HALF)0x6183, (HALF)0xbd3f,
1856 (HALF)0x5578, (HALF)0xde0f, (HALF)0xea8a, (HALF)0xb6a2,
1857 (HALF)0x3c44, (HALF)0xbe99, (HALF)0x3c05, (HALF)0x0e28,
1858 (HALF)0x61e3, (HALF)0xd7cf, (HALF)0x15b6, (HALF)0x40fe,
1859 (HALF)0x967e, (HALF)0x534d, (HALF)0x1046, (HALF)0x3469,
1860 (HALF)0x45bd, (HALF)0xd408, (HALF)0xee4b, (HALF)0xd69c,
1861 (HALF)0xee8d, (HALF)0x557f, (HALF)0x56ba, (HALF)0x51c8,
1862 (HALF)0x51ba, (HALF)0xe6bc, (HALF)0xc173, (HALF)0x587b,
1863 (HALF)0xe379, (HALF)0x1959, (HALF)0x8439, (HALF)0x9282,
1864 (HALF)0x0503, (HALF)0x311e, (HALF)0x9cc2, (HALF)0x7f3c,
1865 (HALF)0x512c, (HALF)0xd426, (HALF)0xb497, (HALF)0xe8b2,
1866 (HALF)0x536a, (HALF)0x9e43, (HALF)0x4cb8, (HALF)0x9f54,
1867 (HALF)0x84c3, (HALF)0xb56f, (HALF)0xbb12, (HALF)0xb82f,
1868 (HALF)0x8549, (HALF)0x6e34, (HALF)0x45
1869 };
1870 STATIC CONST HALF h_rvec16[] = {
1871 (HALF)0x8830, (HALF)0x7b4e, (HALF)0x5db8, (HALF)0x7060, (HALF)0xe4a5,
1872 (HALF)0xa1ab, (HALF)0xbe04, (HALF)0xa70f,
1873 (HALF)0xa8f4, (HALF)0x2bcd, (HALF)0xda9a, (HALF)0xd29a,
1874 (HALF)0x0560, (HALF)0x55ad, (HALF)0x137e, (HALF)0xb367,
1875 (HALF)0x2f1a, (HALF)0xd697, (HALF)0xad45, (HALF)0x809b,
1876 (HALF)0x2454, (HALF)0xb15d, (HALF)0x415f, (HALF)0x0c0d,
1877 (HALF)0x117a, (HALF)0x416e, (HALF)0x9521, (HALF)0xe87a,
1878 (HALF)0x5e1a, (HALF)0x670a, (HALF)0x1772, (HALF)0x53a4,
1879 (HALF)0x5cc1, (HALF)0xfc9c, (HALF)0x45df, (HALF)0xf756,
1880 (HALF)0xd19f, (HALF)0x86f6, (HALF)0xb5bf, (HALF)0x9404,
1881 (HALF)0xd83b, (HALF)0x56e9, (HALF)0x3bc3, (HALF)0xac0f,
1882 (HALF)0x8c4b, (HALF)0xa150, (HALF)0x4977, (HALF)0x4bfd,
1883 (HALF)0x2540, (HALF)0x7192, (HALF)0x4def, (HALF)0xf252,
1884 (HALF)0xa3db, (HALF)0x8a81, (HALF)0x28de, (HALF)0xced8,
1885 (HALF)0xae8f, (HALF)0x7895, (HALF)0x6dcd, (HALF)0x4a4b,
1886 (HALF)0x921a, (HALF)0x973e, (HALF)0x7a07, (HALF)0x9fb2,
1887 (HALF)0xdcb1, (HALF)0xb0d7
1888 };
1889 STATIC CONST HALF h_nvec17[] = {
1890 (HALF)0x2051, (HALF)0x72b7, (HALF)0x4ebf, (HALF)0xedc2,
1891 (HALF)0xa8d1, (HALF)0xe970, (HALF)0xb150, (HALF)0x66c9,
1892 (HALF)0x27f7, (HALF)0xcbb9, (HALF)0xffd9, (HALF)0xb574,
1893 (HALF)0xb249, (HALF)0x4166, (HALF)0x4030, (HALF)0x0fce,
1894 (HALF)0x22ca, (HALF)0xfa69, (HALF)0x14a9, (HALF)0x39cc,
1895 (HALF)0x6e2a, (HALF)0x1439, (HALF)0x4c7f, (HALF)0xaff7,
1896 (HALF)0xa314, (HALF)0xa120, (HALF)0x2700, (HALF)0xe11a,
1897 (HALF)0xad30, (HALF)0x44c9, (HALF)0x8d72, (HALF)0x7f32,
1898 (HALF)0xeaf7, (HALF)0xab2c, (HALF)0xf772, (HALF)0x868f,
1899 (HALF)0xf0b3, (HALF)0x0974, (HALF)0xf0f6, (HALF)0x20e9,
1900 (HALF)0x5b8a, (HALF)0xcd4e, (HALF)0x26e6, (HALF)0x0bde,
1901 (HALF)0xc3ac, (HALF)0x96f6, (HALF)0xc601, (HALF)0x5718,
1902 (HALF)0x7710, (HALF)0x2f11, (HALF)0x876e, (HALF)0x1ab0,
1903 (HALF)0x2c2e, (HALF)0x49ab, (HALF)0x22b9, (HALF)0x7470,
1904 (HALF)0xe4a7, (HALF)0x6e9d, (HALF)0x8f2e, (HALF)0x25e8,
1905 (HALF)0xa00b, (HALF)0xd7d0, (HALF)0x11f6, (HALF)0xc8ff,
1906 (HALF)0xa819, (HALF)0xf50a, (HALF)0x0e9e, (HALF)0xbe53,
1907 (HALF)0xff54, (HALF)0x47b7, (HALF)0x9b46, (HALF)0xe020,
1908 (HALF)0xc5eb, (HALF)0x027e, (HALF)0x3362, (HALF)0x0754,
1909 (HALF)0x9b85, (HALF)0x531e, (HALF)0xb568, (HALF)0x3c23,
1910 (HALF)0xaf7a, (HALF)0x5d07, (HALF)0x8461, (HALF)0xd494,
1911 (HALF)0xa499, (HALF)0x2eb9, (HALF)0xfa0b, (HALF)0x1c71,
1912 (HALF)0xe22b, (HALF)0x7025, (HALF)0x20b9, (HALF)0x4d27,
1913 (HALF)0x8d54, (HALF)0x531b, (HALF)0xfc30, (HALF)0x66e2,
1914 (HALF)0xcafd, (HALF)0xdac7, (HALF)0x953b, (HALF)0x7f29,
1915 (HALF)0x6bff, (HALF)0x9d45, (HALF)0x14bc, (HALF)0x8398,
1916 (HALF)0xfb19, (HALF)0xd0e5, (HALF)0x9f58, (HALF)0x5a6c,
1917 (HALF)0x9dc3, (HALF)0x3a5a, (HALF)0xf28b, (HALF)0x598b,
1918 (HALF)0x1144, (HALF)0xa7a9, (HALF)0x4f76, (HALF)0x6849,
1919 (HALF)0xd7be, (HALF)0xbfed, (HALF)0x4266, (HALF)0x7ca5,
1920 (HALF)0xfaf9, (HALF)0x96e0, (HALF)0x3c0f, (HALF)0x33be,
1921 (HALF)0x040b, (HALF)0xffa3, (HALF)0xeac0, (HALF)0x813a,
1922 (HALF)0x6177
1923 };
1924 STATIC CONST HALF h_rvec17[] = {
1925 (HALF)0x1dac, (HALF)0x22b4, (HALF)0x8005, (HALF)0xd625,
1926 (HALF)0xe0cb, (HALF)0x2aa1, (HALF)0x47b5, (HALF)0x45d1,
1927 (HALF)0x46d9, (HALF)0xbf5c, (HALF)0xdadf, (HALF)0x14c9,
1928 (HALF)0xaec4, (HALF)0x09b0, (HALF)0xbfef, (HALF)0x4286,
1929 (HALF)0xe9d1, (HALF)0xc6f8, (HALF)0x467b, (HALF)0xdd68,
1930 (HALF)0xffb9, (HALF)0x93f4, (HALF)0xeb51, (HALF)0x58f2,
1931 (HALF)0x048f, (HALF)0x2ade, (HALF)0xe6e5, (HALF)0xeaca,
1932 (HALF)0xa807, (HALF)0x8dd2, (HALF)0x8c27, (HALF)0xbcea,
1933 (HALF)0x3281, (HALF)0x02a0, (HALF)0xeb6d, (HALF)0x039a,
1934 (HALF)0x016c, (HALF)0xfa6e, (HALF)0x1b09, (HALF)0x6fda,
1935 (HALF)0x19ed, (HALF)0xea77, (HALF)0x0294, (HALF)0xcf2e,
1936 (HALF)0x4cb9, (HALF)0xa426, (HALF)0x8af1, (HALF)0x4988,
1937 (HALF)0xf0c1, (HALF)0x1b44, (HALF)0xe577, (HALF)0x3cce,
1938 (HALF)0xd170, (HALF)0xdbeb, (HALF)0x1e7e, (HALF)0xf743,
1939 (HALF)0xd584, (HALF)0xfaeb, (HALF)0x4a33, (HALF)0x896e,
1940 (HALF)0xf43c, (HALF)0xe46b, (HALF)0x9a10, (HALF)0x17fe,
1941 (HALF)0x1b51, (HALF)0xb532, (HALF)0xd2a9, (HALF)0xf7e1,
1942 (HALF)0x0a56, (HALF)0xe7fe, (HALF)0x6750, (HALF)0xbd73,
1943 (HALF)0x9a33, (HALF)0xcf02, (HALF)0x99b1, (HALF)0xebee,
1944 (HALF)0xff31, (HALF)0x810f, (HALF)0x8d30, (HALF)0x694c,
1945 (HALF)0x8689, (HALF)0xc8c1, (HALF)0xf4fb, (HALF)0xc2f9,
1946 (HALF)0xfd7f, (HALF)0x5949, (HALF)0xf7b3, (HALF)0x67aa,
1947 (HALF)0x906a, (HALF)0xa82f, (HALF)0xb7b3, (HALF)0x1b84,
1948 (HALF)0x52ee, (HALF)0xeac0, (HALF)0x9345, (HALF)0x1a4e,
1949 (HALF)0x2973, (HALF)0xce3c, (HALF)0x68a7, (HALF)0x4a51,
1950 (HALF)0x51ba, (HALF)0x5c55, (HALF)0xcb26, (HALF)0x77c6,
1951 (HALF)0xa3a6, (HALF)0xfa45, (HALF)0x31e0, (HALF)0x486f,
1952 (HALF)0x7519, (HALF)0xcaf9, (HALF)0x0399, (HALF)0xbe4b,
1953 (HALF)0xc106, (HALF)0x802f, (HALF)0x84da, (HALF)0x5372,
1954 (HALF)0xe167, (HALF)0x20c4, (HALF)0xf329, (HALF)0x2a62,
1955 (HALF)0xfc5b, (HALF)0xc2d2, (HALF)0x5324, (HALF)0xdd66,
1956 (HALF)0xadf1, (HALF)0xc3b8, (HALF)0xaf3b, (HALF)0x0b6e,
1957 (HALF)0x5372
1958 };
1959 STATIC CONST HALF h_nvec18[] = {
1960 (HALF)0x8629, (HALF)0xc8b7, (HALF)0x1b18, (HALF)0x4135, (HALF)0x4ed8,
1961 (HALF)0x28ad, (HALF)0x7df1, (HALF)0xc96f,
1962 (HALF)0xc931, (HALF)0x7cd3, (HALF)0x036a, (HALF)0x0f23,
1963 (HALF)0x7631, (HALF)0xac65, (HALF)0x5812, (HALF)0x6a62,
1964 (HALF)0x4788, (HALF)0x0814, (HALF)0xed62, (HALF)0x8642,
1965 (HALF)0x8a40, (HALF)0x7619, (HALF)0xfd64, (HALF)0x70de,
1966 (HALF)0x673c, (HALF)0x97fb, (HALF)0xddf6, (HALF)0x6f3b,
1967 (HALF)0x2977, (HALF)0x72fe, (HALF)0x82f1, (HALF)0x20ed,
1968 (HALF)0x4fdc, (HALF)0x7e7f, (HALF)0x2d6b, (HALF)0xdc27,
1969 (HALF)0xf317, (HALF)0x6d77, (HALF)0x5b15, (HALF)0x2c59,
1970 (HALF)0x7dd6, (HALF)0x1c3b, (HALF)0x6147, (HALF)0x6e3d,
1971 (HALF)0x640c, (HALF)0x8170, (HALF)0x033f, (HALF)0xb660,
1972 (HALF)0x886b, (HALF)0xf211, (HALF)0x7859, (HALF)0xbcc6,
1973 (HALF)0x4b93, (HALF)0x18ff, (HALF)0xe691, (HALF)0x3068,
1974 (HALF)0xd823, (HALF)0x9db5, (HALF)0xd4ef, (HALF)0x72af,
1975 (HALF)0xcd1a, (HALF)0x5aa4, (HALF)0x3014, (HALF)0xa0a3,
1976 (HALF)0x49e6, (HALF)0x6b83, (HALF)0xe595, (HALF)0x2f4d,
1977 (HALF)0x0384, (HALF)0xd518, (HALF)0x4118, (HALF)0x19d9,
1978 (HALF)0x7534, (HALF)0x369e, (HALF)0x3b18, (HALF)0xce4d,
1979 (HALF)0xf9ee, (HALF)0x119e, (HALF)0x5c25, (HALF)0xe2f4,
1980 (HALF)0xe0a7, (HALF)0xe8ea, (HALF)0x1605, (HALF)0x62e4,
1981 (HALF)0xd2ca, (HALF)0x6346, (HALF)0x625b, (HALF)0x6425,
1982 (HALF)0x33de, (HALF)0x44b0, (HALF)0xe6b3, (HALF)0x1711,
1983 (HALF)0xd02e, (HALF)0xf3a9, (HALF)0xd965, (HALF)0x259c,
1984 (HALF)0x956b, (HALF)0x08aa, (HALF)0x4380, (HALF)0x6ad6,
1985 (HALF)0x0e8e, (HALF)0xe973, (HALF)0x9d28, (HALF)0x539b,
1986 (HALF)0x7950, (HALF)0xe540, (HALF)0x0be4, (HALF)0x8990,
1987 (HALF)0x18f6, (HALF)0xde12, (HALF)0xe52b, (HALF)0x63e1,
1988 (HALF)0x3f4e, (HALF)0x0de0, (HALF)0xa568, (HALF)0x8e21,
1989 (HALF)0x7ee3, (HALF)0x268d, (HALF)0x514e, (HALF)0xafb3,
1990 (HALF)0xefcb, (HALF)0x5378, (HALF)0xc7c6, (HALF)0xfec0,
1991 (HALF)0xb724, (HALF)0xf07c, (HALF)0xb42a, (HALF)0xfb61,
1992 (HALF)0x2a38, (HALF)0x068f
1993 };
1994 STATIC CONST HALF h_rvec18[] = {
1995 (HALF)0x3c63, (HALF)0x35ea, (HALF)0xef97, (HALF)0x8df2,
1996 (HALF)0xafb7, (HALF)0xa2b3, (HALF)0x58f6, (HALF)0x1791,
1997 (HALF)0x0dba, (HALF)0x0492, (HALF)0x077e, (HALF)0xf333,
1998 (HALF)0x4b5a, (HALF)0xf830, (HALF)0xf2ae, (HALF)0x230f,
1999 (HALF)0xf3f0, (HALF)0x84a8, (HALF)0x164e, (HALF)0xadda,
2000 (HALF)0xc944, (HALF)0xc9a1, (HALF)0x02f2, (HALF)0xc705,
2001 (HALF)0xc18f, (HALF)0x41a3, (HALF)0x3254, (HALF)0x09bd,
2002 (HALF)0x65a9, (HALF)0x9736, (HALF)0xc263, (HALF)0x1548,
2003 (HALF)0xd916, (HALF)0x5024, (HALF)0xdde9, (HALF)0x0a3d,
2004 (HALF)0xf1f5, (HALF)0xf2aa, (HALF)0xb92a, (HALF)0x666d,
2005 (HALF)0x5aa5, (HALF)0x3a53, (HALF)0x5775, (HALF)0x49c3,
2006 (HALF)0xa1c4, (HALF)0xc381, (HALF)0x6dbc, (HALF)0xf8d3,
2007 (HALF)0xe870, (HALF)0xe94b, (HALF)0x88a6, (HALF)0x430e,
2008 (HALF)0x9a06, (HALF)0x3721, (HALF)0xdf80, (HALF)0x5109,
2009 (HALF)0xa03f, (HALF)0xe73c, (HALF)0x4541, (HALF)0xf1bb,
2010 (HALF)0x32f3, (HALF)0x3c6f, (HALF)0xfc24, (HALF)0x952c,
2011 (HALF)0x697f, (HALF)0xfbc9, (HALF)0xd472, (HALF)0xc5b8,
2012 (HALF)0xda4b, (HALF)0xcbbe, (HALF)0xdb3b, (HALF)0x7303,
2013 (HALF)0x255e, (HALF)0xa18f, (HALF)0x353d, (HALF)0xe1d3,
2014 (HALF)0x8700, (HALF)0xe8c9, (HALF)0xe8fd, (HALF)0x9e75,
2015 (HALF)0x812e, (HALF)0x3ddd, (HALF)0xf891, (HALF)0x7340,
2016 (HALF)0x69ac, (HALF)0xb1f3, (HALF)0x505b, (HALF)0x764d,
2017 (HALF)0xf51b, (HALF)0xb13e, (HALF)0xaa43, (HALF)0x16c3,
2018 (HALF)0x42c4, (HALF)0x61d0, (HALF)0x0339, (HALF)0x22ac,
2019 (HALF)0x6fd1, (HALF)0x3f30, (HALF)0x6f7f, (HALF)0x4992,
2020 (HALF)0x575c, (HALF)0x2b4c, (HALF)0xb467, (HALF)0x9f3a,
2021 (HALF)0x65af, (HALF)0x5dac, (HALF)0x8dc7, (HALF)0x6277,
2022 (HALF)0x3a89, (HALF)0x9911, (HALF)0x540a, (HALF)0x49c0,
2023 (HALF)0x0ac2, (HALF)0x1df7, (HALF)0x2c5e, (HALF)0x4be1,
2024 (HALF)0x6bdb, (HALF)0xe5d3, (HALF)0x9ff5, (HALF)0x66b9,
2025 (HALF)0xbe89, (HALF)0x5358, (HALF)0x35d7, (HALF)0x4cd8,
2026 (HALF)0xcda8, (HALF)0xf0d5, (HALF)0xc6c3, (HALF)0x1f1a,
2027 (HALF)0x5e92, (HALF)0x0473
2028 };
2029 STATIC CONST HALF h_nvec19[] = {
2030 (HALF)0x9a79, (HALF)0x6b65, (HALF)0xc12d, (HALF)0x0239,
2031 (HALF)0xdf49, (HALF)0xd204, (HALF)0xe0c7, (HALF)0x1d4a,
2032 (HALF)0xf000, (HALF)0x099b, (HALF)0xade8, (HALF)0x6435,
2033 (HALF)0xf029, (HALF)0xdc4a, (HALF)0xe7a2, (HALF)0x2f4e,
2034 (HALF)0xf1e3, (HALF)0xadfc, (HALF)0x8f43, (HALF)0x7335,
2035 (HALF)0xede5, (HALF)0x687e, (HALF)0xcd4d, (HALF)0xb567,
2036 (HALF)0x814f, (HALF)0xc7a7, (HALF)0x624c, (HALF)0xc306,
2037 (HALF)0x80c6, (HALF)0xa82d, (HALF)0x0cd5, (HALF)0x3f39,
2038 (HALF)0xec3f, (HALF)0x7b7d, (HALF)0x1416, (HALF)0x8bdb,
2039 (HALF)0x3a52, (HALF)0x275b, (HALF)0x84fe, (HALF)0x9218,
2040 (HALF)0x02e0, (HALF)0x94ac, (HALF)0xb52c, (HALF)0x62f2,
2041 (HALF)0x92ee, (HALF)0xcdc9, (HALF)0x5eeb, (HALF)0x35e5,
2042 (HALF)0x3fd5, (HALF)0x69a4, (HALF)0xdcfb, (HALF)0x44c1,
2043 (HALF)0x6227, (HALF)0x3cdf, (HALF)0x148f, (HALF)0x23f3,
2044 (HALF)0x8e4c, (HALF)0x4250, (HALF)0x37c3, (HALF)0x95b7,
2045 (HALF)0x831f, (HALF)0x70af, (HALF)0x15c9, (HALF)0x2ee8,
2046 (HALF)0x7251, (HALF)0x47a4, (HALF)0xb2af, (HALF)0x7c11,
2047 (HALF)0x361f, (HALF)0x664c, (HALF)0xa841, (HALF)0xcada,
2048 (HALF)0x172a, (HALF)0x3d97, (HALF)0x2ffb, (HALF)0x8740,
2049 (HALF)0x2ebd, (HALF)0xa0a0, (HALF)0x225f, (HALF)0xb674,
2050 (HALF)0x3fe2, (HALF)0x6559, (HALF)0x98b9, (HALF)0x85f6,
2051 (HALF)0xa7ab, (HALF)0x5ed9, (HALF)0x7371, (HALF)0x6bf3,
2052 (HALF)0x5305, (HALF)0x7da7, (HALF)0x55bf, (HALF)0x0882,
2053 (HALF)0x7684, (HALF)0x8f60, (HALF)0xe57e, (HALF)0x1f3d,
2054 (HALF)0xd501, (HALF)0x118b, (HALF)0x770d, (HALF)0x6833,
2055 (HALF)0x5c51, (HALF)0xd042, (HALF)0xcacb, (HALF)0x2664,
2056 (HALF)0xb920, (HALF)0x9206, (HALF)0x3b8d, (HALF)0xeb90,
2057 (HALF)0x16b4, (HALF)0x0e25, (HALF)0xd841, (HALF)0x36c3,
2058 (HALF)0xcd17, (HALF)0x51d7, (HALF)0xba1e, (HALF)0xe063,
2059 (HALF)0x4a6d, (HALF)0xf28f, (HALF)0xeb0d, (HALF)0x244d,
2060 (HALF)0x0ad4, (HALF)0x2e41, (HALF)0xc315, (HALF)0x0721,
2061 (HALF)0x7654, (HALF)0xdde2, (HALF)0x534b, (HALF)0x2ad6,
2062 (HALF)0x8b25, (HALF)0xd678, (HALF)0xb9e8, (HALF)0xb23b,
2063 (HALF)0x0d7a, (HALF)0x0023
2064 };
2065 STATIC CONST HALF h_rvec19[] = {
2066 (HALF)0xf473, (HALF)0x698e, (HALF)0xa5b7, (HALF)0x3d53,
2067 (HALF)0x8319, (HALF)0x0644, (HALF)0x4445, (HALF)0xd9ad,
2068 (HALF)0xdaa0, (HALF)0x6967, (HALF)0x6240, (HALF)0xa14c,
2069 (HALF)0x7724, (HALF)0x78e7, (HALF)0x2ab7, (HALF)0x63ef,
2070 (HALF)0x2ee2, (HALF)0x8dff, (HALF)0xb424, (HALF)0x662e,
2071 (HALF)0x07d6, (HALF)0xcd93, (HALF)0x6a5d, (HALF)0x0ab0,
2072 (HALF)0xb539, (HALF)0xbdd7, (HALF)0xdd2d, (HALF)0x8621,
2073 (HALF)0xa187, (HALF)0x2bb4, (HALF)0x121d, (HALF)0x1f1e,
2074 (HALF)0xb962, (HALF)0xce8d, (HALF)0xeaf1, (HALF)0xdd3e,
2075 (HALF)0x6ca1, (HALF)0x573b, (HALF)0x0cd6, (HALF)0x6f46,
2076 (HALF)0x4780, (HALF)0x6a8d, (HALF)0xc7e6, (HALF)0xea68,
2077 (HALF)0xeb32, (HALF)0x1148, (HALF)0x44d4, (HALF)0xb43d,
2078 (HALF)0xcb64, (HALF)0xb657, (HALF)0xfdba, (HALF)0x8547,
2079 (HALF)0x3333, (HALF)0x85f7, (HALF)0xa51a, (HALF)0xc1f2,
2080 (HALF)0xee52, (HALF)0x1c05, (HALF)0xc03d, (HALF)0x2847,
2081 (HALF)0x88d8, (HALF)0xfc7a, (HALF)0xd186, (HALF)0x2e8d,
2082 (HALF)0x4683, (HALF)0x5be3, (HALF)0x2ee2, (HALF)0xd43d,
2083 (HALF)0x2bb5, (HALF)0x9f2d, (HALF)0xddea, (HALF)0x89b3,
2084 (HALF)0x4ae0, (HALF)0x7d78, (HALF)0x5e28, (HALF)0xc773,
2085 (HALF)0x8608, (HALF)0x967a, (HALF)0xcf07, (HALF)0xcdcb,
2086 (HALF)0xa423, (HALF)0xfc17, (HALF)0xd053, (HALF)0xd36a,
2087 (HALF)0x8892, (HALF)0xc73d, (HALF)0xc3f4, (HALF)0xa635,
2088 (HALF)0x0cf9, (HALF)0x9b5d, (HALF)0x3fd9, (HALF)0x0ac7,
2089 (HALF)0xfefb, (HALF)0xe801, (HALF)0xffd2, (HALF)0xb31c,
2090 (HALF)0xaa55, (HALF)0xf3ea, (HALF)0xfa23, (HALF)0x0e74,
2091 (HALF)0xb290, (HALF)0x5414, (HALF)0x6101, (HALF)0x6d17,
2092 (HALF)0x93f8, (HALF)0x5229, (HALF)0x3dad, (HALF)0xf829,
2093 (HALF)0x1e82, (HALF)0x713b, (HALF)0xbef0, (HALF)0xb83b,
2094 (HALF)0xcc62, (HALF)0xd001, (HALF)0xda39, (HALF)0x7537,
2095 (HALF)0x158b, (HALF)0x7d80, (HALF)0xc8e2, (HALF)0x9332,
2096 (HALF)0xfa75, (HALF)0x6fa6, (HALF)0x6512, (HALF)0xe21f,
2097 (HALF)0x18b4, (HALF)0x9995, (HALF)0x605b, (HALF)0x2196,
2098 (HALF)0x1798, (HALF)0xe4fc, (HALF)0xe245, (HALF)0x5f21,
2099 (HALF)0xf172, (HALF)0x0008
2100 };
2101 STATIC CONST HALF h_nvec20[] = {
2102 (HALF)0xd081, (HALF)0xc3c1, (HALF)0x2fce, (HALF)0x4d26,
2103 (HALF)0xcc91, (HALF)0x8765, (HALF)0x7f7c, (HALF)0xf372,
2104 (HALF)0x4bbc, (HALF)0xabba, (HALF)0x5801, (HALF)0xe098,
2105 (HALF)0x5c51, (HALF)0xfa36, (HALF)0xb230, (HALF)0xb2a4,
2106 (HALF)0x0a8d, (HALF)0xf443, (HALF)0x98c8, (HALF)0x546b,
2107 (HALF)0x8b26, (HALF)0xd974, (HALF)0xa82f, (HALF)0xe255,
2108 (HALF)0x3e8c, (HALF)0xb00a, (HALF)0x676d, (HALF)0x7069,
2109 (HALF)0xccce, (HALF)0x6233, (HALF)0xa74e, (HALF)0x0299,
2110 (HALF)0xadf9, (HALF)0xd119, (HALF)0xe811, (HALF)0x5273,
2111 (HALF)0xe0bb, (HALF)0xae36, (HALF)0xa486, (HALF)0x32b3,
2112 (HALF)0x21a5, (HALF)0xf049, (HALF)0x28a2, (HALF)0xd31d,
2113 (HALF)0x50de, (HALF)0xda0c, (HALF)0x2b40, (HALF)0x2130,
2114 (HALF)0xe552, (HALF)0xee9d, (HALF)0xac03, (HALF)0x80bb,
2115 (HALF)0x9740, (HALF)0x75f4, (HALF)0xa61b, (HALF)0xc3d0,
2116 (HALF)0xfbf6, (HALF)0x341c, (HALF)0xb5b4, (HALF)0x1615,
2117 (HALF)0x3def, (HALF)0x4e3a, (HALF)0x4dbf, (HALF)0x9573,
2118 (HALF)0x78ac, (HALF)0xe7ab, (HALF)0x5cf9, (HALF)0xffdc,
2119 (HALF)0x6892, (HALF)0xf799, (HALF)0x7ba4, (HALF)0x4740,
2120 (HALF)0x988f, (HALF)0xa78b, (HALF)0x6f42, (HALF)0x2773,
2121 (HALF)0xee50, (HALF)0x5064, (HALF)0x5c56, (HALF)0x6013,
2122 (HALF)0x3283, (HALF)0x1ad7, (HALF)0x4bf7, (HALF)0x2562,
2123 (HALF)0x1419, (HALF)0x2ee2, (HALF)0x5abd, (HALF)0x9319,
2124 (HALF)0x7778, (HALF)0x66b6, (HALF)0xa42d, (HALF)0xb1e0,
2125 (HALF)0xb0f0, (HALF)0x729f, (HALF)0x1864, (HALF)0xd492,
2126 (HALF)0x253f, (HALF)0x2c42, (HALF)0x07a9, (HALF)0x302a,
2127 (HALF)0x1bd4, (HALF)0xbb74, (HALF)0x90ba, (HALF)0x932f,
2128 (HALF)0x4be1, (HALF)0xf335, (HALF)0xd661, (HALF)0x0804,
2129 (HALF)0x9ba1, (HALF)0x010e, (HALF)0x778d, (HALF)0x1a05,
2130 (HALF)0xc833, (HALF)0xa962, (HALF)0x0ee8, (HALF)0xe759,
2131 (HALF)0x03b8, (HALF)0xbe68, (HALF)0x04c1, (HALF)0xc677,
2132 (HALF)0x9660, (HALF)0x56d7, (HALF)0xa3f3, (HALF)0x6066,
2133 (HALF)0x0327, (HALF)0x648b, (HALF)0x5b3a, (HALF)0x267e,
2134 (HALF)0x63a0, (HALF)0xdddc, (HALF)0xe890, (HALF)0x3322,
2135 (HALF)0xd8b1, (HALF)0x20e0, (HALF)0xd2b8, (HALF)0x004f
2136 };
2137 STATIC CONST HALF h_rvec20[] = {
2138 (HALF)0xbd1a, (HALF)0xa048, (HALF)0xdc7b, (HALF)0x95ab,
2139 (HALF)0x7cf8, (HALF)0x98f4, (HALF)0xc98d, (HALF)0x126a,
2140 (HALF)0x85fd, (HALF)0xaebf, (HALF)0x580f, (HALF)0x5650,
2141 (HALF)0xd7dd, (HALF)0x3292, (HALF)0x8377, (HALF)0xf49e,
2142 (HALF)0xed46, (HALF)0x2947, (HALF)0xb26c, (HALF)0xd1a5,
2143 (HALF)0xe6a1, (HALF)0xae14, (HALF)0x5788, (HALF)0x9b1f,
2144 (HALF)0x27b2, (HALF)0x4df7, (HALF)0x5079, (HALF)0xee37,
2145 (HALF)0xc8e4, (HALF)0x131b, (HALF)0x5f53, (HALF)0x294e,
2146 (HALF)0x59bd, (HALF)0x1f57, (HALF)0x8acf, (HALF)0x65d5,
2147 (HALF)0xd3a5, (HALF)0x598e, (HALF)0x61a6, (HALF)0xc393,
2148 (HALF)0xfd7a, (HALF)0xa783, (HALF)0x36a2, (HALF)0x264a,
2149 (HALF)0x2856, (HALF)0x6ca2, (HALF)0x171f, (HALF)0x8ffb,
2150 (HALF)0xea9e, (HALF)0x1d7c, (HALF)0x6fca, (HALF)0xd81d,
2151 (HALF)0x730e, (HALF)0x34ea, (HALF)0x6382, (HALF)0x31f5,
2152 (HALF)0xd9e9, (HALF)0xb39c, (HALF)0x84be, (HALF)0x440e,
2153 (HALF)0x15a1, (HALF)0x4b1d, (HALF)0x75c5, (HALF)0x7bf7,
2154 (HALF)0x4638, (HALF)0xe40f, (HALF)0xf0a7, (HALF)0xe5be,
2155 (HALF)0x8942, (HALF)0x79e5, (HALF)0xe1ba, (HALF)0x881a,
2156 (HALF)0x8372, (HALF)0x01de, (HALF)0x35f8, (HALF)0x14cf,
2157 (HALF)0xb310, (HALF)0xe2d8, (HALF)0x1207, (HALF)0x6696,
2158 (HALF)0x5f91, (HALF)0xde5d, (HALF)0x0849, (HALF)0xe6e7,
2159 (HALF)0x5ac3, (HALF)0x74ec, (HALF)0x4eb1, (HALF)0xe2de,
2160 (HALF)0xdc20, (HALF)0x4a41, (HALF)0xd565, (HALF)0xd306,
2161 (HALF)0x3ff3, (HALF)0xb584, (HALF)0x30d6, (HALF)0x911b,
2162 (HALF)0xd926, (HALF)0x4e9c, (HALF)0xc9ae, (HALF)0x8455,
2163 (HALF)0x8bb5, (HALF)0x6944, (HALF)0x1aad, (HALF)0x0c7b,
2164 (HALF)0xe992, (HALF)0x1da1, (HALF)0x56bd, (HALF)0xc676,
2165 (HALF)0x209e, (HALF)0xc544, (HALF)0x387c, (HALF)0x10ce,
2166 (HALF)0x8df8, (HALF)0xc4e8, (HALF)0xda88, (HALF)0x40e8,
2167 (HALF)0x3028, (HALF)0xbb2c, (HALF)0x4fd9, (HALF)0x4919,
2168 (HALF)0x17ee, (HALF)0xdeef, (HALF)0xc08d, (HALF)0x241b,
2169 (HALF)0xf608, (HALF)0x6fa9, (HALF)0x8b04, (HALF)0x4b0f,
2170 (HALF)0x0da1, (HALF)0xee96, (HALF)0x9293, (HALF)0xa309,
2171 (HALF)0x5fea, (HALF)0x8444, (HALF)0xef01, (HALF)0x0046
2172 };
2173 #else
2174 /\../\ FULL_BITS must be 32 or 64 /\../\ !!!
2175 #endif
2176 /*
2177 * NOTE: set n is found in random_pregen[n-1]
2178 */
2179 STATIC RANDOM random_pregen[BLUM_PREGEN] = {
2180 {1, 0, 7, (HALF)0, (HALF)0x07f,
2181 {(HALF *)h_nvec01, sizeof(h_nvec01)/sizeof(HALF), 0},
2182 {(HALF *)h_rvec01, sizeof(h_rvec01)/sizeof(HALF), 0}},
2183 {1, 0, 7, (HALF)0, (HALF)0x07f,
2184 {(HALF *)h_nvec02, sizeof(h_nvec02)/sizeof(HALF), 0},
2185 {(HALF *)h_rvec02, sizeof(h_rvec02)/sizeof(HALF), 0}},
2186 {1, 0, 7, (HALF)0, (HALF)0x07f,
2187 {(HALF *)h_nvec03, sizeof(h_nvec03)/sizeof(HALF), 0},
2188 {(HALF *)h_rvec03, sizeof(h_rvec03)/sizeof(HALF), 0}},
2189 {1, 0, 7, (HALF)0, (HALF)0x07f,
2190 {(HALF *)h_nvec04, sizeof(h_nvec04)/sizeof(HALF), 0},
2191 {(HALF *)h_rvec04, sizeof(h_rvec04)/sizeof(HALF), 0}},
2192 {1, 0, 8, (HALF)0, (HALF)0x0ff,
2193 {(HALF *)h_nvec05, sizeof(h_nvec05)/sizeof(HALF), 0},
2194 {(HALF *)h_rvec05, sizeof(h_rvec05)/sizeof(HALF), 0}},
2195 {1, 0, 8, (HALF)0, (HALF)0x0ff,
2196 {(HALF *)h_nvec06, sizeof(h_nvec06)/sizeof(HALF), 0},
2197 {(HALF *)h_rvec06, sizeof(h_rvec06)/sizeof(HALF), 0}},
2198 {1, 0, 8, (HALF)0, (HALF)0x0ff,
2199 {(HALF *)h_nvec07, sizeof(h_nvec07)/sizeof(HALF), 0},
2200 {(HALF *)h_rvec07, sizeof(h_rvec07)/sizeof(HALF), 0}},
2201 {1, 0, 8, (HALF)0, (HALF)0x0ff,
2202 {(HALF *)h_nvec08, sizeof(h_nvec08)/sizeof(HALF), 0},
2203 {(HALF *)h_rvec08, sizeof(h_rvec08)/sizeof(HALF), 0}},
2204 {1, 0, 9, (HALF)0, (HALF)0x1ff,
2205 {(HALF *)h_nvec09, sizeof(h_nvec09)/sizeof(HALF), 0},
2206 {(HALF *)h_rvec09, sizeof(h_rvec09)/sizeof(HALF), 0}},
2207 {1, 0, 9, (HALF)0, (HALF)0x1ff,
2208 {(HALF *)h_nvec10, sizeof(h_nvec10)/sizeof(HALF), 0},
2209 {(HALF *)h_rvec10, sizeof(h_rvec10)/sizeof(HALF), 0}},
2210 {1, 0, 9, (HALF)0, (HALF)0x1ff,
2211 {(HALF *)h_nvec11, sizeof(h_nvec11)/sizeof(HALF), 0},
2212 {(HALF *)h_rvec11, sizeof(h_rvec11)/sizeof(HALF), 0}},
2213 {1, 0, 9, (HALF)0, (HALF)0x1ff,
2214 {(HALF *)h_nvec12, sizeof(h_nvec12)/sizeof(HALF), 0},
2215 {(HALF *)h_rvec12, sizeof(h_rvec12)/sizeof(HALF), 0}},
2216 {1, 0, 10, (HALF)0, (HALF)0x3ff,
2217 {(HALF *)h_nvec13, sizeof(h_nvec13)/sizeof(HALF), 0},
2218 {(HALF *)h_rvec13, sizeof(h_rvec13)/sizeof(HALF), 0}},
2219 {1, 0, 10, (HALF)0, (HALF)0x3ff,
2220 {(HALF *)h_nvec14, sizeof(h_nvec14)/sizeof(HALF), 0},
2221 {(HALF *)h_rvec14, sizeof(h_rvec14)/sizeof(HALF), 0}},
2222 {1, 0, 10, (HALF)0, (HALF)0x3ff,
2223 {(HALF *)h_nvec15, sizeof(h_nvec15)/sizeof(HALF), 0},
2224 {(HALF *)h_rvec15, sizeof(h_rvec15)/sizeof(HALF), 0}},
2225 {1, 0, 10, (HALF)0, (HALF)0x3ff,
2226 {(HALF *)h_nvec16, sizeof(h_nvec16)/sizeof(HALF), 0},
2227 {(HALF *)h_rvec16, sizeof(h_rvec16)/sizeof(HALF), 0}},
2228 {1, 0, 11, (HALF)0, (HALF)0x7ff,
2229 {(HALF *)h_nvec17, sizeof(h_nvec17)/sizeof(HALF), 0},
2230 {(HALF *)h_rvec17, sizeof(h_rvec17)/sizeof(HALF), 0}},
2231 {1, 0, 11, (HALF)0, (HALF)0x7ff,
2232 {(HALF *)h_nvec18, sizeof(h_nvec18)/sizeof(HALF), 0},
2233 {(HALF *)h_rvec18, sizeof(h_rvec18)/sizeof(HALF), 0}},
2234 {1, 0, 11, (HALF)0, (HALF)0x7ff,
2235 {(HALF *)h_nvec19, sizeof(h_nvec19)/sizeof(HALF), 0},
2236 {(HALF *)h_rvec19, sizeof(h_rvec19)/sizeof(HALF), 0}},
2237 {1, 0, 11, (HALF)0, (HALF)0x7ff,
2238 {(HALF *)h_nvec20, sizeof(h_nvec20)/sizeof(HALF), 0},
2239 {(HALF *)h_rvec20, sizeof(h_rvec20)/sizeof(HALF), 0}}
2240 };
2241
2242
2243 /*
2244 * forward static declarations
2245 */
2246 S_FUNC void zfree_random(ZVALUE z);
2247
2248
2249 /*
2250 * zsrandom1 - seed the Blum generator 1 arg style
2251 *
2252 * We will seed the Blum generator according 2 argument
2253 * function description described at the top of this file.
2254 * In particular:
2255 *
2256 * given:
2257 * pseed - seed of the generator
2258 * need_ret - TRUE=>malloc return previous Blum state, FALSE=>return NULL
2259 *
2260 * returns:
2261 * previous Blum state
2262 */
2263 RANDOM *
zsrandom1(CONST ZVALUE seed,BOOL need_ret)2264 zsrandom1(CONST ZVALUE seed, BOOL need_ret)
2265 {
2266 RANDOM *ret; /* previous Blum state */
2267 ZVALUE r; /* quadratic residue */
2268 ZVALUE last_r; /* previous quadratic residue */
2269 RANDOM *p_blum; /* malloced RANDOM by randomcopy() */
2270
2271 /*
2272 * initialize state if first call
2273 */
2274 if (!blum.seeded) {
2275 p_blum = randomcopy(&init_blum);
2276 randomfree(&blum);
2277 blum = *p_blum;
2278 free(p_blum);
2279 }
2280
2281 /*
2282 * save the current state to return later, if need_ret says so
2283 */
2284 if (need_ret) {
2285 ret = randomcopy(&blum);
2286 } else {
2287 ret = NULL;
2288 }
2289
2290 /*
2291 * srandom(seed == 0)
2292 *
2293 * If the init arg is TRUE, then restore the initial state and
2294 * modulus of the Blum generator. After this call, the Blum
2295 * generator is restored to its initial state after calc started.
2296 */
2297 if (ziszero(seed)) {
2298
2299 /* set to the default generator state */
2300 p_blum = randomcopy(&init_blum);
2301 randomfree(&blum);
2302 blum = *p_blum;
2303 free(p_blum);
2304
2305 /*
2306 * srandom(seed >= 2^32)
2307 * srandom(seed >= 2^32, newn)
2308 *
2309 * Use seed to compute a new quadratic residue for use with
2310 * the current Blum modulus. We will successively square mod Blum
2311 * modulus until we get a smaller value (modulus wrap).
2312 *
2313 * The Blum modulus will not be changed.
2314 */
2315 } else if (!zisneg(seed) && zge32b(seed)) {
2316
2317 /*
2318 * square the seed mod the Blum modulus until we wrap
2319 */
2320 zcopy(seed, &r);
2321 last_r.v = NULL;
2322 do {
2323 /* free temp storage */
2324 zfree_random(last_r);
2325
2326 /*
2327 * last_r = r;
2328 * r = pmod(last_r, 2, n);
2329 */
2330 zcopy(r, &last_r);
2331 zfree_random(r);
2332 zsquaremod(last_r, blum.n, &r);
2333 } while (zrel(r, last_r) > 0);
2334 zfree_random(blum.r);
2335 blum.r = r;
2336
2337 /* free temp storage */
2338 zfree_random(last_r);
2339
2340 /*
2341 * reserved seed
2342 */
2343 } else {
2344 math_error("srandom seed must be 0 or >= 2^32");
2345 /*NOTREACHED*/
2346 }
2347
2348 /*
2349 * flush the queued up bits
2350 */
2351 blum.bits = 0;
2352 blum.buffer = 0;
2353
2354 /*
2355 * return the previous state
2356 */
2357 return ret;
2358 }
2359
2360
2361 /*
2362 * zsrandom2 - seed the Blum generator 2 arg style
2363 *
2364 * We will seed the Blum generator according 2 argument
2365 * function description described at the top of this file.
2366 * In particular:
2367 *
2368 * given:
2369 * pseed - seed of the generator
2370 * newn - ptr to proposed new n (Blum modulus)
2371 *
2372 * returns:
2373 * previous Blum state
2374 */
2375 RANDOM *
zsrandom2(CONST ZVALUE seed,CONST ZVALUE newn)2376 zsrandom2(CONST ZVALUE seed, CONST ZVALUE newn)
2377 {
2378 RANDOM *ret; /* previous Blum state */
2379 HALF set; /* pre-defined set to use */
2380 FULL nlen; /* length of newn in bits */
2381 RANDOM *p_blum; /* malloced RANDOM by randomcopy() */
2382
2383 /*
2384 * initialize state if first call
2385 */
2386 if (!blum.seeded) {
2387 p_blum = randomcopy(&init_blum);
2388 randomfree(&blum);
2389 blum = *p_blum;
2390 free(p_blum);
2391 }
2392
2393 /*
2394 * save the current state to return later
2395 */
2396 ret = randomcopy(&blum);
2397
2398 /*
2399 * srandom(seed, 0 < newn <= 20)
2400 *
2401 * Set the Blum modulus to one of the pre-defined Blum moduli.
2402 * The new quadratic residue will also be set to one of
2403 * the pre-defined quadratic residues.
2404 */
2405 if (!zisneg(newn) && !zge32b(newn)) {
2406
2407 /*
2408 * preset moduli only if small newn
2409 */
2410 if (ziszero(newn)) {
2411 math_error("srandom newn == 0 reserved for future use");
2412 /*NOTREACHED*/
2413 }
2414 set = (HALF)z1tol(newn);
2415 if (!zistiny(newn) || set > BLUM_PREGEN) {
2416 math_error("srandom small newn must be [1,20]");
2417 /*NOTREACHED*/
2418 }
2419 zfree_random(blum.n);
2420 zcopy(random_pregen[set-1].n, &blum.n);
2421 blum.loglogn = random_pregen[set-1].loglogn;
2422 blum.mask = random_pregen[set-1].mask;
2423
2424 /*
2425 * reset initial seed as well if seed is 0
2426 */
2427 if (ziszero(seed)) {
2428 zfree_random(blum.r);
2429 zcopy(random_pregen[set-1].r, &blum.r);
2430
2431 /*
2432 * Otherwise non-zero seeds are processed as 1 arg calls
2433 */
2434 } else {
2435 zsrandom1(seed, FALSE);
2436 }
2437
2438 /*
2439 * srandom(seed, newn >= 2^32)
2440 *
2441 * Assuming that 'newn' == 3 mod 4, then we will use it as
2442 * the Blum modulus.
2443 *
2444 * We will use the seed arg to compute a new quadratic residue.
2445 * We will successively square it mod Blum modulus until we get
2446 * a smaller value (modulus wrap).
2447 */
2448 } else if (!zisneg(newn)) {
2449
2450 /*
2451 * Blum modulus must be 1 mod 4
2452 */
2453 if (newn.v[0] % 4 != 1) {
2454 math_error("srandom large newn must be 1 mod 4");
2455 /*NOTREACHED*/
2456 }
2457
2458 /*
2459 * For correct Blum moduli, hope they are a product
2460 * of two primes.
2461 */
2462 /* load modulus */
2463 zfree_random(blum.n);
2464 zcopy(newn, &blum.n);
2465
2466 /*
2467 * setup loglogn and mask
2468 *
2469 * If the length if excessive, reduce it down
2470 * so that loglogn is at most BASEB-1.
2471 */
2472 nlen = (FULL)zhighbit(newn);
2473 blum.loglogn = BASEB-1;
2474 if (nlen > 0 && nlen <= TOPHALF) {
2475 for (blum.loglogn=BASEB-1;
2476 ((FULL)1<<blum.loglogn) > nlen && blum.loglogn > 1;
2477 --blum.loglogn) {
2478 }
2479 }
2480 blum.mask = ((HALF)1 << blum.loglogn)-1;
2481
2482 /*
2483 * use default initial seed if seed is 0 and process
2484 * as if this value is given as a 1 arg call
2485 */
2486 if (ziszero(seed)) {
2487 (void) zsrandom1(z_rdefault, FALSE);
2488
2489 /*
2490 * Otherwise non-zero seeds are processed as 1 arg calls
2491 */
2492 } else {
2493 (void) zsrandom1(seed, FALSE);
2494 }
2495
2496 /*
2497 * reserved newn
2498 */
2499 } else {
2500 math_error("srandom newn must be [1,20] or >= 2^32");
2501 /*NOTREACHED*/
2502 }
2503
2504 /*
2505 * flush the queued up bits
2506 */
2507 blum.bits = 0;
2508 blum.buffer = 0;
2509
2510 /*
2511 * return the previous state
2512 */
2513 return ret;
2514 }
2515
2516
2517 /*
2518 * zsrandom4 - seed the Blum generator 4 arg style
2519 *
2520 * We will seed the Blum generator according 2 argument
2521 * function description described at the top of this file.
2522 * In particular:
2523 *
2524 * given:
2525 * pseed - seed of the generator
2526 * ip - initial p search point
2527 * iq - initial q search point
2528 * trials - number of ptests to perform per candidate prime
2529 *
2530 * returns:
2531 * previous Blum state
2532 */
2533 RANDOM *
zsrandom4(CONST ZVALUE seed,CONST ZVALUE ip,CONST ZVALUE iq,long trials)2534 zsrandom4(CONST ZVALUE seed, CONST ZVALUE ip, CONST ZVALUE iq, long trials)
2535 {
2536 RANDOM *ret; /* previous Blum state */
2537 FULL nlen; /* length of n=p*q in bits */
2538 ZVALUE p; /* 1st Blum prime */
2539 ZVALUE q; /* 2nd Blum prime */
2540 RANDOM *p_blum; /* malloced RANDOM by randomcopy() */
2541
2542 /*
2543 * initialize state if first call
2544 */
2545 if (!blum.seeded) {
2546 p_blum = randomcopy(&init_blum);
2547 randomfree(&blum);
2548 blum = *p_blum;
2549 free(p_blum);
2550 }
2551
2552 /*
2553 * save the current state to return later
2554 */
2555 ret = randomcopy(&blum);
2556
2557 /*
2558 * search the 'p' and 'q' Blum prime (3 mod 4) candidates
2559 */
2560 if (!znextcand(ip, trials, _zero_, _three_, _four_, &p)) {
2561 math_error("failed to find 1st Blum prime");
2562 /*NOTREACHED*/
2563 }
2564 if (!znextcand(iq, trials, _zero_, _three_, _four_, &q)) {
2565 math_error("failed to find 2nd Blum prime");
2566 /*NOTREACHED*/
2567 }
2568
2569 /*
2570 * form the Blum modulus
2571 */
2572 zfree_random(blum.n);
2573 zmul(p, q, &blum.n);
2574 /* free temp storage */
2575 zfree_random(p);
2576 zfree_random(q);
2577
2578 /*
2579 * form the loglogn and mask
2580 */
2581 nlen = (FULL)zhighbit(blum.n);
2582 blum.loglogn = BASEB-1;
2583 if (nlen > 0 && nlen <= TOPHALF) {
2584 for (blum.loglogn=BASEB-1;
2585 ((FULL)1<<blum.loglogn) > nlen && blum.loglogn > 1;
2586 --blum.loglogn) {
2587 }
2588 }
2589 blum.mask = ((HALF)1 << blum.loglogn)-1;
2590
2591 /*
2592 * use default initial seed if seed is 0 and process
2593 * as if this value is given as a 1 arg call
2594 */
2595 if (ziszero(seed)) {
2596 (void) zsrandom1(z_rdefault, FALSE);
2597
2598 /*
2599 * Otherwise non-zero seeds are processed as 1 arg calls
2600 */
2601 } else {
2602 (void) zsrandom1(seed, FALSE);
2603 }
2604
2605 /*
2606 * flush the queued up bits
2607 */
2608 blum.bits = 0;
2609 blum.buffer = 0;
2610
2611 /*
2612 * return the previous state
2613 */
2614 return ret;
2615 }
2616
2617
2618 /*
2619 * zsetrandom - set the Blum generator state
2620 *
2621 * given:
2622 * state - the state to copy
2623 *
2624 * In particular:
2625 *
2626 * zsetrandom(pseed) is called by: srandom() and srandom(state)
2627 *
2628 * returns:
2629 * previous RANDOM state
2630 */
2631 RANDOM *
zsetrandom(CONST RANDOM * state)2632 zsetrandom(CONST RANDOM *state)
2633 {
2634 RANDOM *ret; /* previous Blum state */
2635 RANDOM *p_blum; /* malloced RANDOM by randomcopy() */
2636
2637 /*
2638 * initialize state if first call
2639 */
2640 if (!blum.seeded) {
2641 p_blum = randomcopy(&init_blum);
2642 randomfree(&blum);
2643 blum = *p_blum;
2644 free(p_blum);
2645 }
2646
2647 /*
2648 * save the current state to return later
2649 */
2650 ret = randomcopy(&blum);
2651
2652 /*
2653 * load the new state
2654 */
2655 if (state != NULL) {
2656 p_blum = randomcopy(state);
2657 blum = *p_blum;
2658 free(p_blum);
2659 }
2660
2661 /*
2662 * return the previous state
2663 */
2664 return ret;
2665 }
2666
2667
2668 /*
2669 * zrandomskip - skip s bits via the Blum-Blum-Shub generator
2670 *
2671 * given:
2672 * count - number of bits to be skipped
2673 */
2674 void
zrandomskip(long cnt)2675 zrandomskip(long cnt)
2676 {
2677 ZVALUE new_r; /* new quadratic residue */
2678 long loglogn; /* blum.loglogn */
2679 RANDOM *p_blum; /* malloced RANDOM by randomcopy() */
2680
2681 /*
2682 * initialize state if first call
2683 */
2684 if (!blum.seeded) {
2685 p_blum = randomcopy(&init_blum);
2686 randomfree(&blum);
2687 blum = *p_blum;
2688 free(p_blum);
2689 }
2690 loglogn = (long)blum.loglogn;
2691 new_r.len = 0; /* paranoia */
2692 new_r.v = NULL;
2693 new_r.sign = 0;
2694
2695 /*
2696 * skip required bits in the buffer
2697 */
2698 if (blum.bits > 0) {
2699
2700 /*
2701 * depending in if we have too few or too many in the buffer
2702 */
2703 if (blum.bits <= cnt) {
2704
2705 /* too few - just toss the buffer bits */
2706 cnt -= blum.bits;
2707 blum.bits = 0;
2708 blum.buffer = 0;
2709
2710 } else {
2711
2712 /* buffer contains more bits than we need to toss */
2713 blum.buffer >>= cnt;
2714 blum.bits -= cnt;
2715 return; /* skip need satisfied */
2716 }
2717 }
2718
2719 /*
2720 * skip loglogn bits at a time
2721 */
2722 while (cnt >= loglogn) {
2723
2724 /* turn the Blum-Blum-Shub crank */
2725 zsquaremod(blum.r, blum.n, &new_r);
2726 zfree_random(blum.r);
2727 blum.r = new_r;
2728 cnt -= blum.loglogn;
2729 }
2730
2731 /*
2732 * skip the final bits
2733 */
2734 if (cnt > 0) {
2735
2736 /* turn the Blum-Blum-Shub crank */
2737 zsquaremod(blum.r, blum.n, &new_r);
2738 zfree_random(blum.r);
2739 blum.r = new_r;
2740
2741 /* fill the buffer with the unused bits */
2742 blum.bits = loglogn - cnt;
2743 blum.buffer = (blum.r.v[0] & lowhalf[blum.bits]);
2744 }
2745 return;
2746 }
2747
2748
2749 /*
2750 * zrandom - crank the Blum-Blum-Shub generator for some bits
2751 *
2752 * We will load the ZVALUE with random bits starting from the
2753 * most significant and ending with the lowest bit in the
2754 * least significant HALF.
2755 *
2756 * given:
2757 * count - number of bits required
2758 * res - where to place the random bits as ZVALUE
2759 */
2760 void
zrandom(long cnt,ZVALUE * res)2761 zrandom(long cnt, ZVALUE *res)
2762 {
2763 BITSTR dest; /* destination bit string */
2764 int loglogn; /* blum.loglogn */
2765 HALF mask; /* mask for bottom loglogn bits */
2766 ZVALUE new_r; /* new quadratic residue */
2767 RANDOM *p_blum; /* malloced RANDOM by randomcopy() */
2768 int t; /* temp shift value */
2769
2770 /*
2771 * firewall
2772 */
2773 if (cnt <= 0) {
2774 if (cnt == 0) {
2775 /* zero length random number is always 0 */
2776 itoz(0, res);
2777 return;
2778 } else {
2779 math_error("negative zrandom bit count");
2780 /*NOTREACHED*/
2781 }
2782 #if LONG_BITS > 32
2783 } else if (cnt > (1L<<31)) {
2784 math_error("huge random count in internal zrandom function");
2785 /*NOTREACHED*/
2786 #endif
2787 }
2788
2789 /*
2790 * initialize state if first call
2791 */
2792 if (!blum.seeded) {
2793 p_blum = randomcopy(&init_blum);
2794 randomfree(&blum);
2795 blum = *p_blum;
2796 free(p_blum);
2797 }
2798 loglogn = blum.loglogn;
2799 mask = blum.mask;
2800 new_r.len = 0; /* paranoia */
2801 new_r.v = NULL;
2802 new_r.sign = 0;
2803
2804 /*
2805 * allocate storage
2806 */
2807 res->len = (LEN)((cnt+BASEB-1)/BASEB);
2808 res->v = alloc((LEN)((cnt+BASEB-1)/BASEB));
2809
2810 /*
2811 * dest bit string
2812 */
2813 dest.len = (int)cnt;
2814 dest.loc = res->v + (((cnt+BASEB-1)/BASEB)-1);
2815 dest.bit = (int)((cnt-1) % BASEB);
2816 *dest.loc = 0;
2817
2818 /*
2819 * load from buffer first
2820 */
2821 if (blum.bits > 0) {
2822
2823 /*
2824 * If we need only part of the buffer, use
2825 * the top bits and keep the bottom in place.
2826 * If we need exactly all of the buffer,
2827 * process it as a partial buffer fill.
2828 */
2829 if (dest.len <= blum.bits) {
2830
2831 /* load part of the buffer */
2832 *dest.loc = (blum.buffer >> (blum.bits-dest.len));
2833
2834 /* update buffer */
2835 blum.buffer &= ((1 << (blum.bits-dest.len))-1);
2836 blum.bits -= dest.len;
2837
2838 /* cleanup */
2839 res->sign = 0;
2840 ztrim(res);
2841
2842 /* we are done now */
2843 return;
2844 }
2845
2846 /*
2847 * Otherwise we need all of the buffer and then some ...
2848 *
2849 * dest.len > blum.bits
2850 *
2851 * NOTE: We use = instead of |= as this will ensure that
2852 * bit bits above dest.bit are set to 0.
2853 */
2854 if (dest.bit >= blum.bits) {
2855 /* copy all of buffer into upper element */
2856 *dest.loc = (blum.buffer << (dest.bit+1-blum.bits));
2857 dest.bit -= blum.bits;
2858 } else {
2859 /* copy buffer into upper and next element */
2860 t = blum.bits-(dest.bit+1);
2861 *dest.loc-- = (blum.buffer >> t);
2862 dest.bit = BASEB-t-1;
2863 *dest.loc = ((blum.buffer&lowhalf[t]) << (dest.bit+1));
2864 }
2865 dest.len -= blum.bits;
2866 }
2867
2868 /*
2869 * Crank the generator up until, but not including, the
2870 * time when we will write into the least significant bit.
2871 *
2872 * In this loop we know that we have exactly blum.loglogn bits
2873 * to load.
2874 */
2875 while (dest.len > loglogn) {
2876
2877 /*
2878 * turn the Blum-Blum-Shub crank
2879 */
2880 zsquaremod(blum.r, blum.n, &new_r);
2881 zfree_random(blum.r);
2882 blum.r = new_r;
2883 /* peal off the bottom loglogn bits */
2884 blum.buffer = (blum.r.v[0] & mask);
2885
2886 /*
2887 * load the loglogn bits into dest.loc starting at bit dest.bit
2888 */
2889 if (dest.bit >= loglogn) {
2890 /* copy all of buffer into upper element */
2891 *dest.loc |= (blum.buffer << (dest.bit+1-loglogn));
2892 dest.bit -= loglogn;
2893 } else {
2894 /* copy buffer into upper and next element */
2895 t = loglogn-(dest.bit+1);
2896 *dest.loc-- |= (blum.buffer >> t);
2897 dest.bit = BASEB-t-1;
2898 *dest.loc = ((blum.buffer&lowhalf[t]) << (dest.bit+1));
2899 }
2900 dest.len -= loglogn;
2901 }
2902
2903 /*
2904 * We have a full or less than a full crank (loglogn bits) left
2905 * to generate and load into the least significant bits.
2906 *
2907 * If we have any bits left over, we will save them in the
2908 * buffer for use by the next call.
2909 */
2910 /* turn the Blum-Blum-Shub crank */
2911 zsquaremod(blum.r, blum.n, &new_r);
2912 zfree_random(blum.r);
2913 blum.r = new_r;
2914 /* peal off the bottom loglogn bits */
2915 blum.buffer = (blum.r.v[0] & mask);
2916 blum.bits = loglogn;
2917
2918 /*
2919 * load dest.len bits into the lowest order bits
2920 */
2921 *dest.loc |= (blum.buffer >> (loglogn - dest.len));
2922
2923 /*
2924 * leave any unused bits in the buffer for next time
2925 */
2926 blum.buffer &= lowhalf[loglogn - dest.len];
2927 blum.bits -= dest.len;
2928
2929 /*
2930 * cleanup
2931 */
2932 res->sign = 0;
2933 ztrim(res);
2934 }
2935
2936
2937 /*
2938 * zrandomrange - generate a Blum-Blum-Shub random value in [low, beyond)
2939 *
2940 * given:
2941 * low - low value of range
2942 * beyond - beyond end of range
2943 * res - where to place the random bits as ZVALUE
2944 */
2945 void
zrandomrange(CONST ZVALUE low,CONST ZVALUE beyond,ZVALUE * res)2946 zrandomrange(CONST ZVALUE low, CONST ZVALUE beyond, ZVALUE *res)
2947 {
2948 ZVALUE range; /* beyond-low */
2949 ZVALUE rval; /* random value [0, 2^bitlen) */
2950 ZVALUE rangem1; /* range - 1 */
2951 long bitlen; /* smallest power of 2 >= diff */
2952
2953 /*
2954 * firewall
2955 */
2956 if (zrel(low, beyond) >= 0) {
2957 math_error("srand low range >= beyond range");
2958 /*NOTREACHED*/
2959 }
2960
2961 /*
2962 * determine the size of the random number needed
2963 */
2964 zsub(beyond, low, &range);
2965 if (zisone(range)) {
2966 zfree_random(range);
2967 zcopy(low, res);
2968 return;
2969 }
2970 zsub(range, _one_, &rangem1);
2971 bitlen = 1+zhighbit(rangem1);
2972 zfree_random(rangem1);
2973
2974 /*
2975 * generate a random value between [0, diff)
2976 *
2977 * We will not fall into the trap of thinking that we can simply take
2978 * a value mod 'range'. Consider the case where 'range' is '80'
2979 * and we are given pseudo-random numbers [0,100). If we took them
2980 * mod 80, then the numbers [0,20) would be produced more frequently
2981 * because the numbers [81,100) mod 80 wrap back into [0,20).
2982 */
2983 rval.v = NULL;
2984 do {
2985 if (rval.v != NULL) {
2986 zfree_random(rval);
2987 }
2988 zrandom(bitlen, &rval);
2989 } while (zrel(rval, range) >= 0);
2990
2991 /*
2992 * add in low value to produce the range [0+low, diff+low)
2993 * which is the range [low, beyond)
2994 */
2995 zadd(rval, low, res);
2996 zfree_random(rval);
2997 zfree_random(range);
2998 }
2999
3000
3001 /*
3002 * irandom - generate a Blum-Blum-Shub random long in the range [0, s)
3003 *
3004 * given:
3005 * s - limit of the range
3006 *
3007 * returns:
3008 * random long in the range [0, s)
3009 */
3010 long
irandom(long s)3011 irandom(long s)
3012 {
3013 ZVALUE z1, z2;
3014 long res;
3015
3016 if (s <= 0) {
3017 math_error("Non-positive argument for irandom()");
3018 /*NOTREACHED*/
3019 }
3020 if (s == 1)
3021 return 0;
3022 itoz(s, &z1);
3023 zrandomrange(_zero_, z1, &z2);
3024 res = ztoi(z2);
3025 zfree_random(z1);
3026 zfree_random(z2);
3027 return res;
3028 }
3029
3030
3031 /*
3032 * randomcopy - make a copy of a Blum state
3033 *
3034 * given:
3035 * state - the state to copy
3036 *
3037 * returns:
3038 * a malloced copy of the state
3039 */
3040 RANDOM *
randomcopy(CONST RANDOM * state)3041 randomcopy(CONST RANDOM *state)
3042 {
3043 RANDOM *ret; /* return copy of state */
3044
3045 /*
3046 * malloc state
3047 */
3048 ret = (RANDOM *)malloc(sizeof(RANDOM));
3049 if (ret == NULL) {
3050 math_error("can't allocate RANDOM state");
3051 /*NOTREACHED*/
3052 }
3053
3054 /*
3055 * clone data
3056 */
3057 *ret = *state;
3058 if (state->r.v == NULL) {
3059 ret->r.v = NULL;
3060 } else {
3061 if (state->r.v == h_rdefvec) {
3062 ret->r.v = state->r.v;
3063 } else {
3064 zcopy(state->r, &ret->r);
3065 }
3066 }
3067 if (state->n.v == NULL) {
3068 ret->n.v = NULL;
3069 } else {
3070 if (state->n.v == h_ndefvec) {
3071 ret->n.v = state->n.v;
3072 } else {
3073 zcopy(state->n, &ret->n);
3074 }
3075 }
3076
3077 /*
3078 * return copy
3079 */
3080 return ret;
3081 }
3082
3083
3084 /*
3085 * randomfree - free a Blum state
3086 *
3087 * We avoid freeing the pre-compiled states as they were
3088 * never malloced in the first place.
3089 *
3090 * given:
3091 * state - the state to free
3092 */
3093 void
randomfree(RANDOM * state)3094 randomfree(RANDOM *state)
3095 {
3096 /* free the values */
3097 zfree_random(state->n);
3098 zfree_random(state->r);
3099
3100 /* free it if it is not pre-defined */
3101 state->seeded = 0;
3102 state->bits = 0; /* paranoia */
3103 state->buffer = 0;
3104
3105 /* free it if it is not pre-defined */
3106 if (state != &blum) {
3107 free(state);
3108 }
3109 }
3110
3111
3112 /*
3113 * randomcmp - compare two Blum states
3114 *
3115 * given:
3116 * s1 - first state to compare
3117 * s2 - second state to compare
3118 *
3119 * return:
3120 * TRUE if states differ
3121 */
3122 BOOL
randomcmp(CONST RANDOM * s1,CONST RANDOM * s2)3123 randomcmp(CONST RANDOM *s1, CONST RANDOM *s2)
3124 {
3125 /*
3126 * assume uninitialized state == the default seeded state
3127 */
3128 if (!s1->seeded) {
3129 if (!s2->seeded) {
3130 /* uninitialized == uninitialized */
3131 return FALSE;
3132 } else {
3133 /* uninitialized only equals default state */
3134 return randomcmp(s2, &init_blum);
3135 }
3136 } else if (!s2->seeded) {
3137 /* uninitialized only equals default state */
3138 return randomcmp(s1, &init_blum);
3139 }
3140
3141 /*
3142 * compare operating mask parameters
3143 */
3144 if ((s1->loglogn != s2->loglogn) || (s1->mask != s2->mask)) {
3145 return TRUE;
3146 }
3147
3148 /*
3149 * compare bit buffer
3150 */
3151 if ((s1->bits != s2->bits) || (s1->buffer != s2->buffer)) {
3152 return TRUE;
3153 }
3154
3155 /*
3156 * compare quadratic residues and moduli
3157 */
3158 return (zcmp(s1->r, s2->r) && zcmp(s1->n, s2->n));
3159 }
3160
3161
3162 /*
3163 * randomprint - print a Blum state
3164 *
3165 * given:
3166 * state - state to print
3167 * flags - print flags passed from printvalue() in value.c
3168 */
3169 /*ARGSUSED*/
3170 void
randomprint(CONST RANDOM * UNUSED (state),int UNUSED (flags))3171 randomprint(CONST RANDOM *UNUSED(state), int UNUSED(flags))
3172 {
3173 math_str("RANDOM state");
3174 }
3175
3176
3177 /*
3178 * random_libcalc_cleanup - cleanup code for final libcalc_call_me_last() call
3179 *
3180 * This call is needed only by libcalc_call_me_last() to help clean up any
3181 * unneeded storage.
3182 *
3183 * Do not call this function directly! Let libcalc_call_me_last() do it.
3184 */
3185 void
random_libcalc_cleanup(void)3186 random_libcalc_cleanup(void)
3187 {
3188 /* free our state - let zfree_random protect the default state */
3189 randomfree(&blum);
3190 return;
3191 }
3192
3193
3194 /*
3195 * zfree_random - perform a zfree if we are not trying to free static data
3196 *
3197 * given:
3198 * z the ZVALUE to zfree(z) if not pointing to static data
3199 */
3200 S_FUNC void
zfree_random(ZVALUE z)3201 zfree_random(ZVALUE z)
3202 {
3203 if (z.v != NULL &&
3204 z.v != h_ndefvec && z.v != h_rdefvec && z.v != h_rdefvec_2 &&
3205 z.v != h_nvec01 && z.v != h_rvec01 &&
3206 z.v != h_nvec02 && z.v != h_rvec02 &&
3207 z.v != h_nvec03 && z.v != h_rvec03 &&
3208 z.v != h_nvec04 && z.v != h_rvec04 &&
3209 z.v != h_nvec05 && z.v != h_rvec05 &&
3210 z.v != h_nvec06 && z.v != h_rvec06 &&
3211 z.v != h_nvec07 && z.v != h_rvec07 &&
3212 z.v != h_nvec08 && z.v != h_rvec08 &&
3213 z.v != h_nvec09 && z.v != h_rvec09 &&
3214 z.v != h_nvec10 && z.v != h_rvec10 &&
3215 z.v != h_nvec11 && z.v != h_rvec11 &&
3216 z.v != h_nvec12 && z.v != h_rvec12 &&
3217 z.v != h_nvec13 && z.v != h_rvec13 &&
3218 z.v != h_nvec14 && z.v != h_rvec14 &&
3219 z.v != h_nvec15 && z.v != h_rvec15 &&
3220 z.v != h_nvec16 && z.v != h_rvec16 &&
3221 z.v != h_nvec17 && z.v != h_rvec17 &&
3222 z.v != h_nvec18 && z.v != h_rvec18 &&
3223 z.v != h_nvec19 && z.v != h_rvec19 &&
3224 z.v != h_nvec20 && z.v != h_rvec20) {
3225 zfree(z);
3226 }
3227 return;
3228 }
3229