| Geant4 Cross Reference |
1 // 1 //
2 // -*- C++ -*- 2 // -*- C++ -*-
3 // 3 //
4 // ------------------------------------------- 4 // -----------------------------------------------------------------------
5 // HEP Random 5 // HEP Random
6 // --- RanluxppEngine -- 6 // --- RanluxppEngine ---
7 // helper implementation f 7 // helper implementation file
8 // ------------------------------------------- 8 // -----------------------------------------------------------------------
9 9
10 #ifndef RANLUXPP_RANLUX_LCG_H 10 #ifndef RANLUXPP_RANLUX_LCG_H
11 #define RANLUXPP_RANLUX_LCG_H 11 #define RANLUXPP_RANLUX_LCG_H
12 12
13 #include "helpers.h" 13 #include "helpers.h"
14 14
15 #include <cstdint> 15 #include <cstdint>
16 16
17 /// Convert RANLUX numbers to an LCG state 17 /// Convert RANLUX numbers to an LCG state
18 /// 18 ///
19 /// \param[in] ranlux the RANLUX numbers as 57 19 /// \param[in] ranlux the RANLUX numbers as 576 bits
20 /// \param[out] lcg the 576 bits of the LCG st 20 /// \param[out] lcg the 576 bits of the LCG state, smaller than m
21 /// \param[in] c the carry bit of the RANLUX s 21 /// \param[in] c the carry bit of the RANLUX state
22 /// 22 ///
23 /// \f$ m = 2^{576} - 2^{240} + 1 \f$ 23 /// \f$ m = 2^{576} - 2^{240} + 1 \f$
24 static void to_lcg(const uint64_t *ranlux, uns 24 static void to_lcg(const uint64_t *ranlux, unsigned c, uint64_t *lcg) {
25 unsigned carry = 0; 25 unsigned carry = 0;
26 // Subtract the final 240 bits. 26 // Subtract the final 240 bits.
27 for (int i = 0; i < 9; i++) { 27 for (int i = 0; i < 9; i++) {
28 uint64_t ranlux_i = ranlux[i]; 28 uint64_t ranlux_i = ranlux[i];
29 uint64_t lcg_i = sub_overflow(ranlux_i, ca 29 uint64_t lcg_i = sub_overflow(ranlux_i, carry, carry);
30 30
31 uint64_t bits = 0; 31 uint64_t bits = 0;
32 if (i < 4) { 32 if (i < 4) {
33 bits += ranlux[i + 5] >> 16; 33 bits += ranlux[i + 5] >> 16;
34 if (i < 3) { 34 if (i < 3) {
35 bits += ranlux[i + 6] << 48; 35 bits += ranlux[i + 6] << 48;
36 } 36 }
37 } 37 }
38 lcg_i = sub_carry(lcg_i, bits, carry); 38 lcg_i = sub_carry(lcg_i, bits, carry);
39 lcg[i] = lcg_i; 39 lcg[i] = lcg_i;
40 } 40 }
41 41
42 // Add and propagate the carry bit. 42 // Add and propagate the carry bit.
43 for (int i = 0; i < 9; i++) { 43 for (int i = 0; i < 9; i++) {
44 lcg[i] = add_overflow(lcg[i], c, c); 44 lcg[i] = add_overflow(lcg[i], c, c);
45 } 45 }
46 } 46 }
47 47
48 /// Convert an LCG state to RANLUX numbers 48 /// Convert an LCG state to RANLUX numbers
49 /// 49 ///
50 /// \param[in] lcg the 576 bits of the LCG sta 50 /// \param[in] lcg the 576 bits of the LCG state, must be smaller than m
51 /// \param[out] ranlux the RANLUX numbers as 5 51 /// \param[out] ranlux the RANLUX numbers as 576 bits
52 /// \param[out] c the carry bit of the RANLUX 52 /// \param[out] c the carry bit of the RANLUX state
53 /// 53 ///
54 /// \f$ m = 2^{576} - 2^{240} + 1 \f$ 54 /// \f$ m = 2^{576} - 2^{240} + 1 \f$
55 static void to_ranlux(const uint64_t *lcg, uin 55 static void to_ranlux(const uint64_t *lcg, uint64_t *ranlux, unsigned &c_out) {
56 uint64_t r[9] = {0}; 56 uint64_t r[9] = {0};
57 int64_t c = compute_r(lcg, r); 57 int64_t c = compute_r(lcg, r);
58 58
59 // ranlux = t1 + t2 + c 59 // ranlux = t1 + t2 + c
60 unsigned carry = 0; 60 unsigned carry = 0;
61 for (int i = 0; i < 9; i++) { 61 for (int i = 0; i < 9; i++) {
62 uint64_t in_i = lcg[i]; 62 uint64_t in_i = lcg[i];
63 uint64_t tmp_i = add_overflow(in_i, carry, 63 uint64_t tmp_i = add_overflow(in_i, carry, carry);
64 64
65 uint64_t bits = 0; 65 uint64_t bits = 0;
66 if (i < 4) { 66 if (i < 4) {
67 bits += lcg[i + 5] >> 16; 67 bits += lcg[i + 5] >> 16;
68 if (i < 3) { 68 if (i < 3) {
69 bits += lcg[i + 6] << 48; 69 bits += lcg[i + 6] << 48;
70 } 70 }
71 } 71 }
72 tmp_i = add_carry(tmp_i, bits, carry); 72 tmp_i = add_carry(tmp_i, bits, carry);
73 ranlux[i] = tmp_i; 73 ranlux[i] = tmp_i;
74 } 74 }
75 75
76 // If c = -1, we need to add it to all compo 76 // If c = -1, we need to add it to all components.
77 int64_t c1 = c >> 1; 77 int64_t c1 = c >> 1;
78 ranlux[0] = add_overflow(ranlux[0], c, carry 78 ranlux[0] = add_overflow(ranlux[0], c, carry);
79 for (int i = 1; i < 9; i++) { 79 for (int i = 1; i < 9; i++) {
80 uint64_t ranlux_i = ranlux[i]; 80 uint64_t ranlux_i = ranlux[i];
81 ranlux_i = add_overflow(ranlux_i, carry, c 81 ranlux_i = add_overflow(ranlux_i, carry, carry);
82 ranlux_i = add_carry(ranlux_i, c1, carry); 82 ranlux_i = add_carry(ranlux_i, c1, carry);
83 } 83 }
84 84
85 c_out = carry; 85 c_out = carry;
86 } 86 }
87 87
88 #endif 88 #endif
89 89