Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
ranlux_lcg.h
Go to the documentation of this file.
1// @(#)root/mathcore:$Id$
2// Author: Jonas Hahnfeld 05/2021
3
4/*************************************************************************
5 * Copyright (C) 1995-2021, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12#ifndef RANLUXPP_RANLUX_LCG_H
13#define RANLUXPP_RANLUX_LCG_H
14
15#include "helpers.h"
16
17#include <cstdint>
18
19/// Convert RANLUX numbers to an LCG state
20///
21/// \param[in] ranlux the RANLUX numbers as 576 bits
22/// \param[out] lcg the 576 bits of the LCG state, smaller than m
23/// \param[in] c the carry bit of the RANLUX state
24///
25/// \f$ m = 2^{576} - 2^{240} + 1 \f$
26static void to_lcg(const uint64_t *ranlux, unsigned c, uint64_t *lcg)
27{
28 unsigned carry = 0;
29 // Subtract the final 240 bits.
30 for (int i = 0; i < 9; i++) {
31 uint64_t ranlux_i = ranlux[i];
32 uint64_t lcg_i = sub_overflow(ranlux_i, carry, carry);
33
34 uint64_t bits = 0;
35 if (i < 4) {
36 bits += ranlux[i + 5] >> 16;
37 if (i < 3) {
38 bits += ranlux[i + 6] << 48;
39 }
40 }
41 lcg_i = sub_carry(lcg_i, bits, carry);
42 lcg[i] = lcg_i;
43 }
44
45 // Add and propagate the carry bit.
46 for (int i = 0; i < 9; i++) {
47 lcg[i] = add_overflow(lcg[i], c, c);
48 }
49}
50
51/// Convert an LCG state to RANLUX numbers
52///
53/// \param[in] lcg the 576 bits of the LCG state, must be smaller than m
54/// \param[out] ranlux the RANLUX numbers as 576 bits
55/// \param[out] c the carry bit of the RANLUX state
56///
57/// \f$ m = 2^{576} - 2^{240} + 1 \f$
58static void to_ranlux(const uint64_t *lcg, uint64_t *ranlux, unsigned &c_out)
59{
60 uint64_t r[9] = {0};
61 int64_t c = compute_r(lcg, r);
62
63 // ranlux = t1 + t2 + c
64 unsigned carry = 0;
65 for (int i = 0; i < 9; i++) {
66 uint64_t in_i = lcg[i];
67 uint64_t tmp_i = add_overflow(in_i, carry, carry);
68
69 uint64_t bits = 0;
70 if (i < 4) {
71 bits += lcg[i + 5] >> 16;
72 if (i < 3) {
73 bits += lcg[i + 6] << 48;
74 }
75 }
76 tmp_i = add_carry(tmp_i, bits, carry);
77 ranlux[i] = tmp_i;
78 }
79
80 // If c = -1, we need to add it to all components.
81 int64_t c1 = c >> 1;
82 ranlux[0] = add_overflow(ranlux[0], c, carry);
83 for (int i = 1; i < 9; i++) {
84 uint64_t ranlux_i = ranlux[i];
85 ranlux_i = add_overflow(ranlux_i, carry, carry);
86 ranlux_i = add_carry(ranlux_i, c1, carry);
87 }
88
89 c_out = carry;
90}
91
92#endif
ROOT::R::TRInterface & r
Definition Object.C:4
#define c(i)
Definition RSha256.hxx:101
static uint64_t sub_overflow(uint64_t a, uint64_t b, unsigned &overflow)
Compute a - b and set overflow accordingly.
Definition helpers.h:37
static uint64_t sub_carry(uint64_t a, uint64_t b, unsigned &carry)
Compute a - b and increment carry if there was an overflow.
Definition helpers.h:45
static uint64_t add_carry(uint64_t a, uint64_t b, unsigned &carry)
Compute a + b and increment carry if there was an overflow.
Definition helpers.h:26
static int64_t compute_r(const uint64_t *upper, uint64_t *r)
Update r = r - (t1 + t2) + (t3 + t2) * b ** 10.
Definition helpers.h:62
static uint64_t add_overflow(uint64_t a, uint64_t b, unsigned &overflow)
Compute a + b and set overflow accordingly.
Definition helpers.h:18
return c1
Definition legend1.C:41
static void to_lcg(const uint64_t *ranlux, unsigned c, uint64_t *lcg)
Convert RANLUX numbers to an LCG state.
Definition ranlux_lcg.h:26
static void to_ranlux(const uint64_t *lcg, uint64_t *ranlux, unsigned &c_out)
Convert an LCG state to RANLUX numbers.
Definition ranlux_lcg.h:58