SFMT.c revision b1941c615023cab9baf0a78a28df1e3b4972434f
1/* 2 * This file derives from SFMT 1.3.3 3 * (http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/index.html), which was 4 * released under the terms of the following license: 5 * 6 * Copyright (c) 2006,2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima 7 * University. All rights reserved. 8 * 9 * Redistribution and use in source and binary forms, with or without 10 * modification, are permitted provided that the following conditions are 11 * met: 12 * 13 * * Redistributions of source code must retain the above copyright 14 * notice, this list of conditions and the following disclaimer. 15 * * Redistributions in binary form must reproduce the above 16 * copyright notice, this list of conditions and the following 17 * disclaimer in the documentation and/or other materials provided 18 * with the distribution. 19 * * Neither the name of the Hiroshima University nor the names of 20 * its contributors may be used to endorse or promote products 21 * derived from this software without specific prior written 22 * permission. 23 * 24 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 25 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 26 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 27 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 28 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 29 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 30 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 31 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 32 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 33 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 34 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 35 */ 36/** 37 * @file SFMT.c 38 * @brief SIMD oriented Fast Mersenne Twister(SFMT) 39 * 40 * @author Mutsuo Saito (Hiroshima University) 41 * @author Makoto Matsumoto (Hiroshima University) 42 * 43 * Copyright (C) 2006,2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima 44 * University. All rights reserved. 45 * 46 * The new BSD License is applied to this software, see LICENSE.txt 47 */ 48#define SFMT_C_ 49#include "test/jemalloc_test.h" 50#include "test/SFMT-params.h" 51 52#if defined(__BIG_ENDIAN__) && !defined(__amd64) && !defined(BIG_ENDIAN64) 53#define BIG_ENDIAN64 1 54#endif 55#if defined(HAVE_ALTIVEC) && !defined(BIG_ENDIAN64) 56#define BIG_ENDIAN64 1 57#endif 58#if defined(ONLY64) && !defined(BIG_ENDIAN64) 59 #if defined(__GNUC__) 60 #error "-DONLY64 must be specified with -DBIG_ENDIAN64" 61 #endif 62#undef ONLY64 63#endif 64/*------------------------------------------------------ 65 128-bit SIMD data type for Altivec, SSE2 or standard C 66 ------------------------------------------------------*/ 67#if defined(HAVE_ALTIVEC) 68/** 128-bit data structure */ 69union W128_T { 70 vector unsigned int s; 71 uint32_t u[4]; 72}; 73/** 128-bit data type */ 74typedef union W128_T w128_t; 75 76#elif defined(HAVE_SSE2) 77/** 128-bit data structure */ 78union W128_T { 79 __m128i si; 80 uint32_t u[4]; 81}; 82/** 128-bit data type */ 83typedef union W128_T w128_t; 84 85#else 86 87/** 128-bit data structure */ 88struct W128_T { 89 uint32_t u[4]; 90}; 91/** 128-bit data type */ 92typedef struct W128_T w128_t; 93 94#endif 95 96struct sfmt_s { 97 /** the 128-bit internal state array */ 98 w128_t sfmt[N]; 99 /** index counter to the 32-bit internal state array */ 100 int idx; 101 /** a flag: it is 0 if and only if the internal state is not yet 102 * initialized. */ 103 int initialized; 104}; 105 106/*-------------------------------------- 107 FILE GLOBAL VARIABLES 108 internal state, index counter and flag 109 --------------------------------------*/ 110 111/** a parity check vector which certificate the period of 2^{MEXP} */ 112static uint32_t parity[4] = {PARITY1, PARITY2, PARITY3, PARITY4}; 113 114/*---------------- 115 STATIC FUNCTIONS 116 ----------------*/ 117inline static int idxof(int i); 118#if (!defined(HAVE_ALTIVEC)) && (!defined(HAVE_SSE2)) 119inline static void rshift128(w128_t *out, w128_t const *in, int shift); 120inline static void lshift128(w128_t *out, w128_t const *in, int shift); 121#endif 122inline static void gen_rand_all(sfmt_t *ctx); 123inline static void gen_rand_array(sfmt_t *ctx, w128_t *array, int size); 124inline static uint32_t func1(uint32_t x); 125inline static uint32_t func2(uint32_t x); 126static void period_certification(sfmt_t *ctx); 127#if defined(BIG_ENDIAN64) && !defined(ONLY64) 128inline static void swap(w128_t *array, int size); 129#endif 130 131#if defined(HAVE_ALTIVEC) 132 #include "test/SFMT-alti.h" 133#elif defined(HAVE_SSE2) 134 #include "test/SFMT-sse2.h" 135#endif 136 137/** 138 * This function simulate a 64-bit index of LITTLE ENDIAN 139 * in BIG ENDIAN machine. 140 */ 141#ifdef ONLY64 142inline static int idxof(int i) { 143 return i ^ 1; 144} 145#else 146inline static int idxof(int i) { 147 return i; 148} 149#endif 150/** 151 * This function simulates SIMD 128-bit right shift by the standard C. 152 * The 128-bit integer given in in is shifted by (shift * 8) bits. 153 * This function simulates the LITTLE ENDIAN SIMD. 154 * @param out the output of this function 155 * @param in the 128-bit data to be shifted 156 * @param shift the shift value 157 */ 158#if (!defined(HAVE_ALTIVEC)) && (!defined(HAVE_SSE2)) 159#ifdef ONLY64 160inline static void rshift128(w128_t *out, w128_t const *in, int shift) { 161 uint64_t th, tl, oh, ol; 162 163 th = ((uint64_t)in->u[2] << 32) | ((uint64_t)in->u[3]); 164 tl = ((uint64_t)in->u[0] << 32) | ((uint64_t)in->u[1]); 165 166 oh = th >> (shift * 8); 167 ol = tl >> (shift * 8); 168 ol |= th << (64 - shift * 8); 169 out->u[0] = (uint32_t)(ol >> 32); 170 out->u[1] = (uint32_t)ol; 171 out->u[2] = (uint32_t)(oh >> 32); 172 out->u[3] = (uint32_t)oh; 173} 174#else 175inline static void rshift128(w128_t *out, w128_t const *in, int shift) { 176 uint64_t th, tl, oh, ol; 177 178 th = ((uint64_t)in->u[3] << 32) | ((uint64_t)in->u[2]); 179 tl = ((uint64_t)in->u[1] << 32) | ((uint64_t)in->u[0]); 180 181 oh = th >> (shift * 8); 182 ol = tl >> (shift * 8); 183 ol |= th << (64 - shift * 8); 184 out->u[1] = (uint32_t)(ol >> 32); 185 out->u[0] = (uint32_t)ol; 186 out->u[3] = (uint32_t)(oh >> 32); 187 out->u[2] = (uint32_t)oh; 188} 189#endif 190/** 191 * This function simulates SIMD 128-bit left shift by the standard C. 192 * The 128-bit integer given in in is shifted by (shift * 8) bits. 193 * This function simulates the LITTLE ENDIAN SIMD. 194 * @param out the output of this function 195 * @param in the 128-bit data to be shifted 196 * @param shift the shift value 197 */ 198#ifdef ONLY64 199inline static void lshift128(w128_t *out, w128_t const *in, int shift) { 200 uint64_t th, tl, oh, ol; 201 202 th = ((uint64_t)in->u[2] << 32) | ((uint64_t)in->u[3]); 203 tl = ((uint64_t)in->u[0] << 32) | ((uint64_t)in->u[1]); 204 205 oh = th << (shift * 8); 206 ol = tl << (shift * 8); 207 oh |= tl >> (64 - shift * 8); 208 out->u[0] = (uint32_t)(ol >> 32); 209 out->u[1] = (uint32_t)ol; 210 out->u[2] = (uint32_t)(oh >> 32); 211 out->u[3] = (uint32_t)oh; 212} 213#else 214inline static void lshift128(w128_t *out, w128_t const *in, int shift) { 215 uint64_t th, tl, oh, ol; 216 217 th = ((uint64_t)in->u[3] << 32) | ((uint64_t)in->u[2]); 218 tl = ((uint64_t)in->u[1] << 32) | ((uint64_t)in->u[0]); 219 220 oh = th << (shift * 8); 221 ol = tl << (shift * 8); 222 oh |= tl >> (64 - shift * 8); 223 out->u[1] = (uint32_t)(ol >> 32); 224 out->u[0] = (uint32_t)ol; 225 out->u[3] = (uint32_t)(oh >> 32); 226 out->u[2] = (uint32_t)oh; 227} 228#endif 229#endif 230 231/** 232 * This function represents the recursion formula. 233 * @param r output 234 * @param a a 128-bit part of the internal state array 235 * @param b a 128-bit part of the internal state array 236 * @param c a 128-bit part of the internal state array 237 * @param d a 128-bit part of the internal state array 238 */ 239#if (!defined(HAVE_ALTIVEC)) && (!defined(HAVE_SSE2)) 240#ifdef ONLY64 241inline static void do_recursion(w128_t *r, w128_t *a, w128_t *b, w128_t *c, 242 w128_t *d) { 243 w128_t x; 244 w128_t y; 245 246 lshift128(&x, a, SL2); 247 rshift128(&y, c, SR2); 248 r->u[0] = a->u[0] ^ x.u[0] ^ ((b->u[0] >> SR1) & MSK2) ^ y.u[0] 249 ^ (d->u[0] << SL1); 250 r->u[1] = a->u[1] ^ x.u[1] ^ ((b->u[1] >> SR1) & MSK1) ^ y.u[1] 251 ^ (d->u[1] << SL1); 252 r->u[2] = a->u[2] ^ x.u[2] ^ ((b->u[2] >> SR1) & MSK4) ^ y.u[2] 253 ^ (d->u[2] << SL1); 254 r->u[3] = a->u[3] ^ x.u[3] ^ ((b->u[3] >> SR1) & MSK3) ^ y.u[3] 255 ^ (d->u[3] << SL1); 256} 257#else 258inline static void do_recursion(w128_t *r, w128_t *a, w128_t *b, w128_t *c, 259 w128_t *d) { 260 w128_t x; 261 w128_t y; 262 263 lshift128(&x, a, SL2); 264 rshift128(&y, c, SR2); 265 r->u[0] = a->u[0] ^ x.u[0] ^ ((b->u[0] >> SR1) & MSK1) ^ y.u[0] 266 ^ (d->u[0] << SL1); 267 r->u[1] = a->u[1] ^ x.u[1] ^ ((b->u[1] >> SR1) & MSK2) ^ y.u[1] 268 ^ (d->u[1] << SL1); 269 r->u[2] = a->u[2] ^ x.u[2] ^ ((b->u[2] >> SR1) & MSK3) ^ y.u[2] 270 ^ (d->u[2] << SL1); 271 r->u[3] = a->u[3] ^ x.u[3] ^ ((b->u[3] >> SR1) & MSK4) ^ y.u[3] 272 ^ (d->u[3] << SL1); 273} 274#endif 275#endif 276 277#if (!defined(HAVE_ALTIVEC)) && (!defined(HAVE_SSE2)) 278/** 279 * This function fills the internal state array with pseudorandom 280 * integers. 281 */ 282inline static void gen_rand_all(sfmt_t *ctx) { 283 int i; 284 w128_t *r1, *r2; 285 286 r1 = &ctx->sfmt[N - 2]; 287 r2 = &ctx->sfmt[N - 1]; 288 for (i = 0; i < N - POS1; i++) { 289 do_recursion(&ctx->sfmt[i], &ctx->sfmt[i], &ctx->sfmt[i + POS1], r1, 290 r2); 291 r1 = r2; 292 r2 = &ctx->sfmt[i]; 293 } 294 for (; i < N; i++) { 295 do_recursion(&ctx->sfmt[i], &ctx->sfmt[i], &ctx->sfmt[i + POS1 - N], r1, 296 r2); 297 r1 = r2; 298 r2 = &ctx->sfmt[i]; 299 } 300} 301 302/** 303 * This function fills the user-specified array with pseudorandom 304 * integers. 305 * 306 * @param array an 128-bit array to be filled by pseudorandom numbers. 307 * @param size number of 128-bit pseudorandom numbers to be generated. 308 */ 309inline static void gen_rand_array(sfmt_t *ctx, w128_t *array, int size) { 310 int i, j; 311 w128_t *r1, *r2; 312 313 r1 = &ctx->sfmt[N - 2]; 314 r2 = &ctx->sfmt[N - 1]; 315 for (i = 0; i < N - POS1; i++) { 316 do_recursion(&array[i], &ctx->sfmt[i], &ctx->sfmt[i + POS1], r1, r2); 317 r1 = r2; 318 r2 = &array[i]; 319 } 320 for (; i < N; i++) { 321 do_recursion(&array[i], &ctx->sfmt[i], &array[i + POS1 - N], r1, r2); 322 r1 = r2; 323 r2 = &array[i]; 324 } 325 for (; i < size - N; i++) { 326 do_recursion(&array[i], &array[i - N], &array[i + POS1 - N], r1, r2); 327 r1 = r2; 328 r2 = &array[i]; 329 } 330 for (j = 0; j < 2 * N - size; j++) { 331 ctx->sfmt[j] = array[j + size - N]; 332 } 333 for (; i < size; i++, j++) { 334 do_recursion(&array[i], &array[i - N], &array[i + POS1 - N], r1, r2); 335 r1 = r2; 336 r2 = &array[i]; 337 ctx->sfmt[j] = array[i]; 338 } 339} 340#endif 341 342#if defined(BIG_ENDIAN64) && !defined(ONLY64) && !defined(HAVE_ALTIVEC) 343inline static void swap(w128_t *array, int size) { 344 int i; 345 uint32_t x, y; 346 347 for (i = 0; i < size; i++) { 348 x = array[i].u[0]; 349 y = array[i].u[2]; 350 array[i].u[0] = array[i].u[1]; 351 array[i].u[2] = array[i].u[3]; 352 array[i].u[1] = x; 353 array[i].u[3] = y; 354 } 355} 356#endif 357/** 358 * This function represents a function used in the initialization 359 * by init_by_array 360 * @param x 32-bit integer 361 * @return 32-bit integer 362 */ 363static uint32_t func1(uint32_t x) { 364 return (x ^ (x >> 27)) * (uint32_t)1664525UL; 365} 366 367/** 368 * This function represents a function used in the initialization 369 * by init_by_array 370 * @param x 32-bit integer 371 * @return 32-bit integer 372 */ 373static uint32_t func2(uint32_t x) { 374 return (x ^ (x >> 27)) * (uint32_t)1566083941UL; 375} 376 377/** 378 * This function certificate the period of 2^{MEXP} 379 */ 380static void period_certification(sfmt_t *ctx) { 381 int inner = 0; 382 int i, j; 383 uint32_t work; 384 uint32_t *psfmt32 = &ctx->sfmt[0].u[0]; 385 386 for (i = 0; i < 4; i++) 387 inner ^= psfmt32[idxof(i)] & parity[i]; 388 for (i = 16; i > 0; i >>= 1) 389 inner ^= inner >> i; 390 inner &= 1; 391 /* check OK */ 392 if (inner == 1) { 393 return; 394 } 395 /* check NG, and modification */ 396 for (i = 0; i < 4; i++) { 397 work = 1; 398 for (j = 0; j < 32; j++) { 399 if ((work & parity[i]) != 0) { 400 psfmt32[idxof(i)] ^= work; 401 return; 402 } 403 work = work << 1; 404 } 405 } 406} 407 408/*---------------- 409 PUBLIC FUNCTIONS 410 ----------------*/ 411/** 412 * This function returns the identification string. 413 * The string shows the word size, the Mersenne exponent, 414 * and all parameters of this generator. 415 */ 416const char *get_idstring(void) { 417 return IDSTR; 418} 419 420/** 421 * This function returns the minimum size of array used for \b 422 * fill_array32() function. 423 * @return minimum size of array used for fill_array32() function. 424 */ 425int get_min_array_size32(void) { 426 return N32; 427} 428 429/** 430 * This function returns the minimum size of array used for \b 431 * fill_array64() function. 432 * @return minimum size of array used for fill_array64() function. 433 */ 434int get_min_array_size64(void) { 435 return N64; 436} 437 438#ifndef ONLY64 439/** 440 * This function generates and returns 32-bit pseudorandom number. 441 * init_gen_rand or init_by_array must be called before this function. 442 * @return 32-bit pseudorandom number 443 */ 444uint32_t gen_rand32(sfmt_t *ctx) { 445 uint32_t r; 446 uint32_t *psfmt32 = &ctx->sfmt[0].u[0]; 447 448 assert(ctx->initialized); 449 if (ctx->idx >= N32) { 450 gen_rand_all(ctx); 451 ctx->idx = 0; 452 } 453 r = psfmt32[ctx->idx++]; 454 return r; 455} 456 457/* Generate a random integer in [0..limit). */ 458uint32_t gen_rand32_range(sfmt_t *ctx, uint32_t limit) { 459 uint32_t ret, above; 460 461 above = 0xffffffffU - (0xffffffffU % limit); 462 while (1) { 463 ret = gen_rand32(ctx); 464 if (ret < above) { 465 ret %= limit; 466 break; 467 } 468 } 469 return ret; 470} 471#endif 472/** 473 * This function generates and returns 64-bit pseudorandom number. 474 * init_gen_rand or init_by_array must be called before this function. 475 * The function gen_rand64 should not be called after gen_rand32, 476 * unless an initialization is again executed. 477 * @return 64-bit pseudorandom number 478 */ 479uint64_t gen_rand64(sfmt_t *ctx) { 480#if defined(BIG_ENDIAN64) && !defined(ONLY64) 481 uint32_t r1, r2; 482 uint32_t *psfmt32 = &ctx->sfmt[0].u[0]; 483#else 484 uint64_t r; 485 uint64_t *psfmt64 = (uint64_t *)&ctx->sfmt[0].u[0]; 486#endif 487 488 assert(ctx->initialized); 489 assert(ctx->idx % 2 == 0); 490 491 if (ctx->idx >= N32) { 492 gen_rand_all(ctx); 493 ctx->idx = 0; 494 } 495#if defined(BIG_ENDIAN64) && !defined(ONLY64) 496 r1 = psfmt32[ctx->idx]; 497 r2 = psfmt32[ctx->idx + 1]; 498 ctx->idx += 2; 499 return ((uint64_t)r2 << 32) | r1; 500#else 501 r = psfmt64[ctx->idx / 2]; 502 ctx->idx += 2; 503 return r; 504#endif 505} 506 507/* Generate a random integer in [0..limit). */ 508uint64_t gen_rand64_range(sfmt_t *ctx, uint64_t limit) { 509 uint64_t ret, above; 510 511 above = 0xffffffffffffffffLLU - (0xffffffffffffffffLLU % limit); 512 while (1) { 513 ret = gen_rand64(ctx); 514 if (ret < above) { 515 ret %= limit; 516 break; 517 } 518 } 519 return ret; 520} 521 522#ifndef ONLY64 523/** 524 * This function generates pseudorandom 32-bit integers in the 525 * specified array[] by one call. The number of pseudorandom integers 526 * is specified by the argument size, which must be at least 624 and a 527 * multiple of four. The generation by this function is much faster 528 * than the following gen_rand function. 529 * 530 * For initialization, init_gen_rand or init_by_array must be called 531 * before the first call of this function. This function can not be 532 * used after calling gen_rand function, without initialization. 533 * 534 * @param array an array where pseudorandom 32-bit integers are filled 535 * by this function. The pointer to the array must be \b "aligned" 536 * (namely, must be a multiple of 16) in the SIMD version, since it 537 * refers to the address of a 128-bit integer. In the standard C 538 * version, the pointer is arbitrary. 539 * 540 * @param size the number of 32-bit pseudorandom integers to be 541 * generated. size must be a multiple of 4, and greater than or equal 542 * to (MEXP / 128 + 1) * 4. 543 * 544 * @note \b memalign or \b posix_memalign is available to get aligned 545 * memory. Mac OSX doesn't have these functions, but \b malloc of OSX 546 * returns the pointer to the aligned memory block. 547 */ 548void fill_array32(sfmt_t *ctx, uint32_t *array, int size) { 549 assert(ctx->initialized); 550 assert(ctx->idx == N32); 551 assert(size % 4 == 0); 552 assert(size >= N32); 553 554 gen_rand_array(ctx, (w128_t *)array, size / 4); 555 ctx->idx = N32; 556} 557#endif 558 559/** 560 * This function generates pseudorandom 64-bit integers in the 561 * specified array[] by one call. The number of pseudorandom integers 562 * is specified by the argument size, which must be at least 312 and a 563 * multiple of two. The generation by this function is much faster 564 * than the following gen_rand function. 565 * 566 * For initialization, init_gen_rand or init_by_array must be called 567 * before the first call of this function. This function can not be 568 * used after calling gen_rand function, without initialization. 569 * 570 * @param array an array where pseudorandom 64-bit integers are filled 571 * by this function. The pointer to the array must be "aligned" 572 * (namely, must be a multiple of 16) in the SIMD version, since it 573 * refers to the address of a 128-bit integer. In the standard C 574 * version, the pointer is arbitrary. 575 * 576 * @param size the number of 64-bit pseudorandom integers to be 577 * generated. size must be a multiple of 2, and greater than or equal 578 * to (MEXP / 128 + 1) * 2 579 * 580 * @note \b memalign or \b posix_memalign is available to get aligned 581 * memory. Mac OSX doesn't have these functions, but \b malloc of OSX 582 * returns the pointer to the aligned memory block. 583 */ 584void fill_array64(sfmt_t *ctx, uint64_t *array, int size) { 585 assert(ctx->initialized); 586 assert(ctx->idx == N32); 587 assert(size % 2 == 0); 588 assert(size >= N64); 589 590 gen_rand_array(ctx, (w128_t *)array, size / 2); 591 ctx->idx = N32; 592 593#if defined(BIG_ENDIAN64) && !defined(ONLY64) 594 swap((w128_t *)array, size /2); 595#endif 596} 597 598/** 599 * This function initializes the internal state array with a 32-bit 600 * integer seed. 601 * 602 * @param seed a 32-bit integer used as the seed. 603 */ 604sfmt_t *init_gen_rand(uint32_t seed) { 605 sfmt_t *ctx; 606 int i; 607 uint32_t *psfmt32; 608 609 if (posix_memalign((void **)&ctx, sizeof(w128_t), sizeof(sfmt_t)) != 0) { 610 return NULL; 611 } 612 psfmt32 = &ctx->sfmt[0].u[0]; 613 614 psfmt32[idxof(0)] = seed; 615 for (i = 1; i < N32; i++) { 616 psfmt32[idxof(i)] = 1812433253UL * (psfmt32[idxof(i - 1)] 617 ^ (psfmt32[idxof(i - 1)] >> 30)) 618 + i; 619 } 620 ctx->idx = N32; 621 period_certification(ctx); 622 ctx->initialized = 1; 623 624 return ctx; 625} 626 627/** 628 * This function initializes the internal state array, 629 * with an array of 32-bit integers used as the seeds 630 * @param init_key the array of 32-bit integers, used as a seed. 631 * @param key_length the length of init_key. 632 */ 633sfmt_t *init_by_array(uint32_t *init_key, int key_length) { 634 sfmt_t *ctx; 635 int i, j, count; 636 uint32_t r; 637 int lag; 638 int mid; 639 int size = N * 4; 640 uint32_t *psfmt32; 641 642 if (posix_memalign((void **)&ctx, sizeof(w128_t), sizeof(sfmt_t)) != 0) { 643 return NULL; 644 } 645 psfmt32 = &ctx->sfmt[0].u[0]; 646 647 if (size >= 623) { 648 lag = 11; 649 } else if (size >= 68) { 650 lag = 7; 651 } else if (size >= 39) { 652 lag = 5; 653 } else { 654 lag = 3; 655 } 656 mid = (size - lag) / 2; 657 658 memset(ctx->sfmt, 0x8b, sizeof(ctx->sfmt)); 659 if (key_length + 1 > N32) { 660 count = key_length + 1; 661 } else { 662 count = N32; 663 } 664 r = func1(psfmt32[idxof(0)] ^ psfmt32[idxof(mid)] 665 ^ psfmt32[idxof(N32 - 1)]); 666 psfmt32[idxof(mid)] += r; 667 r += key_length; 668 psfmt32[idxof(mid + lag)] += r; 669 psfmt32[idxof(0)] = r; 670 671 count--; 672 for (i = 1, j = 0; (j < count) && (j < key_length); j++) { 673 r = func1(psfmt32[idxof(i)] ^ psfmt32[idxof((i + mid) % N32)] 674 ^ psfmt32[idxof((i + N32 - 1) % N32)]); 675 psfmt32[idxof((i + mid) % N32)] += r; 676 r += init_key[j] + i; 677 psfmt32[idxof((i + mid + lag) % N32)] += r; 678 psfmt32[idxof(i)] = r; 679 i = (i + 1) % N32; 680 } 681 for (; j < count; j++) { 682 r = func1(psfmt32[idxof(i)] ^ psfmt32[idxof((i + mid) % N32)] 683 ^ psfmt32[idxof((i + N32 - 1) % N32)]); 684 psfmt32[idxof((i + mid) % N32)] += r; 685 r += i; 686 psfmt32[idxof((i + mid + lag) % N32)] += r; 687 psfmt32[idxof(i)] = r; 688 i = (i + 1) % N32; 689 } 690 for (j = 0; j < N32; j++) { 691 r = func2(psfmt32[idxof(i)] + psfmt32[idxof((i + mid) % N32)] 692 + psfmt32[idxof((i + N32 - 1) % N32)]); 693 psfmt32[idxof((i + mid) % N32)] ^= r; 694 r -= i; 695 psfmt32[idxof((i + mid + lag) % N32)] ^= r; 696 psfmt32[idxof(i)] = r; 697 i = (i + 1) % N32; 698 } 699 700 ctx->idx = N32; 701 period_certification(ctx); 702 ctx->initialized = 1; 703 704 return ctx; 705} 706 707void fini_gen_rand(sfmt_t *ctx) { 708 assert(ctx != NULL); 709 710 ctx->initialized = 0; 711 free(ctx); 712} 713