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