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