c_cpp/etc/calc/zrandom.c

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