1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkMersenneTwister_Private.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13 =========================================================================*/
14
15 /* Many thanks to M. Matsumoto, T. Nishimura and M. Saito for the following */
16 /* implementation of their algorithm, the Mersenne Twister, taken from */
17 /* http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/DC/dc.html */
18
19 /* Dynamic Creation (DC) of Mersenne Twister generators */
20 /* */
21 /* Reference: */
22 /* Makoto Matsumoto and Takuji Nishimura, */
23 /* "Dynamic Creation of Pseudorandom Number Generators", */
24 /* Monte Carlo and Quasi-Monte Carlo Methods 1998, */
25 /* Springer, 2000, pp 56--69. */
26
27 /*
28 Copyright (C) 2001-2009 Makoto Matsumoto and Takuji Nishimura.
29 Copyright (C) 2009 Mutsuo Saito
30 All rights reserved.
31
32 Redistribution and use in source and binary forms, with or without
33 modification, are permitted provided that the following conditions are
34 met:
35
36 * Redistributions of source code must retain the above copyright
37 notice, this list of conditions and the following disclaimer.
38 * Redistributions in binary form must reproduce the above
39 copyright notice, this list of conditions and the following
40 disclaimer in the documentation and/or other materials provided
41 with the distribution.
42
43 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
44 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
45 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
46 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
47 OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
48 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
49 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
50 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
51 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
52 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
53 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
54 */
55
56 /* A C-program for MT19937: Integer version (1999/10/28) */
57 /* genrand() generates one pseudorandom unsigned integer (32bit) */
58 /* which is uniformly distributed among 0 to 2^32-1 for each */
59 /* call. sgenrand(seed) sets initial values to the working area */
60 /* of 624 words. Before genrand(), sgenrand(seed) must be */
61 /* called once. (seed is any 32-bit integer.) */
62 /* Coded by Takuji Nishimura, considering the suggestions by */
63 /* Topher Cooper and Marc Rieffel in July-Aug. 1997. */
64
65 /* When you use this, send an email to: matumoto@math.keio.ac.jp */
66 /* with an appropriate reference to your work. */
67
68 /* REFERENCE */
69 /* M. Matsumoto and T. Nishimura, */
70 /* "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform */
71 /* Pseudo-Random Number Generator", */
72 /* ACM Transactions on Modeling and Computer Simulation, */
73 /* Vol. 8, No. 1, January 1998, pp 3--30. */
74
75 #include <cstdio>
76 #include <cstdlib>
77 #include <cstring>
78
79 #include "vtkType.h"
80 typedef vtkTypeUInt32 uint32_t;
81 typedef vtkTypeUInt8 uint8_t;
82
83 #ifndef UINT32_C
84 #ifndef UINT32_MAX
85 #define UINT32_MAX 0xffffffff
86 #endif
87 #define UINT32_C(x) ((x) + (UINT32_MAX - UINT32_MAX))
88 #endif
89
90 // dc.h
91
92 typedef struct _mt_struct {
93 uint32_t aaa;
94 int mm,nn,rr,ww;
95 uint32_t wmask,umask,lmask;
96 int shift0, shift1, shiftB, shiftC;
97 uint32_t maskB, maskC;
98 int i;
99 uint32_t *state;
100 }mt_struct;
101
102 /* old interface */
103 void init_dc(uint32_t seed);
104 mt_struct *get_mt_parameter(int w, int p);
105 mt_struct *get_mt_parameter_id(int w, int p, int id);
106 mt_struct **get_mt_parameters(int w, int p, int max_id, int *count);
107
108 /* new interface */
109 mt_struct *get_mt_parameter_st(int w, int p, uint32_t seed);
110 mt_struct *get_mt_parameter_id_st(int w, int p, int id, uint32_t seed);
111 mt_struct **get_mt_parameters_st(int w, int p, int start_id, int max_id,
112 uint32_t seed, int *count);
113 /* common */
114 void free_mt_struct(mt_struct *mts);
115 void free_mt_struct_array(mt_struct **mtss, int count);
116 void sgenrand_mt(uint32_t seed, mt_struct *mts);
117 uint32_t genrand_mt(mt_struct *mts);
118
119 // mt19937.h
120
121 #define N 624
122
123 typedef struct _ORG_STATE {
124 uint32_t mt[N];
125 int mti;
126 } _org_state;
127
128 void _sgenrand_dc(_org_state *st, uint32_t seed);
129 uint32_t _genrand_dc(_org_state *st);
130
131 // dci.h
132
133 #define NOT_REJECTED 1
134 #define REJECTED 0
135 #define REDU 0
136 #define IRRED 1
137 #define NONREDU 1
138
139 extern _org_state global_mt19937;
140 typedef struct Polynomial_t {int *x; int deg;} Polynomial;
141
142 typedef struct PRESCR_T {
143 int sizeofA; /* parameter size */
144 uint32_t **modlist;
145 Polynomial **preModPolys;
146 } prescr_t;
147
148 typedef struct CHECK32_T {
149 uint32_t upper_mask;
150 uint32_t lower_mask;
151 uint32_t word_mask;
152 } check32_t;
153
154 typedef struct EQDEG_T {
155 uint32_t bitmask[32];
156 uint32_t mask_b;
157 uint32_t mask_c;
158 uint32_t upper_v_bits;
159 int shift_0;
160 int shift_1;
161 int shift_s;
162 int shift_t;
163 int mmm;
164 int nnn;
165 int rrr;
166 int www;
167 uint32_t aaa[2];
168 uint32_t gupper_mask; /** most significant (WWW - RRR) bits **/
169 uint32_t glower_mask; /** least significant RRR bits **/
170 uint32_t greal_mask; /** upper WWW bitmask **/
171 int ggap; /** difference between machine wordsize and dest wordsize **/
172 int gcur_maxlengs[32]; /** for optimize_v_hard **/
173 uint32_t gmax_b, gmax_c;
174 } eqdeg_t;
175
176 int _prescreening_dc(prescr_t *pre, uint32_t aaa);
177 void _InitPrescreening_dc(prescr_t *pre, int m, int n, int r, int w);
178 void _EndPrescreening_dc(prescr_t *pre);
179 int _CheckPeriod_dc(check32_t *ck, _org_state *st,
180 uint32_t a, int m, int n, int r, int w);
181 void _get_tempering_parameter_dc(mt_struct *mts);
182 void _get_tempering_parameter_hard_dc(mt_struct *mts);
183 void _InitCheck32_dc(check32_t *ck, int r, int w);
184
185 // eqdeg.c
186
187 /**************************************/
188 #define SSS 7
189 #define TTT 15
190 /* #define S00 11 */
191 #define S00 12
192 #define S01 18
193 /**************************************/
194
195 /** for get_tempering_parameter_hard **/
196 #define LIMIT_V_BEST_OPT 15
197 /**************************************/
198
199 #define WORD_LEN 32
200 #define MIN_INFINITE (-2147483647-1)
201
202 typedef struct Vector_t{
203 uint32_t *cf; /* fraction part */ // status
204 int start; /* beginning of fraction part */ // idx
205 int count; /* maximum (degree) */
206 uint32_t next; /* (bp) rm (shifted&bitmasked) at the maximum degree */
207 } Vector;
208
209 typedef struct mask_node{
210 uint32_t b,c;
211 int v,leng;
212 struct mask_node *next;
213 } MaskNode;
214
trnstmp(eqdeg_t * eq,uint32_t tmp)215 static inline uint32_t trnstmp(eqdeg_t *eq, uint32_t tmp) {
216 tmp ^= (tmp >> eq->shift_0) & eq->greal_mask;
217 return tmp;
218 }
219
masktmp(eqdeg_t * eq,uint32_t tmp)220 static inline uint32_t masktmp(eqdeg_t *eq, uint32_t tmp) {
221 tmp ^= (tmp << eq->shift_s) & eq->mask_b;
222 tmp ^= (tmp << eq->shift_t) & eq->mask_c;
223 return tmp;
224 }
225
lsb(eqdeg_t * eq,uint32_t x)226 static inline uint32_t lsb(eqdeg_t *eq, uint32_t x) {
227 return (x >> eq->ggap) & 1;
228 }
229
230 static const uint8_t pivot_calc_tbl[256] = {
231 0, 8, 7, 8, 6, 8, 7, 8, 5, 8, 7, 8, 6, 8, 7, 8,
232 4, 8, 7, 8, 6, 8, 7, 8, 5, 8, 7, 8, 6, 8, 7, 8,
233 3, 8, 7, 8, 6, 8, 7, 8, 5, 8, 7, 8, 6, 8, 7, 8,
234 4, 8, 7, 8, 6, 8, 7, 8, 5, 8, 7, 8, 6, 8, 7, 8,
235 2, 8, 7, 8, 6, 8, 7, 8, 5, 8, 7, 8, 6, 8, 7, 8,
236 4, 8, 7, 8, 6, 8, 7, 8, 5, 8, 7, 8, 6, 8, 7, 8,
237 3, 8, 7, 8, 6, 8, 7, 8, 5, 8, 7, 8, 6, 8, 7, 8,
238 4, 8, 7, 8, 6, 8, 7, 8, 5, 8, 7, 8, 6, 8, 7, 8,
239 1, 8, 7, 8, 6, 8, 7, 8, 5, 8, 7, 8, 6, 8, 7, 8,
240 4, 8, 7, 8, 6, 8, 7, 8, 5, 8, 7, 8, 6, 8, 7, 8,
241 3, 8, 7, 8, 6, 8, 7, 8, 5, 8, 7, 8, 6, 8, 7, 8,
242 4, 8, 7, 8, 6, 8, 7, 8, 5, 8, 7, 8, 6, 8, 7, 8,
243 2, 8, 7, 8, 6, 8, 7, 8, 5, 8, 7, 8, 6, 8, 7, 8,
244 4, 8, 7, 8, 6, 8, 7, 8, 5, 8, 7, 8, 6, 8, 7, 8,
245 3, 8, 7, 8, 6, 8, 7, 8, 5, 8, 7, 8, 6, 8, 7, 8,
246 4, 8, 7, 8, 6, 8, 7, 8, 5, 8, 7, 8, 6, 8, 7, 8,
247 };
248
249 static int calc_pivot(uint32_t v);
250 static int push_stack(eqdeg_t *eq, uint32_t b, uint32_t c,
251 int v, uint32_t *bbb, uint32_t *ccc);
252 static int push_mask(eqdeg_t * eq, int l, int v,
253 uint32_t b, uint32_t c, uint32_t *bbb, uint32_t *ccc);
254 static int pivot_reduction(eqdeg_t *eq, int v);
255 static void init_tempering(eqdeg_t *eq, mt_struct *mts);
256 static void freeVector_t( Vector *v );
257 static void free_lattice( Vector **lattice, int v);
258 static void add(int nnn, Vector *u, Vector *v);
259 static void optimize_v(eqdeg_t *eq, uint32_t b, uint32_t c, int v);
260 static MaskNode *optimize_v_hard(eqdeg_t *eq, int v, MaskNode *prev);
261 static Vector *newVector_t(int nnn);
262 static Vector **make_lattice(eqdeg_t *eq, int v);
263 static void delete_MaskNodes(MaskNode *head);
264 static MaskNode *delete_lower_MaskNodes(MaskNode *head, int l);
265 static MaskNode *cons_MaskNode(MaskNode *head, uint32_t b, uint32_t c, int leng);
266 /* static void count_MaskNodes(MaskNode *head); */
267 static void next_state(eqdeg_t *eq, Vector *v, int *count);
268
_get_tempering_parameter_dc(mt_struct * mts)269 void _get_tempering_parameter_dc(mt_struct *mts)
270 {
271 eqdeg_t eq;
272 init_tempering(&eq, mts);
273 optimize_v(&eq, 0, 0, 0);
274 mts->shift0 = eq.shift_0;
275 mts->shift1 = eq.shift_1;
276 mts->shiftB = eq.shift_s;
277 mts->shiftC = eq.shift_t;
278 mts->maskB = eq.mask_b >> eq.ggap;
279 mts->maskC = eq.mask_c >> eq.ggap;
280 }
281
_get_tempering_parameter_hard_dc(mt_struct * mts)282 void _get_tempering_parameter_hard_dc(mt_struct *mts)
283 {
284 int i;
285 MaskNode mn0, *cur, *next;
286 eqdeg_t eq;
287
288 init_tempering(&eq, mts);
289
290 for (i=0; i<eq.www; i++)
291 eq.gcur_maxlengs[i] = -1;
292
293 mn0.b = mn0.c = mn0.leng = 0;
294 mn0.next = nullptr;
295
296 cur = &mn0;
297 for (i=0; i<LIMIT_V_BEST_OPT; i++) {
298 next = optimize_v_hard(&eq, i, cur);
299 if (i > 0)
300 delete_MaskNodes(cur);
301 cur = next;
302 }
303 delete_MaskNodes(cur);
304
305 optimize_v(&eq, eq.gmax_b, eq.gmax_c,i);
306 mts->shift0 = eq.shift_0;
307 mts->shift1 = eq.shift_1;
308 mts->shiftB = eq.shift_s;
309 mts->shiftC = eq.shift_t;
310 mts->maskB = eq.mask_b >> eq.ggap;
311 mts->maskC = eq.mask_c >> eq.ggap;
312
313 /* show_distrib(mts); */
314 }
315
calc_pivot(uint32_t v)316 static int calc_pivot(uint32_t v) {
317 int p1, p2, p3, p4;
318
319 p1 = pivot_calc_tbl[v & 0xff];
320 if (p1) {
321 return p1 + 24 - 1;
322 }
323 p2 = pivot_calc_tbl[(v >> 8) & 0xff];
324 if (p2) {
325 return p2 + 16 - 1;
326 }
327 p3 = pivot_calc_tbl[(v >> 16) & 0xff];
328 if (p3) {
329 return p3 + 8 - 1;
330 }
331 p4 = pivot_calc_tbl[(v >> 24) & 0xff];
332 if (p4) {
333 return p4 - 1;
334 }
335 return -1;
336 }
337
is_zero(int size,Vector * v)338 static int is_zero(int size, Vector *v) {
339 if (v->cf[0] != 0) {
340 return 0;
341 } else {
342 return (memcmp(v->cf, v->cf + 1, sizeof(uint32_t) * (size - 1)) == 0);
343 }
344 }
345
init_tempering(eqdeg_t * eq,mt_struct * mts)346 static void init_tempering(eqdeg_t *eq, mt_struct *mts)
347 {
348 int i;
349
350 eq->mmm = mts->mm;
351 eq->nnn = mts->nn;
352 eq->rrr = mts->rr;
353 eq->www = mts->ww;
354 eq->shift_0 = S00;
355 eq->shift_1 = S01;
356 eq->shift_s = SSS;
357 eq->shift_t = TTT;
358 eq->ggap = WORD_LEN - eq->www;
359 /* bits are filled in mts->aaa from MSB */
360 eq->aaa[0] = 0; eq->aaa[1] = (mts->aaa) << eq->ggap;
361
362
363 for( i=0; i<WORD_LEN; i++)
364 eq->bitmask[i] = UINT32_C(0x80000000) >> i;
365
366 for( i=0, eq->glower_mask=0; i<eq->rrr; i++)
367 eq->glower_mask = (eq->glower_mask<<1)| 0x1;
368
369 eq->gupper_mask = ~eq->glower_mask;
370 eq->gupper_mask <<= eq->ggap;
371 eq->glower_mask <<= eq->ggap;
372
373 eq->greal_mask = (eq->gupper_mask | eq->glower_mask);
374 }
375
376 /* (v-1) bitmasks of b,c */
optimize_v_hard(eqdeg_t * eq,int v,MaskNode * prev_masks)377 static MaskNode *optimize_v_hard(eqdeg_t *eq, int v, MaskNode *prev_masks)
378 {
379 int i, ll, t;
380 uint32_t bbb[8], ccc[8];
381 MaskNode *cur_masks;
382
383 cur_masks = nullptr;
384
385 while (prev_masks != nullptr) {
386
387 ll = push_stack(eq, prev_masks->b,prev_masks->c,v,bbb,ccc);
388
389 for (i=0; i<ll; ++i) {
390 eq->mask_b = bbb[i];
391 eq->mask_c = ccc[i];
392 t = pivot_reduction(eq, v+1);
393 if (t >= eq->gcur_maxlengs[v]) {
394 eq->gcur_maxlengs[v] = t;
395 eq->gmax_b = eq->mask_b;
396 eq->gmax_c = eq->mask_c;
397 cur_masks = cons_MaskNode(cur_masks, eq->mask_b, eq->mask_c, t);
398 }
399 }
400 prev_masks = prev_masks->next;
401 }
402
403 cur_masks = delete_lower_MaskNodes(cur_masks, eq->gcur_maxlengs[v]);
404
405 return cur_masks;
406 }
407
408
409 /* (v-1) bitmasks of b,c */
optimize_v(eqdeg_t * eq,uint32_t b,uint32_t c,int v)410 static void optimize_v(eqdeg_t *eq, uint32_t b, uint32_t c, int v)
411 {
412 int i, max_len, max_i, ll, t;
413 uint32_t bbb[8], ccc[8];
414
415 ll = push_stack(eq, b,c,v,bbb,ccc);
416
417 max_len = max_i = 0;
418 if (ll > 1) {
419 for (i=0; i<ll; ++i) {
420 eq->mask_b = bbb[i];
421 eq->mask_c = ccc[i];
422 t = pivot_reduction(eq, v+1);
423 if (t > max_len) {
424 max_len = t;
425 max_i = i;
426 }
427 }
428 }
429
430 if ( v >= eq->www-1 ) {
431 eq->mask_b = bbb[max_i];
432 eq->mask_c = ccc[max_i];
433 return;
434 }
435
436 optimize_v(eq, bbb[max_i], ccc[max_i], v+1);
437 }
438
push_stack(eqdeg_t * eq,uint32_t b,uint32_t c,int v,uint32_t * bbb,uint32_t * ccc)439 static int push_stack(eqdeg_t *eq, uint32_t b, uint32_t c, int v,
440 uint32_t *bbb, uint32_t *ccc)
441 {
442 int i, ll, ncv;
443 uint32_t cv_buf[2];
444
445 ll = 0;
446
447 if( (v+eq->shift_t) < eq->www ){
448 ncv = 2; cv_buf[0] = c | eq->bitmask[v]; cv_buf[1] = c;
449 }
450 else {
451 ncv = 1; cv_buf[0] = c;
452 }
453
454 for( i=0; i<ncv; ++i)
455 ll += push_mask(eq, ll, v, b, cv_buf[i], bbb, ccc);
456
457 return ll;
458 }
459
push_mask(eqdeg_t * eq,int l,int v,uint32_t b,uint32_t c,uint32_t * bbb,uint32_t * ccc)460 static int push_mask(eqdeg_t *eq, int l, int v, uint32_t b, uint32_t c,
461 uint32_t *bbb, uint32_t *ccc)
462 {
463 int i, j, k, nbv, nbvt;
464 uint32_t bmask, bv_buf[2], bvt_buf[2];
465
466 k = l;
467 if( (eq->shift_s+v) >= eq->www ){
468 nbv = 1; bv_buf[0] = 0;
469 }
470 else if( (v>=eq->shift_t) && (c&eq->bitmask[v-eq->shift_t] ) ){
471 nbv = 1; bv_buf[0] = b&eq->bitmask[v];
472 }
473 else {
474 nbv = 2; bv_buf[0] = eq->bitmask[v]; bv_buf[1] = 0;
475 }
476
477 if( ((v+eq->shift_t+eq->shift_s) < eq->www) && (c&eq->bitmask[v]) ){
478 nbvt = 2; bvt_buf[0] = eq->bitmask[v+eq->shift_t]; bvt_buf[1] = 0;
479 }
480 else {
481 nbvt = 1; bvt_buf[0] = 0;
482 }
483
484 bmask = eq->bitmask[v];
485 if( (v+eq->shift_t) < eq->www )
486 bmask |= eq->bitmask[v+eq->shift_t];
487 bmask = ~bmask;
488 for( i=0; i<nbvt; ++i){
489 for( j=0; j<nbv; ++j){
490 bbb[k] = (b&bmask) | bv_buf[j] | bvt_buf[i];
491 ccc[k] = c;
492 ++k;
493 }
494 }
495
496 return k-l;
497 }
498
499
500 /**********************************/
501 /**** subroutines for lattice ****/
502 /**********************************/
pivot_reduction(eqdeg_t * eq,int v)503 static int pivot_reduction(eqdeg_t *eq, int v)
504 {
505 Vector **lattice, *ltmp;
506 int i;
507 int pivot;
508 int count;
509 int min;
510
511 eq->upper_v_bits = 0;
512 for( i=0; i<v; i++) {
513 eq->upper_v_bits |= eq->bitmask[i];
514 }
515
516 lattice = make_lattice(eq, v );
517
518 for (;;) {
519 pivot = calc_pivot(lattice[v]->next);
520 if (lattice[pivot]->count < lattice[v]->count) {
521 ltmp = lattice[pivot];
522 lattice[pivot] = lattice[v];
523 lattice[v] = ltmp;
524 }
525 add(eq->nnn, lattice[v], lattice[pivot]);
526 if (lattice[v]->next == 0) {
527 count = 0;
528 next_state(eq, lattice[v], &count);
529 if (lattice[v]->next == 0) {
530 if (is_zero(eq->nnn, lattice[v])) {
531 break;
532 }
533 while (lattice[v]->next == 0) {
534 count++;
535 next_state(eq, lattice[v], &count);
536 if (count > eq->nnn * (eq->www-1) - eq->rrr) {
537 break;
538 }
539 }
540 if (lattice[v]->next == 0) {
541 break;
542 }
543 }
544 }
545 }
546
547 min = lattice[0]->count;
548 for (i = 1; i < v; i++) {
549 if (min > lattice[i]->count) {
550 min = lattice[i]->count;
551 }
552 }
553 free_lattice( lattice, v );
554 return min;
555 }
556
557
558
559
560 /********************************/
561 /** allocate momory for Vector **/
562 /********************************/
newVector_t(int nnn)563 static Vector *newVector_t(int nnn)
564 {
565 Vector *v;
566
567 v = (Vector *)malloc( sizeof( Vector ) );
568 if( v == nullptr ){
569 printf("malloc error in \"newVector_t()\"\n");
570 exit(1);
571 }
572
573 v->cf = (uint32_t *)calloc( nnn, sizeof( uint32_t ) );
574 if( v->cf == nullptr ){
575 printf("calloc error in \"newVector_t()\"\n");
576 exit(1);
577 }
578
579 v->start = 0;
580
581 return v;
582 }
583
584
585 /************************************************/
586 /* frees *v which was allocated by newVector_t() */
587 /************************************************/
freeVector_t(Vector * v)588 static void freeVector_t( Vector *v )
589 {
590 if( nullptr != v->cf ) free( v->cf );
591 if( nullptr != v ) free( v );
592 }
593
free_lattice(Vector ** lattice,int v)594 static void free_lattice( Vector **lattice, int v)
595 {
596 int i;
597
598 for( i=0; i<=v; i++)
599 freeVector_t( lattice[i] );
600 free( lattice );
601 }
602
603 /* adds v to u (then u will change) */
add(int nnn,Vector * u,Vector * v)604 static void add(int nnn, Vector *u, Vector *v)
605 {
606 int i;
607 int diff = (v->start - u->start + nnn) % nnn;
608 for (i = 0; i < nnn - diff; i++) {
609 u->cf[i] ^= v->cf[i + diff];
610 }
611 diff = diff - nnn;
612 for (; i < nnn; i++) {
613 u->cf[i] ^= v->cf[i + diff];
614 }
615 u->next ^= v->next;
616 }
617
618 /* makes a initial lattice */
make_lattice(eqdeg_t * eq,int v)619 static Vector **make_lattice(eqdeg_t *eq, int v)
620 {
621 int i;
622 int count;
623 Vector **lattice, *bottom;
624
625 lattice = (Vector **)malloc( (v+1) * sizeof( Vector *) );
626 if( nullptr == lattice ){
627 printf("malloc error in \"make_lattice\"\n");
628 exit(1);
629 }
630
631 for( i=0; i<v; i++){ /* from 0th row to v-1-th row */
632 lattice[i] = newVector_t(eq->nnn);
633 lattice[i]->next = eq->bitmask[i];
634 lattice[i]->start = 0;
635 lattice[i]->count = 0;
636 }
637
638 bottom = newVector_t(eq->nnn); /* last row */
639 for(i=0; i< eq->nnn; i++) {
640 bottom->cf[i] = 0;
641 }
642 bottom->cf[eq->nnn -1] = 0xc0000000 & eq->greal_mask;
643 bottom->start = 0;
644 bottom->count = 0;
645 count = 0;
646 do {
647 next_state(eq, bottom, &count);
648 } while (bottom->next == 0);
649 // degree_of_vector(eq, top );
650 lattice[v] = bottom;
651
652 return lattice;
653 }
654
next_state(eqdeg_t * eq,Vector * v,int * count)655 static void next_state(eqdeg_t *eq, Vector *v, int *count) {
656 uint32_t tmp;
657
658 do {
659 tmp = ( v->cf[v->start] & eq->gupper_mask )
660 | ( v->cf[(v->start + 1) % eq->nnn] & eq->glower_mask );
661 v->cf[v->start] = v->cf[(v->start + eq->mmm) % eq->nnn]
662 ^ ( (tmp>>1) ^ eq->aaa[lsb(eq, tmp)] );
663 v->cf[v->start] &= eq->greal_mask;
664 tmp = v->cf[v->start];
665 v->start = (v->start + 1) % eq->nnn;
666 v->count++;
667 tmp = trnstmp(eq, tmp);
668 tmp = masktmp(eq, tmp);
669 v->next = tmp & eq->upper_v_bits;
670 (*count)++;
671 if (*count > eq->nnn * (eq->www-1) - eq->rrr) {
672 break;
673 }
674 } while (v->next == 0);
675 }
676
677 /***********/
cons_MaskNode(MaskNode * head,uint32_t b,uint32_t c,int leng)678 static MaskNode *cons_MaskNode(MaskNode *head, uint32_t b, uint32_t c, int leng)
679 {
680 MaskNode *t;
681
682 t = (MaskNode*)malloc(sizeof(MaskNode));
683 if (t == nullptr) {
684 printf("malloc error in \"cons_MaskNode\"\n");
685 exit(1);
686 }
687
688 t->b = b;
689 t->c = c;
690 t->leng = leng;
691 t->next = head;
692
693 return t;
694 }
695
delete_MaskNodes(MaskNode * head)696 static void delete_MaskNodes(MaskNode *head)
697 {
698 MaskNode *t;
699
700 while(head != nullptr) {
701 t = head->next;
702 free(head);
703 head = t;
704 }
705 }
706
delete_lower_MaskNodes(MaskNode * head,int l)707 static MaskNode *delete_lower_MaskNodes(MaskNode *head, int l)
708 {
709 MaskNode *s, *t, *tail;
710
711 s = head;
712 for(;;) { /* heading */
713 if (s == nullptr)
714 return nullptr;
715 if (s->leng >= l)
716 break;
717 t = s->next;
718 free(s);
719 s = t;
720 }
721
722 head = tail = s;
723
724 while (head != nullptr) {
725 t = head->next;
726 if (head->leng < l) {
727 free(head);
728 head = nullptr;
729 }
730 else {
731 tail->next = head;
732 tail = head;
733 }
734 head = t;
735 }
736
737 tail->next = nullptr;
738 return s;
739 }
740
741 // mt19937.c
742
743 /* Period parameters */
744 /* #define N 624 */
745 #define M 397
746 #define MATRIX_A UINT32_C(0x9908b0df) /* constant vector a */
747 #define UPPER_MASK UINT32_C(0x80000000) /* most significant w-r bits */
748 #define LOWER_MASK UINT32_C(0x7fffffff) /* least significant r bits */
749
750 /* Tempering parameters */
751 #define TEMPERING_MASK_B UINT32_C(0x9d2c5680)
752 #define TEMPERING_MASK_C UINT32_C(0xefc60000)
753 #define TEMPERING_SHIFT_U(y) (y >> 11)
754 #define TEMPERING_SHIFT_S(y) (y << 7)
755 #define TEMPERING_SHIFT_T(y) (y << 15)
756 #define TEMPERING_SHIFT_L(y) (y >> 18)
757
758 /* Initializing the array with a seed */
_sgenrand_dc(_org_state * st,uint32_t seed)759 void _sgenrand_dc(_org_state *st, uint32_t seed)
760 {
761 int i;
762
763 for (i=0;i<N;i++) {
764 st->mt[i] = seed;
765 seed = (UINT32_C(1812433253) * (seed ^ (seed >> 30))) + i + 1;
766 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
767 /* In the previous versions, MSBs of the seed affect */
768 /* only MSBs of the array mt[]. */
769 }
770 st->mti = N;
771 }
772
773
_genrand_dc(_org_state * st)774 uint32_t _genrand_dc(_org_state *st)
775 {
776 uint32_t y;
777 static const uint32_t mag01[2]={0x0, MATRIX_A};
778 /* mag01[x] = x * MATRIX_A for x=0,1 */
779
780 if (st->mti >= N) { /* generate N words at one time */
781 int kk;
782
783 for (kk=0;kk<N-M;kk++) {
784 y = (st->mt[kk]&UPPER_MASK)|(st->mt[kk+1]&LOWER_MASK);
785 st->mt[kk] = st->mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1];
786 }
787 for (;kk<N-1;kk++) {
788 y = (st->mt[kk]&UPPER_MASK)|(st->mt[kk+1]&LOWER_MASK);
789 st->mt[kk] = st->mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1];
790 }
791 y = (st->mt[N-1]&UPPER_MASK)|(st->mt[0]&LOWER_MASK);
792 st->mt[N-1] = st->mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];
793
794 st->mti = 0;
795 }
796
797 y = st->mt[st->mti++];
798 y ^= TEMPERING_SHIFT_U(y);
799 y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
800 y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
801 y ^= TEMPERING_SHIFT_L(y);
802
803 return y;
804 }
805
806 // seive.c
807
808 #define WORDLEN 32
809 #define LSB 0x1
810 #define MAX_SEARCH 10000
811
812
813 _org_state global_mt19937;
814 /*******************************************************************/
815 static uint32_t nextA(_org_state *org, int w);
816 static uint32_t nextA_id(_org_state *org, int w, int id, int idw);
817 static void make_masks(int r, int w, mt_struct *mts);
818 static int get_irred_param(check32_t *ck, prescr_t *pre, _org_state *org,
819 mt_struct *mts,int id, int idw);
820 static mt_struct *alloc_mt_struct(int n);
821 static mt_struct *init_mt_search(check32_t *ck, prescr_t *pre, int w, int p);
822 static void end_mt_search(prescr_t *pre);
823 static void copy_params_of_mt_struct(mt_struct *src, mt_struct *dst);
824 static int proper_mersenne_exponent(int p);
825 /*******************************************************************/
826
827 /* When idw==0, id is not embedded into "a" */
828 #define FOUND 1
829 #define NOT_FOUND 0
get_irred_param(check32_t * ck,prescr_t * pre,_org_state * org,mt_struct * mts,int id,int idw)830 static int get_irred_param(check32_t *ck, prescr_t *pre, _org_state *org,
831 mt_struct *mts, int id, int idw)
832 {
833 int i;
834 uint32_t a;
835
836 for (i=0; i<MAX_SEARCH; i++) {
837 if (idw == 0)
838 a = nextA(org, mts->ww);
839 else
840 a = nextA_id(org, mts->ww, id, idw);
841 if (NOT_REJECTED == _prescreening_dc(pre, a) ) {
842 if (IRRED
843 == _CheckPeriod_dc(ck, org, a,mts->mm,mts->nn,mts->rr,mts->ww)) {
844 mts->aaa = a;
845 break;
846 }
847 }
848 }
849
850 if (MAX_SEARCH == i) return NOT_FOUND;
851 return FOUND;
852 }
853
854
nextA(_org_state * org,int w)855 static uint32_t nextA(_org_state *org, int w)
856 {
857 uint32_t x, word_mask;
858
859 word_mask = 0xFFFFFFFF;
860 word_mask <<= WORDLEN - w;
861 word_mask >>= WORDLEN - w;
862
863 x = _genrand_dc(org);
864 x &= word_mask;
865 x |= (LSB << (w-1));
866
867 return x;
868 }
869
nextA_id(_org_state * org,int w,int id,int idw)870 static uint32_t nextA_id(_org_state *org, int w, int id, int idw)
871 {
872 uint32_t x, word_mask;
873
874 word_mask = 0xFFFFFFFF;
875 word_mask <<= WORDLEN - w;
876 word_mask >>= WORDLEN - w;
877 word_mask >>= idw;
878 word_mask <<= idw;
879
880 x = _genrand_dc(org);
881 x &= word_mask;
882 x |= (LSB << (w-1));
883 x |= (uint32_t)id; /* embedding id */
884
885 return x;
886 }
887
make_masks(int r,int w,mt_struct * mts)888 static void make_masks(int r, int w, mt_struct *mts)
889 {
890 int i;
891 uint32_t ut, wm, um, lm;
892
893 wm = 0xFFFFFFFF;
894 wm >>= (WORDLEN - w);
895
896 ut = 0;
897 for (i=0; i<r; i++) {
898 ut <<= 1;
899 ut |= LSB;
900 }
901
902 lm = ut;
903 um = (~ut) & wm;
904
905 mts->wmask = wm;
906 mts->umask = um;
907 mts->lmask = lm;
908 }
909
init_mt_search(check32_t * ck,prescr_t * pre,int w,int p)910 static mt_struct *init_mt_search(check32_t *ck, prescr_t *pre, int w, int p)
911 {
912 int n, m, r;
913 mt_struct *mts;
914
915 if ( (w>32) || (w<31) ) {
916 printf ("Sorry, currently only w = 32 or 31 is allowded.\n");
917 return nullptr;
918 }
919
920 if ( !proper_mersenne_exponent(p) ) {
921 if (p<521) {
922 printf ("\"p\" is too small.\n");
923 return nullptr;
924 }
925 else if (p>44497){
926 printf ("\"p\" is too large.\n");
927 return nullptr;
928 }
929 else {
930 printf ("\"p\" is not a Mersenne exponent.\n");
931 return nullptr;
932 }
933 }
934
935 n = p/w + 1; /* since p is Mersenne Exponent, w never divids p */
936 mts = alloc_mt_struct(n);
937 if (nullptr == mts) return nullptr;
938
939 m = n/2;
940 if (m < 2) m = n-1;
941 r = n * w - p;
942
943 make_masks(r, w, mts);
944 _InitPrescreening_dc(pre, m, n, r, w);
945 _InitCheck32_dc(ck, r, w);
946
947 mts->mm = m;
948 mts->nn = n;
949 mts->rr = r;
950 mts->ww = w;
951
952 return mts;
953 }
954
end_mt_search(prescr_t * pre)955 static void end_mt_search(prescr_t *pre)
956 {
957 _EndPrescreening_dc(pre);
958 }
959
960 /*
961 w -- word size
962 p -- Mersenne Exponent
963 seed -- seed for original mt19937 to generate parameter.
964 */
get_mt_parameter_st(int w,int p,uint32_t seed)965 mt_struct *get_mt_parameter_st(int w, int p, uint32_t seed)
966 {
967 mt_struct *mts;
968 prescr_t pre;
969 _org_state org;
970 check32_t ck;
971
972 _sgenrand_dc(&org, seed);
973 mts = init_mt_search(&ck, &pre, w, p);
974 if (mts == nullptr) return nullptr;
975
976 if ( NOT_FOUND == get_irred_param(&ck, &pre, &org, mts,0,0) ) {
977 free_mt_struct(mts);
978 return nullptr;
979 }
980 _get_tempering_parameter_hard_dc(mts);
981 end_mt_search(&pre);
982
983 return mts;
984 }
985
986 /*
987 w -- word size
988 p -- Mersenne Exponent
989 */
get_mt_parameter(int w,int p)990 mt_struct *get_mt_parameter(int w, int p)
991 {
992 mt_struct *mts;
993 prescr_t pre;
994 check32_t ck;
995
996 mts = init_mt_search(&ck, &pre, w, p);
997 if (mts == nullptr) return nullptr;
998
999 if ( NOT_FOUND == get_irred_param(&ck, &pre, &global_mt19937, mts,0,0) ) {
1000 free_mt_struct(mts);
1001 return nullptr;
1002 }
1003 _get_tempering_parameter_hard_dc(mts);
1004 end_mt_search(&pre);
1005
1006 return mts;
1007 }
1008
1009 /*
1010 w -- word size
1011 p -- Mersenne Exponent
1012 */
1013 #if 0
1014 mt_struct *get_mt_parameter_opt_temper(int w, int p, uint32_t seed)
1015 {
1016 mt_struct *mts;
1017 prescr_t pre;
1018 _org_state org;
1019 check32_t ck;
1020
1021 _sgenrand_dc(&org, seed);
1022 mts = init_mt_search(&ck, &pre, w, p);
1023 if (mts == nullptr) return nullptr;
1024
1025 if ( NOT_FOUND == get_irred_param(&ck, &pre, &org, mts,0,0) ) {
1026 free_mt_struct(mts);
1027 return nullptr;
1028 }
1029 _get_tempering_parameter_hard_dc(mts);
1030 end_mt_search(&pre);
1031
1032 return mts;
1033 }
1034 #endif
1035 /*
1036 w -- word size
1037 p -- Mersenne Exponent
1038 */
1039 #define DEFAULT_ID_SIZE 16
1040 /* id <= 0xffff */
get_mt_parameter_id_st(int w,int p,int id,uint32_t seed)1041 mt_struct *get_mt_parameter_id_st(int w, int p, int id, uint32_t seed)
1042 {
1043 mt_struct *mts;
1044 prescr_t pre;
1045 _org_state org;
1046 check32_t ck;
1047
1048 _sgenrand_dc(&org, seed);
1049 if (id > 0xffff) {
1050 printf("\"id\" must be less than 65536\n");
1051 return nullptr;
1052 }
1053 if (id < 0) {
1054 printf("\"id\" must be positive\n");
1055 return nullptr;
1056 }
1057
1058 mts = init_mt_search(&ck, &pre, w, p);
1059 if (mts == nullptr) return nullptr;
1060
1061 if ( NOT_FOUND == get_irred_param(&ck, &pre, &org,
1062 mts, id, DEFAULT_ID_SIZE) ) {
1063 free_mt_struct(mts);
1064 return nullptr;
1065 }
1066 _get_tempering_parameter_hard_dc(mts);
1067 end_mt_search(&pre);
1068
1069 return mts;
1070 }
1071
get_mt_parameter_id(int w,int p,int id)1072 mt_struct *get_mt_parameter_id(int w, int p, int id)
1073 {
1074 mt_struct *mts;
1075 prescr_t pre;
1076 check32_t ck;
1077
1078 if (id > 0xffff) {
1079 printf("\"id\" must be less than 65536\n");
1080 return nullptr;
1081 }
1082 if (id < 0) {
1083 printf("\"id\" must be positive\n");
1084 return nullptr;
1085 }
1086
1087 mts = init_mt_search(&ck, &pre, w, p);
1088 if (mts == nullptr) return nullptr;
1089
1090 if ( NOT_FOUND == get_irred_param(&ck, &pre, &global_mt19937,
1091 mts, id, DEFAULT_ID_SIZE) ) {
1092 free_mt_struct(mts);
1093 return nullptr;
1094 }
1095 _get_tempering_parameter_hard_dc(mts);
1096 end_mt_search(&pre);
1097
1098 return mts;
1099 }
1100
get_mt_parameters_st(int w,int p,int start_id,int max_id,uint32_t seed,int * count)1101 mt_struct **get_mt_parameters_st(int w, int p, int start_id,
1102 int max_id, uint32_t seed, int *count)
1103 {
1104 mt_struct **mtss, *template_mts;
1105 int i;
1106 prescr_t pre;
1107 _org_state org;
1108 check32_t ck;
1109
1110 if ((start_id > max_id) || (max_id > 0xffff) || (start_id < 0)) {
1111 printf("\"id\" error\n");
1112 return nullptr;
1113 }
1114
1115 _sgenrand_dc(&org, seed);
1116 mtss = (mt_struct**)malloc(sizeof(mt_struct*)*(max_id-start_id+1));
1117 if (nullptr == mtss) return nullptr;
1118
1119 template_mts = init_mt_search(&ck, &pre, w, p);
1120 if (template_mts == nullptr) {
1121 free(mtss);
1122 return nullptr;
1123 }
1124 *count = 0;
1125 for (i=0; i<=max_id-start_id; i++) {
1126 mtss[i] = alloc_mt_struct(template_mts->nn);
1127 if (nullptr == mtss[i]) {
1128 break;
1129 }
1130
1131 copy_params_of_mt_struct(template_mts, mtss[i]);
1132
1133 if ( NOT_FOUND == get_irred_param(&ck, &pre, &org, mtss[i],
1134 i+start_id,DEFAULT_ID_SIZE) ) {
1135 free_mt_struct(mtss[i]);
1136 break;
1137 }
1138 _get_tempering_parameter_hard_dc(mtss[i]);
1139 ++(*count);
1140 }
1141
1142 free_mt_struct(template_mts);
1143 end_mt_search(&pre);
1144 if (*count > 0) {
1145 return mtss;
1146 } else {
1147 free(mtss);
1148 return nullptr;
1149 }
1150 }
1151
get_mt_parameters(int w,int p,int max_id,int * count)1152 mt_struct **get_mt_parameters(int w, int p, int max_id, int *count)
1153 {
1154 mt_struct **mtss, *template_mts;
1155 int i;
1156 prescr_t pre;
1157 check32_t ck;
1158 int start_id = 0;
1159
1160 if ((start_id > max_id) || (max_id > 0xffff) || (start_id < 0)) {
1161 printf("\"id\" error\n");
1162 return nullptr;
1163 }
1164
1165 mtss = (mt_struct**)malloc(sizeof(mt_struct*)*(max_id-start_id+1));
1166 if (nullptr == mtss) return nullptr;
1167
1168 template_mts = init_mt_search(&ck, &pre, w, p);
1169 if (template_mts == nullptr) {
1170 free(mtss);
1171 return nullptr;
1172 }
1173 *count = 0;
1174 for (i=0; i<=max_id-start_id; i++) {
1175 mtss[i] = alloc_mt_struct(template_mts->nn);
1176 if (nullptr == mtss[i]) {
1177 break;
1178 }
1179
1180 copy_params_of_mt_struct(template_mts, mtss[i]);
1181
1182 if ( NOT_FOUND == get_irred_param(&ck, &pre, &global_mt19937, mtss[i],
1183 i+start_id,DEFAULT_ID_SIZE) ) {
1184 free_mt_struct(mtss[i]);
1185 break;
1186 }
1187 _get_tempering_parameter_hard_dc(mtss[i]);
1188 ++(*count);
1189 }
1190
1191 free_mt_struct(template_mts);
1192 end_mt_search(&pre);
1193 if (*count > 0) {
1194 return mtss;
1195 } else {
1196 free(mtss);
1197 return nullptr;
1198 }
1199 }
1200
1201 /* n : sizeof state vector */
alloc_mt_struct(int n)1202 static mt_struct *alloc_mt_struct(int n)
1203 {
1204 mt_struct *mts;
1205
1206 mts = (mt_struct*)malloc(sizeof(mt_struct));
1207 if (nullptr == mts) return nullptr;
1208 mts->state = (uint32_t*)malloc(n*sizeof(uint32_t));
1209 if (nullptr == mts->state) {
1210 free(mts);
1211 return nullptr;
1212 }
1213
1214 return mts;
1215 }
1216
free_mt_struct(mt_struct * mts)1217 void free_mt_struct(mt_struct *mts)
1218 {
1219 free(mts->state);
1220 free(mts);
1221 }
1222
free_mt_struct_array(mt_struct ** mtss,int count)1223 void free_mt_struct_array(mt_struct **mtss, int count)
1224 {
1225 int i;
1226
1227 if (mtss == nullptr) {
1228 return;
1229 }
1230 for (i=0; i < count; i++) {
1231 free_mt_struct(mtss[i]);
1232 }
1233 free(mtss);
1234 }
1235
copy_params_of_mt_struct(mt_struct * src,mt_struct * dst)1236 static void copy_params_of_mt_struct(mt_struct *src, mt_struct *dst)
1237 {
1238 dst->nn = src->nn;
1239 dst->mm = src->mm;
1240 dst->rr = src->rr;
1241 dst->ww = src->ww;
1242 dst->wmask = src->wmask;
1243 dst->umask = src->umask;
1244 dst->lmask = src->lmask;
1245 }
1246
proper_mersenne_exponent(int p)1247 static int proper_mersenne_exponent(int p)
1248 {
1249 switch(p) {
1250 case 521:
1251 case 607:
1252 case 1279:
1253 case 2203:
1254 case 2281:
1255 case 3217:
1256 case 4253:
1257 case 4423:
1258 case 9689:
1259 case 9941:
1260 case 11213:
1261 case 19937:
1262 case 21701:
1263 case 23209:
1264 case 44497:
1265 return 1;
1266 default:
1267 return 0;
1268 }
1269 }
1270
1271 // prescr.c
1272
1273 #define LIMIT_IRRED_DEG 31
1274 #define NIRREDPOLY 127
1275 #define MAX_IRRED_DEG 9
1276
1277 /* list of irreducible polynomials whose degrees are less than 10 */
1278 static const int irredpolylist[NIRREDPOLY][MAX_IRRED_DEG+1] = {
1279 {0,1,0,0,0,0,0,0,0,0,},{1,1,0,0,0,0,0,0,0,0,},{1,1,1,0,0,0,0,0,0,0,},
1280 {1,1,0,1,0,0,0,0,0,0,},{1,0,1,1,0,0,0,0,0,0,},{1,1,0,0,1,0,0,0,0,0,},
1281 {1,0,0,1,1,0,0,0,0,0,},{1,1,1,1,1,0,0,0,0,0,},{1,0,1,0,0,1,0,0,0,0,},
1282 {1,0,0,1,0,1,0,0,0,0,},{1,1,1,1,0,1,0,0,0,0,},{1,1,1,0,1,1,0,0,0,0,},
1283 {1,1,0,1,1,1,0,0,0,0,},{1,0,1,1,1,1,0,0,0,0,},{1,1,0,0,0,0,1,0,0,0,},
1284 {1,0,0,1,0,0,1,0,0,0,},{1,1,1,0,1,0,1,0,0,0,},{1,1,0,1,1,0,1,0,0,0,},
1285 {1,0,0,0,0,1,1,0,0,0,},{1,1,1,0,0,1,1,0,0,0,},{1,0,1,1,0,1,1,0,0,0,},
1286 {1,1,0,0,1,1,1,0,0,0,},{1,0,1,0,1,1,1,0,0,0,},{1,1,0,0,0,0,0,1,0,0,},
1287 {1,0,0,1,0,0,0,1,0,0,},{1,1,1,1,0,0,0,1,0,0,},{1,0,0,0,1,0,0,1,0,0,},
1288 {1,0,1,1,1,0,0,1,0,0,},{1,1,1,0,0,1,0,1,0,0,},{1,1,0,1,0,1,0,1,0,0,},
1289 {1,0,0,1,1,1,0,1,0,0,},{1,1,1,1,1,1,0,1,0,0,},{1,0,0,0,0,0,1,1,0,0,},
1290 {1,1,0,1,0,0,1,1,0,0,},{1,1,0,0,1,0,1,1,0,0,},{1,0,1,0,1,0,1,1,0,0,},
1291 {1,0,1,0,0,1,1,1,0,0,},{1,1,1,1,0,1,1,1,0,0,},{1,0,0,0,1,1,1,1,0,0,},
1292 {1,1,1,0,1,1,1,1,0,0,},{1,0,1,1,1,1,1,1,0,0,},{1,1,0,1,1,0,0,0,1,0,},
1293 {1,0,1,1,1,0,0,0,1,0,},{1,1,0,1,0,1,0,0,1,0,},{1,0,1,1,0,1,0,0,1,0,},
1294 {1,0,0,1,1,1,0,0,1,0,},{1,1,1,1,1,1,0,0,1,0,},{1,0,1,1,0,0,1,0,1,0,},
1295 {1,1,1,1,1,0,1,0,1,0,},{1,1,0,0,0,1,1,0,1,0,},{1,0,1,0,0,1,1,0,1,0,},
1296 {1,0,0,1,0,1,1,0,1,0,},{1,0,0,0,1,1,1,0,1,0,},{1,1,1,0,1,1,1,0,1,0,},
1297 {1,1,0,1,1,1,1,0,1,0,},{1,1,1,0,0,0,0,1,1,0,},{1,1,0,1,0,0,0,1,1,0,},
1298 {1,0,1,1,0,0,0,1,1,0,},{1,1,1,1,1,0,0,1,1,0,},{1,1,0,0,0,1,0,1,1,0,},
1299 {1,0,0,1,0,1,0,1,1,0,},{1,0,0,0,1,1,0,1,1,0,},{1,0,1,1,1,1,0,1,1,0,},
1300 {1,1,0,0,0,0,1,1,1,0,},{1,1,1,1,0,0,1,1,1,0,},{1,1,1,0,1,0,1,1,1,0,},
1301 {1,0,1,1,1,0,1,1,1,0,},{1,1,1,0,0,1,1,1,1,0,},{1,1,0,0,1,1,1,1,1,0,},
1302 {1,0,1,0,1,1,1,1,1,0,},{1,0,0,1,1,1,1,1,1,0,},{1,1,0,0,0,0,0,0,0,1,},
1303 {1,0,0,0,1,0,0,0,0,1,},{1,1,1,0,1,0,0,0,0,1,},{1,1,0,1,1,0,0,0,0,1,},
1304 {1,0,0,0,0,1,0,0,0,1,},{1,0,1,1,0,1,0,0,0,1,},{1,1,0,0,1,1,0,0,0,1,},
1305 {1,1,0,1,0,0,1,0,0,1,},{1,0,0,1,1,0,1,0,0,1,},{1,1,1,1,1,0,1,0,0,1,},
1306 {1,0,1,0,0,1,1,0,0,1,},{1,0,0,1,0,1,1,0,0,1,},{1,1,1,1,0,1,1,0,0,1,},
1307 {1,1,1,0,1,1,1,0,0,1,},{1,0,1,1,1,1,1,0,0,1,},{1,1,1,0,0,0,0,1,0,1,},
1308 {1,0,1,0,1,0,0,1,0,1,},{1,0,0,1,1,0,0,1,0,1,},{1,1,0,0,0,1,0,1,0,1,},
1309 {1,0,1,0,0,1,0,1,0,1,},{1,1,1,1,0,1,0,1,0,1,},{1,1,1,0,1,1,0,1,0,1,},
1310 {1,0,1,1,1,1,0,1,0,1,},{1,1,1,1,0,0,1,1,0,1,},{1,0,0,0,1,0,1,1,0,1,},
1311 {1,1,0,1,1,0,1,1,0,1,},{1,0,1,0,1,1,1,1,0,1,},{1,0,0,1,1,1,1,1,0,1,},
1312 {1,0,0,0,0,0,0,0,1,1,},{1,1,0,0,1,0,0,0,1,1,},{1,0,1,0,1,0,0,0,1,1,},
1313 {1,1,1,1,1,0,0,0,1,1,},{1,1,0,0,0,1,0,0,1,1,},{1,0,0,0,1,1,0,0,1,1,},
1314 {1,1,0,1,1,1,0,0,1,1,},{1,0,0,1,0,0,1,0,1,1,},{1,1,1,1,0,0,1,0,1,1,},
1315 {1,1,0,1,1,0,1,0,1,1,},{1,0,0,0,0,1,1,0,1,1,},{1,1,0,1,0,1,1,0,1,1,},
1316 {1,0,1,1,0,1,1,0,1,1,},{1,1,0,0,1,1,1,0,1,1,},{1,1,1,1,1,1,1,0,1,1,},
1317 {1,0,1,0,0,0,0,1,1,1,},{1,1,1,1,0,0,0,1,1,1,},{1,0,0,0,0,1,0,1,1,1,},
1318 {1,0,1,0,1,1,0,1,1,1,},{1,0,0,1,1,1,0,1,1,1,},{1,1,1,0,0,0,1,1,1,1,},
1319 {1,1,0,1,0,0,1,1,1,1,},{1,0,1,1,0,0,1,1,1,1,},{1,0,1,0,1,0,1,1,1,1,},
1320 {1,0,0,1,1,0,1,1,1,1,},{1,1,0,0,0,1,1,1,1,1,},{1,0,0,1,0,1,1,1,1,1,},
1321 {1,1,0,1,1,1,1,1,1,1,},
1322 };
1323
1324 static void MakepreModPolys(prescr_t *pre, int mm, int nn, int rr, int ww);
1325 static Polynomial *make_tntm( int n, int m);
1326 static Polynomial *PolynomialDup(Polynomial *pl);
1327 static void PolynomialMod(Polynomial *wara, const Polynomial *waru);
1328 static Polynomial *PolynomialMult(Polynomial *p0, Polynomial *p1);
1329 static void FreePoly( Polynomial *p);
1330 static Polynomial *NewPoly(int degree);
1331 static int IsReducible(prescr_t *pre, uint32_t aaa, uint32_t *polylist);
1332 static uint32_t word2bit(Polynomial *pl);
1333 static void makemodlist(prescr_t *pre, Polynomial *pl, int nPoly);
1334 static void NextIrredPoly(Polynomial *pl, int nth);
1335
1336
1337 /*************************************************/
1338 /*************************************************/
_prescreening_dc(prescr_t * pre,uint32_t aaa)1339 int _prescreening_dc(prescr_t *pre, uint32_t aaa)
1340 {
1341
1342 int i;
1343
1344 for (i=0; i<NIRREDPOLY; i++) {
1345 if (IsReducible(pre, aaa,pre->modlist[i])==REDU)
1346 return REJECTED;
1347 }
1348 return NOT_REJECTED;
1349 }
1350
_InitPrescreening_dc(prescr_t * pre,int m,int n,int r,int w)1351 void _InitPrescreening_dc(prescr_t *pre, int m, int n, int r, int w)
1352 {
1353 int i;
1354 Polynomial *pl;
1355
1356 pre->sizeofA = w;
1357
1358 pre->preModPolys = (Polynomial **)malloc(
1359 (pre->sizeofA+1)*(sizeof(Polynomial*)));
1360 if (nullptr == pre->preModPolys) {
1361 printf ("malloc error in \"InitPrescreening\"\n");
1362 exit(1);
1363 }
1364 MakepreModPolys(pre, m,n,r,w);
1365
1366 pre->modlist = (uint32_t**)malloc(NIRREDPOLY * sizeof(uint32_t*));
1367 if (nullptr == pre->modlist) {
1368 printf ("malloc error in \"InitPrescreening()\"\n");
1369 exit(1);
1370 }
1371 for (i=0; i<NIRREDPOLY; i++) {
1372 pre->modlist[i]
1373 = (uint32_t*)malloc( (pre->sizeofA + 1) * (sizeof(uint32_t)) );
1374 if (nullptr == pre->modlist[i]) {
1375 printf ("malloc error in \"InitPrescreening()\"\n");
1376 exit(1);
1377 }
1378 }
1379
1380
1381 for (i=0; i<NIRREDPOLY; i++) {
1382 pl = NewPoly(MAX_IRRED_DEG);
1383 NextIrredPoly(pl,i);
1384 makemodlist(pre, pl, i);
1385 FreePoly(pl);
1386 }
1387
1388 for (i=pre->sizeofA; i>=0; i--)
1389 FreePoly(pre->preModPolys[i]);
1390 free(pre->preModPolys);
1391
1392 }
1393
_EndPrescreening_dc(prescr_t * pre)1394 void _EndPrescreening_dc(prescr_t *pre)
1395 {
1396 int i;
1397
1398 for (i=0; i<NIRREDPOLY; i++)
1399 free(pre->modlist[i]);
1400 free(pre->modlist);
1401 }
1402
1403 /*************************************************/
1404 /****** static functions ******/
1405 /*************************************************/
1406
NextIrredPoly(Polynomial * pl,int nth)1407 void NextIrredPoly(Polynomial *pl, int nth)
1408 {
1409 int i, max_deg;
1410
1411 for (max_deg=0,i=0; i<=MAX_IRRED_DEG; i++) {
1412 if ( irredpolylist[nth][i] )
1413 max_deg = i;
1414 pl->x[i] = irredpolylist[nth][i];
1415 }
1416
1417 pl->deg = max_deg;
1418
1419 }
1420
makemodlist(prescr_t * pre,Polynomial * pl,int nPoly)1421 static void makemodlist(prescr_t *pre, Polynomial *pl, int nPoly)
1422 {
1423 Polynomial *tmpPl;
1424 int i;
1425
1426 for (i=0; i<=pre->sizeofA; i++) {
1427 tmpPl = PolynomialDup(pre->preModPolys[i]);
1428 PolynomialMod(tmpPl,pl);
1429 pre->modlist[nPoly][i] = word2bit(tmpPl);
1430 FreePoly(tmpPl);
1431 }
1432 }
1433
1434 /* Pack Polynomial into a word */
word2bit(Polynomial * pl)1435 static uint32_t word2bit(Polynomial *pl)
1436 {
1437 int i;
1438 uint32_t bx;
1439
1440 bx = 0;
1441 for (i=pl->deg; i>0; i--) {
1442 if (pl->x[i]) bx |= 0x1;
1443 bx <<= 1;
1444 }
1445 if (pl->x[0]) bx |= 0x1;
1446
1447 return bx;
1448 }
1449
1450 /* REDU -- reducible */
1451 /* aaa = (a_{w-1}a_{w-2}...a_1a_0 */
IsReducible(prescr_t * pre,uint32_t aaa,uint32_t * polylist)1452 static int IsReducible(prescr_t *pre, uint32_t aaa, uint32_t *polylist)
1453 {
1454 int i;
1455 uint32_t x;
1456
1457 x = polylist[pre->sizeofA];
1458 for (i=pre->sizeofA-1; i>=0; i--) {
1459 if (aaa&0x1)
1460 x ^= polylist[i];
1461 aaa >>= 1;
1462 }
1463
1464 if ( x == 0 ) return REDU;
1465 else return NONREDU;
1466 }
1467
1468
1469 /***********************************/
1470 /** functions for polynomial **/
1471 /***********************************/
NewPoly(int degree)1472 static Polynomial *NewPoly(int degree)
1473 {
1474 Polynomial *p;
1475
1476 p = (Polynomial *)calloc( 1, sizeof(Polynomial));
1477 if( p==nullptr ){
1478 printf("calloc error in \"NewPoly()\"\n");
1479 exit(1);
1480 }
1481 p->deg = degree;
1482
1483 if (degree < 0) {
1484 p->x = nullptr;
1485 return p;
1486 }
1487
1488 p->x = (int *)calloc( degree + 1, sizeof(int));
1489 if( p->x == nullptr ){
1490 printf("calloc error\n");
1491 exit(1);
1492 }
1493
1494 return p;
1495 }
1496
FreePoly(Polynomial * p)1497 static void FreePoly( Polynomial *p)
1498 {
1499 if (p->x != nullptr)
1500 free( p->x );
1501 free( p );
1502 }
1503
1504
1505 /** multiplication **/
PolynomialMult(Polynomial * p0,Polynomial * p1)1506 static Polynomial *PolynomialMult(Polynomial *p0,Polynomial *p1)
1507 {
1508 int i, j;
1509 Polynomial *p;
1510
1511 /* if either p0 or p1 is 0, return 0 */
1512 if ( (p0->deg < 0) || (p1->deg < 0) ) {
1513 p = NewPoly(-1);
1514 return p;
1515 }
1516
1517 p = NewPoly(p0->deg + p1->deg);
1518 for( i=0; i<=p1->deg; i++){
1519 if( p1->x[i] ){
1520 for( j=0; j<=p0->deg; j++){
1521 p->x[i+j] ^= p0->x[j];
1522 }
1523 }
1524 }
1525
1526 return p;
1527 }
1528
1529 /** wara mod waru **/
1530 /** the result is stored in wara ********/
PolynomialMod(Polynomial * wara,const Polynomial * waru)1531 static void PolynomialMod( Polynomial *wara, const Polynomial *waru)
1532 {
1533 int i;
1534 int deg_diff;
1535
1536 while( wara->deg >= waru->deg ){
1537 deg_diff = wara->deg - waru->deg;
1538 for( i=0; i<=waru->deg; i++){
1539 wara->x[ i+deg_diff ] ^= waru->x[i];
1540 }
1541
1542 for( i=wara->deg; i>=0; i--){
1543 if( wara->x[i] ) break;
1544 }
1545 wara->deg=i;
1546
1547 }
1548 }
1549
PolynomialDup(Polynomial * pl)1550 static Polynomial *PolynomialDup(Polynomial *pl)
1551 {
1552 Polynomial *pt;
1553 int i;
1554
1555 pt = NewPoly(pl->deg);
1556 for (i=pl->deg; i>=0; i--)
1557 pt->x[i] = pl->x[i];
1558
1559 return pt;
1560 }
1561
1562 /** make the polynomial "t**n + t**m" **/
make_tntm(int n,int m)1563 static Polynomial *make_tntm( int n, int m)
1564 {
1565 Polynomial *p;
1566
1567 p = NewPoly(n);
1568 p->x[n] = p->x[m] = 1;
1569
1570 return p;
1571 }
1572
MakepreModPolys(prescr_t * pre,int mm,int nn,int rr,int ww)1573 static void MakepreModPolys(prescr_t *pre, int mm, int nn, int rr, int ww)
1574 {
1575 Polynomial *t, *t0, *t1, *s, *s0, *s1;
1576 int i,j;
1577
1578 j = 0;
1579 t = NewPoly(0);
1580 t->deg = 0;
1581 t->x[0] = 1;
1582 pre->preModPolys[j++] = t;
1583
1584 t = make_tntm (nn, mm);
1585 t0 = make_tntm (nn, mm);
1586 s = make_tntm (nn-1, mm-1);
1587
1588 for( i=1; i<(ww - rr); i++){
1589 pre->preModPolys[j++] = PolynomialDup(t0);
1590 t1 = t0;
1591 t0 = PolynomialMult(t0, t);
1592 FreePoly(t1);
1593 }
1594
1595 pre->preModPolys[j++] = PolynomialDup(t0);
1596
1597 s0 =PolynomialMult( t0, s);
1598 FreePoly(t0); FreePoly(t);
1599 for( i=(rr-2); i>=0; i--){
1600 pre->preModPolys[j++] = PolynomialDup(s0);
1601 s1 = s0;
1602 s0 = PolynomialMult( s0, s);
1603 FreePoly(s1);
1604 }
1605
1606 pre->preModPolys[j++] = PolynomialDup(s0);
1607
1608 FreePoly(s0); FreePoly(s);
1609 }
1610
1611 // init.c
1612
init_dc(uint32_t seed)1613 void init_dc(uint32_t seed)
1614 {
1615 _sgenrand_dc(&global_mt19937, seed);
1616 }
1617
1618 // genmtrand.c
1619
1620 #define SHIFT1 18
1621
sgenrand_mt(uint32_t seed,mt_struct * mts)1622 void sgenrand_mt(uint32_t seed, mt_struct *mts)
1623 {
1624 int i;
1625
1626 for (i=0; i<mts->nn; i++) {
1627 mts->state[i] = seed;
1628 seed = (UINT32_C(1812433253) * (seed ^ (seed >> 30))) + i + 1;
1629 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
1630 /* In the previous versions, MSBs of the seed affect */
1631 /* only MSBs of the array mt[]. */
1632 }
1633 mts->i = mts->nn;
1634
1635 for (i=0; i<mts->nn; i++)
1636 mts->state[i] &= mts->wmask;
1637 }
1638
genrand_mt(mt_struct * mts)1639 uint32_t genrand_mt(mt_struct *mts)
1640 {
1641 uint32_t *st, uuu, lll, aa, x;
1642 int k,n,m,lim;
1643
1644 if ( mts->i >= mts->nn ) {
1645 n = mts->nn; m = mts->mm;
1646 aa = mts->aaa;
1647 st = mts->state;
1648 uuu = mts->umask; lll = mts->lmask;
1649
1650 lim = n - m;
1651 for (k=0; k<lim; k++) {
1652 x = (st[k]&uuu)|(st[k+1]&lll);
1653 st[k] = st[k+m] ^ (x>>1) ^ (x&1U ? aa : 0U);
1654 }
1655 lim = n - 1;
1656 for (; k<lim; k++) {
1657 x = (st[k]&uuu)|(st[k+1]&lll);
1658 st[k] = st[k+m-n] ^ (x>>1) ^ (x&1U ? aa : 0U);
1659 }
1660 x = (st[n-1]&uuu)|(st[0]&lll);
1661 st[n-1] = st[m-1] ^ (x>>1) ^ (x&1U ? aa : 0U);
1662 mts->i=0;
1663 }
1664
1665 x = mts->state[mts->i];
1666 mts->i += 1;
1667 x ^= x >> mts->shift0;
1668 x ^= (x << mts->shiftB) & mts->maskB;
1669 x ^= (x << mts->shiftC) & mts->maskC;
1670 x ^= x >> mts->shift1;
1671
1672 return x;
1673 }
1674
1675 // check32.c
1676
1677 #define LSB 0x1
1678 #define WORDLEN 32
1679
_InitCheck32_dc(check32_t * ck,int r,int w)1680 void _InitCheck32_dc(check32_t *ck, int r, int w)
1681 {
1682 int i;
1683
1684 /* word_mask (least significant w bits) */
1685 ck->word_mask = 0xFFFFFFFF;
1686 ck->word_mask <<= WORDLEN - w;
1687 ck->word_mask >>= WORDLEN - w;
1688 /* lower_mask (least significant r bits) */
1689 for (ck->lower_mask=0,i=0; i<r; ++i) {
1690 ck->lower_mask <<= 1;
1691 ck->lower_mask |= LSB;
1692 }
1693 /* upper_mask (most significant (w-r) bits */
1694 ck->upper_mask = (~ck->lower_mask) & ck->word_mask;
1695 }
1696
_CheckPeriod_dc(check32_t * ck,_org_state * st,uint32_t a,int m,int n,int r,int w)1697 int _CheckPeriod_dc(check32_t *ck, _org_state *st,
1698 uint32_t a, int m, int n, int r, int w)
1699 {
1700 int i, j, p, pp;
1701 uint32_t y, *x, *init, mat[2];
1702
1703
1704 p = n*w-r;
1705 x = (uint32_t*) malloc (2*p*sizeof(uint32_t));
1706 if (nullptr==x) {
1707 printf("malloc error in \"_CheckPeriod_dc()\"\n");
1708 exit(1);
1709 }
1710
1711 init = (uint32_t*) malloc (n*sizeof(uint32_t));
1712 if (nullptr==init) {
1713 printf("malloc error \"_CheckPeriod_dc()\"\n");
1714 free(x);
1715 exit(1);
1716 }
1717
1718 /* set initial values */
1719 for (i=0; i<n; ++i)
1720 x[i] = init[i] = (ck->word_mask & _genrand_dc(st));
1721 /* it is better that LSBs of x[2] and x[3] are different */
1722 if ( (x[2]&LSB) == (x[3]&LSB) ) {
1723 x[3] ^= 1;
1724 init[3] ^= 1;
1725 }
1726
1727 pp = 2*p-n;
1728 mat[0] = 0; mat[1] = a;
1729 for (j=0; j<p; ++j) {
1730
1731 /* generate */
1732 for (i=0; i<pp; ++i){
1733 y = (x[i]&ck->upper_mask) | (x[i+1]&ck->lower_mask);
1734 x[i+n] = x[i+m] ^ ( (y>>1) ^ mat[y&LSB] );
1735 }
1736
1737 /* pick up odd subscript elements */
1738 for (i=2; i<=p; ++i)
1739 x[i] = x[(i<<1)-1];
1740
1741 /* reverse generate */
1742 for (i=p-n; i>=0; --i) {
1743 y = x[i+n] ^ x[i+m] ^ mat[ x[i+1]&LSB ];
1744 y <<=1; y |= x[i+1]&LSB;
1745
1746 x[i+1] = (x[i+1]&ck->upper_mask) | (y&ck->lower_mask);
1747 x[i] = (y&ck->upper_mask) | (x[i]&ck->lower_mask);
1748 }
1749
1750 }
1751
1752 if ((x[0]&ck->upper_mask)==(init[0]&ck->upper_mask)) {
1753 for (i=1; i<n; ++i) {
1754 if (x[i] != init[i])
1755 break;
1756 }
1757 if (i==n) {
1758 free(x); free(init);
1759 return IRRED;
1760 }
1761 }
1762
1763
1764 free(x); free(init);
1765 return REDU;
1766 }
1767
1768 #undef N
1769 #undef NOT_REJECTED
1770 #undef REJECTED
1771 #undef REDU
1772 #undef IRRED
1773 #undef NONREDU
1774 #undef SSS
1775 #undef TTT
1776 #undef S00
1777 #undef S01
1778 #undef LIMIT_V_BEST_OPT
1779 #undef WORD_LEN
1780 #undef MIN_INFINITE
1781 #undef M
1782 #undef MATRIX_A
1783 #undef UPPER_MASK
1784 #undef LOWER_MASK
1785 #undef TEMPERING_MASK_B
1786 #undef TEMPERING_MASK_C
1787 #undef TEMPERING_SHIFT_U
1788 #undef TEMPERING_SHIFT_S
1789 #undef TEMPERING_SHIFT_T
1790 #undef TEMPERING_SHIFT_L
1791 #undef LSB
1792 #undef MAX_SEARCH
1793 #undef FOUND
1794 #undef NOT_FOUND
1795 #undef DEFAULT_ID_SIZE
1796 #undef LIMIT_IRRED_DEG
1797 #undef NIRREDPOLY
1798 #undef MAX_IRRED_DEG
1799 #undef SHIFT1
1800 #undef LSB
1801 #undef WORDLEN
1802