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