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