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