Geant4  10.01.p03
G4Exp.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 // The basic idea is to exploit Pade polynomials.
35 // A lot of ideas were inspired by the cephes math library
36 // (by Stephen L. Moshier moshier@na-net.ornl.gov) as well as actual code.
37 // The Cephes library can be found here: http://www.netlib.org/cephes/
38 // Code and algorithms for G4Exp have been extracted and adapted for Geant4
39 // from the original implementation in the VDT mathematical library
40 // (https://svnweb.cern.ch/trac/vdt), version 0.3.7.
41 
42 // Original implementation created on: Jun 23, 2012
43 // Author: Danilo Piparo, Thomas Hauth, Vincenzo Innocente
44 //
45 // --------------------------------------------------------------------
46 /*
47  * VDT is free software: you can redistribute it and/or modify
48  * it under the terms of the GNU Lesser Public License as published by
49  * the Free Software Foundation, either version 3 of the License, or
50  * (at your option) any later version.
51  *
52  * This program is distributed in the hope that it will be useful,
53  * but WITHOUT ANY WARRANTY; without even the implied warranty of
54  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
55  * GNU Lesser Public License for more details.
56  *
57  * You should have received a copy of the GNU Lesser Public License
58  * along with this program. If not, see <http://www.gnu.org/licenses/>.
59  */
60 // --------------------------------------------------------------------
61 #ifndef G4Exp_h
62 #define G4Exp_h 1
63 
64 #ifdef WIN32
65 
66  #define G4Exp std::exp
67 
68 #else
69 
70 #include <limits>
71 #include <stdint.h>
72 #include "G4Types.hh"
73 
74 namespace G4ExpConsts
75 {
76  const G4double EXP_LIMIT = 708;
77 
78  const G4double PX1exp = 1.26177193074810590878E-4;
79  const G4double PX2exp = 3.02994407707441961300E-2;
80  const G4double PX3exp = 9.99999999999999999910E-1;
81  const G4double QX1exp = 3.00198505138664455042E-6;
82  const G4double QX2exp = 2.52448340349684104192E-3;
83  const G4double QX3exp = 2.27265548208155028766E-1;
84  const G4double QX4exp = 2.00000000000000000009E0;
85 
86  const G4double LOG2E = 1.4426950408889634073599; // 1/log(2)
87 
88  const G4float MAXLOGF = 88.72283905206835f;
89  const G4float MINLOGF = -88.f;
90 
91  const G4float C1F = 0.693359375f;
92  const G4float C2F = -2.12194440e-4f;
93 
94  const G4float PX1expf = 1.9875691500E-4f;
95  const G4float PX2expf =1.3981999507E-3f;
96  const G4float PX3expf =8.3334519073E-3f;
97  const G4float PX4expf =4.1665795894E-2f;
98  const G4float PX5expf =1.6666665459E-1f;
99  const G4float PX6expf =5.0000001201E-1f;
100 
101  const G4float LOG2EF = 1.44269504088896341f;
102 
103  //----------------------------------------------------------------------------
104  // Used to switch between different type of interpretations of the data
105  // (64 bits)
106  //
107  union ieee754
108  {
109  ieee754 () {};
110  ieee754 (G4double thed) {d=thed;};
111  ieee754 (uint64_t thell) {ll=thell;};
112  ieee754 (G4float thef) {f[0]=thef;};
113  ieee754 (uint32_t thei) {i[0]=thei;};
114  G4double d;
115  G4float f[2];
116  uint32_t i[2];
117  uint64_t ll;
118  uint16_t s[4];
119  };
120 
121  //----------------------------------------------------------------------------
122  // Converts an unsigned long long to a double
123  //
124  inline G4double uint642dp(uint64_t ll)
125  {
126  ieee754 tmp;
127  tmp.ll=ll;
128  return tmp.d;
129  }
130 
131  //----------------------------------------------------------------------------
132  // Converts an int to a float
133  //
135  {
136  ieee754 tmp;
137  tmp.i[0]=x;
138  return tmp.f[0];
139  }
140 
141  //----------------------------------------------------------------------------
142  // Converts a float to an int
143  //
144  inline uint32_t sp2uint32(G4float x)
145  {
146  ieee754 tmp;
147  tmp.f[0]=x;
148  return tmp.i[0];
149  }
150 
151  //----------------------------------------------------------------------------
157  inline G4double fpfloor(const G4double x)
158  {
159  // no problem since exp is defined between -708 and 708. Int is enough for it!
160  int32_t ret = int32_t (x);
161  ret-=(sp2uint32(x)>>31);
162  return ret;
163  }
164 
165  //----------------------------------------------------------------------------
171  inline G4float fpfloor(const G4float x)
172  {
173  int32_t ret = int32_t (x);
174  ret-=(sp2uint32(x)>>31);
175  return ret;
176  }
177 }
178 
179 // Exp double precision --------------------------------------------------------
180 
181 
183 inline G4double G4Exp(G4double initial_x)
184 {
185  G4double x = initial_x;
187 
188  const int32_t n = int32_t(px);
189 
190  x -= px * 6.93145751953125E-1;
191  x -= px * 1.42860682030941723212E-6;
192 
193  const G4double xx = x * x;
194 
195  // px = x * P(x**2).
196  px = G4ExpConsts::PX1exp;
197  px *= xx;
198  px += G4ExpConsts::PX2exp;
199  px *= xx;
200  px += G4ExpConsts::PX3exp;
201  px *= x;
202 
203  // Evaluate Q(x**2).
205  qx *= xx;
206  qx += G4ExpConsts::QX2exp;
207  qx *= xx;
208  qx += G4ExpConsts::QX3exp;
209  qx *= xx;
210  qx += G4ExpConsts::QX4exp;
211 
212  // e**x = 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
213  x = px / (qx - px);
214  x = 1.0 + 2.0 * x;
215 
216  // Build 2^n in double.
217  x *= G4ExpConsts::uint642dp(( ((uint64_t)n) +1023)<<52);
218 
219  if (initial_x > G4ExpConsts::EXP_LIMIT)
220  x = std::numeric_limits<G4double>::infinity();
221  if (initial_x < -G4ExpConsts::EXP_LIMIT)
222  x = 0.;
223 
224  return x;
225 }
226 
227 // Exp single precision --------------------------------------------------------
228 
230 inline G4float G4Expf(G4float initial_x)
231 {
232  G4float x = initial_x;
233 
234  G4float z = G4ExpConsts::fpfloor( G4ExpConsts::LOG2EF * x +0.5f ); /* floor() truncates toward -infinity. */
235 
236  x -= z * G4ExpConsts::C1F;
237  x -= z * G4ExpConsts::C2F;
238  const int32_t n = int32_t ( z );
239 
240  const G4float x2 = x * x;
241 
242  z = x*G4ExpConsts::PX1expf;
244  z *= x;
246  z *= x;
248  z *= x;
250  z *= x;
252  z *= x2;
253  z += x + 1.0f;
254 
255  /* multiply by power of 2 */
256  z *= G4ExpConsts::uint322sp((n+0x7f)<<23);
257 
258  if (initial_x > G4ExpConsts::MAXLOGF) z=std::numeric_limits<G4float>::infinity();
259  if (initial_x < G4ExpConsts::MINLOGF) z=0.f;
260 
261  return z;
262 }
263 
264 //------------------------------------------------------------------------------
265 
266 void expv(const uint32_t size, G4double const * __restrict__ iarray, G4double* __restrict__ oarray);
267 void G4Expv(const uint32_t size, G4double const * __restrict__ iarray, G4double* __restrict__ oarray);
268 void expfv(const uint32_t size, G4float const * __restrict__ iarray, G4float* __restrict__ oarray);
269 void G4Expfv(const uint32_t size, G4float const * __restrict__ iarray, G4float* __restrict__ oarray);
270 
271 #endif /* WIN32 */
272 
273 #endif
const G4double QX1exp
Definition: G4Exp.hh:81
ieee754(G4float thef)
Definition: G4Exp.hh:112
G4double fpfloor(const G4double x)
A vectorisable floor implementation, not only triggered by fast-math.
Definition: G4Exp.hh:157
uint32_t sp2uint32(G4float x)
Definition: G4Exp.hh:144
void expfv(const uint32_t size, G4float const *__restrict__ iarray, G4float *__restrict__ oarray)
const G4double EXP_LIMIT
Definition: G4Exp.hh:76
G4double z
Definition: TRTMaterials.hh:39
void G4Expfv(const uint32_t size, G4float const *__restrict__ iarray, G4float *__restrict__ oarray)
float G4float
Definition: G4Types.hh:77
const G4float PX1expf
Definition: G4Exp.hh:94
const G4double LOG2E
Definition: G4Exp.hh:86
const G4double QX2exp
Definition: G4Exp.hh:82
const G4double PX3exp
Definition: G4Exp.hh:80
const G4float C1F
Definition: G4Exp.hh:91
const G4float PX4expf
Definition: G4Exp.hh:97
const G4double PX1exp
Definition: G4Exp.hh:78
G4float f[2]
Definition: G4Exp.hh:115
int G4int
Definition: G4Types.hh:78
const G4float PX2expf
Definition: G4Exp.hh:95
const G4double QX4exp
Definition: G4Exp.hh:84
const G4float MAXLOGF
Definition: G4Exp.hh:88
G4float G4Expf(G4float initial_x)
Exponential Function single precision.
Definition: G4Exp.hh:230
const G4float PX3expf
Definition: G4Exp.hh:96
G4float uint322sp(G4int x)
Definition: G4Exp.hh:134
const G4float C2F
Definition: G4Exp.hh:92
const G4int n
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
const G4float LOG2EF
Definition: G4Exp.hh:101
const G4float PX6expf
Definition: G4Exp.hh:99
ieee754(uint64_t thell)
Definition: G4Exp.hh:111
const G4double QX3exp
Definition: G4Exp.hh:83
uint16_t s[4]
Definition: G4Exp.hh:118
ieee754(uint32_t thei)
Definition: G4Exp.hh:113
void G4Expv(const uint32_t size, G4double const *__restrict__ iarray, G4double *__restrict__ oarray)
const G4float MINLOGF
Definition: G4Exp.hh:89
ieee754(G4double thed)
Definition: G4Exp.hh:110
G4double uint642dp(uint64_t ll)
Definition: G4Exp.hh:124
double G4double
Definition: G4Types.hh:76
uint32_t i[2]
Definition: G4Exp.hh:116
void expv(const uint32_t size, G4double const *__restrict__ iarray, G4double *__restrict__ oarray)
const G4double PX2exp
Definition: G4Exp.hh:79
const G4float PX5expf
Definition: G4Exp.hh:98