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