1 /**
2 * @file SFMT.c
3 * @brief SIMD oriented Fast Mersenne Twister(SFMT)
4 *
5 * @author Mutsuo Saito (Hiroshima University)
6 * @author Makoto Matsumoto (Hiroshima University)
7 *
8 * Copyright (C) 2006, 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima
9 * University.
10 * Copyright (C) 2012 Mutsuo Saito, Makoto Matsumoto, Hiroshima
11 * University and The University of Tokyo.
12 * All rights reserved.
13 *
14 * The 3-clause BSD License is applied to this software:
15 *
16 * * * * * *
17 * Redistribution and use in source and binary forms, with or without
18 * modification, are permitted provided that the following conditions are met:
19 *
20 * - Redistributions of source code must retain the above copyright notice,
21 * this list of conditions and the following disclaimer.
22 * - Redistributions in binary form must reproduce the above copyright notice,
23 * this list of conditions and the following disclaimer in the documentation
24 * and/or other materials provided with the distribution.
25 * - Neither the names of Hiroshima University, The University of Tokyo nor the
26 * names of its contributors may be used to endorse or promote products
27 * derived from this software without specific prior written permission.
28 *
29 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
30 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
31 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
32 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
33 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
34 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
35 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
36 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
37 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
38 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
39 * POSSIBILITY OF SUCH DAMAGE.
40 * * * * * *
41 */
42
43 #if defined(__cplusplus)
44 extern "C" {
45 #endif
46
47 #include <string.h>
48 #include <assert.h>
49 #include "SFMT.h"
50
51 #ifndef __LP64__
52 inline static void do_recursion(w128_t * r, w128_t * a, w128_t * b,
53 w128_t * c, w128_t * d);
54 #endif
55
56 #ifndef __LP64__
57 inline static void rshift128(w128_t *out, w128_t const *in, int shift);
58 inline static void lshift128(w128_t *out, w128_t const *in, int shift);
59
60 /**
61 * This function simulates SIMD 128-bit right shift by the standard C.
62 * The 128-bit integer given in in is shifted by (shift * 8) bits.
63 * This function simulates the LITTLE ENDIAN SIMD.
64 * @param out the output of this function
65 * @param in the 128-bit data to be shifted
66 * @param shift the shift value
67 */
rshift128(w128_t * out,w128_t const * in,int shift)68 inline static void rshift128(w128_t *out, w128_t const *in, int shift)
69 {
70 uint64_t th, tl, oh, ol;
71
72 th = ((uint64_t)in->u[3] << 32) | ((uint64_t)in->u[2]);
73 tl = ((uint64_t)in->u[1] << 32) | ((uint64_t)in->u[0]);
74
75 oh = th >> (shift * 8);
76 ol = tl >> (shift * 8);
77 ol |= th << (64 - shift * 8);
78 out->u[1] = (uint32_t)(ol >> 32);
79 out->u[0] = (uint32_t)ol;
80 out->u[3] = (uint32_t)(oh >> 32);
81 out->u[2] = (uint32_t)oh;
82 }
83
84 /**
85 * This function simulates SIMD 128-bit left shift by the standard C.
86 * The 128-bit integer given in in is shifted by (shift * 8) bits.
87 * This function simulates the LITTLE ENDIAN SIMD.
88 * @param out the output of this function
89 * @param in the 128-bit data to be shifted
90 * @param shift the shift value
91 */
lshift128(w128_t * out,w128_t const * in,int shift)92 inline static void lshift128(w128_t *out, w128_t const *in, int shift)
93 {
94 uint64_t th, tl, oh, ol;
95
96 th = ((uint64_t)in->u[3] << 32) | ((uint64_t)in->u[2]);
97 tl = ((uint64_t)in->u[1] << 32) | ((uint64_t)in->u[0]);
98
99 oh = th << (shift * 8);
100 ol = tl << (shift * 8);
101 oh |= tl >> (64 - shift * 8);
102 out->u[1] = (uint32_t)(ol >> 32);
103 out->u[0] = (uint32_t)ol;
104 out->u[3] = (uint32_t)(oh >> 32);
105 out->u[2] = (uint32_t)oh;
106 }
107
108 /**
109 * This function represents the recursion formula.
110 * @param r output
111 * @param a a 128-bit part of the internal state array
112 * @param b a 128-bit part of the internal state array
113 * @param c a 128-bit part of the internal state array
114 * @param d a 128-bit part of the internal state array
115 */
do_recursion(w128_t * r,w128_t * a,w128_t * b,w128_t * c,w128_t * d)116 inline static void do_recursion(w128_t *r, w128_t *a, w128_t *b,
117 w128_t *c, w128_t *d)
118 {
119 w128_t x;
120 w128_t y;
121
122 lshift128(&x, a, SFMT_SL2);
123 rshift128(&y, c, SFMT_SR2);
124 r->u[0] = a->u[0] ^ x.u[0] ^ ((b->u[0] >> SFMT_SR1) & SFMT_MSK1)
125 ^ y.u[0] ^ (d->u[0] << SFMT_SL1);
126 r->u[1] = a->u[1] ^ x.u[1] ^ ((b->u[1] >> SFMT_SR1) & SFMT_MSK2)
127 ^ y.u[1] ^ (d->u[1] << SFMT_SL1);
128 r->u[2] = a->u[2] ^ x.u[2] ^ ((b->u[2] >> SFMT_SR1) & SFMT_MSK3)
129 ^ y.u[2] ^ (d->u[2] << SFMT_SL1);
130 r->u[3] = a->u[3] ^ x.u[3] ^ ((b->u[3] >> SFMT_SR1) & SFMT_MSK4)
131 ^ y.u[3] ^ (d->u[3] << SFMT_SL1);
132 }
133 #endif
134
135 /**
136 * parameters used by sse2.
137 */
138 #ifdef __LP64__
139 static const w128_t sse2_param_mask = {{SFMT_MSK1, SFMT_MSK2,
140 SFMT_MSK3, SFMT_MSK4}};
141 #endif
142 /*----------------
143 STATIC FUNCTIONS
144 ----------------*/
145 inline static int idxof(int i);
146 inline static void gen_rand_array(sfmt_t * sfmt, w128_t *array, int size);
147 inline static uint32_t func1(uint32_t x);
148 inline static uint32_t func2(uint32_t x);
149 static void period_certification(sfmt_t * sfmt);
150
151 #ifdef __LP64__
152 inline static void mm_recursion(__m128i * r, __m128i a, __m128i b,
153 __m128i c, __m128i d);
154
155 /**
156 * This function represents the recursion formula.
157 * @param r an output
158 * @param a a 128-bit part of the interal state array
159 * @param b a 128-bit part of the interal state array
160 * @param c a 128-bit part of the interal state array
161 * @param d a 128-bit part of the interal state array
162 */
mm_recursion(__m128i * r,__m128i a,__m128i b,__m128i c,__m128i d)163 inline static void mm_recursion(__m128i * r, __m128i a, __m128i b,
164 __m128i c, __m128i d)
165 {
166 __m128i v, x, y, z;
167
168 y = _mm_srli_epi32(b, SFMT_SR1);
169 z = _mm_srli_si128(c, SFMT_SR2);
170 v = _mm_slli_epi32(d, SFMT_SL1);
171 z = _mm_xor_si128(z, a);
172 z = _mm_xor_si128(z, v);
173 x = _mm_slli_si128(a, SFMT_SL2);
174 y = _mm_and_si128(y, sse2_param_mask.si);
175 z = _mm_xor_si128(z, x);
176 z = _mm_xor_si128(z, y);
177 *r = z;
178 }
179
180 /**
181 * This function fills the internal state array with pseudorandom
182 * integers.
183 * @param sfmt SFMT internal state
184 */
sfmt_gen_rand_all(sfmt_t * sfmt)185 void sfmt_gen_rand_all(sfmt_t * sfmt) {
186 int i;
187 __m128i r1, r2;
188 w128_t * pstate = sfmt->state;
189
190 r1 = pstate[SFMT_N - 2].si;
191 r2 = pstate[SFMT_N - 1].si;
192 for (i = 0; i < SFMT_N - SFMT_POS1; i++) {
193 mm_recursion(&pstate[i].si, pstate[i].si,
194 pstate[i + SFMT_POS1].si, r1, r2);
195 r1 = r2;
196 r2 = pstate[i].si;
197 }
198 for (; i < SFMT_N; i++) {
199 mm_recursion(&pstate[i].si, pstate[i].si,
200 pstate[i + SFMT_POS1 - SFMT_N].si,
201 r1, r2);
202 r1 = r2;
203 r2 = pstate[i].si;
204 }
205 }
206
207 /**
208 * This function fills the user-specified array with pseudorandom
209 * integers.
210 * @param sfmt SFMT internal state.
211 * @param array an 128-bit array to be filled by pseudorandom numbers.
212 * @param size number of 128-bit pseudorandom numbers to be generated.
213 */
gen_rand_array(sfmt_t * sfmt,w128_t * array,int size)214 static void gen_rand_array(sfmt_t * sfmt, w128_t * array, int size)
215 {
216 int i, j;
217 __m128i r1, r2;
218 w128_t * pstate = sfmt->state;
219
220 r1 = pstate[SFMT_N - 2].si;
221 r2 = pstate[SFMT_N - 1].si;
222 for (i = 0; i < SFMT_N - SFMT_POS1; i++) {
223 mm_recursion(&array[i].si, pstate[i].si,
224 pstate[i + SFMT_POS1].si, r1, r2);
225 r1 = r2;
226 r2 = array[i].si;
227 }
228 for (; i < SFMT_N; i++) {
229 mm_recursion(&array[i].si, pstate[i].si,
230 array[i + SFMT_POS1 - SFMT_N].si, r1, r2);
231 r1 = r2;
232 r2 = array[i].si;
233 }
234 for (; i < size - SFMT_N; i++) {
235 mm_recursion(&array[i].si, array[i - SFMT_N].si,
236 array[i + SFMT_POS1 - SFMT_N].si, r1, r2);
237 r1 = r2;
238 r2 = array[i].si;
239 }
240 for (j = 0; j < 2 * SFMT_N - size; j++) {
241 pstate[j] = array[j + size - SFMT_N];
242 }
243 for (; i < size; i++, j++) {
244 mm_recursion(&array[i].si, array[i - SFMT_N].si,
245 array[i + SFMT_POS1 - SFMT_N].si, r1, r2);
246 r1 = r2;
247 r2 = array[i].si;
248 pstate[j] = array[i];
249 }
250 }
251
252 #endif
253
254 /**
255 * This function simulate a 64-bit index of LITTLE ENDIAN
256 * in BIG ENDIAN machine.
257 */
idxof(int i)258 inline static int idxof(int i) {
259 return i;
260 }
261
262 #ifndef __LP64__
263 /**
264 * This function fills the user-specified array with pseudorandom
265 * integers.
266 *
267 * @param sfmt SFMT internal state
268 * @param array an 128-bit array to be filled by pseudorandom numbers.
269 * @param size number of 128-bit pseudorandom numbers to be generated.
270 */
gen_rand_array(sfmt_t * sfmt,w128_t * array,int size)271 inline static void gen_rand_array(sfmt_t * sfmt, w128_t *array, int size) {
272 int i, j;
273 w128_t *r1, *r2;
274
275 r1 = &sfmt->state[SFMT_N - 2];
276 r2 = &sfmt->state[SFMT_N - 1];
277 for (i = 0; i < SFMT_N - SFMT_POS1; i++) {
278 do_recursion(&array[i], &sfmt->state[i], &sfmt->state[i + SFMT_POS1], r1, r2);
279 r1 = r2;
280 r2 = &array[i];
281 }
282 for (; i < SFMT_N; i++) {
283 do_recursion(&array[i], &sfmt->state[i],
284 &array[i + SFMT_POS1 - SFMT_N], r1, r2);
285 r1 = r2;
286 r2 = &array[i];
287 }
288 for (; i < size - SFMT_N; i++) {
289 do_recursion(&array[i], &array[i - SFMT_N],
290 &array[i + SFMT_POS1 - SFMT_N], r1, r2);
291 r1 = r2;
292 r2 = &array[i];
293 }
294 for (j = 0; j < 2 * SFMT_N - size; j++) {
295 sfmt->state[j] = array[j + size - SFMT_N];
296 }
297 for (; i < size; i++, j++) {
298 do_recursion(&array[i], &array[i - SFMT_N],
299 &array[i + SFMT_POS1 - SFMT_N], r1, r2);
300 r1 = r2;
301 r2 = &array[i];
302 sfmt->state[j] = array[i];
303 }
304 }
305 #endif
306
307 /**
308 * This function represents a function used in the initialization
309 * by init_by_array
310 * @param x 32-bit integer
311 * @return 32-bit integer
312 */
func1(uint32_t x)313 static uint32_t func1(uint32_t x) {
314 return (x ^ (x >> 27)) * (uint32_t)1664525UL;
315 }
316
317 /**
318 * This function represents a function used in the initialization
319 * by init_by_array
320 * @param x 32-bit integer
321 * @return 32-bit integer
322 */
func2(uint32_t x)323 static uint32_t func2(uint32_t x) {
324 return (x ^ (x >> 27)) * (uint32_t)1566083941UL;
325 }
326
327 /**
328 * This function certificate the period of 2^{MEXP}
329 * @param sfmt SFMT internal state
330 */
period_certification(sfmt_t * sfmt)331 static void period_certification(sfmt_t * sfmt) {
332 int inner = 0;
333 int i, j;
334 uint32_t work;
335 uint32_t *psfmt32 = &sfmt->state[0].u[0];
336 const uint32_t parity[4] = {SFMT_PARITY1, SFMT_PARITY2,
337 SFMT_PARITY3, SFMT_PARITY4};
338
339 for (i = 0; i < 4; i++)
340 inner ^= psfmt32[idxof(i)] & parity[i];
341 for (i = 16; i > 0; i >>= 1)
342 inner ^= inner >> i;
343 inner &= 1;
344 /* check OK */
345 if (inner == 1) {
346 return;
347 }
348 /* check NG, and modification */
349 for (i = 0; i < 4; i++) {
350 work = 1;
351 for (j = 0; j < 32; j++) {
352 if ((work & parity[i]) != 0) {
353 psfmt32[idxof(i)] ^= work;
354 return;
355 }
356 work = work << 1;
357 }
358 }
359 }
360
361 /*----------------
362 PUBLIC FUNCTIONS
363 ----------------*/
364 #define UNUSED_VARIABLE(x) (void)(x)
365 /**
366 * This function returns the identification string.
367 * The string shows the word size, the Mersenne exponent,
368 * and all parameters of this generator.
369 * @param sfmt SFMT internal state
370 */
sfmt_get_idstring(sfmt_t * sfmt)371 const char *sfmt_get_idstring(sfmt_t * sfmt) {
372 UNUSED_VARIABLE(sfmt);
373 return SFMT_IDSTR;
374 }
375
376 /**
377 * This function returns the minimum size of array used for \b
378 * fill_array32() function.
379 * @param sfmt SFMT internal state
380 * @return minimum size of array used for fill_array32() function.
381 */
sfmt_get_min_array_size32(sfmt_t * sfmt)382 int sfmt_get_min_array_size32(sfmt_t * sfmt) {
383 UNUSED_VARIABLE(sfmt);
384 return SFMT_N32;
385 }
386
387 /**
388 * This function returns the minimum size of array used for \b
389 * fill_array64() function.
390 * @param sfmt SFMT internal state
391 * @return minimum size of array used for fill_array64() function.
392 */
sfmt_get_min_array_size64(sfmt_t * sfmt)393 int sfmt_get_min_array_size64(sfmt_t * sfmt) {
394 UNUSED_VARIABLE(sfmt);
395 return SFMT_N64;
396 }
397
398 #ifndef __LP64__
399 /**
400 * This function fills the internal state array with pseudorandom
401 * integers.
402 * @param sfmt SFMT internal state
403 */
sfmt_gen_rand_all(sfmt_t * sfmt)404 void sfmt_gen_rand_all(sfmt_t * sfmt) {
405 int i;
406 w128_t *r1, *r2;
407
408 r1 = &sfmt->state[SFMT_N - 2];
409 r2 = &sfmt->state[SFMT_N - 1];
410 for (i = 0; i < SFMT_N - SFMT_POS1; i++) {
411 do_recursion(&sfmt->state[i], &sfmt->state[i],
412 &sfmt->state[i + SFMT_POS1], r1, r2);
413 r1 = r2;
414 r2 = &sfmt->state[i];
415 }
416 for (; i < SFMT_N; i++) {
417 do_recursion(&sfmt->state[i], &sfmt->state[i],
418 &sfmt->state[i + SFMT_POS1 - SFMT_N], r1, r2);
419 r1 = r2;
420 r2 = &sfmt->state[i];
421 }
422 }
423 #endif
424
425 /**
426 * This function generates pseudorandom 32-bit integers in the
427 * specified array[] by one call. The number of pseudorandom integers
428 * is specified by the argument size, which must be at least 624 and a
429 * multiple of four. The generation by this function is much faster
430 * than the following gen_rand function.
431 *
432 * For initialization, init_gen_rand or init_by_array must be called
433 * before the first call of this function. This function can not be
434 * used after calling gen_rand function, without initialization.
435 *
436 * @param sfmt SFMT internal state
437 * @param array an array where pseudorandom 32-bit integers are filled
438 * by this function. The pointer to the array must be \b "aligned"
439 * (namely, must be a multiple of 16) in the SIMD version, since it
440 * refers to the address of a 128-bit integer. In the standard C
441 * version, the pointer is arbitrary.
442 *
443 * @param size the number of 32-bit pseudorandom integers to be
444 * generated. size must be a multiple of 4, and greater than or equal
445 * to (MEXP / 128 + 1) * 4.
446 *
447 * @note \b memalign or \b posix_memalign is available to get aligned
448 * memory. Mac OSX doesn't have these functions, but \b malloc of OSX
449 * returns the pointer to the aligned memory block.
450 */
sfmt_fill_array32(sfmt_t * sfmt,uint32_t * array,int size)451 void sfmt_fill_array32(sfmt_t * sfmt, uint32_t *array, int size) {
452 assert(sfmt->idx == SFMT_N32);
453 assert(size % 4 == 0);
454 assert(size >= SFMT_N32);
455
456 gen_rand_array(sfmt, (w128_t *)array, size / 4);
457 sfmt->idx = SFMT_N32;
458 }
459
460 /**
461 * This function generates pseudorandom 64-bit integers in the
462 * specified array[] by one call. The number of pseudorandom integers
463 * is specified by the argument size, which must be at least 312 and a
464 * multiple of two. The generation by this function is much faster
465 * than the following gen_rand function.
466 *
467 * @param sfmt SFMT internal state
468 * For initialization, init_gen_rand or init_by_array must be called
469 * before the first call of this function. This function can not be
470 * used after calling gen_rand function, without initialization.
471 *
472 * @param array an array where pseudorandom 64-bit integers are filled
473 * by this function. The pointer to the array must be "aligned"
474 * (namely, must be a multiple of 16) in the SIMD version, since it
475 * refers to the address of a 128-bit integer. In the standard C
476 * version, the pointer is arbitrary.
477 *
478 * @param size the number of 64-bit pseudorandom integers to be
479 * generated. size must be a multiple of 2, and greater than or equal
480 * to (MEXP / 128 + 1) * 2
481 *
482 * @note \b memalign or \b posix_memalign is available to get aligned
483 * memory. Mac OSX doesn't have these functions, but \b malloc of OSX
484 * returns the pointer to the aligned memory block.
485 */
sfmt_fill_array64(sfmt_t * sfmt,uint64_t * array,int size)486 void sfmt_fill_array64(sfmt_t * sfmt, uint64_t *array, int size) {
487 assert(sfmt->idx == SFMT_N32);
488 assert(size % 2 == 0);
489 assert(size >= SFMT_N64);
490
491 gen_rand_array(sfmt, (w128_t *)array, size / 2);
492 sfmt->idx = SFMT_N32;
493
494 }
495
496 /**
497 * This function initializes the internal state array with a 32-bit
498 * integer seed.
499 *
500 * @param sfmt SFMT internal state
501 * @param seed a 32-bit integer used as the seed.
502 */
sfmt_init_gen_rand(sfmt_t * sfmt,uint32_t seed)503 void sfmt_init_gen_rand(sfmt_t * sfmt, uint32_t seed) {
504 int i;
505
506 uint32_t *psfmt32 = &sfmt->state[0].u[0];
507
508 psfmt32[idxof(0)] = seed;
509 for (i = 1; i < SFMT_N32; i++) {
510 psfmt32[idxof(i)] = 1812433253UL * (psfmt32[idxof(i - 1)]
511 ^ (psfmt32[idxof(i - 1)] >> 30))
512 + i;
513 }
514 sfmt->idx = SFMT_N32;
515 period_certification(sfmt);
516 }
517
518 /**
519 * This function initializes the internal state array,
520 * with an array of 32-bit integers used as the seeds
521 * @param sfmt SFMT internal state
522 * @param init_key the array of 32-bit integers, used as a seed.
523 * @param key_length the length of init_key.
524 */
sfmt_init_by_array(sfmt_t * sfmt,uint32_t * init_key,int key_length)525 void sfmt_init_by_array(sfmt_t * sfmt, uint32_t *init_key, int key_length) {
526 int i, j, count;
527 uint32_t r;
528 int lag;
529 int mid;
530 int size = SFMT_N * 4;
531 uint32_t *psfmt32 = &sfmt->state[0].u[0];
532
533 if (size >= 623) {
534 lag = 11;
535 } else if (size >= 68) {
536 lag = 7;
537 } else if (size >= 39) {
538 lag = 5;
539 } else {
540 lag = 3;
541 }
542 mid = (size - lag) / 2;
543
544 memset(sfmt, 0x8b, sizeof(sfmt_t));
545 if (key_length + 1 > SFMT_N32) {
546 count = key_length + 1;
547 } else {
548 count = SFMT_N32;
549 }
550 r = func1(psfmt32[idxof(0)] ^ psfmt32[idxof(mid)]
551 ^ psfmt32[idxof(SFMT_N32 - 1)]);
552 psfmt32[idxof(mid)] += r;
553 r += key_length;
554 psfmt32[idxof(mid + lag)] += r;
555 psfmt32[idxof(0)] = r;
556
557 count--;
558 for (i = 1, j = 0; (j < count) && (j < key_length); j++) {
559 r = func1(psfmt32[idxof(i)] ^ psfmt32[idxof((i + mid) % SFMT_N32)]
560 ^ psfmt32[idxof((i + SFMT_N32 - 1) % SFMT_N32)]);
561 psfmt32[idxof((i + mid) % SFMT_N32)] += r;
562 r += init_key[j] + i;
563 psfmt32[idxof((i + mid + lag) % SFMT_N32)] += r;
564 psfmt32[idxof(i)] = r;
565 i = (i + 1) % SFMT_N32;
566 }
567 for (; j < count; j++) {
568 r = func1(psfmt32[idxof(i)] ^ psfmt32[idxof((i + mid) % SFMT_N32)]
569 ^ psfmt32[idxof((i + SFMT_N32 - 1) % SFMT_N32)]);
570 psfmt32[idxof((i + mid) % SFMT_N32)] += r;
571 r += i;
572 psfmt32[idxof((i + mid + lag) % SFMT_N32)] += r;
573 psfmt32[idxof(i)] = r;
574 i = (i + 1) % SFMT_N32;
575 }
576 for (j = 0; j < SFMT_N32; j++) {
577 r = func2(psfmt32[idxof(i)] + psfmt32[idxof((i + mid) % SFMT_N32)]
578 + psfmt32[idxof((i + SFMT_N32 - 1) % SFMT_N32)]);
579 psfmt32[idxof((i + mid) % SFMT_N32)] ^= r;
580 r -= i;
581 psfmt32[idxof((i + mid + lag) % SFMT_N32)] ^= r;
582 psfmt32[idxof(i)] = r;
583 i = (i + 1) % SFMT_N32;
584 }
585
586 sfmt->idx = SFMT_N32;
587 period_certification(sfmt);
588 }
589 #if defined(__cplusplus)
590 }
591 #endif
592