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