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