1 /*------------------------------------------------------
2 Maximum likelihood estimation
3 of migration rate and effectice population size
4 using a Metropolis-Hastings Monte Carlo algorithm
5 -------------------------------------------------------
6 R A N D O M G E N E R A T O R R O U T I N E S
7
8 creates options structures,
9 reads options from parmfile if present
10
11 prints options,
12 and finally helps to destroy itself.
13
14 Peter Beerli 1996, Seattle
15 beerli@fsu.edu
16
17 Copyright 1996-2002 Peter Beerli and Joseph Felsenstein, Seattle WA
18 Copyright 2003-2007 Peter Beerli, Tallahassee FL
19
20 Permission is hereby granted, free of charge, to any person obtaining
21 a copy of this software and associated documentation files (the
22 "Software"), to deal in the Software without restriction, including
23 without limitation the rights to use, copy, modify, merge, publish,
24 distribute, sublicense, and/or sell copies of the Software, and to
25 permit persons to whom the Software is furnished to do so, subject
26 to the following conditions:
27
28 The above copyright notice and this permission notice shall be
29 included in all copies or substantial portions of the Software.
30
31 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
32 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
33 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
34 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR
35 ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
36 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
37 WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
38
39
40 $Id: random.c 2169 2013-08-24 19:02:04Z beerli $
41 -------------------------------------------------------*/
42 /* \file random.c */
43 #include <stdlib.h>
44 #include "sighandler.h"
45 #include "migration.h"
46 #include "heating.h"
47 #include "random.h"
48 #include "migrate_mpi.h"
49 #ifdef WIN32
50 #include <sys/timeb.h>
51 #else /* WIN32 */
52 #include <sys/time.h>
53 #endif /* WIN32 */
54
55
56 #ifdef MERSENNE_TWISTER
57 #include "SFMT.c"
58 #endif
59
60 #ifdef DMALLOC_FUNC_CHECK
61 #include <dmalloc.h>
62 #endif
63
64
65 #ifdef GRANDCENTRAL
66 //#include <dispatch/dispatch.h>
67 #endif
68
69
70 #ifdef PTHREADS
71 extern tpool_t heating_pool;
72 #endif
73
74 /* prototypes ----------------------------------------- */
75 void getseed (option_fmt * options);
76 void swap_ptr (long **ptr, long **newptr);
77
78 MYREAL randum (void);
79 #ifdef PTHREADS
80 MYREAL randum_thread (void);
81 #endif
82 long random_integer(long low, long high);
83
84 #ifdef QUASIRANDOM
85 #define MYINDEX 0
86
87 //char[100] generator;
88
89 #ifdef CUD
90 int nextn=1, kk=1,ii=0,jj=0;
setup_cud(int thiskk)91 void setup_cud(int thiskk)
92 {
93 strncpy(generator,"Quasi-random number generator: Completely Uniform Distributed numbers", 79);
94 kk = thiskk;
95 }
96
get_cud()97 MYINLINE MYREAL get_cud()
98 {
99 MYREAL rr;
100 MYREAL value;
101 const MYREAL magicnumbers[]={2.,2.,4.,15.,65.,315.,1586.,8036.,40448.,200401.,972536.,
102 4609846.,21310545.,96017492., 421606654.,1804551131. };
103 //double ps[]={2.,3.,5.,7.,11.,13.,17.,19.,23.,29.,31.,37.,41.,47.,53.,59.,61.,67.,71.,73.,79.,83.,89.,97.,101.};
104 const MYREAL lps[] =
105 {0.69314718055994530942, 1.0986122886681096914,
106 1.6094379124341003746, 1.9459101490553133051,
107 2.3978952727983705441, 2.5649493574615367361,
108 2.8332133440562160802, 2.9444389791664404600,
109 3.1354942159291496908, 3.3672958299864740272,
110 3.4339872044851462459, 3.6109179126442244444,
111 3.7135720667043078039, 3.8501476017100585868,
112 3.9702919135521218341, 4.0775374439057194506,
113 4.1108738641733112488, 4.2046926193909660597,
114 4.2626798770413154213, 4.2904594411483911291,
115 4.3694478524670214942, 4.4188406077965979235,
116 4.4886363697321398383, 4.5747109785033828221,
117 4.6151205168412594509} ;
118 //rr=kk*log(ps[jj]);
119 #ifdef PTHREADS
120
121 if ((pthread_mutex_lock (&(heating_pool->random_lock))) != 0)
122 error ("pthread_mutex_lock failed in get_cud()");
123 #endif
124
125 rr = kk * lps[jj];
126
127 if (jj >= ii)
128 {
129 kk++;
130 jj=0;
131 }
132 else
133 jj++;
134
135 if(kk >= magicnumbers[ii])
136 {
137 kk = 1;
138 ii++;
139 }
140 value = rr - floor(rr);
141 #ifdef PTHREADS
142
143 if ((pthread_mutex_unlock (&(heating_pool->random_lock))) != 0)
144 error ("pthread_mutex_unlock failed in get_cud()");
145 #endif
146
147 return value;
148 //printf("%f\n",rr);
149 }
150
get_quasi()151 MYINLINE MYREAL get_quasi()
152 {
153 return get_cud();
154 }
155
156 #endif
157
158 #ifdef KOROBOV
159
160 unsigned int x0=137,ii=0; nextn=1;
161
setup_korobov()162 void setup_korobov()
163 {
164 strncpy(generator,"Quasi-random number generator: Korobov sequence", 79);
165 }
166
get_korobov()167 MYREAL get_korobov()
168 {
169 unsigned long a =17364, m=65521;
170 unsigned long xn;
171 double x;
172
173 xn=(a*x0)%m;
174 x=(xn+0.0)/(m+0.0);
175 x0=xn; nextn++;
176 if (nextn>m){ ii++; nextn=2+ii;}
177 return x;
178 }
179
180
get_quasi()181 MYINLINE MYREAL get_quasi()
182 {
183 return get_korobov();
184 }
185
186 #endif
187
188 #ifdef WEYL
189 static int s, nextn=8;
190
setup_weyl()191 void setup_weyl()
192 {
193 strncpy(generator,"Quasi-random number generator: Weyl sequence", 79);
194 }
195
196
get_weyl()197 MYINLINE MYREAL get_weyl()
198 {
199 double r=sqrt(7.0),rr;
200 int next;
201
202 next=nextn*nextn;
203 rr=r*next;
204
205 // rr= (r*nextn-floor(r*nextn))*nextn;
206 rr=rr-floor(rr);
207
208 nextn++;
209 return rr;
210 }
211
get_quasi()212 MYINLINE MYREAL get_quasi()
213 {
214 return get_weyl();
215 }
216 #endif /*WEYL*/
217
218 #ifdef HALTON
219
220 #define MAX_D 500
221 static int s, nextn=8;
222
223 double e, quasi[100];
224 static int prime[MAX_D];
225 static double iprime[MAX_D];
226 static int primroots[][10]={{1, 2, 3, 3, 8, 11,12,14, 7,18},
227 {12,13,17,18,29,14,18,43,41,44},
228 {40,30,47,65,71,28,40,60,79,89},
229 {56,50,52,61,108,56,66,63,60,66},
230 {104,76,111,142,71,154,118,84,127,142},
231 {84,105,186,178,188,152,165,159,103,205},
232 {166,173,188,181,91,233,210,217,153,212},
233 };
234
235 static int warnockOpt[]={1, 2, 2, 5, 3, 7, 3, 10, 18, 11,
236 17, 5, 17, 26, 40, 14, 40, 44, 12, 31,
237 45, 70,8, 38, 82, 8, 12, 38, 47, 70,
238 29, 57, 97, 110,32, 48, 84, 124,155,26,
239 69, 83, 157,171, 8, 22, 112,205, 15, 31,
240 61, 105,127,212,12, 57, 109,133,179,210,
241 231,34, 161,199,222,255,59, 120,218,237,
242 278,341,54, 110,176,218,280,369,17, 97,
243 193,221,331,350,419,21, 85, 173,221,243,
244 288,424,45, 78, 173,213,288,426,455,138,
245 };
246
247 int gohalt(double *,double *, double *);
248 double get_halton();
249 int primes();
250 int power(int, int, int);
251 int inhalt(int dimen, int atmost, double tiny, double *quasi);
power(int a,int b,int m)252 int power(int a, int b, int m)
253 { int i,c=1;
254 for(i=0;i<b;i++)
255 c=(c*a)%m;
256 return c;
257 }
258
setup_halton()259 void setup_halton()
260 {
261 strncpy(generator,"Quasi-random number generator: Halton sequence", 79);
262 }
263
get_halton()264 double get_halton()
265 { double wq[100],dq[100];
266 gohalt(quasi,dq,wq);
267 return wq[MYINDEX];
268 }
269
inhalt(int dimen,int atmost,double tiny,double * quasi)270 int inhalt(int dimen, int atmost, double tiny, double *quasi)
271 {
272 double delta, f;
273 int i,m;
274
275 // check dimen
276 primes();
277
278 s=dimen;
279 if (s<1||s>1000)
280 return(-1);
281
282 // compute and check tolerance
283
284 e=0.9*(1.0/(atmost*prime[s-1])-10.0*tiny);
285 delta=100*tiny*(double)(atmost+1)*log10((double)atmost);
286 if (delta>=0.09*(e-10.0*tiny))
287 return(-2);
288
289 // now compute first vector
290
291 m=1;
292 for (i=0;i<s;i++)
293 {
294 iprime[i]=1.0/iprime[i];
295 quasi[i]=iprime[i];
296 m=i*prime[i];
297 }
298
299 printf("largest prime=%d, %f \n",prime[s-1], quasi[1]);
300
301 nextn=2;
302
303 return 0;
304 }
305
gohalt(double * quasi,double * dq,double * wq)306 int gohalt(double *quasi, double *dq, double *wq)
307 {
308 int i, j, k, ytemp[40],xtemp[40], ztemp, ktemp, ltemp, mtemp;
309 double r;
310 double t,f,g,h;
311
312 // generate quasi one compoment at a time using radix prime[k] for
313 // component k
314
315
316 for (i=0;i<s;i++)
317 {
318 t=iprime[i];
319 f=1.0-quasi[i];
320 g=1.0;
321 h=t;
322 while ((f-h)<e)
323 // this checks whether q+h>1-e
324 {
325 g=h;
326 h*=t;
327 }
328 quasi[i]=g+h-f;
329 }
330
331 for(i=0;i<s;i++)
332 {
333 k=0; mtemp=nextn;
334 ltemp=prime[i];
335
336 while(mtemp!=0){
337 ytemp[k]=mtemp%ltemp;
338 mtemp=mtemp/ltemp;
339 k++;
340 }
341
342 //generating Optimal primitive root
343 for(j=0;j<k;j++)
344 {
345 // xtemp[j] = (ytemp[j]*power(primroots[i/10][i%10], nextn%ltemp, ltemp))%ltemp;
346 if(j>=1)
347 xtemp[j] =(warnockOpt[i]*power(primroots[i/10][i%10], ytemp[j],
348 ltemp)+ytemp[j-1])%ltemp;
349 else xtemp[j] =(warnockOpt[i]*power(primroots[i/10][i%10], ytemp[j],
350 ltemp))%ltemp;
351 xtemp[j] -= ytemp[j];
352 }
353
354 dq[i]=0;t=iprime[i];
355 for(j=0;j<k; j++)
356 {
357 dq[i] += xtemp[j]*t;
358 t *= iprime[i];
359 }
360
361 dq[i] += quasi[i];
362
363
364 // generating Warnock Optimal sequences
365 for(j=0;j<k;j++)
366 {
367 if(j>=1)
368 xtemp[j]= (ytemp[j]*power(warnockOpt[i],i+1,ltemp)+ytemp[j-1])%ltemp;
369 else
370 xtemp[j]= (ytemp[j]*power(warnockOpt[i],i+1,ltemp))%ltemp;
371
372 xtemp[j] -= ytemp[j];
373 }
374
375 wq[i]=0;t=iprime[i];
376 for(j=0;j<k; j++)
377 {
378 wq[i] += xtemp[j]*t;
379 t *= iprime[i];
380 }
381
382 wq[i] += quasi[i];
383 }
384
385 nextn++;
386 return(0);
387 }
388
primes()389 int primes()
390 {
391 int i, j, a[MAX_D+1];
392 for (a[1] = 0, i = 2; i <= MAX_D; i++)
393 a[i] = 1;
394 for (i = 2; i <= MAX_D/2; i++)
395 for (j = 2; j <= MAX_D/i; j++)
396 a[i*j] = 0;
397 for (i = 1, j = 0; i <= MAX_D; i++){
398 if (a[i]){
399 prime[j] =i;
400 iprime[j]=i;
401 j++;
402 }
403 }
404 return j;
405 }
406
get_quasi()407 MYINLINE MYREAL get_quasi()
408 {
409 return get_halton();
410 }
411
412 #endif /*HALTON */
413
414
415 #endif
416
417 #ifdef MERSENNE_TWISTER
setup_mersennetwister()418 void setup_mersennetwister()
419 {
420 strncpy(generator,"Pseudo-random number generator: Mersenne twister", 79);
421 }
422 #endif
423
424 /// random integer
425 //=============================================
random_integer(long low,long high)426 MYINLINE long random_integer(long low, long high)
427 {
428 //Math.floor(Math.random()*(N-M+1))%(N-M+1)+M
429 long r;
430 long tt = high-low+1;
431 r = low + (long) (RANDUM() * tt);
432 if(r<low || r>high)
433 {
434 // warning("Random %li integer is out of bounds (%li, %li)\n",r, low, high);
435 if(r>high)
436 return high;
437 else
438 return low;
439 // error("Stop because this should never happen");
440 }
441 return r;
442 }
443
444
445 ///
446 /// better random_seed function taken from a discussion on
447 /// http://sourceware.org/ml/gsl-discuss/2004-q1/msg00071.html
448 /// March 30, 2007, Robert G. Brown (http://www.phy.duke.edu/~rgb/)
449 /// at Duke University Dept. of Physics, Box 90305 Durham, N.C. 27708-030
450 /// suggested the code
451 /// changed December 2012 because on som cluster machines little entropy
452 /// is generated thus the machine stalls (AMD?)
random_seed()453 unsigned long int random_seed()
454 {
455
456 unsigned int seed;
457 struct timeval tv;
458 FILE *devrandom;
459 size_t retval;
460 if ((devrandom = fopen("/dev/urandom","r")) == NULL)
461 {
462 gettimeofday(&tv,0);
463 seed = tv.tv_sec + tv.tv_usec;
464 }
465 else
466 {
467 retval = fread(&seed,sizeof(seed),1,devrandom);
468 fclose(devrandom);
469 if(seed == 0)
470 {
471 gettimeofday(&tv,0);
472 seed = tv.tv_sec + tv.tv_usec;
473 }
474 }
475 return(seed);
476 }
477
set_seed(long autoseed,unsigned long * inseed)478 void set_seed(long autoseed, unsigned long * inseed)
479 {
480 unsigned long timeseed = 0;
481 switch (autoseed)
482 {
483 case AUTO:
484 timeseed = random_seed();
485 #ifdef LCG
486 switch (timeseed % 4)
487 {
488 case 0:
489 ++timeseed; break;
490 case 2:
491 ++timeseed;
492 break;
493 case 1:
494 case 3:
495 break;
496 }
497 #endif
498 *inseed = (unsigned long) labs((long)timeseed);
499 break;
500 case NOAUTO:
501 case NOAUTOSELF:
502 break;
503 default:
504 error ("Error: Seed value not defined");
505 break;
506 }
507 }
508
init_lcg(option_fmt * options)509 void init_lcg(option_fmt *options)
510 {
511 long i;
512 for (i = 0; i <= 2; i++)
513 seed[i] = 0;
514 i = 0;
515 do
516 {
517 seed[i] = options->inseed & 2047;
518 printf("RANDOM SEEDING: i=%li seed[i]=%li inseed=%li\n",i,seed[i], options->inseed);
519 options->inseed /= 2048;
520 i++;
521 }
522 while (options->inseed != 0);
523 if(i>3)
524 error("lcg random number generator is broken! random.c:507");
525 }
526
getseed_lcg(option_fmt * options)527 void getseed_lcg (option_fmt *options)
528 {
529 strncpy(generator,"Pseudo-random number generator: Least Congruental Generator", 79);
530 set_seed(options->autoseed, &options->inseed);
531 options->saveseed = options->inseed;
532 init_lcg(options);
533 if(options->randomsubset>0 && options->randomsubsetseed > 0)
534 {
535 #ifdef WIN32
536 srand(options->randomsubsetseed);
537 #else
538 srand48(options->randomsubsetseed);
539 #endif
540 }
541 }
542
getseed_mt(option_fmt * options)543 void getseed_mt (option_fmt *options)
544 {
545 strncpy(generator,"Pseudo-random number generator: Mersenne-Twister", 79);
546 set_seed(options->autoseed, &options->inseed);
547 options->saveseed = options->inseed;
548 init_gen_rand((unsigned long)options->inseed);
549 if(options->randomsubset>0 && options->randomsubsetseed > 0)
550 {
551 #ifdef WIN32
552 srand(options->randomsubsetseed);
553 #else
554 srand48(options->randomsubsetseed);
555 #endif
556 }
557 }
558 ///
559 /// set up the quasi random number material
560 /// this function is dependent on running the standard MT or LCG random number generator first.
getseed_quasi(option_fmt * options)561 void getseed_quasi (option_fmt *options)
562 {
563 #ifdef CUD
564 setup_cud((int) RANDINT (1, MAXLONG-1));
565 #endif
566 #ifdef WEYL
567 setup_weyl();
568 #endif
569 #ifdef KOROBOV
570 setup_korobov();
571 #endif
572 #ifdef HALTON
573 setup_halton();
574 #endif
575
576 }
577
578
579 #ifdef MPI
580 /// Random number seed function for MPI usage, the master generates numcpu random numbers
581 /// and distributes them to the workers.
getseed(option_fmt * options)582 void getseed(option_fmt *options)
583 {
584 //long test = -1;
585 long i;
586 long worker;
587 long *newseeds;
588 long newseed;
589 boolean found=FALSE;
590 MPI_Status status;
591 if (myID==MASTER)
592 {
593 #ifdef MERSENNE_TWISTER
594 getseed_mt(options);
595 #else
596 getseed_lcg(options);
597 #endif
598 if(options->progress)
599 printf("Master random number seed: %li\n", options->inseed);
600 newseeds = (long*) calloc(numcpu,sizeof(long));
601 for (worker=1;worker<numcpu;worker++)
602 {
603 newseed = RANDINT(0,MAXLONG-1);
604 for (i=1;i<worker;i++)
605 {
606 if (newseed == newseeds[i])
607 {
608 found=TRUE;
609 break;
610 }
611 }
612 if (found)
613 worker--;
614 else
615 newseeds[worker] = newseed;
616 }
617 for (worker=1;worker<numcpu;worker++)
618 {
619 newseed = newseeds[worker];
620 MYMPISEND(&newseed,(MYINT) ONE, MPI_LONG, (MYINT) worker, 0, comm_world);
621 }
622 }
623 else /*worker*/
624 {
625 MYMPIRECV (&options->inseed, (MYINT) ONE, MPI_LONG, MPI_ANY_SOURCE, MPI_ANY_TAG, comm_world, &status);
626 #ifdef MERSENNE_TWISTER
627 getseed_mt(options);
628 #else
629 getseed_lcg(options);
630 #endif
631 #ifdef QUASIRANDOM
632 getseed_quasi(options);
633 #endif
634 #ifdef DEBUG
635 char nowstr[LINESIZE];
636 get_time (nowstr, "%c");
637 printf("%i> RANDOM NUMBER SEED:%li at time %s\n",myID,options->inseed,nowstr);
638 #endif
639 }
640 }
641 #if 0 /*does not execute this section old getseed function*/
642 /// Random number seed function for MPI usage, the main seed gets distributed to
643 /// all nodes and then the nodes pick a random number from the master seed for
644 /// their own seed. given that the nodes typically analyze different loci
645 /// it will be unlikely that even when picked by the off chance the same seed on two
646 /// nodes we get the same stream of events. the node ID gets the randomnumber[#ID] from the stream.
647 void getseed(option_fmt *options)
648 {
649 long test = -1;
650 long i;
651 if (myID==MASTER)
652 {
653 #ifdef MERSENNE_TWISTER
654 getseed_mt(options);
655 #else
656 getseed_lcg(options);
657 #endif
658 }
659 #ifdef DEBUG
660 char nowstr[LINESIZE];
661 get_time (nowstr, "%c");
662 printf("%i> random number seed broadcast start %s\n",myID, nowstr);
663 #endif
664 MYMPIBCAST (&options->inseed, 1, MPI_LONG, MASTER, comm_world);
665 #ifdef DEBUG
666 get_time (nowstr, "%c");
667 printf("%i> random number seed broadcast finished %s\n",myID, nowstr);
668 #endif
669 if(myID != MASTER)
670 {
671 // sets all worker to same seed as master
672 #ifdef MERSENNE_TWISTER
673 getseed_mt(options);
674 #else
675 getseed_lcg(options);
676 #endif
677 #ifdef DEBUG_MPI
678 printf("WARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNING\n");
679 printf("WARNING all random number seeds on the nodes are identical to the master WARNING\n");
680 printf("WARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNING\n");
681 return;
682 #endif
683 // gets a random value using the master seed
684 test = RANDINT (1, MAXLONG-1);
685 options->inseed = test;
686 // reset the master seed to the own seed.
687 #ifdef MERSENNE_TWISTER
688 getseed_mt(options);
689 #else
690 getseed_lcd(options);
691 #endif
692 #ifdef QUASIRANDOM
693 getseed_quasi(options);
694 #endif
695 #ifdef DEBUG
696 char nowstr[LINESIZE];
697 get_time (nowstr, "%c");
698 printf("%i> Random number seed: %li at time %s\n",myID,test,nowstr);
699 #endif
700 }
701 else
702 {
703 printf("%i> Master random number seed: %li\n",myID,options->inseed);
704 }
705 }
706 #endif /*0*/
707 #else
getseed(option_fmt * options)708 void getseed(option_fmt *options)
709 {
710 #ifdef MERSENNE_TWISTER
711 getseed_mt(options);
712 #else
713 getseed_lcg(options);
714 #endif
715 #ifdef QUASIRANDOM
716 getseed_quasi(options);
717 #endif
718 }
719 #endif
720
721
722 void
swap_ptr(long ** ptr,long ** newptr)723 swap_ptr (long **ptr, long **newptr)
724 {
725 long *temp;
726 temp = *ptr;
727 *ptr = *newptr;
728 *newptr = temp;
729 }
730
731 #ifdef PTHREADS
732 MYREAL
randum_thread(void)733 randum_thread (void)
734 /* thread save random number generator */
735 {
736 MYREAL value;
737 if (heating_pool != NULL && (pthread_mutex_lock (&(heating_pool->random_lock))) != 0)
738 error ("pthread_mutex_lock failed in random_thread()");
739 #ifdef MERSENNE_TWISTER
740 value = genrand_res53();
741 #else
742 newseed[0] = 1549 * seed[0];
743 newseed[1] = newseed[0] / 2048;
744 newseed[0] &= 2047;
745 newseed[2] = newseed[1] / 2048;
746 newseed[1] &= 2047;
747 newseed[1] += 1549 * seed[1] + 812 * seed[0];
748 newseed[2] += newseed[1] / 2048;
749 newseed[1] &= 2047;
750 newseed[2] += 1549 * seed[2] + 812 * seed[1];
751 swap_ptr (&newseed, &seed);
752 seed[2] &= 1023;
753 value = (((seed[0] / 2048.0 + seed[1]) / 2048.0 + seed[2]) / 1024.0);
754 #endif
755 if (heating_pool != NULL && (pthread_mutex_unlock (&(heating_pool->random_lock))) != 0)
756 error ("pthread_mutex_unlock failed in randum_thread()");
757 return value;
758 } /* threaded randum */
759 #else
760 MYREAL
randum(void)761 randum (void)
762 /* Non thread-safe random number generator (faster) */
763 {
764 #ifdef MERSENNE_TWISTER
765 #ifdef MESS
766 MYREAL val = genrand_res53();
767 MYREAL d = fabs(val-0.7489040);
768 fprintf(stdout,"%i> R:%f\n",myID,val);
769 return val;
770 #else
771 #ifdef GRANDCENTRAL
772 // MYREAL r;
773 // dispatch_queue_t queue = dispatch_get_global_queue(DISPATCH_QUEUE_PRIORITY_DEFAULT, 0);
774 //dispatch_async(queue,
775 // ^{
776 MYREAL r = genrand_res53();
777 // }
778 //);
779 //fprintf(stdout,"%i> R:%f\n", myID, r);
780 return r;
781 #else
782 MYREAL r = genrand_res53();
783 return r;
784 #endif
785 #endif
786 #else
787 MYREAL value;
788 newseed[0] = 1549 * seed[0];
789 newseed[1] = newseed[0] / 2048;
790 newseed[0] &= 2047;
791 newseed[2] = newseed[1] / 2048;
792 newseed[1] &= 2047;
793 newseed[1] += 1549 * seed[1] + 812 * seed[0];
794 newseed[2] += newseed[1] / 2048;
795 newseed[1] &= 2047;
796 newseed[2] += 1549 * seed[2] + 812 * seed[1];
797 swap_ptr (&newseed, &seed);
798 seed[2] &= 1023;
799 value = (((seed[0] / 2048.0 + seed[1]) / 2048.0 + seed[2]) / 1024.0);
800 return value;
801 #endif /*mersenne_twister*/
802 } /* randum */
803 #endif /*phtreads*/
804