Coverage Report

Created: 2019-07-03 22:50

/home/liu/buildslave/linux-x64-runtests/build/lib/3rdparty/dsfmt/dSFMT.c
Line
Count
Source (jump to first uncovered line)
1
/**
2
 * @file dSFMT.c
3
 * @brief double precision SIMD-oriented Fast Mersenne Twister (dSFMT)
4
 * based on IEEE 754 format.
5
 *
6
 * @author Mutsuo Saito (Hiroshima University)
7
 * @author Makoto Matsumoto (Hiroshima University)
8
 *
9
 * Copyright (C) 2007,2008 Mutsuo Saito, Makoto Matsumoto and Hiroshima
10
 * University. All rights reserved.
11
 *
12
 * The new BSD License is applied to this software, see LICENSE.txt
13
 */
14
#include <stdio.h>
15
#include <string.h>
16
#include <stdlib.h>
17
#include "dSFMT-params.h"
18
#include "dSFMT-common.h"
19
20
#if defined(__cplusplus)
21
extern "C" {
22
#endif
23
24
/** dsfmt internal state vector */
25
dsfmt_t dsfmt_global_data;
26
/** dsfmt mexp for check */
27
static const int dsfmt_mexp = DSFMT_MEXP;
28
29
/*----------------
30
  STATIC FUNCTIONS
31
  ----------------*/
32
inline static uint32_t ini_func1(uint32_t x);
33
inline static uint32_t ini_func2(uint32_t x);
34
inline static void gen_rand_array_c1o2(dsfmt_t *dsfmt, dw128_t *array,
35
               int size);
36
inline static void gen_rand_array_c0o1(dsfmt_t *dsfmt, dw128_t *array,
37
               int size);
38
inline static void gen_rand_array_o0c1(dsfmt_t *dsfmt, dw128_t *array,
39
               int size);
40
inline static void gen_rand_array_o0o1(dsfmt_t *dsfmt, dw128_t *array,
41
               int size);
42
inline static int idxof(int i);
43
static void initial_mask(dsfmt_t *dsfmt);
44
static void period_certification(dsfmt_t *dsfmt);
45
46
#if defined(HAVE_SSE2)
47
/** 1 in 64bit for sse2 */
48
static const union X128I_T sse2_int_one = {{1, 1}};
49
/** 2.0 double for sse2 */
50
static const union X128D_T sse2_double_two = {{2.0, 2.0}};
51
/** -1.0 double for sse2 */
52
static const union X128D_T sse2_double_m_one = {{-1.0, -1.0}};
53
#endif
54
55
/**
56
 * This function simulate a 32-bit array index overlapped to 64-bit
57
 * array of LITTLE ENDIAN in BIG ENDIAN machine.
58
 */
59
#if defined(DSFMT_BIG_ENDIAN)
60
inline static int idxof(int i) {
61
    return i ^ 1;
62
}
63
#else
64
255k
inline static int idxof(int i) {
65
255k
    return i;
66
255k
}
67
#endif
68
69
#if defined(HAVE_SSE2)
70
/**
71
 * This function converts the double precision floating point numbers which
72
 * distribute uniformly in the range [1, 2) to those which distribute uniformly
73
 * in the range [0, 1).
74
 * @param w 128bit stracture of double precision floating point numbers (I/O)
75
 */
76
0
inline static void convert_c0o1(dw128_t *w) {
77
0
    w->sd = _mm_add_pd(w->sd, sse2_double_m_one.d128);
78
0
}
79
80
/**
81
 * This function converts the double precision floating point numbers which
82
 * distribute uniformly in the range [1, 2) to those which distribute uniformly
83
 * in the range (0, 1].
84
 * @param w 128bit stracture of double precision floating point numbers (I/O)
85
 */
86
0
inline static void convert_o0c1(dw128_t *w) {
87
0
    w->sd = _mm_sub_pd(sse2_double_two.d128, w->sd);
88
0
}
89
90
/**
91
 * This function converts the double precision floating point numbers which
92
 * distribute uniformly in the range [1, 2) to those which distribute uniformly
93
 * in the range (0, 1).
94
 * @param w 128bit stracture of double precision floating point numbers (I/O)
95
 */
96
0
inline static void convert_o0o1(dw128_t *w) {
97
0
    w->si = _mm_or_si128(w->si, sse2_int_one.i128);
98
0
    w->sd = _mm_add_pd(w->sd, sse2_double_m_one.d128);
99
0
}
100
#else /* standard C and altivec */
101
/**
102
 * This function converts the double precision floating point numbers which
103
 * distribute uniformly in the range [1, 2) to those which distribute uniformly
104
 * in the range [0, 1).
105
 * @param w 128bit stracture of double precision floating point numbers (I/O)
106
 */
107
inline static void convert_c0o1(dw128_t *w) {
108
    w->d[0] -= 1.0;
109
    w->d[1] -= 1.0;
110
}
111
112
/**
113
 * This function converts the double precision floating point numbers which
114
 * distribute uniformly in the range [1, 2) to those which distribute uniformly
115
 * in the range (0, 1].
116
 * @param w 128bit stracture of double precision floating point numbers (I/O)
117
 */
118
inline static void convert_o0c1(dw128_t *w) {
119
    w->d[0] = 2.0 - w->d[0];
120
    w->d[1] = 2.0 - w->d[1];
121
}
122
123
/**
124
 * This function converts the double precision floating point numbers which
125
 * distribute uniformly in the range [1, 2) to those which distribute uniformly
126
 * in the range (0, 1).
127
 * @param w 128bit stracture of double precision floating point numbers (I/O)
128
 */
129
inline static void convert_o0o1(dw128_t *w) {
130
    w->u[0] |= 1;
131
    w->u[1] |= 1;
132
    w->d[0] -= 1.0;
133
    w->d[1] -= 1.0;
134
}
135
#endif
136
137
/**
138
 * This function fills the user-specified array with double precision
139
 * floating point pseudorandom numbers of the IEEE 754 format.
140
 * @param dsfmt dsfmt state vector.
141
 * @param array an 128-bit array to be filled by pseudorandom numbers.
142
 * @param size number of 128-bit pseudorandom numbers to be generated.
143
 */
144
inline static void gen_rand_array_c1o2(dsfmt_t *dsfmt, dw128_t *array,
145
0
               int size) {
146
0
    int i, j;
147
0
    dw128_t lung;
148
0
149
0
    lung = dsfmt->status[DSFMT_N];
150
0
    do_recursion(&array[0], &dsfmt->status[0], &dsfmt->status[DSFMT_POS1],
151
0
     &lung);
152
0
    for (i = 1; i < DSFMT_N - DSFMT_POS1; i++) {
153
0
  do_recursion(&array[i], &dsfmt->status[i],
154
0
         &dsfmt->status[i + DSFMT_POS1], &lung);
155
0
    }
156
0
    for (; i < DSFMT_N; i++) {
157
0
  do_recursion(&array[i], &dsfmt->status[i],
158
0
         &array[i + DSFMT_POS1 - DSFMT_N], &lung);
159
0
    }
160
0
    for (; i < size - DSFMT_N; i++) {
161
0
  do_recursion(&array[i], &array[i - DSFMT_N],
162
0
         &array[i + DSFMT_POS1 - DSFMT_N], &lung);
163
0
    }
164
0
    for (j = 0; j < 2 * DSFMT_N - size; j++) {
165
0
  dsfmt->status[j] = array[j + size - DSFMT_N];
166
0
    }
167
0
    for (; i < size; i++, j++) {
168
0
  do_recursion(&array[i], &array[i - DSFMT_N],
169
0
         &array[i + DSFMT_POS1 - DSFMT_N], &lung);
170
0
  dsfmt->status[j] = array[i];
171
0
    }
172
0
    dsfmt->status[DSFMT_N] = lung;
173
0
}
174
175
/**
176
 * This function fills the user-specified array with double precision
177
 * floating point pseudorandom numbers of the IEEE 754 format.
178
 * @param dsfmt dsfmt state vector.
179
 * @param array an 128-bit array to be filled by pseudorandom numbers.
180
 * @param size number of 128-bit pseudorandom numbers to be generated.
181
 */
182
inline static void gen_rand_array_c0o1(dsfmt_t *dsfmt, dw128_t *array,
183
0
               int size) {
184
0
    int i, j;
185
0
    dw128_t lung;
186
0
187
0
    lung = dsfmt->status[DSFMT_N];
188
0
    do_recursion(&array[0], &dsfmt->status[0], &dsfmt->status[DSFMT_POS1],
189
0
     &lung);
190
0
    for (i = 1; i < DSFMT_N - DSFMT_POS1; i++) {
191
0
  do_recursion(&array[i], &dsfmt->status[i],
192
0
         &dsfmt->status[i + DSFMT_POS1], &lung);
193
0
    }
194
0
    for (; i < DSFMT_N; i++) {
195
0
  do_recursion(&array[i], &dsfmt->status[i],
196
0
         &array[i + DSFMT_POS1 - DSFMT_N], &lung);
197
0
    }
198
0
    for (; i < size - DSFMT_N; i++) {
199
0
  do_recursion(&array[i], &array[i - DSFMT_N],
200
0
         &array[i + DSFMT_POS1 - DSFMT_N], &lung);
201
0
  convert_c0o1(&array[i - DSFMT_N]);
202
0
    }
203
0
    for (j = 0; j < 2 * DSFMT_N - size; j++) {
204
0
  dsfmt->status[j] = array[j + size - DSFMT_N];
205
0
    }
206
0
    for (; i < size; i++, j++) {
207
0
  do_recursion(&array[i], &array[i - DSFMT_N],
208
0
         &array[i + DSFMT_POS1 - DSFMT_N], &lung);
209
0
  dsfmt->status[j] = array[i];
210
0
  convert_c0o1(&array[i - DSFMT_N]);
211
0
    }
212
0
    for (i = size - DSFMT_N; i < size; i++) {
213
0
  convert_c0o1(&array[i]);
214
0
    }
215
0
    dsfmt->status[DSFMT_N] = lung;
216
0
}
217
218
/**
219
 * This function fills the user-specified array with double precision
220
 * floating point pseudorandom numbers of the IEEE 754 format.
221
 * @param dsfmt dsfmt state vector.
222
 * @param array an 128-bit array to be filled by pseudorandom numbers.
223
 * @param size number of 128-bit pseudorandom numbers to be generated.
224
 */
225
inline static void gen_rand_array_o0o1(dsfmt_t *dsfmt, dw128_t *array,
226
0
               int size) {
227
0
    int i, j;
228
0
    dw128_t lung;
229
0
230
0
    lung = dsfmt->status[DSFMT_N];
231
0
    do_recursion(&array[0], &dsfmt->status[0], &dsfmt->status[DSFMT_POS1],
232
0
     &lung);
233
0
    for (i = 1; i < DSFMT_N - DSFMT_POS1; i++) {
234
0
  do_recursion(&array[i], &dsfmt->status[i],
235
0
         &dsfmt->status[i + DSFMT_POS1], &lung);
236
0
    }
237
0
    for (; i < DSFMT_N; i++) {
238
0
  do_recursion(&array[i], &dsfmt->status[i],
239
0
         &array[i + DSFMT_POS1 - DSFMT_N], &lung);
240
0
    }
241
0
    for (; i < size - DSFMT_N; i++) {
242
0
  do_recursion(&array[i], &array[i - DSFMT_N],
243
0
         &array[i + DSFMT_POS1 - DSFMT_N], &lung);
244
0
  convert_o0o1(&array[i - DSFMT_N]);
245
0
    }
246
0
    for (j = 0; j < 2 * DSFMT_N - size; j++) {
247
0
  dsfmt->status[j] = array[j + size - DSFMT_N];
248
0
    }
249
0
    for (; i < size; i++, j++) {
250
0
  do_recursion(&array[i], &array[i - DSFMT_N],
251
0
         &array[i + DSFMT_POS1 - DSFMT_N], &lung);
252
0
  dsfmt->status[j] = array[i];
253
0
  convert_o0o1(&array[i - DSFMT_N]);
254
0
    }
255
0
    for (i = size - DSFMT_N; i < size; i++) {
256
0
  convert_o0o1(&array[i]);
257
0
    }
258
0
    dsfmt->status[DSFMT_N] = lung;
259
0
}
260
261
/**
262
 * This function fills the user-specified array with double precision
263
 * floating point pseudorandom numbers of the IEEE 754 format.
264
 * @param dsfmt dsfmt state vector.
265
 * @param array an 128-bit array to be filled by pseudorandom numbers.
266
 * @param size number of 128-bit pseudorandom numbers to be generated.
267
 */
268
inline static void gen_rand_array_o0c1(dsfmt_t *dsfmt, dw128_t *array,
269
0
               int size) {
270
0
    int i, j;
271
0
    dw128_t lung;
272
0
273
0
    lung = dsfmt->status[DSFMT_N];
274
0
    do_recursion(&array[0], &dsfmt->status[0], &dsfmt->status[DSFMT_POS1],
275
0
     &lung);
276
0
    for (i = 1; i < DSFMT_N - DSFMT_POS1; i++) {
277
0
  do_recursion(&array[i], &dsfmt->status[i],
278
0
         &dsfmt->status[i + DSFMT_POS1], &lung);
279
0
    }
280
0
    for (; i < DSFMT_N; i++) {
281
0
  do_recursion(&array[i], &dsfmt->status[i],
282
0
         &array[i + DSFMT_POS1 - DSFMT_N], &lung);
283
0
    }
284
0
    for (; i < size - DSFMT_N; i++) {
285
0
  do_recursion(&array[i], &array[i - DSFMT_N],
286
0
         &array[i + DSFMT_POS1 - DSFMT_N], &lung);
287
0
  convert_o0c1(&array[i - DSFMT_N]);
288
0
    }
289
0
    for (j = 0; j < 2 * DSFMT_N - size; j++) {
290
0
  dsfmt->status[j] = array[j + size - DSFMT_N];
291
0
    }
292
0
    for (; i < size; i++, j++) {
293
0
  do_recursion(&array[i], &array[i - DSFMT_N],
294
0
         &array[i + DSFMT_POS1 - DSFMT_N], &lung);
295
0
  dsfmt->status[j] = array[i];
296
0
  convert_o0c1(&array[i - DSFMT_N]);
297
0
    }
298
0
    for (i = size - DSFMT_N; i < size; i++) {
299
0
  convert_o0c1(&array[i]);
300
0
    }
301
0
    dsfmt->status[DSFMT_N] = lung;
302
0
}
303
304
/**
305
 * This function represents a function used in the initialization
306
 * by init_by_array
307
 * @param x 32-bit integer
308
 * @return 32-bit integer
309
 */
310
0
static uint32_t ini_func1(uint32_t x) {
311
0
    return (x ^ (x >> 27)) * (uint32_t)1664525UL;
312
0
}
313
314
/**
315
 * This function represents a function used in the initialization
316
 * by init_by_array
317
 * @param x 32-bit integer
318
 * @return 32-bit integer
319
 */
320
0
static uint32_t ini_func2(uint32_t x) {
321
0
    return (x ^ (x >> 27)) * (uint32_t)1566083941UL;
322
0
}
323
324
/**
325
 * This function initializes the internal state array to fit the IEEE
326
 * 754 format.
327
 * @param dsfmt dsfmt state vector.
328
 */
329
111
static void initial_mask(dsfmt_t *dsfmt) {
330
111
    int i;
331
111
    uint64_t *psfmt;
332
111
333
111
    psfmt = &dsfmt->status[0].u[0];
334
42.5k
    for (i = 0; i < DSFMT_N * 2; 
i++42.4k
) {
335
42.4k
        psfmt[i] = (psfmt[i] & DSFMT_LOW_MASK) | DSFMT_HIGH_CONST;
336
42.4k
    }
337
111
}
338
339
/**
340
 * This function certificate the period of 2^{SFMT_MEXP}-1.
341
 * @param dsfmt dsfmt state vector.
342
 */
343
111
static void period_certification(dsfmt_t *dsfmt) {
344
111
    uint64_t pcv[2] = {DSFMT_PCV1, DSFMT_PCV2};
345
111
    uint64_t tmp[2];
346
111
    uint64_t inner;
347
111
    int i;
348
#if (DSFMT_PCV2 & 1) != 1
349
    int j;
350
    uint64_t work;
351
#endif
352
353
111
    tmp[0] = (dsfmt->status[DSFMT_N].u[0] ^ DSFMT_FIX1);
354
111
    tmp[1] = (dsfmt->status[DSFMT_N].u[1] ^ DSFMT_FIX2);
355
111
356
111
    inner = tmp[0] & pcv[0];
357
111
    inner ^= tmp[1] & pcv[1];
358
777
    for (i = 32; i > 0; 
i >>= 1666
) {
359
666
        inner ^= inner >> i;
360
666
    }
361
111
    inner &= 1;
362
111
    /* check OK */
363
111
    if (inner == 1) {
364
29
  return;
365
29
    }
366
82
    /* check NG, and modification */
367
82
#if (DSFMT_PCV2 & 1) == 1
368
82
    dsfmt->status[DSFMT_N].u[1] ^= 1;
369
#else
370
    for (i = 1; i >= 0; i--) {
371
  work = 1;
372
  for (j = 0; j < 64; j++) {
373
      if ((work & pcv[i]) != 0) {
374
    dsfmt->status[DSFMT_N].u[i] ^= work;
375
    return;
376
      }
377
      work = work << 1;
378
  }
379
    }
380
#endif
381
    return;
382
82
}
383
384
/*----------------
385
  PUBLIC FUNCTIONS
386
  ----------------*/
387
/**
388
 * This function returns the identification string.  The string shows
389
 * the Mersenne exponent, and all parameters of this generator.
390
 * @return id string.
391
 */
392
0
const char *dsfmt_get_idstring(void) {
393
0
    return DSFMT_IDSTR;
394
0
}
395
396
/**
397
 * This function returns the minimum size of array used for \b
398
 * fill_array functions.
399
 * @return minimum size of array used for fill_array functions.
400
 */
401
0
int dsfmt_get_min_array_size(void) {
402
0
    return DSFMT_N64;
403
0
}
404
405
/**
406
 * This function fills the internal state array with double precision
407
 * floating point pseudorandom numbers of the IEEE 754 format.
408
 * @param dsfmt dsfmt state vector.
409
 */
410
83.4k
void dsfmt_gen_rand_all(dsfmt_t *dsfmt) {
411
83.4k
    int i;
412
83.4k
    dw128_t lung;
413
83.4k
414
83.4k
    lung = dsfmt->status[DSFMT_N];
415
83.4k
    do_recursion(&dsfmt->status[0], &dsfmt->status[0],
416
83.4k
     &dsfmt->status[DSFMT_POS1], &lung);
417
6.17M
    for (i = 1; i < DSFMT_N - DSFMT_POS1; 
i++6.09M
) {
418
6.09M
  do_recursion(&dsfmt->status[i], &dsfmt->status[i],
419
6.09M
         &dsfmt->status[i + DSFMT_POS1], &lung);
420
6.09M
    }
421
9.85M
    for (; i < DSFMT_N; 
i++9.76M
) {
422
9.76M
  do_recursion(&dsfmt->status[i], &dsfmt->status[i],
423
9.76M
         &dsfmt->status[i + DSFMT_POS1 - DSFMT_N], &lung);
424
9.76M
    }
425
83.4k
    dsfmt->status[DSFMT_N] = lung;
426
83.4k
}
427
428
/**
429
 * This function generates double precision floating point
430
 * pseudorandom numbers which distribute in the range [1, 2) to the
431
 * specified array[] by one call. The number of pseudorandom numbers
432
 * is specified by the argument \b size, which must be at least (SFMT_MEXP
433
 * / 128) * 2 and a multiple of two.  The function
434
 * get_min_array_size() returns this minimum size.  The generation by
435
 * this function is much faster than the following fill_array_xxx functions.
436
 *
437
 * For initialization, init_gen_rand() or init_by_array() must be called
438
 * before the first call of this function. This function can not be
439
 * used after calling genrand_xxx functions, without initialization.
440
 *
441
 * @param dsfmt dsfmt state vector.
442
 * @param array an array where pseudorandom numbers are filled
443
 * by this function.  The pointer to the array must be "aligned"
444
 * (namely, must be a multiple of 16) in the SIMD version, since it
445
 * refers to the address of a 128-bit integer.  In the standard C
446
 * version, the pointer is arbitrary.
447
 *
448
 * @param size the number of 64-bit pseudorandom integers to be
449
 * generated.  size must be a multiple of 2, and greater than or equal
450
 * to (SFMT_MEXP / 128) * 2.
451
 *
452
 * @note \b memalign or \b posix_memalign is available to get aligned
453
 * memory. Mac OSX doesn't have these functions, but \b malloc of OSX
454
 * returns the pointer to the aligned memory block.
455
 */
456
void dsfmt_fill_array_close1_open2(dsfmt_t *dsfmt, double array[], int size) {
457
    assert(size % 2 == 0);
458
    assert(size >= DSFMT_N64);
459
    gen_rand_array_c1o2(dsfmt, (dw128_t *)array, size / 2);
460
}
461
462
/**
463
 * This function generates double precision floating point
464
 * pseudorandom numbers which distribute in the range (0, 1] to the
465
 * specified array[] by one call. This function is the same as
466
 * fill_array_close1_open2() except the distribution range.
467
 *
468
 * @param dsfmt dsfmt state vector.
469
 * @param array an array where pseudorandom numbers are filled
470
 * by this function.
471
 * @param size the number of pseudorandom numbers to be generated.
472
 * see also \sa fill_array_close1_open2()
473
 */
474
void dsfmt_fill_array_open_close(dsfmt_t *dsfmt, double array[], int size) {
475
    assert(size % 2 == 0);
476
    assert(size >= DSFMT_N64);
477
    gen_rand_array_o0c1(dsfmt, (dw128_t *)array, size / 2);
478
}
479
480
/**
481
 * This function generates double precision floating point
482
 * pseudorandom numbers which distribute in the range [0, 1) to the
483
 * specified array[] by one call. This function is the same as
484
 * fill_array_close1_open2() except the distribution range.
485
 *
486
 * @param array an array where pseudorandom numbers are filled
487
 * by this function.
488
 * @param dsfmt dsfmt state vector.
489
 * @param size the number of pseudorandom numbers to be generated.
490
 * see also \sa fill_array_close1_open2()
491
 */
492
void dsfmt_fill_array_close_open(dsfmt_t *dsfmt, double array[], int size) {
493
    assert(size % 2 == 0);
494
    assert(size >= DSFMT_N64);
495
    gen_rand_array_c0o1(dsfmt, (dw128_t *)array, size / 2);
496
}
497
498
/**
499
 * This function generates double precision floating point
500
 * pseudorandom numbers which distribute in the range (0, 1) to the
501
 * specified array[] by one call. This function is the same as
502
 * fill_array_close1_open2() except the distribution range.
503
 *
504
 * @param dsfmt dsfmt state vector.
505
 * @param array an array where pseudorandom numbers are filled
506
 * by this function.
507
 * @param size the number of pseudorandom numbers to be generated.
508
 * see also \sa fill_array_close1_open2()
509
 */
510
void dsfmt_fill_array_open_open(dsfmt_t *dsfmt, double array[], int size) {
511
    assert(size % 2 == 0);
512
    assert(size >= DSFMT_N64);
513
    gen_rand_array_o0o1(dsfmt, (dw128_t *)array, size / 2);
514
}
515
516
#if defined(__INTEL_COMPILER)
517
#  pragma warning(disable:981)
518
#endif
519
/**
520
 * This function initializes the internal state array with a 32-bit
521
 * integer seed.
522
 * @param dsfmt dsfmt state vector.
523
 * @param seed a 32-bit integer used as the seed.
524
 * @param mexp caller's mersenne expornent
525
 */
526
111
void dsfmt_chk_init_gen_rand(dsfmt_t *dsfmt, uint32_t seed, int mexp) {
527
111
    int i;
528
111
    uint32_t *psfmt;
529
111
530
111
    /* make sure caller program is compiled with the same MEXP */
531
111
    if (mexp != dsfmt_mexp) {
532
0
  fprintf(stderr, "DSFMT_MEXP doesn't match with dSFMT.c\n");
533
0
  exit(1);
534
0
    }
535
111
    psfmt = &dsfmt->status[0].u32[0];
536
111
    psfmt[idxof(0)] = seed;
537
85.2k
    for (i = 1; i < (DSFMT_N + 1) * 4; 
i++85.1k
) {
538
85.1k
        psfmt[idxof(i)] = 1812433253UL
539
85.1k
      * (psfmt[idxof(i - 1)] ^ (psfmt[idxof(i - 1)] >> 30)) + i;
540
85.1k
    }
541
111
    initial_mask(dsfmt);
542
111
    period_certification(dsfmt);
543
111
    dsfmt->idx = DSFMT_N64;
544
111
}
545
546
/**
547
 * This function initializes the internal state array,
548
 * with an array of 32-bit integers used as the seeds
549
 * @param dsfmt dsfmt state vector.
550
 * @param init_key the array of 32-bit integers, used as a seed.
551
 * @param key_length the length of init_key.
552
 * @param mexp caller's mersenne expornent
553
 */
554
void dsfmt_chk_init_by_array(dsfmt_t *dsfmt, uint32_t init_key[],
555
0
           int key_length, int mexp) {
556
0
    int i, j, count;
557
0
    uint32_t r;
558
0
    uint32_t *psfmt32;
559
0
    int lag;
560
0
    int mid;
561
0
    int size = (DSFMT_N + 1) * 4;  /* pulmonary */
562
0
563
0
    /* make sure caller program is compiled with the same MEXP */
564
0
    if (mexp != dsfmt_mexp) {
565
0
  fprintf(stderr, "DSFMT_MEXP doesn't match with dSFMT.c\n");
566
0
  exit(1);
567
0
    }
568
0
    if (size >= 623) {
569
0
  lag = 11;
570
0
    } else if (size >= 68) {
571
0
  lag = 7;
572
0
    } else if (size >= 39) {
573
0
  lag = 5;
574
0
    } else {
575
0
  lag = 3;
576
0
    }
577
0
    mid = (size - lag) / 2;
578
0
579
0
    psfmt32 = &dsfmt->status[0].u32[0];
580
0
    memset(dsfmt->status, 0x8b, sizeof(dsfmt->status));
581
0
    if (key_length + 1 > size) {
582
0
  count = key_length + 1;
583
0
    } else {
584
0
  count = size;
585
0
    }
586
0
    r = ini_func1(psfmt32[idxof(0)] ^ psfmt32[idxof(mid % size)]
587
0
      ^ psfmt32[idxof((size - 1) % size)]);
588
0
    psfmt32[idxof(mid % size)] += r;
589
0
    r += key_length;
590
0
    psfmt32[idxof((mid + lag) % size)] += r;
591
0
    psfmt32[idxof(0)] = r;
592
0
    count--;
593
0
    for (i = 1, j = 0; (j < count) && (j < key_length); j++) {
594
0
  r = ini_func1(psfmt32[idxof(i)]
595
0
          ^ psfmt32[idxof((i + mid) % size)]
596
0
          ^ psfmt32[idxof((i + size - 1) % size)]);
597
0
  psfmt32[idxof((i + mid) % size)] += r;
598
0
  r += init_key[j] + i;
599
0
  psfmt32[idxof((i + mid + lag) % size)] += r;
600
0
  psfmt32[idxof(i)] = r;
601
0
  i = (i + 1) % size;
602
0
    }
603
0
    for (; j < count; j++) {
604
0
  r = ini_func1(psfmt32[idxof(i)]
605
0
          ^ psfmt32[idxof((i + mid) % size)]
606
0
          ^ psfmt32[idxof((i + size - 1) % size)]);
607
0
  psfmt32[idxof((i + mid) % size)] += r;
608
0
  r += i;
609
0
  psfmt32[idxof((i + mid + lag) % size)] += r;
610
0
  psfmt32[idxof(i)] = r;
611
0
  i = (i + 1) % size;
612
0
    }
613
0
    for (j = 0; j < size; j++) {
614
0
  r = ini_func2(psfmt32[idxof(i)]
615
0
          + psfmt32[idxof((i + mid) % size)]
616
0
          + psfmt32[idxof((i + size - 1) % size)]);
617
0
  psfmt32[idxof((i + mid) % size)] ^= r;
618
0
  r -= i;
619
0
  psfmt32[idxof((i + mid + lag) % size)] ^= r;
620
0
  psfmt32[idxof(i)] = r;
621
0
  i = (i + 1) % size;
622
0
    }
623
0
    initial_mask(dsfmt);
624
0
    period_certification(dsfmt);
625
0
    dsfmt->idx = DSFMT_N64;
626
0
}
627
#if defined(__INTEL_COMPILER)
628
#  pragma warning(default:981)
629
#endif
630
631
#if defined(__cplusplus)
632
}
633
#endif