Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/global/management/include/G4Log.hh

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 
 68 #  include <cstdint>
 69 #  include <limits>
 70 
 71 // local namespace for the constants/functions which are necessary only here
 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.70710678118654752440;
 79   const G4float MAXNUMF = 3.4028234663852885981170418348451692544e38f;
 80 
 81   //----------------------------------------------------------------------------
 82   // Used to switch between different type of interpretations of the data
 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.01875663804580931796E-4;
102     const G4double PX2log = 4.97494994976747001425E-1;
103     const G4double PX3log = 4.70579119878881725854E0;
104     const G4double PX4log = 1.44989225341610930846E1;
105     const G4double PX5log = 1.79368678507819816313E1;
106     const G4double PX6log = 7.70838733755885391666E0;
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.12873587189167450590E1;
125     const G4double QX2log = 4.52279145837532221105E1;
126     const G4double QX3log = 8.29875266912776603211E1;
127     const G4double QX4log = 7.11544750618563894466E1;
128     const G4double QX5log = 2.31251620126765340583E1;
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 long
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 double
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 exponent is a double.
185   inline G4double getMantExponent(const G4double x, G4double& fe)
186   {
187     uint64_t n = dp2uint64(x);
188 
189     // Shift to the right up to the beginning of the exponent.
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 contains more than 11 bits (32)
194     int32_t e =
195       (int32_t)le;  // This is important since sums on uint64_t do not vectorise
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 exponent equal to 1022
201     // which means *2, see the above +1.
202     const uint64_t p05 = 0x3FE0000000000000ULL;  // dp2uint64(0.5);
203     n |= p05;
204 
205     return uint642dp(n);
206   }
207 
208   //----------------------------------------------------------------------------
209   /// Like frexp but vectorising and the exponent is a float.
210   inline G4float getMantExponentf(const G4float x, G4float& fe)
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;  // //sp2uint32(0.5);
218     n &= 0x807fffff;                   // ~0x7f800000;
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(x);
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>::infinity();
259   if(original_x < G4LogConsts::LOG_LOWER_LIMIT)  // THIS IS NAN!
260     res = -std::numeric_limits<G4double>::quiet_NaN();
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_LIMIT)
331     res = std::numeric_limits<G4float>::infinity();
332   if(original_x < G4LogConsts::LOGF_LOWER_LIMIT)
333     res = -std::numeric_limits<G4float>::quiet_NaN();
334 
335   return res;
336 }
337 
338 //------------------------------------------------------------------------------
339 
340 void logv(const uint32_t size, G4double const* __restrict__ iarray,
341           G4double* __restrict__ oarray);
342 void G4Logv(const uint32_t size, G4double const* __restrict__ iarray,
343             G4double* __restrict__ oarray);
344 void logfv(const uint32_t size, G4float const* __restrict__ iarray,
345            G4float* __restrict__ oarray);
346 void G4Logfv(const uint32_t size, G4float const* __restrict__ iarray,
347              G4float* __restrict__ oarray);
348 
349 #endif /* WIN32 */
350 
351 #endif /* LOG_H_ */
352