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