| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4Log 26 // G4Log
27 // 27 //
28 // Class description: 28 // Class description:
29 // 29 //
30 // The basic idea is to exploit Pade polynomia 30 // The basic idea is to exploit Pade polynomials.
31 // A lot of ideas were inspired by the cephes 31 // A lot of ideas were inspired by the cephes math library
32 // (by Stephen L. Moshier moshier@na-net.ornl. 32 // (by Stephen L. Moshier moshier@na-net.ornl.gov) as well as actual code.
33 // The Cephes library can be found here: http 33 // The Cephes library can be found here: http://www.netlib.org/cephes/
34 // Code and algorithms for G4Exp have been ext 34 // Code and algorithms for G4Exp have been extracted and adapted for Geant4
35 // from the original implementation in the VDT 35 // from the original implementation in the VDT mathematical library
36 // (https://svnweb.cern.ch/trac/vdt), version 36 // (https://svnweb.cern.ch/trac/vdt), version 0.3.7.
37 37
38 // Original implementation created on: Jun 23, 38 // Original implementation created on: Jun 23, 2012
39 // Author: Danilo Piparo, Thomas Hauth, V 39 // Author: Danilo Piparo, Thomas Hauth, Vincenzo Innocente
40 // 40 //
41 // ------------------------------------------- 41 // --------------------------------------------------------------------
42 /* 42 /*
43 * VDT is free software: you can redistribute 43 * VDT is free software: you can redistribute it and/or modify
44 * it under the terms of the GNU Lesser Public 44 * it under the terms of the GNU Lesser Public License as published by
45 * the Free Software Foundation, either versio 45 * the Free Software Foundation, either version 3 of the License, or
46 * (at your option) any later version. 46 * (at your option) any later version.
47 * 47 *
48 * This program is distributed in the hope tha 48 * This program is distributed in the hope that it will be useful,
49 * but WITHOUT ANY WARRANTY; without even the 49 * but WITHOUT ANY WARRANTY; without even the implied warranty of
50 * MERCHANTABILITY or FITNESS FOR A PARTICULAR 50 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
51 * GNU Lesser Public License for more details. 51 * GNU Lesser Public License for more details.
52 * 52 *
53 * You should have received a copy of the GNU 53 * You should have received a copy of the GNU Lesser Public License
54 * along with this program. If not, see <http 54 * along with this program. If not, see <http://www.gnu.org/licenses/>.
55 */ 55 */
56 // ------------------------------------------- 56 // --------------------------------------------------------------------
57 #ifndef G4Log_hh 57 #ifndef G4Log_hh
58 #define G4Log_hh 1 58 #define G4Log_hh 1
59 59
60 #ifdef WIN32 60 #ifdef WIN32
61 61
62 # define G4Log std::log 62 # define G4Log std::log
63 63
64 #else 64 #else
65 65
66 # include "G4Types.hh" 66 # include "G4Types.hh"
67 <<
68 # include <cstdint> <<
69 # include <limits> 67 # include <limits>
>> 68 # include <stdint.h>
70 69
71 // local namespace for the constants/functions 70 // local namespace for the constants/functions which are necessary only here
72 // 71 //
73 namespace G4LogConsts 72 namespace G4LogConsts
74 { 73 {
75 const G4double LOG_UPPER_LIMIT = 1e307; 74 const G4double LOG_UPPER_LIMIT = 1e307;
76 const G4double LOG_LOWER_LIMIT = 0; 75 const G4double LOG_LOWER_LIMIT = 0;
77 76
78 const G4double SQRTH = 0.707106781186547524 77 const G4double SQRTH = 0.70710678118654752440;
79 const G4float MAXNUMF = 3.402823466385288598 78 const G4float MAXNUMF = 3.4028234663852885981170418348451692544e38f;
80 79
81 //------------------------------------------ 80 //----------------------------------------------------------------------------
82 // Used to switch between different type of 81 // Used to switch between different type of interpretations of the data
83 // (64 bits) 82 // (64 bits)
84 // 83 //
85 union ieee754 84 union ieee754
86 { 85 {
87 ieee754()= default; << 86 ieee754(){};
88 ieee754(G4double thed) { d = thed; }; 87 ieee754(G4double thed) { d = thed; };
89 ieee754(uint64_t thell) { ll = thell; }; 88 ieee754(uint64_t thell) { ll = thell; };
90 ieee754(G4float thef) { f[0] = thef; }; 89 ieee754(G4float thef) { f[0] = thef; };
91 ieee754(uint32_t thei) { i[0] = thei; }; 90 ieee754(uint32_t thei) { i[0] = thei; };
92 G4double d; 91 G4double d;
93 G4float f[2]; 92 G4float f[2];
94 uint32_t i[2]; 93 uint32_t i[2];
95 uint64_t ll; 94 uint64_t ll;
96 uint16_t s[4]; 95 uint16_t s[4];
97 }; 96 };
98 97
99 inline G4double get_log_px(const G4double x) 98 inline G4double get_log_px(const G4double x)
100 { 99 {
101 const G4double PX1log = 1.0187566380458093 100 const G4double PX1log = 1.01875663804580931796E-4;
102 const G4double PX2log = 4.9749499497674700 101 const G4double PX2log = 4.97494994976747001425E-1;
103 const G4double PX3log = 4.7057911987888172 102 const G4double PX3log = 4.70579119878881725854E0;
104 const G4double PX4log = 1.4498922534161093 103 const G4double PX4log = 1.44989225341610930846E1;
105 const G4double PX5log = 1.7936867850781981 104 const G4double PX5log = 1.79368678507819816313E1;
106 const G4double PX6log = 7.7083873375588539 105 const G4double PX6log = 7.70838733755885391666E0;
107 106
108 G4double px = PX1log; 107 G4double px = PX1log;
109 px *= x; 108 px *= x;
110 px += PX2log; 109 px += PX2log;
111 px *= x; 110 px *= x;
112 px += PX3log; 111 px += PX3log;
113 px *= x; 112 px *= x;
114 px += PX4log; 113 px += PX4log;
115 px *= x; 114 px *= x;
116 px += PX5log; 115 px += PX5log;
117 px *= x; 116 px *= x;
118 px += PX6log; 117 px += PX6log;
119 return px; 118 return px;
120 } 119 }
121 120
122 inline G4double get_log_qx(const G4double x) 121 inline G4double get_log_qx(const G4double x)
123 { 122 {
124 const G4double QX1log = 1.1287358718916745 123 const G4double QX1log = 1.12873587189167450590E1;
125 const G4double QX2log = 4.5227914583753222 124 const G4double QX2log = 4.52279145837532221105E1;
126 const G4double QX3log = 8.2987526691277660 125 const G4double QX3log = 8.29875266912776603211E1;
127 const G4double QX4log = 7.1154475061856389 126 const G4double QX4log = 7.11544750618563894466E1;
128 const G4double QX5log = 2.3125162012676534 127 const G4double QX5log = 2.31251620126765340583E1;
129 128
130 G4double qx = x; 129 G4double qx = x;
131 qx += QX1log; 130 qx += QX1log;
132 qx *= x; 131 qx *= x;
133 qx += QX2log; 132 qx += QX2log;
134 qx *= x; 133 qx *= x;
135 qx += QX3log; 134 qx += QX3log;
136 qx *= x; 135 qx *= x;
137 qx += QX4log; 136 qx += QX4log;
138 qx *= x; 137 qx *= x;
139 qx += QX5log; 138 qx += QX5log;
140 return qx; 139 return qx;
141 } 140 }
142 141
143 //------------------------------------------ 142 //----------------------------------------------------------------------------
144 // Converts a double to an unsigned long lon 143 // Converts a double to an unsigned long long
145 // 144 //
146 inline uint64_t dp2uint64(G4double x) 145 inline uint64_t dp2uint64(G4double x)
147 { 146 {
148 ieee754 tmp; 147 ieee754 tmp;
149 tmp.d = x; 148 tmp.d = x;
150 return tmp.ll; 149 return tmp.ll;
151 } 150 }
152 151
153 //------------------------------------------ 152 //----------------------------------------------------------------------------
154 // Converts an unsigned long long to a doubl 153 // Converts an unsigned long long to a double
155 // 154 //
156 inline G4double uint642dp(uint64_t ll) 155 inline G4double uint642dp(uint64_t ll)
157 { 156 {
158 ieee754 tmp; 157 ieee754 tmp;
159 tmp.ll = ll; 158 tmp.ll = ll;
160 return tmp.d; 159 return tmp.d;
161 } 160 }
162 161
163 //------------------------------------------ 162 //----------------------------------------------------------------------------
164 // Converts an int to a float 163 // Converts an int to a float
165 // 164 //
166 inline G4float uint322sp(G4int x) 165 inline G4float uint322sp(G4int x)
167 { 166 {
168 ieee754 tmp; 167 ieee754 tmp;
169 tmp.i[0] = x; 168 tmp.i[0] = x;
170 return tmp.f[0]; 169 return tmp.f[0];
171 } 170 }
172 171
173 //------------------------------------------ 172 //----------------------------------------------------------------------------
174 // Converts a float to an int 173 // Converts a float to an int
175 // 174 //
176 inline uint32_t sp2uint32(G4float x) 175 inline uint32_t sp2uint32(G4float x)
177 { 176 {
178 ieee754 tmp; 177 ieee754 tmp;
179 tmp.f[0] = x; 178 tmp.f[0] = x;
180 return tmp.i[0]; 179 return tmp.i[0];
181 } 180 }
182 181
183 //------------------------------------------ 182 //----------------------------------------------------------------------------
184 /// Like frexp but vectorising and the expon 183 /// Like frexp but vectorising and the exponent is a double.
185 inline G4double getMantExponent(const G4doub 184 inline G4double getMantExponent(const G4double x, G4double& fe)
186 { 185 {
187 uint64_t n = dp2uint64(x); 186 uint64_t n = dp2uint64(x);
188 187
189 // Shift to the right up to the beginning 188 // Shift to the right up to the beginning of the exponent.
190 // Then with a mask, cut off the sign bit 189 // Then with a mask, cut off the sign bit
191 uint64_t le = (n >> 52); 190 uint64_t le = (n >> 52);
192 191
193 // chop the head of the number: an int con 192 // chop the head of the number: an int contains more than 11 bits (32)
194 int32_t e = 193 int32_t e =
195 (int32_t)le; // This is important since << 194 le; // This is important since sums on uint64_t do not vectorise
196 fe = e - 1023; 195 fe = e - 1023;
197 196
198 // This puts to 11 zeroes the exponent 197 // This puts to 11 zeroes the exponent
199 n &= 0x800FFFFFFFFFFFFFULL; 198 n &= 0x800FFFFFFFFFFFFFULL;
200 // build a mask which is 0.5, i.e. an expo 199 // build a mask which is 0.5, i.e. an exponent equal to 1022
201 // which means *2, see the above +1. 200 // which means *2, see the above +1.
202 const uint64_t p05 = 0x3FE0000000000000ULL 201 const uint64_t p05 = 0x3FE0000000000000ULL; // dp2uint64(0.5);
203 n |= p05; 202 n |= p05;
204 203
205 return uint642dp(n); 204 return uint642dp(n);
206 } 205 }
207 206
208 //------------------------------------------ 207 //----------------------------------------------------------------------------
209 /// Like frexp but vectorising and the expon 208 /// Like frexp but vectorising and the exponent is a float.
210 inline G4float getMantExponentf(const G4floa 209 inline G4float getMantExponentf(const G4float x, G4float& fe)
211 { 210 {
212 uint32_t n = sp2uint32(x); 211 uint32_t n = sp2uint32(x);
213 int32_t e = (n >> 23) - 127; 212 int32_t e = (n >> 23) - 127;
214 fe = e; 213 fe = e;
215 214
216 // fractional part 215 // fractional part
217 const uint32_t p05f = 0x3f000000; // //sp 216 const uint32_t p05f = 0x3f000000; // //sp2uint32(0.5);
218 n &= 0x807fffff; // ~0x7 217 n &= 0x807fffff; // ~0x7f800000;
219 n |= p05f; 218 n |= p05f;
220 219
221 return uint322sp(n); 220 return uint322sp(n);
222 } 221 }
223 } // namespace G4LogConsts 222 } // namespace G4LogConsts
224 223
225 // Log double precision ---------------------- 224 // Log double precision --------------------------------------------------------
226 225
227 inline G4double G4Log(G4double x) 226 inline G4double G4Log(G4double x)
228 { 227 {
229 const G4double original_x = x; 228 const G4double original_x = x;
230 229
231 /* separate mantissa from exponent */ 230 /* separate mantissa from exponent */
232 G4double fe; 231 G4double fe;
233 x = G4LogConsts::getMantExponent(x, fe); 232 x = G4LogConsts::getMantExponent(x, fe);
234 233
235 // blending 234 // blending
236 x > G4LogConsts::SQRTH ? fe += 1. : x += x; 235 x > G4LogConsts::SQRTH ? fe += 1. : x += x;
237 x -= 1.0; 236 x -= 1.0;
238 237
239 /* rational form */ 238 /* rational form */
240 G4double px = G4LogConsts::get_log_px(x); 239 G4double px = G4LogConsts::get_log_px(x);
241 240
242 // for the final formula 241 // for the final formula
243 const G4double x2 = x * x; 242 const G4double x2 = x * x;
244 px *= x; 243 px *= x;
245 px *= x2; 244 px *= x2;
246 245
247 const G4double qx = G4LogConsts::get_log_qx( 246 const G4double qx = G4LogConsts::get_log_qx(x);
248 247
249 G4double res = px / qx; 248 G4double res = px / qx;
250 249
251 res -= fe * 2.121944400546905827679e-4; 250 res -= fe * 2.121944400546905827679e-4;
252 res -= 0.5 * x2; 251 res -= 0.5 * x2;
253 252
254 res = x + res; 253 res = x + res;
255 res += fe * 0.693359375; 254 res += fe * 0.693359375;
256 255
257 if(original_x > G4LogConsts::LOG_UPPER_LIMIT 256 if(original_x > G4LogConsts::LOG_UPPER_LIMIT)
258 res = std::numeric_limits<G4double>::infin 257 res = std::numeric_limits<G4double>::infinity();
259 if(original_x < G4LogConsts::LOG_LOWER_LIMIT 258 if(original_x < G4LogConsts::LOG_LOWER_LIMIT) // THIS IS NAN!
260 res = -std::numeric_limits<G4double>::quie 259 res = -std::numeric_limits<G4double>::quiet_NaN();
261 260
262 return res; 261 return res;
263 } 262 }
264 263
265 // Log single precision ---------------------- 264 // Log single precision --------------------------------------------------------
266 265
267 namespace G4LogConsts 266 namespace G4LogConsts
268 { 267 {
269 const G4float LOGF_UPPER_LIMIT = MAXNUMF; 268 const G4float LOGF_UPPER_LIMIT = MAXNUMF;
270 const G4float LOGF_LOWER_LIMIT = 0; 269 const G4float LOGF_LOWER_LIMIT = 0;
271 270
272 const G4float PX1logf = 7.0376836292E-2f; 271 const G4float PX1logf = 7.0376836292E-2f;
273 const G4float PX2logf = -1.1514610310E-1f; 272 const G4float PX2logf = -1.1514610310E-1f;
274 const G4float PX3logf = 1.1676998740E-1f; 273 const G4float PX3logf = 1.1676998740E-1f;
275 const G4float PX4logf = -1.2420140846E-1f; 274 const G4float PX4logf = -1.2420140846E-1f;
276 const G4float PX5logf = 1.4249322787E-1f; 275 const G4float PX5logf = 1.4249322787E-1f;
277 const G4float PX6logf = -1.6668057665E-1f; 276 const G4float PX6logf = -1.6668057665E-1f;
278 const G4float PX7logf = 2.0000714765E-1f; 277 const G4float PX7logf = 2.0000714765E-1f;
279 const G4float PX8logf = -2.4999993993E-1f; 278 const G4float PX8logf = -2.4999993993E-1f;
280 const G4float PX9logf = 3.3333331174E-1f; 279 const G4float PX9logf = 3.3333331174E-1f;
281 280
282 inline G4float get_log_poly(const G4float x) 281 inline G4float get_log_poly(const G4float x)
283 { 282 {
284 G4float y = x * PX1logf; 283 G4float y = x * PX1logf;
285 y += PX2logf; 284 y += PX2logf;
286 y *= x; 285 y *= x;
287 y += PX3logf; 286 y += PX3logf;
288 y *= x; 287 y *= x;
289 y += PX4logf; 288 y += PX4logf;
290 y *= x; 289 y *= x;
291 y += PX5logf; 290 y += PX5logf;
292 y *= x; 291 y *= x;
293 y += PX6logf; 292 y += PX6logf;
294 y *= x; 293 y *= x;
295 y += PX7logf; 294 y += PX7logf;
296 y *= x; 295 y *= x;
297 y += PX8logf; 296 y += PX8logf;
298 y *= x; 297 y *= x;
299 y += PX9logf; 298 y += PX9logf;
300 return y; 299 return y;
301 } 300 }
302 301
303 const G4float SQRTHF = 0.707106781186547524f 302 const G4float SQRTHF = 0.707106781186547524f;
304 } // namespace G4LogConsts 303 } // namespace G4LogConsts
305 304
306 // Log single precision ---------------------- 305 // Log single precision --------------------------------------------------------
307 306
308 inline G4float G4Logf(G4float x) 307 inline G4float G4Logf(G4float x)
309 { 308 {
310 const G4float original_x = x; 309 const G4float original_x = x;
311 310
312 G4float fe; 311 G4float fe;
313 x = G4LogConsts::getMantExponentf(x, fe); 312 x = G4LogConsts::getMantExponentf(x, fe);
314 313
315 x > G4LogConsts::SQRTHF ? fe += 1.f : x += x 314 x > G4LogConsts::SQRTHF ? fe += 1.f : x += x;
316 x -= 1.0f; 315 x -= 1.0f;
317 316
318 const G4float x2 = x * x; 317 const G4float x2 = x * x;
319 318
320 G4float res = G4LogConsts::get_log_poly(x); 319 G4float res = G4LogConsts::get_log_poly(x);
321 res *= x2 * x; 320 res *= x2 * x;
322 321
323 res += -2.12194440e-4f * fe; 322 res += -2.12194440e-4f * fe;
324 res += -0.5f * x2; 323 res += -0.5f * x2;
325 324
326 res = x + res; 325 res = x + res;
327 326
328 res += 0.693359375f * fe; 327 res += 0.693359375f * fe;
329 328
330 if(original_x > G4LogConsts::LOGF_UPPER_LIMI 329 if(original_x > G4LogConsts::LOGF_UPPER_LIMIT)
331 res = std::numeric_limits<G4float>::infini 330 res = std::numeric_limits<G4float>::infinity();
332 if(original_x < G4LogConsts::LOGF_LOWER_LIMI 331 if(original_x < G4LogConsts::LOGF_LOWER_LIMIT)
333 res = -std::numeric_limits<G4float>::quiet 332 res = -std::numeric_limits<G4float>::quiet_NaN();
334 333
335 return res; 334 return res;
336 } 335 }
337 336
338 //-------------------------------------------- 337 //------------------------------------------------------------------------------
339 338
340 void logv(const uint32_t size, G4double const* 339 void logv(const uint32_t size, G4double const* __restrict__ iarray,
341 G4double* __restrict__ oarray); 340 G4double* __restrict__ oarray);
342 void G4Logv(const uint32_t size, G4double cons 341 void G4Logv(const uint32_t size, G4double const* __restrict__ iarray,
343 G4double* __restrict__ oarray); 342 G4double* __restrict__ oarray);
344 void logfv(const uint32_t size, G4float const* 343 void logfv(const uint32_t size, G4float const* __restrict__ iarray,
345 G4float* __restrict__ oarray); 344 G4float* __restrict__ oarray);
346 void G4Logfv(const uint32_t size, G4float cons 345 void G4Logfv(const uint32_t size, G4float const* __restrict__ iarray,
347 G4float* __restrict__ oarray); 346 G4float* __restrict__ oarray);
348 347
349 #endif /* WIN32 */ 348 #endif /* WIN32 */
350 349
351 #endif /* LOG_H_ */ 350 #endif /* LOG_H_ */
352 351