1 /*************************************************************************/
2 /*************************************************************************/
3 /* Parallel 48 bit Linear Congruential Generator */
4 /* */
5 /* Author: Ashok Srinivasan, */
6 /* NCSA, University of Illinois, Urbana-Champaign */
7 /* E-Mail: ashoks@ncsa.uiuc.edu */
8 /* */
9 /* Disclaimer: NCSA expressly disclaims any and all warranties, expressed*/
10 /* or implied, concerning the enclosed software. The intent in sharing */
11 /* this software is to promote the productive interchange of ideas */
12 /* throughout the research community. All software is furnished on an */
13 /* "as is" basis. No further updates to this software should be */
14 /* expected. Although this may occur, no commitment exists. The authors */
15 /* certainly invite your comments as well as the reporting of any bugs. */
16 /* NCSA cannot commit that any or all bugs will be fixed. */
17 /*************************************************************************/
18 /*************************************************************************/
19
20 /* This is version 0.2, created 13 April 1998 */
21
22 #include <stdio.h>
23 #include <string.h>
24 #include <stdlib.h>
25 #include "memory.h"
26 #include "primes-lcg.h"
27 #include "interface.h"
28 #include <limits.h>
29 #define NDEBUG
30 #include <assert.h>
31 #include "store.h"
32
33 #define VERSION "00"
34 #define GENTYPE VERSION "48 bit Linear Congruential Generator with Prime Addend"
35
36 #if LONG_MAX > 2147483647L /* 32 bit integer */
37 #if LONG_MAX > 35184372088831L /* 46 bit integer */
38 #if LONG_MAX >= 9223372036854775807L /* 64 bit integer */
39 #define LONG_SPRNG
40 #define LONG64 long
41 #define store_long64 store_long
42 #define store_long64array store_longarray
43 #define load_long64 load_long
44 #define load_long64array load_longarray
45
46 #define INIT_SEED 0x2bc68cfe166dL
47 #define MSB_SET 0x3ff0000000000000L
48 #define LSB48 0xffffffffffffL
49 #define AN1 0xdadf0ac00001L
50 #define AN2 0xfefd7a400001L
51 #define AN3 0x6417b5c00001L
52 #define AN4 0xcf9f72c00001L
53 #define AN5 0xbdf07b400001L
54 #define AN6 0xf33747c00001L
55 #define AN7 0xcbe632c00001L
56 #define PMULT1 0xa42c22700000L
57 #define PMULT2 0xfa858cb00000L
58 #define PMULT3 0xd0c4ef00000L
59 #define PMULT4 0xc3cc8e300000L
60 #define PMULT5 0x11bdbe700000L
61 #define PMULT6 0xb0f0e9f00000L
62 #define PMULT7 0x6407de700000L
63 #define MULT1 0x2875a2e7b175L /* 44485709377909 */
64 #define MULT2 0x5deece66dL /* 1575931494 */
65 #define MULT3 0x3eac44605265L /* 68909602460261 */
66 #define MULT4 0x275b38eb4bbdL /* 4327250451645 */
67 #define MULT5 0x1ee1429cc9f5L /* 33952834046453 */
68 #define MULT6 0x739a9cb08605L /* 127107890972165 */
69 #define MULT7 0x3228d7cc25f5L /* 55151000561141 */
70 #endif
71 #endif
72 #endif
73
74 #if !defined(LONG_SPRNG) && defined(_LONG_LONG)
75 #define LONG64 long long
76 #define store_long64 store_longlong
77 #define store_long64array store_longlongarray
78 #define load_long64 load_longlong
79 #define load_long64array load_longlongarray
80
81 #define INIT_SEED 0x2bc68cfe166dLL
82 #define MSB_SET 0x3ff0000000000000LL
83 #define LSB48 0xffffffffffffLL
84 #define AN1 0xdadf0ac00001LL
85 #define AN2 0xfefd7a400001LL
86 #define AN3 0x6417b5c00001LL
87 #define AN4 0xcf9f72c00001LL
88 #define AN5 0xbdf07b400001LL
89 #define AN6 0xf33747c00001LL
90 #define AN7 0xcbe632c00001LL
91 #define PMULT1 0xa42c22700000LL
92 #define PMULT2 0xfa858cb00000LL
93 #define PMULT3 0xd0c4ef00000LL
94 #define PMULT4 0x11bdbe700000LL
95 #define PMULT5 0xc3cc8e300000LL
96 #define PMULT6 0xb0f0e9f00000LL
97 #define PMULT7 0x6407de700000LL
98 #define MULT1 0x2875a2e7b175LL
99 #define MULT2 0x5deece66dLL
100 #define MULT3 0x3eac44605265LL
101 #define MULT4 0x1ee1429cc9f5LL
102 #define MULT5 0x275b38eb4bbdLL
103 #define MULT6 0x739a9cb08605LL
104 #define MULT7 0x3228d7cc25f5LL
105 #endif
106
107 #define TWO_M24 5.96046447753906234e-8
108 #define TWO_M48 3.5527136788005008323e-15
109
110 #include "multiply.h"
111
112 #ifdef LittleEndian
113 #define MSB 1
114 #else
115 #define MSB 0
116 #endif
117 #define LSB (1-MSB)
118
119 #define LCGRUNUP 29
120
121 int MAX_STREAMS = 1<<19;
122
123 #ifndef TOOMANY
124 #define TOOMANY "generator has branched maximum number of times;\nindependence of streams cannot be guranteed\n"
125 #endif
126
127 #ifdef LONG64
128 struct rngen
129 {
130 unsigned LONG64 seed;
131 int init_seed;
132 int prime;
133 int prime_position;
134 int prime_next;
135 char *gentype;
136 int parameter;
137 unsigned LONG64 multiplier;
138 };
139
140 unsigned LONG64 mults[] = {MULT1,MULT2,MULT3,MULT4,MULT5,MULT6,MULT7};
141 unsigned LONG64 multiplier=0;
142
143 #else
144 struct rngen
145 {
146 int seed[2];
147 int init_seed;
148 int prime;
149 int prime_position;
150 int prime_next;
151 char *gentype;
152 int parameter;
153 int *multiplier;
154 };
155
156 int mults[][4] = {{0x175,0xe7b,0x5a2,0x287},{0x66d,0xece,0x5de,0x000},
157 {0x265,0x605,0xc44,0x3ea},{0x9f5,0x9cc,0x142,0x1ee},
158 {0xbbd,0xeb4,0xb38,0x275},{0x605,0xb08,0xa9c,0x739},
159 {0x5f5,0xcc2,0x8d7,0x322}};
160 int *multiplier=NULL;
161 #endif
162
163 #define NPARAMS 7
164
165 int NGENS=0;
166
167 void plus ANSI_ARGS((int *a, int *b, int *c));
168 void mult ANSI_ARGS((int *a, int *b, int *c, int size));
169 void advance_seed ANSI_ARGS((struct rngen *gen));
170 double get_rn_dbl ANSI_ARGS((int *gen));
171
172 #ifdef __STDC__
bit_reverse(int n)173 int bit_reverse(int n)
174 #else
175 int bit_reverse(n)
176 int n;
177 #endif
178 {
179 int i=31, rev=0;
180
181 for(i=30; i>=0; i--)
182 {
183 rev |= (n&1)<<i;
184 n >>= 1;
185 }
186
187 return rev;
188 }
189
190
191 #ifdef __STDC__
errprint(char * level,char * routine,char * error)192 void errprint(char *level, char *routine, char *error)
193 #else
194 void errprint(level, routine, error)
195 char *level,*routine,*error;
196 #endif
197 {
198 #ifdef CRAY
199 #pragma _CRI guard 63
200 #endif
201 fprintf(stderr,"%s from %s: %s\n",level,routine,error);
202 #ifdef CRAY
203 #pragma _CRI endguard 63
204 #endif
205 }
206
207
208 #ifdef __STDC__
strxncmp(char * s1,char * s2,int n)209 int strxncmp(char *s1, char *s2, int n)
210 #else
211 int strxncmp(s1, s2, n)
212 char *s1, *s2;
213 int n;
214 #endif
215 {
216 int i;
217
218 for(i=0; i<n; i++)
219 if(s1[i] != s2[i])
220 return s1[i]-s2[i];
221
222 return 0; /* First n characters of strings are equal. */
223 }
224
225
226
227 #ifdef __STDC__
init_rng(int gennum,int total_gen,int seed,int mult)228 int *init_rng( int gennum, int total_gen, int seed, int mult)
229 #else
230 int *init_rng(gennum,total_gen,seed,mult)
231 int gennum,mult,seed,total_gen;
232 #endif
233 {
234 /* gives back one generator (node gennum) with updated spawning */
235 /* info; should be called total_gen times, with different value */
236 /* of gennum in [0,total_gen) each call */
237 struct rngen *genptr;
238 int i;
239
240 if (total_gen <= 0) /* check if total_gen is valid */
241 {
242 total_gen = 1;
243 errprint("WARNING","init_rng","Total_gen <= 0. Default value of 1 used for total_gen");
244 }
245
246 if (gennum >= MAX_STREAMS) /* check if gen_num is valid */
247 fprintf(stderr,"WARNING - init_rng: gennum: %d > maximum number of independent streams: %d\n\tIndependence of streams cannot be guranteed.\n",
248 gennum, MAX_STREAMS);
249
250 if (gennum < 0 || gennum >= total_gen) /* check if gen_num is valid */
251 {
252 errprint("ERROR","init_rng","gennum out of range. ");
253 return (int *) NULL;
254 }
255
256 seed &= 0x7fffffff; /* Only 31 LSB of seed considered */
257
258 if (mult < 0 || mult >= NPARAMS)
259 {
260 errprint("WARNING","init_rng","multiplier not valid. Using Default param");
261 mult = 0;
262 }
263
264 #ifdef LONG64
265 if(multiplier == 0)
266 multiplier = mults[mult];
267 /* else
268 if(multiplier != mults[mult])
269 fprintf(stderr,"WARNING: init_rng_d: Proposed multiplier does not agree with current multiplier.\n\t Independence of streams is not guaranteed\n");*/
270 #else
271 if(multiplier == NULL)
272 multiplier = mults[mult];
273 /*else
274 if(strxncmp((char *) multiplier,(char *) mults[mult],4*sizeof(int)) != 0)
275 fprintf(stderr,"WARNING: init_rng_d: Proposed multiplier does not agree with current multiplier.\n\t Independence of streams is not guaranteed\n");*/
276 #endif
277
278 genptr = (struct rngen *) mymalloc(1*sizeof(struct rngen));
279 if(genptr == NULL)
280 return NULL;
281
282 genptr->gentype = GENTYPE;
283 genptr->init_seed = seed;
284 getprime(1, &(genptr->prime), gennum);
285 genptr->prime_position = gennum;
286 genptr->prime_next = total_gen;
287 genptr->parameter = mult;
288
289 #ifdef LONG64
290 genptr->seed = INIT_SEED; /* initialize generator */
291 genptr->seed ^= ((unsigned LONG64) seed)<<16;
292 genptr->multiplier = mults[mult];
293 if (genptr->prime == 0)
294 genptr->seed |= 1;
295 #else
296 genptr->seed[1] = 16651885^((seed<<16)&(0xff0000));/* initialize generator */
297 genptr->seed[0] = 2868876^((seed>>8)&(0xffffff));
298 genptr->multiplier = mults[mult];
299 if (genptr->prime == 0)
300 genptr->seed[1] |= 1;
301 #endif
302
303 for(i=0; i<LCGRUNUP*(genptr->prime_position); i++)
304 get_rn_dbl( (int *) genptr);
305
306 NGENS++;
307
308 return (int *) genptr;
309 }
310
311
312
313
314 /* On machines with 32 bit integers, */
315 /* the Cray's 48 bit integer math is duplicated by breaking the problem into*/
316 /* steps. The algorithm is linear congruential. M is the multiplier and S*/
317 /* is the current seed. The 31 High order bits out of the 48 bits are
318 returned*/
319 #ifdef __STDC__
get_rn_int(int * igenptr)320 int get_rn_int(int *igenptr)
321 #else
322 int get_rn_int(igenptr)
323 int *igenptr;
324 #endif
325 {
326 struct rngen *genptr = (struct rngen *) igenptr;
327
328 #ifdef LONG64
329 multiply(genptr);
330
331 return ((unsigned LONG64) genptr->seed) >> 17;
332 #else
333 int s[4], res[4];
334 multiply(genptr,genptr->multiplier,s,res);
335
336
337 return (genptr->seed[0]<<7) | ((unsigned int) genptr->seed[1] >> 17) ;
338 #endif
339
340
341 }
342
343 #ifdef __STDC__
get_rn_flt(int * igenptr)344 float get_rn_flt(int *igenptr)
345 #else
346 float get_rn_flt(igenptr)
347 int *igenptr;
348 #endif
349 {
350 return (float) get_rn_dbl(igenptr);
351 } /* get_rn_float */
352
353
354 #ifdef __STDC__
get_rn_dbl(int * igenptr)355 double get_rn_dbl(int *igenptr)
356 #else
357 double get_rn_dbl(igenptr)
358 int *igenptr;
359 #endif
360 {
361 struct rngen *genptr = (struct rngen *) igenptr;
362
363 #ifdef LONG64
364 double temp[1];
365 unsigned LONG64 *ltemp;
366
367 temp[0] = 0.0;
368 multiply(genptr);
369 /* Add defined(O2K) || defined(SGI) if optimization level is -O2
370 or lower */
371 #if defined(CONVEX) || defined(GENERIC)
372 ltemp = (unsigned LONG64 *) temp;
373 *ltemp = (genptr->seed<<4) | MSB_SET;
374
375 return temp[0] - (double) 1.0;
376 #else
377 return genptr->seed*3.5527136788005008e-15;
378 #endif
379
380 #else
381
382 static double equiv[1];
383 #define iran ((int *)equiv)
384 #define ran (equiv)
385
386 int expo, s[4], res[4];
387
388
389 multiply(genptr,genptr->multiplier,s,res);
390
391 #if defined(HP) || defined(SUN) || defined(SOLARIS) || defined(GENERIC)
392 expo = 1072693248;
393
394 /*IEEE mantissa is 52 bits. We have only 48 bits, so we shift our result 4*/
395 /* bits left. 32-(24+4) = 4 bits are still blank in the lower word, so we*/
396 /* grab the low 4 bits of seedhi to fill these. */
397 iran[LSB] = genptr->seed[1] << 4 | genptr->seed[0] << 28;
398
399 /* place the remaining (24-4)=20 bits of seedhi in bits 20-0 of ran. */
400 /* Expo occupies bits 30-20. Bit 31 (sign) is always zero. */
401 iran[MSB] = expo | genptr->seed[0] >> 4;
402
403 return (*ran - (double) 1.0);
404 #else
405 return genptr->seed[0]*TWO_M24 + genptr->seed[1]*TWO_M48;
406 #endif
407
408 #undef ran
409 #undef iran
410 #endif
411 } /* get_rn_dbl */
412
413
414
415 /*************************************************************************/
416 /*************************************************************************/
417 /* SPAWN_RNG: spawns new generators */
418 /*************************************************************************/
419 /*************************************************************************/
420
421 #ifdef __STDC__
spawn_rng(int * igenptr,int nspawned,int *** newgens,int checkid)422 int spawn_rng(int *igenptr, int nspawned, int ***newgens, int checkid)
423 #else
424 int spawn_rng(igenptr,nspawned, newgens, checkid)
425 int *igenptr,nspawned, ***newgens, checkid;
426 #endif
427 {
428 struct rngen **genptr, *tempptr = (struct rngen *) igenptr;
429 int i, j;
430
431 if (nspawned <= 0) /* check if nspawned is valid */
432 {
433 nspawned = 1;
434 errprint("WARNING","spawn_rng","nspawned <= 0. Default value of 1 used for nspawned");
435 }
436
437 genptr = (struct rngen **) mymalloc(nspawned*sizeof(struct rngen *));
438 if(genptr == NULL)
439 {
440 *newgens = NULL;
441 return 0;
442 }
443
444 for(i=0; i<nspawned; i++)
445 {
446 genptr[i] = (struct rngen *) mymalloc(sizeof(struct rngen));
447 if(genptr[i] == NULL)
448 {
449 nspawned = i;
450 break;
451 }
452
453 genptr[i]->init_seed = tempptr->init_seed;
454 genptr[i]->prime_position = tempptr->prime_position +
455 tempptr->prime_next*(i+1);
456 if(genptr[i]->prime_position > MAXPRIMEOFFSET)
457 {
458 fprintf(stderr,"WARNING - spawn_rng: gennum: %d > maximum number of independent streams: %d\n\tIndependence of streams cannot be guranteed.\n",
459 genptr[i]->prime_position, MAX_STREAMS);
460 genptr[i]->prime_position %= MAXPRIMEOFFSET;
461 }
462
463 genptr[i]->prime_next = (nspawned+1)*tempptr->prime_next;
464 getprime(1, &(genptr[i]->prime), genptr[i]->prime_position);
465 genptr[i]->multiplier = tempptr->multiplier;
466 genptr[i]->parameter = tempptr->parameter;
467 genptr[i]->gentype = tempptr->gentype;
468
469 #ifdef LONG64
470 genptr[i]->seed = INIT_SEED; /* initialize generator */
471 genptr[i]->seed ^= ((unsigned LONG64) tempptr->init_seed)<<16;
472
473 if (genptr[i]->prime == 0)
474 genptr[i]->seed |= 1;
475 #else
476 genptr[i]->seed[1] = 16651885^((tempptr->init_seed<<16)&(0xff0000));
477 genptr[i]->seed[0] = 2868876^((tempptr->init_seed>>8)&(0xffffff));
478 if (genptr[i]->prime == 0)
479 genptr[i]->seed[1] |= 1;
480 #endif
481
482 if(genptr[i]->prime_position > MAXPRIMEOFFSET)
483 advance_seed(genptr[i]); /* advance lcg 10^6 steps from initial seed */
484
485 for(j=0; j<LCGRUNUP*(genptr[i]->prime_position); j++)
486 get_rn_dbl( (int *) genptr[i]);
487 }
488 tempptr->prime_next = (nspawned+1)*tempptr->prime_next;
489
490 NGENS += nspawned;
491
492 *newgens = (int **) genptr;
493
494 if(checkid != 0)
495 for(i=0; i<nspawned; i++)
496 if(addID(( int *) genptr[i]) == NULL)
497 return i;
498
499 return nspawned;
500 }
501
502
503
504 /*Compute a + b. a and b are positive 4 digit integers */
505 /* in base 2^12, modulo 2^48 */
506 #ifdef __STDC__
plus(int * a,int * b,int * result)507 void plus(int *a, int *b, int *result)
508 #else
509 void plus(a,b,result)
510 int *a, *b, *result;
511 #endif
512 {
513 int temp[5];
514 int i;
515
516 temp[4] = 0;
517
518 for(i=0; i<4; i++)
519 temp[i] = a[i] + b[i];
520
521 for(i=1; i<5; i++)
522 {
523 temp[i] += temp[i-1]>>12;
524 temp[i-1] &= 4095;
525 }
526
527 for(i=0; i<4 ; i++)
528 result[i] = temp[i];
529 }
530
531
532
533 /*multiply two 4 digit numbers in base 2^12 and return 'size' lowest order */
534 /* digits*/
535 #ifdef __STDC__
mult(int * a,int * b,int * c,int size)536 void mult(int *a, int *b, int *c, int size)
537 #else
538 void mult(a,b,c,size)
539 int *a, *b, *c, size;
540 #endif
541 {
542 int temp[8];
543 int i, j;
544
545 for(i=0; i<8; i++)
546 temp[i] = 0;
547
548 for(i=0; i<4; i++)
549 for(j=0; j<4; j++)
550 temp[i+j] += a[i]*b[j];
551
552 for(i=1; i<8; i++)
553 {
554 temp[i] += temp[i-1]>>12;
555 temp[i-1] &= 4095;
556 }
557
558 for(i=0; i<size && i<8 ; i++)
559 c[i] = temp[i];
560 }
561
562
563 #ifdef __STDC__
advance_seed(struct rngen * gen)564 void advance_seed(struct rngen *gen)
565 #else
566 void advance_seed(gen)
567 struct rngen *gen;
568 #endif
569 {
570 #ifdef LONG64
571 int i, found;
572 unsigned LONG64 an, pmult;
573
574 for(i=0, found=0; i<NPARAMS; i++)
575 if (gen->multiplier == mults[i])
576 {
577 found = 1;
578 break;
579 }
580 if(found == 0)
581 {
582 fprintf(stderr,"WARNING: advance_seed: multiplier not acceptable.\n");
583 return ;
584 }
585
586 /* a^n, n = 10^6 and pmult = (a^n-1)/(a-1), n = 10^6 */
587 switch(i)
588 {
589 case 0 :
590 an = AN1;
591 pmult = PMULT1;
592 break;
593 case 1 :
594 an = AN2;
595 pmult = PMULT2;
596 break;
597 case 2 :
598 an = AN3;
599 pmult = PMULT3;
600 break;
601 case 3 :
602 an = AN4;
603 pmult = PMULT4;
604 break;
605 case 4 :
606 an = AN5;
607 pmult = PMULT5;
608 break;
609 case 5 :
610 an = AN6;
611 pmult = PMULT6;
612 break;
613 case 6 :
614 an = AN7;
615 pmult = PMULT7;
616 break;
617 default:
618 fprintf(stderr,"WARNING: advance_seed parameters for multiplier %d not available\n", i);
619 return;
620 }
621
622 gen->seed = gen->seed*an + pmult*gen->prime;
623 gen->seed &= LSB48;
624
625 #else
626 int an[4], pmult[4], x0, x1, temp[4],temp2[4], i, found;
627
628 for(i=0, found=0; i<NPARAMS; i++)
629 if (strxncmp((char *) gen->multiplier, (char *) (mults[i]), 4*sizeof(int))
630 == 0)
631 {
632 found = 1;
633 break;
634 }
635 if(found == 0)
636 {
637 fprintf(stderr,"WARNING: advance_seed: multiplier not acceptable.\n");
638 return ;
639 }
640
641 /* a^n, n = 10^6 and pmult = (a^n-1)/(a-1), n = 10^6 */
642 switch(i)
643 {
644 case 0 :
645 an[0] = 0x001; an[1] = 0xc00; an[2] = 0xf0a; an[3] = 0xdad;
646 pmult[0] = 0x000; pmult[1] = 0x700; pmult[2] = 0xc22; pmult[3] = 0xa42;
647 break;
648 case 1 :
649 an[0] = 0x001; an[1] = 0x400; an[2] = 0xd7a; an[3] = 0xfef;
650 pmult[0] = 0x000; pmult[1] = 0xb00; pmult[2] = 0x58c; pmult[3] = 0xfa8;
651 break;
652 case 2 :
653 an[0] = 0x001; an[1] = 0xc00; an[2] = 0x7b5; an[3] = 0x641;
654 pmult[0] = 0x000; pmult[1] = 0xf00; pmult[2] = 0xc4e; pmult[3] = 0x0d0;
655 break;
656 case 3 :
657 an[0] = 0x001; an[1] = 0xc00; an[2] = 0xf72; an[3] = 0xcf9;
658 pmult[0] = 0x000; pmult[1] = 0x700; pmult[2] = 0xdbe; pmult[3] = 0x11b;
659 break;
660 case 4 :
661 an[0] = 0x001; an[1] = 0x400; an[2] = 0x07b; an[3] = 0xbdf;
662 pmult[0] = 0x000; pmult[1] = 0x300; pmult[2] = 0xc8e; pmult[3] = 0xc3c;
663 break;
664 case 5 :
665 an[0] = 0x001; an[1] = 0xc00; an[2] = 0x747; an[3] = 0xf33;
666 pmult[0] = 0x000; pmult[1] = 0xf00; pmult[2] = 0x0e9; pmult[3] = 0xb0f;
667 break;
668 case 6 :
669 an[0] = 0x001; an[1] = 0xc00; an[2] = 0x632; an[3] = 0xcbe;
670 pmult[0] = 0x000; pmult[1] = 0x700; pmult[2] = 0x7de; pmult[3] = 0x640;
671 break;
672 default:
673 fprintf(stderr,"WARNING: advance_seed parameters for multiplier %d not available\n", i);
674 return;
675 }
676
677 x0 = gen->seed[0]; x1 = gen->seed[1];
678
679 temp[0] = x1&4095; temp[1] = (x1>>12)&4095; temp[2] = x0&4095; /* seed */
680 temp[3] = (x0>>12)&4095;
681
682 temp2[0] = gen->prime%(1<<12); temp2[1] = (gen->prime>>12)%(1<<12);
683 temp2[2] = (gen->prime>>24)%(1<<12); temp2[3] = 0; /* prime */
684
685 mult(temp,an,temp,4);
686 mult(temp2,pmult,temp2,4);
687
688 plus(temp,temp2,temp);
689
690 gen->seed[1] = (temp[1]<<12) + temp[0];
691 gen->seed[0] = (temp[3]<<12) + temp[2];
692 #endif
693 }
694
695
696 #ifdef __STDC__
free_rng(int * genptr)697 int free_rng(int *genptr)
698 #else
699 int free_rng(genptr)
700 int *genptr;
701 #endif
702 {
703 struct rngen *q;
704
705 q = (struct rngen *) genptr;
706 free(q);
707
708 NGENS--;
709
710 return NGENS;
711 }
712
713
714 #ifdef __STDC__
pack_rng(int * genptr,char ** buffer)715 int pack_rng( int *genptr, char **buffer)
716 #else
717 int pack_rng(genptr,buffer)
718 int *genptr;
719 char **buffer;
720 #endif
721 {
722 unsigned char *p, *initp;
723 unsigned int size, m[2], i;
724 struct rngen *q;
725
726 q = (struct rngen *) genptr;
727 size = 5*4 /*sizeof(int)*/ + 2*8/*sizeof(unsigned LONG64)*/
728 + strlen(q->gentype)+1;
729 /* The new load/store routines make using sizeof unnecessary. Infact, */
730 /* using sizeof could be erroneous. */
731 initp = p = (unsigned char *) mymalloc(size);
732 if(p == NULL)
733 {
734 *buffer = NULL;
735 return 0;
736 }
737
738
739 strcpy((char *)p,q->gentype);
740 p += strlen(q->gentype)+1;
741 #ifdef LONG64
742 p += store_long64(q->seed,8,p);
743 p += store_int(q->init_seed,4,p);
744 p += store_int(q->prime,4,p);
745 p += store_int(q->prime_position,4,p);
746 p += store_int(q->prime_next,4,p);
747 p += store_int(q->parameter,4,p);
748 p += store_long64(q->multiplier,8,p);
749 #else
750 m[0] = q->seed[0]>>8;m[1] = q->seed[0]<<24 | q->seed[1];
751 p += store_intarray(m,2,4,p);
752 p += store_int(q->init_seed,4,p);
753 p += store_int(q->prime,4,p);
754 p += store_int(q->prime_position,4,p);
755 p += store_int(q->prime_next,4,p);
756 p += store_int(q->parameter,4,p);
757 /* The following is done since the multiplier is stored in */
758 /* pieces of 12 bits each */
759 m[1] = q->multiplier[2]&0xff; m[1] <<= 24;
760 m[1] |= q->multiplier[1]<<12; m[1] |= q->multiplier[0];
761 m[0] = (q->multiplier[3]<<4) | (q->multiplier[2]>>8);
762
763 p += store_intarray(m,2,4,p);
764 #endif
765 *buffer = (char *) initp;
766
767 assert(p-initp == size);
768
769 return p-initp;
770 }
771
772
773
774 #ifdef __STDC__
get_seed_rng(int * gen)775 int get_seed_rng(int *gen)
776 #else
777 int get_seed_rng(gen)
778 int *gen;
779 #endif
780 {
781 return ((struct rngen *) gen)->init_seed;
782 }
783
784
785 #ifdef __STDC__
unpack_rng(char * packed)786 int *unpack_rng( char *packed)
787 #else
788 int *unpack_rng(packed)
789 char *packed;
790 #endif
791 {
792 struct rngen *q;
793 unsigned char *p;
794 unsigned int m[4], m2[2], i;
795
796 p = (unsigned char *) packed;
797 q = (struct rngen *) mymalloc(sizeof(struct rngen));
798 if(q == NULL)
799 return NULL;
800
801 if(strcmp((char *)p,GENTYPE) != 0)
802 {
803 fprintf(stderr,"ERROR: Unpacked ' %.24s ' instead of ' %s '\n",
804 p, GENTYPE);
805 return NULL;
806 }
807 else
808 q->gentype = GENTYPE;
809 p += strlen(q->gentype)+1;
810
811 #ifdef LONG64
812 p += load_long64(p,8,&q->seed);
813 p += load_int(p,4,(unsigned int *)&q->init_seed);
814 p += load_int(p,4,(unsigned int *)&q->prime);
815 p += load_int(p,4,(unsigned int *)&q->prime_position);
816 p += load_int(p,4,(unsigned int *)&q->prime_next);
817 p += load_int(p,4,(unsigned int *)&q->parameter);
818 p += load_long64(p,8,&q->multiplier);
819
820 #else
821 p += load_intarray(p,2,4,m2);
822 q->seed[1] = m2[1]&0xffffff; q->seed[0] = m2[1]>>24 | m2[0]<<8;
823 p += load_int(p,4,(unsigned int *)&q->init_seed);
824 p += load_int(p,4,(unsigned int *)&q->prime);
825 p += load_int(p,4,(unsigned int *)&q->prime_position);
826 p += load_int(p,4,(unsigned int *)&q->prime_next);
827 p += load_int(p,4,(unsigned int *)&q->parameter);
828 p += load_intarray(p,2,4,m2);
829 /* The following is done since the multiplier is stored in */
830 /* pieces of 12 bits each */
831 m[0] = m2[1]&0xfff; m[1] = (m2[1]&0xfff000)>>12;
832 m[2] = (m2[1]&0xff000000)>>24 | (m2[0]&0xf)<<8;
833 m[3] = (m2[0]&0xfff0)>>4;
834
835 #endif
836
837
838 if(q->parameter < 0 || q->parameter >= NPARAMS)
839 {
840 fprintf(stderr,"ERROR: Unpacked parameters not acceptable.\n");
841 free(q);
842 return NULL;
843 }
844
845 q->multiplier = mults[q->parameter];
846
847 NGENS++;
848
849 return (int *) q;
850 }
851
852
853
854 #ifdef __STDC__
print_rng(int * igen)855 int print_rng( int *igen)
856 #else
857 int print_rng(igen)
858 int *igen;
859 #endif
860 {
861 struct rngen *gen;
862
863 printf("\n%s\n", GENTYPE+2);
864
865 gen = (struct rngen *) igen;
866 printf("\n \tseed = %d, stream_number = %d\tparameter = %d\n\n", gen->init_seed, gen->prime_position, gen->parameter);
867
868 return 1;
869 }
870
871
872 /* HAS ;-) */
873 #if 0
874 #include "simple_.h"
875 #include "../fwrap_.h"
876 #endif
877