| Geant4 Cross Reference |
1 //
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- RanluxppEngine ---
7 // helper implementation file
8 // -----------------------------------------------------------------------
9
10 #ifndef RANLUXPP_HELPERS_H
11 #define RANLUXPP_HELPERS_H
12
13 #include <cstdint>
14
15 /// Compute `a + b` and set `overflow` accordingly.
16 static inline uint64_t add_overflow(uint64_t a, uint64_t b,
17 unsigned &overflow) {
18 uint64_t add = a + b;
19 overflow = (add < a);
20 return add;
21 }
22
23 /// Compute `a + b` and increment `carry` if there was an overflow
24 static inline uint64_t add_carry(uint64_t a, uint64_t b, unsigned &carry) {
25 unsigned overflow;
26 uint64_t add = add_overflow(a, b, overflow);
27 // Do NOT branch on overflow to avoid jumping code, just add 0 if there was
28 // no overflow.
29 carry += overflow;
30 return add;
31 }
32
33 /// Compute `a - b` and set `overflow` accordingly
34 static inline uint64_t sub_overflow(uint64_t a, uint64_t b,
35 unsigned &overflow) {
36 uint64_t sub = a - b;
37 overflow = (sub > a);
38 return sub;
39 }
40
41 /// Compute `a - b` and increment `carry` if there was an overflow
42 static inline uint64_t sub_carry(uint64_t a, uint64_t b, unsigned &carry) {
43 unsigned overflow;
44 uint64_t sub = sub_overflow(a, b, overflow);
45 // Do NOT branch on overflow to avoid jumping code, just add 0 if there was
46 // no overflow.
47 carry += overflow;
48 return sub;
49 }
50
51 /// Update r = r - (t1 + t2) + (t3 + t2) * b ** 10
52 ///
53 /// This function also yields cbar = floor(r / m) as its return value (int64_t
54 /// because the value can be -1). With an initial value of r = t0, this can
55 /// be used for computing the remainder after division by m (see the function
56 /// mod_m in mulmod.h). The function to_ranlux passes r = 0 and uses only the
57 /// return value to obtain the decimal expansion after divison by m.
58 static inline int64_t compute_r(const uint64_t *upper, uint64_t *r) {
59 // Subtract t1 (24 * 24 = 576 bits)
60 unsigned carry = 0;
61 for (int i = 0; i < 9; i++) {
62 uint64_t r_i = r[i];
63 r_i = sub_overflow(r_i, carry, carry);
64
65 uint64_t t1_i = upper[i];
66 r_i = sub_carry(r_i, t1_i, carry);
67 r[i] = r_i;
68 }
69 int64_t c = -((int64_t)carry);
70
71 // Subtract t2 (only 240 bits, so need to extend)
72 carry = 0;
73 for (int i = 0; i < 9; i++) {
74 uint64_t r_i = r[i];
75 r_i = sub_overflow(r_i, carry, carry);
76
77 uint64_t t2_bits = 0;
78 if (i < 4) {
79 t2_bits += upper[i + 5] >> 16;
80 if (i < 3) {
81 t2_bits += upper[i + 6] << 48;
82 }
83 }
84 r_i = sub_carry(r_i, t2_bits, carry);
85 r[i] = r_i;
86 }
87 c -= carry;
88
89 // r += (t3 + t2) * 2 ** 240
90 carry = 0;
91 {
92 uint64_t r_3 = r[3];
93 // 16 upper bits
94 uint64_t t2_bits = (upper[5] >> 16) << 48;
95 uint64_t t3_bits = (upper[0] << 48);
96
97 r_3 = add_carry(r_3, t2_bits, carry);
98 r_3 = add_carry(r_3, t3_bits, carry);
99
100 r[3] = r_3;
101 }
102 for (int i = 0; i < 3; i++) {
103 uint64_t r_i = r[i + 4];
104 r_i = add_overflow(r_i, carry, carry);
105
106 uint64_t t2_bits = (upper[5 + i] >> 32) + (upper[6 + i] << 32);
107 uint64_t t3_bits = (upper[i] >> 16) + (upper[1 + i] << 48);
108
109 r_i = add_carry(r_i, t2_bits, carry);
110 r_i = add_carry(r_i, t3_bits, carry);
111
112 r[i + 4] = r_i;
113 }
114 {
115 uint64_t r_7 = r[7];
116 r_7 = add_overflow(r_7, carry, carry);
117
118 uint64_t t2_bits = (upper[8] >> 32);
119 uint64_t t3_bits = (upper[3] >> 16) + (upper[4] << 48);
120
121 r_7 = add_carry(r_7, t2_bits, carry);
122 r_7 = add_carry(r_7, t3_bits, carry);
123
124 r[7] = r_7;
125 }
126 {
127 uint64_t r_8 = r[8];
128 r_8 = add_overflow(r_8, carry, carry);
129
130 uint64_t t3_bits = (upper[4] >> 16) + (upper[5] << 48);
131
132 r_8 = add_carry(r_8, t3_bits, carry);
133
134 r[8] = r_8;
135 }
136 c += carry;
137
138 // c = floor(r / 2 ** 576) has been computed along the way via the carry
139 // flags. Now if c = 0 and the value currently stored in r is greater or
140 // equal to m, we need cbar = 1 and subtract m, otherwise cbar = c. The
141 // value currently in r is greater or equal to m, if and only if one of
142 // the last 240 bits is set and the upper bits are all set.
143 bool greater_m = r[0] | r[1] | r[2] | (r[3] & 0x0000ffffffffffff);
144 greater_m &= (r[3] >> 48) == 0xffff;
145 for (int i = 4; i < 9; i++) {
146 greater_m &= (r[i] == UINT64_MAX);
147 }
148 return c + (c == 0 && greater_m);
149 }
150
151 #endif
152