BOSS 6.6.4.p03
BESIII Offline Software System
Loading...
Searching...
No Matches
ranlxd.c
Go to the documentation of this file.
1
2/*******************************************************************************
3*
4* file ranlxd.c
5*
6* Random number generator "ranlxd". See the notes
7*
8* "User's guide for ranlxs and ranlxd v3.0" (May 2001)
9*
10* "Algorithms used in ranlux v3.0" (May 2001)
11*
12* for a detailed description
13*
14* The externally accessible functions are
15*
16* void ranlxd(double r[],int n)
17* Computes the next n single-precision random numbers and
18* assigns them to the elements r[0],...,r[n-1] of the array r[]
19*
20* void rlxd_init(int level,int seed)
21* Initialization of the generator
22*
23* int rlxd_size(void)
24* Returns the number of integers required to save the state of
25* the generator
26*
27* void rlxd_get(int state[])
28* Extracts the current state of the generator and stores the
29* information in the array state[N] where N>=rlxd_size()
30*
31* void rlxd_reset(int state[])
32* Resets the generator to the state defined by the array state[N]
33*
34* Version: 3.0
35* Author: Martin Luescher <[email protected]>
36* Date: 06.05.2001
37*
38*******************************************************************************/
39
40#include <limits.h>
41#include <float.h>
42#include <stdlib.h>
43#include <stdio.h>
44#include <math.h>
45
46#if ((defined SSE)||(defined SSE2))
47
48typedef struct
49{
50 float c1,c2,c3,c4;
51} vec_t __attribute__ ((aligned (16)));
52
53typedef struct
54{
55 vec_t c1,c2;
56} dble_vec_t __attribute__ ((aligned (16)));
57
58static int init=0,pr,prm,ir,jr,is,is_old,next[96];
59static vec_t one,one_bit,carry;
60
61static union
62{
63 dble_vec_t vec[12];
64 float num[96];
65} x __attribute__ ((aligned (16)));
66
67#define STEP(pi,pj) \
68 __asm__ __volatile__ ("movaps %2, %%xmm4 \n\t" \
69 "movaps %%xmm2, %%xmm3 \n\t" \
70 "subps %0, %%xmm4 \n\t" \
71 "movaps %%xmm1, %%xmm5 \n\t" \
72 "cmpps $0x6, %%xmm4, %%xmm2 \n\t" \
73 "andps %%xmm2, %%xmm5 \n\t" \
74 "subps %%xmm3, %%xmm4 \n\t" \
75 "andps %%xmm0, %%xmm2 \n\t" \
76 "addps %%xmm4, %%xmm5 \n\t" \
77 "movaps %%xmm5, %0 \n\t" \
78 "movaps %3, %%xmm6 \n\t" \
79 "movaps %%xmm2, %%xmm3 \n\t" \
80 "subps %1, %%xmm6 \n\t" \
81 "movaps %%xmm1, %%xmm7 \n\t" \
82 "cmpps $0x6, %%xmm6, %%xmm2 \n\t" \
83 "andps %%xmm2, %%xmm7 \n\t" \
84 "subps %%xmm3, %%xmm6 \n\t" \
85 "andps %%xmm0, %%xmm2 \n\t" \
86 "addps %%xmm6, %%xmm7 \n\t" \
87 "movaps %%xmm7, %1" \
88 : \
89 "+m" ((*pi).c1), \
90 "+m" ((*pi).c2) \
91 : \
92 "m" ((*pj).c1), \
93 "m" ((*pj).c2))
94
95
96static void error(int no)
97{
98 switch(no)
99 {
100 case 1:
101 printf("Error in subroutine rlxd_init\n");
102 printf("Bad choice of luxury level (should be 1 or 2)\n");
103 break;
104 case 2:
105 printf("Error in subroutine rlxd_init\n");
106 printf("Bad choice of seed (should be between 1 and 2^31-1)\n");
107 break;
108 case 3:
109 printf("Error in rlxd_get\n");
110 printf("Undefined state (ranlxd is not initialized\n");
111 break;
112 case 5:
113 printf("Error in rlxd_reset\n");
114 printf("Unexpected input data\n");
115 break;
116 }
117 printf("Program aborted\n");
118 exit(0);
119}
120
121
122static void update()
123{
124 int k,kmax;
125 dble_vec_t *pmin,*pmax,*pi,*pj;
126
127 kmax=pr;
128 pmin=&x.vec[0];
129 pmax=pmin+12;
130 pi=&x.vec[ir];
131 pj=&x.vec[jr];
132
133 __asm__ __volatile__ ("movaps %0, %%xmm0 \n\t"
134 "movaps %1, %%xmm1 \n\t"
135 "movaps %2, %%xmm2"
136 :
137 :
138 "m" (one_bit),
139 "m" (one),
140 "m" (carry));
141
142 for (k=0;k<kmax;k++)
143 {
144 STEP(pi,pj);
145 pi+=1;
146 pj+=1;
147 if (pi==pmax)
148 pi=pmin;
149 if (pj==pmax)
150 pj=pmin;
151 }
152
153 __asm__ __volatile__ ("movaps %%xmm2, %0"
154 :
155 :
156 "m" (carry));
157
158 ir+=prm;
159 jr+=prm;
160 if (ir>=12)
161 ir-=12;
162 if (jr>=12)
163 jr-=12;
164 is=8*ir;
165 is_old=is;
166}
167
168
169static void define_constants()
170{
171 int k;
172 float b;
173
174 one.c1=1.0f;
175 one.c2=1.0f;
176 one.c3=1.0f;
177 one.c4=1.0f;
178
179 b=(float)(ldexp(1.0,-24));
180 one_bit.c1=b;
181 one_bit.c2=b;
182 one_bit.c3=b;
183 one_bit.c4=b;
184
185 for (k=0;k<96;k++)
186 {
187 next[k]=(k+1)%96;
188 if ((k%4)==3)
189 next[k]=(k+5)%96;
190 }
191}
192
193
194void rlxd_init(int level,int seed)
195{
196 int i,k,l;
197 int ibit,jbit,xbit[31];
198 int ix,iy;
199
200 define_constants();
201
202 if (level==1)
203 pr=202;
204 else if (level==2)
205 pr=397;
206 else
207 error(1);
208
209 i=seed;
210
211 for (k=0;k<31;k++)
212 {
213 xbit[k]=i%2;
214 i/=2;
215 }
216
217 if ((seed<=0)||(i!=0))
218 error(2);
219
220 ibit=0;
221 jbit=18;
222
223 for (i=0;i<4;i++)
224 {
225 for (k=0;k<24;k++)
226 {
227 ix=0;
228
229 for (l=0;l<24;l++)
230 {
231 iy=xbit[ibit];
232 ix=2*ix+iy;
233
234 xbit[ibit]=(xbit[ibit]+xbit[jbit])%2;
235 ibit=(ibit+1)%31;
236 jbit=(jbit+1)%31;
237 }
238
239 if ((k%4)!=i)
240 ix=16777215-ix;
241
242 x.num[4*k+i]=(float)(ldexp((double)(ix),-24));
243 }
244 }
245
246 carry.c1=0.0f;
247 carry.c2=0.0f;
248 carry.c3=0.0f;
249 carry.c4=0.0f;
250
251 ir=0;
252 jr=7;
253 is=91;
254 is_old=0;
255 prm=pr%12;
256 init=1;
257}
258
259
260void ranlxd(double r[],int n)
261{
262 int k;
263
264 if (init==0)
265 rlxd_init(1,1);
266
267 for (k=0;k<n;k++)
268 {
269 is=next[is];
270 if (is==is_old)
271 update();
272 r[k]=(double)(x.num[is+4])+(double)(one_bit.c1*x.num[is]);
273 }
274}
275
276
277int rlxd_size(void)
278{
279 return(105);
280}
281
282
283void rlxd_get(int state[])
284{
285 int k;
286 float base;
287
288 if (init==0)
289 error(3);
290
291 base=(float)(ldexp(1.0,24));
292 state[0]=rlxd_size();
293
294 for (k=0;k<96;k++)
295 state[k+1]=(int)(base*x.num[k]);
296
297 state[97]=(int)(base*carry.c1);
298 state[98]=(int)(base*carry.c2);
299 state[99]=(int)(base*carry.c3);
300 state[100]=(int)(base*carry.c4);
301
302 state[101]=pr;
303 state[102]=ir;
304 state[103]=jr;
305 state[104]=is;
306}
307
308
309void rlxd_reset(int state[])
310{
311 int k;
312
313 define_constants();
314
315 if (state[0]!=rlxd_size())
316 error(5);
317
318 for (k=0;k<96;k++)
319 {
320 if ((state[k+1]<0)||(state[k+1]>=167777216))
321 error(5);
322
323 x.num[k]=(float)(ldexp((double)(state[k+1]),-24));
324 }
325
326 if (((state[97]!=0)&&(state[97]!=1))||
327 ((state[98]!=0)&&(state[98]!=1))||
328 ((state[99]!=0)&&(state[99]!=1))||
329 ((state[100]!=0)&&(state[100]!=1)))
330 error(5);
331
332 carry.c1=(float)(ldexp((double)(state[97]),-24));
333 carry.c2=(float)(ldexp((double)(state[98]),-24));
334 carry.c3=(float)(ldexp((double)(state[99]),-24));
335 carry.c4=(float)(ldexp((double)(state[100]),-24));
336
337 pr=state[101];
338 ir=state[102];
339 jr=state[103];
340 is=state[104];
341 is_old=8*ir;
342 prm=pr%12;
343 init=1;
344
345 if (((pr!=202)&&(pr!=397))||
346 (ir<0)||(ir>11)||(jr<0)||(jr>11)||(jr!=((ir+7)%12))||
347 (is<0)||(is>91))
348 error(5);
349}
350
351#else
352
353#define BASE 0x1000000
354#define MASK 0xffffff
355
356typedef struct
357{
358 int c1,c2,c3,c4;
359} vec_t;
360
361typedef struct
362{
363 vec_t c1,c2;
364} dble_vec_t;
365
366static int init=0,pr,prm,ir,jr,is,is_old,next[96];
367static double one_bit;
368static vec_t carry;
369
370static union
371{
372 dble_vec_t vec[12];
373 int num[96];
374} x;
375
376#define STEP(pi,pj) \
377 d=(*pj).c1.c1-(*pi).c1.c1-carry.c1; \
378 (*pi).c2.c1+=(d<0); \
379 d+=BASE; \
380 (*pi).c1.c1=d&MASK; \
381 d=(*pj).c1.c2-(*pi).c1.c2-carry.c2; \
382 (*pi).c2.c2+=(d<0); \
383 d+=BASE; \
384 (*pi).c1.c2=d&MASK; \
385 d=(*pj).c1.c3-(*pi).c1.c3-carry.c3; \
386 (*pi).c2.c3+=(d<0); \
387 d+=BASE; \
388 (*pi).c1.c3=d&MASK; \
389 d=(*pj).c1.c4-(*pi).c1.c4-carry.c4; \
390 (*pi).c2.c4+=(d<0); \
391 d+=BASE; \
392 (*pi).c1.c4=d&MASK; \
393 d=(*pj).c2.c1-(*pi).c2.c1; \
394 carry.c1=(d<0); \
395 d+=BASE; \
396 (*pi).c2.c1=d&MASK; \
397 d=(*pj).c2.c2-(*pi).c2.c2; \
398 carry.c2=(d<0); \
399 d+=BASE; \
400 (*pi).c2.c2=d&MASK; \
401 d=(*pj).c2.c3-(*pi).c2.c3; \
402 carry.c3=(d<0); \
403 d+=BASE; \
404 (*pi).c2.c3=d&MASK; \
405 d=(*pj).c2.c4-(*pi).c2.c4; \
406 carry.c4=(d<0); \
407 d+=BASE; \
408 (*pi).c2.c4=d&MASK
409
410
411static void error(int no)
412{
413 switch(no)
414 {
415 case 0:
416 printf("Error in rlxd_init\n");
417 printf("Arithmetic on this machine is not suitable for ranlxd\n");
418 break;
419 case 1:
420 printf("Error in subroutine rlxd_init\n");
421 printf("Bad choice of luxury level (should be 1 or 2)\n");
422 break;
423 case 2:
424 printf("Error in subroutine rlxd_init\n");
425 printf("Bad choice of seed (should be between 1 and 2^31-1)\n");
426 break;
427 case 3:
428 printf("Error in rlxd_get\n");
429 printf("Undefined state (ranlxd is not initialized)\n");
430 break;
431 case 4:
432 printf("Error in rlxd_reset\n");
433 printf("Arithmetic on this machine is not suitable for ranlxd\n");
434 break;
435 case 5:
436 printf("Error in rlxd_reset\n");
437 printf("Unexpected input data\n");
438 break;
439 }
440 printf("Program aborted\n");
441 exit(0);
442}
443
444
445static void update()
446{
447 int k,kmax,d;
448 dble_vec_t *pmin,*pmax,*pi,*pj;
449
450 kmax=pr;
451 pmin=&x.vec[0];
452 pmax=pmin+12;
453 pi=&x.vec[ir];
454 pj=&x.vec[jr];
455
456 for (k=0;k<kmax;k++)
457 {
458 STEP(pi,pj);
459 pi+=1;
460 pj+=1;
461 if (pi==pmax)
462 pi=pmin;
463 if (pj==pmax)
464 pj=pmin;
465 }
466
467 ir+=prm;
468 jr+=prm;
469 if (ir>=12)
470 ir-=12;
471 if (jr>=12)
472 jr-=12;
473 is=8*ir;
474 is_old=is;
475}
476
477
478static void define_constants()
479{
480 int k;
481
482 one_bit=ldexp(1.0,-24);
483
484 for (k=0;k<96;k++)
485 {
486 next[k]=(k+1)%96;
487 if ((k%4)==3)
488 next[k]=(k+5)%96;
489 }
490}
491
492
493void rlxd_init(int level,int seed)
494{
495 int i,k,l;
496 int ibit,jbit,xbit[31];
497 int ix,iy;
498
499 if ((INT_MAX<2147483647)||(FLT_RADIX!=2)||(FLT_MANT_DIG<24)||
500 (DBL_MANT_DIG<48))
501 error(0);
502
503 define_constants();
504
505 if (level==1)
506 pr=202;
507 else if (level==2)
508 pr=397;
509 else
510 error(1);
511
512 i=seed;
513
514 for (k=0;k<31;k++)
515 {
516 xbit[k]=i%2;
517 i/=2;
518 }
519
520 if ((seed<=0)||(i!=0))
521 error(2);
522
523 ibit=0;
524 jbit=18;
525
526 for (i=0;i<4;i++)
527 {
528 for (k=0;k<24;k++)
529 {
530 ix=0;
531
532 for (l=0;l<24;l++)
533 {
534 iy=xbit[ibit];
535 ix=2*ix+iy;
536
537 xbit[ibit]=(xbit[ibit]+xbit[jbit])%2;
538 ibit=(ibit+1)%31;
539 jbit=(jbit+1)%31;
540 }
541
542 if ((k%4)!=i)
543 ix=16777215-ix;
544
545 x.num[4*k+i]=ix;
546 }
547 }
548
549 carry.c1=0;
550 carry.c2=0;
551 carry.c3=0;
552 carry.c4=0;
553
554 ir=0;
555 jr=7;
556 is=91;
557 is_old=0;
558 prm=pr%12;
559 init=1;
560}
561
562
563void ranlxd(double r[],int n)
564{
565 int k;
566
567 if (init==0)
568 rlxd_init(0,1);
569
570 for (k=0;k<n;k++)
571 {
572 is=next[is];
573 if (is==is_old)
574 update();
575 r[k]=one_bit*((double)(x.num[is+4])+one_bit*(double)(x.num[is]));
576 }
577}
578
579
580int rlxd_size(void)
581{
582 return(105);
583}
584
585
586void rlxd_get(int state[])
587{
588 int k;
589
590 if (init==0)
591 error(3);
592
593 state[0]=rlxd_size();
594
595 for (k=0;k<96;k++)
596 state[k+1]=x.num[k];
597
598 state[97]=carry.c1;
599 state[98]=carry.c2;
600 state[99]=carry.c3;
601 state[100]=carry.c4;
602
603 state[101]=pr;
604 state[102]=ir;
605 state[103]=jr;
606 state[104]=is;
607}
608
609
610void rlxd_reset(int state[])
611{
612 int k;
613
614 if ((INT_MAX<2147483647)||(FLT_RADIX!=2)||(FLT_MANT_DIG<24)||
615 (DBL_MANT_DIG<48))
616 error(4);
617
618 define_constants();
619
620 if (state[0]!=rlxd_size())
621 error(5);
622
623 for (k=0;k<96;k++)
624 {
625 if ((state[k+1]<0)||(state[k+1]>=167777216))
626 error(5);
627
628 x.num[k]=state[k+1];
629 }
630
631 if (((state[97]!=0)&&(state[97]!=1))||
632 ((state[98]!=0)&&(state[98]!=1))||
633 ((state[99]!=0)&&(state[99]!=1))||
634 ((state[100]!=0)&&(state[100]!=1)))
635 error(5);
636
637 carry.c1=state[97];
638 carry.c2=state[98];
639 carry.c3=state[99];
640 carry.c4=state[100];
641
642 pr=state[101];
643 ir=state[102];
644 jr=state[103];
645 is=state[104];
646 is_old=8*ir;
647 prm=pr%12;
648 init=1;
649
650 if (((pr!=202)&&(pr!=397))||
651 (ir<0)||(ir>11)||(jr<0)||(jr>11)||(jr!=((ir+7)%12))||
652 (is<0)||(is>91))
653 error(5);
654}
655
656#endif
657
n
const Int_t n
Definition: DataBase/tau_mode.c:65
one
const DifNumber one
ers::error
@ error
Definition: Core.h:24
rlxd_get
void rlxd_get(int state[])
Definition: ranlxd.c:586
rlxd_reset
void rlxd_reset(int state[])
Definition: ranlxd.c:610
vec
dble_vec_t vec[12]
Definition: ranlxd.c:372
rlxd_size
int rlxd_size(void)
Definition: ranlxd.c:580
STEP
#define STEP(pi, pj)
Definition: ranlxd.c:376
rlxd_init
void rlxd_init(int level, int seed)
Definition: ranlxd.c:493
num
int num[96]
Definition: ranlxd.c:373
ranlxd
void ranlxd(double r[], int n)
Definition: ranlxd.c:563
dble_vec_t
Definition: ranlxd.c:362
dble_vec_t::c1
vec_t c1
Definition: ranlxd.c:363
vec_t
Definition: ranlxd.c:357
vec_t::c2
int c2
Definition: ranlxd.c:358
vec_t::c1
int c1
Definition: ranlxd.c:358
vec_t::c3
int c3
Definition: ranlxd.c:358
vec_t::c4
int c4
Definition: ranlxd.c:358
c1
TCanvas * c1
Definition: tau_mode.c:75
pi
const float pi
Definition: vector3.h:133