| 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 // G4Pow << 26 // $Id: G4Pow.hh,v 1.2 2009/07/03 11:35:07 vnivanch Exp $
>> 27 // GEANT4 tag $Name: geant4-09-03 $
>> 28 //
>> 29 //
>> 30 // -------------------------------------------------------------------
>> 31 //
>> 32 // Class G4Pow
27 // 33 //
28 // Class description: 34 // Class description:
29 // 35 //
30 // Utility singleton class for the fast comput 36 // Utility singleton class for the fast computation of log and pow
31 // functions. Integer argument should be in th << 37 // functions
32 // check is performed inside these methods for << 38
33 // For factorial integer argument should be in << 39 // Author: Vladimir Ivanchenko
34 // Computations with double arguments are fast << 40 //
35 // 0.002-511.5 for all functions except expone << 41 // Creation date: 23.05.2009
36 // for the interval 0-84.4, standard library i << 42 // -------------------------------------------------------------------
37 << 43
38 // Author: Vladimir Ivanchenko, 23.05.2009 << 44 #ifndef G4Pow_h
39 // ------------------------------------------- << 45 #define G4Pow_h 1
40 #ifndef G4Pow_hh <<
41 #define G4Pow_hh 1 <<
42 46
43 #include "G4DataVector.hh" <<
44 #include "G4Exp.hh" <<
45 #include "G4Log.hh" <<
46 #include "globals.hh" 47 #include "globals.hh"
>> 48 #include "G4DataVector.hh"
47 49
48 class G4Pow 50 class G4Pow
49 { 51 {
50 public: <<
51 static G4Pow* GetInstance(); <<
52 ~G4Pow() = default; <<
53 <<
54 // Fast computation of Z^1/3 <<
55 // <<
56 inline G4double Z13(G4int Z) const; <<
57 G4double A13(G4double A) const; <<
58 <<
59 // Fast computation of Z^2/3 <<
60 // <<
61 inline G4double Z23(G4int Z) const; <<
62 inline G4double A23(G4double A) const; <<
63 <<
64 // Fast computation of log(Z) <<
65 // <<
66 inline G4double logZ(G4int Z) const; <<
67 inline G4double logA(G4double A) const; <<
68 inline G4double logX(G4double x) const; <<
69 <<
70 // Fast computation of log10(Z) <<
71 // <<
72 inline G4double log10Z(G4int Z) const; <<
73 inline G4double log10A(G4double A) const; <<
74 <<
75 // Fast computation of exp(X) <<
76 // <<
77 inline G4double expA(G4double A) const; <<
78 <<
79 // Fast computation of pow(Z,X) <<
80 // <<
81 inline G4double powZ(G4int Z, G4double y) co <<
82 inline G4double powA(G4double A, G4double y) <<
83 G4double powN(G4double x, G4int n) const; <<
84 <<
85 // Fast factorial <<
86 // <<
87 inline G4double factorial(G4int Z) const; <<
88 inline G4double logfactorial(G4int Z) const; <<
89 <<
90 private: <<
91 G4Pow(); <<
92 <<
93 G4double A13Low(const G4double, const G4bool <<
94 G4double A13High(const G4double, const G4boo <<
95 <<
96 inline G4double logBase(G4double x) const; <<
97 <<
98 static G4Pow* fpInstance; <<
99 <<
100 const G4double onethird = 1.0 / 3.0; <<
101 const G4int max2 = 5; <<
102 <<
103 G4double maxA; <<
104 G4double maxLowA; <<
105 G4double maxA2; <<
106 G4double maxAexp; <<
107 <<
108 G4DataVector ener; <<
109 G4DataVector logen; <<
110 G4DataVector pz13; <<
111 G4DataVector lowa13; <<
112 G4DataVector lz; <<
113 G4DataVector lz2; <<
114 G4DataVector fexp; <<
115 G4DataVector fact; <<
116 G4DataVector logfact; <<
117 }; <<
118 52
119 // ----------------------------- << 53 public:
120 // Inline methods implementation <<
121 // ----------------------------- <<
122 54
123 inline G4double G4Pow::Z13(G4int Z) const { re << 55 static G4Pow* GetInstance();
124 56
125 inline G4double G4Pow::Z23(G4int Z) const << 57 // Fast computation of Z^1/3
126 { << 58 //
127 G4double x = Z13(Z); << 59 inline G4double Z13(G4int Z);
128 return x * x; << 60 inline G4double A13(G4double A);
129 } << 61
>> 62 // Fast computation of Z^2/3
>> 63 //
>> 64 inline G4double Z23(G4int Z);
>> 65 inline G4double A23(G4double A);
>> 66
>> 67 // Fast computation of log(Z)
>> 68 //
>> 69 inline G4double logZ(G4int Z);
>> 70 inline G4double logA(G4double A);
>> 71
>> 72 // Fast computation of log10(Z)
>> 73 //
>> 74 inline G4double log10Z(G4int Z);
>> 75 inline G4double log10A(G4double A);
>> 76
>> 77 // Fast computation of pow(Z,X)
>> 78 //
>> 79 inline G4double powZ(G4int Z, G4double y);
>> 80 inline G4double powA(G4double A, G4double y);
130 81
131 inline G4double G4Pow::A23(G4double A) const << 82 // Fast factorial
>> 83 //
>> 84 inline G4double factorial(G4int Z);
>> 85
>> 86 private:
>> 87
>> 88 G4Pow();
>> 89 ~G4Pow();
>> 90
>> 91 private:
>> 92
>> 93 static G4Pow* fInstance;
>> 94
>> 95 G4double onethird;
>> 96
>> 97 G4DataVector pz13;
>> 98 G4DataVector lz;
>> 99 G4DataVector fact;
>> 100 };
>> 101
>> 102 // -------------------------------------------------------------------
>> 103
>> 104 inline G4double G4Pow::Z13(G4int Z)
132 { 105 {
133 G4double x = A13(A); << 106 return pz13[Z];
134 return x * x; <<
135 } 107 }
136 108
137 inline G4double G4Pow::logZ(G4int Z) const { r << 109 inline G4double G4Pow::A13(G4double A)
138 <<
139 inline G4double G4Pow::logBase(G4double a) con <<
140 { 110 {
>> 111 const G4int maxZ = 256;
>> 112
141 G4double res; 113 G4double res;
142 if(a <= maxA2) << 114 G4int i = G4int(A + 0.5);
>> 115 if((i >= 1) && (i < maxZ))
143 { 116 {
144 G4int i = G4int(max2 * (a - 1) + 0.5); << 117 G4double x = (1.0 - A/G4double(i))*onethird;
145 if(i > max2) << 118 res = pz13[i]*(1.0 + x - x*x*(1.0 - 1.66666666*x));
146 { <<
147 i = max2; <<
148 } <<
149 G4double x = a / (G4double(i) / max2 + 1) <<
150 res = lz2[i] + x * (1.0 - (0.5 - on <<
151 } <<
152 else if(a <= maxA) <<
153 { <<
154 G4int i = G4int(a + 0.5); <<
155 G4double x = a / G4double(i) - 1; <<
156 res = lz[i] + x * (1.0 - (0.5 - one <<
157 } 119 }
158 else 120 else
159 { 121 {
160 res = G4Log(a); << 122 res = std::pow(A, onethird);
161 } 123 }
162 return res; 124 return res;
163 } 125 }
164 126
165 inline G4double G4Pow::logA(G4double A) const << 127 inline G4double G4Pow::Z23(G4int Z)
>> 128 {
>> 129 G4double x = Z13(Z);
>> 130 return x*x;
>> 131 }
>> 132
>> 133 inline G4double G4Pow::A23(G4double A)
166 { 134 {
167 return (1.0 <= A ? logBase(A) : -logBase(1. << 135 G4double x = A13(A);
>> 136 return x*x;
168 } 137 }
169 138
170 inline G4double G4Pow::logX(G4double x) const << 139 inline G4double G4Pow::logZ(G4int Z)
171 { 140 {
172 G4double res = 0.0; << 141 return lz[Z];
173 G4double a = (1.0 <= x) ? x : 1.0 / x; << 142 }
174 143
175 if(a <= maxA) << 144 inline G4double G4Pow::logA(G4double A)
176 { << 145 {
177 res = logBase(a); << 146 const G4int maxZ = 256;
178 } << 147
179 else if(a <= ener[2]) << 148 G4double res;
180 { << 149 G4int i = G4int(A + 0.5);
181 res = logen[1] + logBase(a / ener[1]); << 150 if((i >= 1) && (i < maxZ))
182 } <<
183 else if(a <= ener[3]) <<
184 { 151 {
185 res = logen[2] + logBase(a / ener[2]); << 152 G4double x = 1.0 - A/G4double(i);
>> 153 res = lz[i] + x*(1.0 - 0.5*x + onethird*x*x);
186 } 154 }
187 else 155 else
188 { 156 {
189 res = G4Log(a); << 157 res = std::log(A);
190 } <<
191 <<
192 if(1.0 > x) <<
193 { <<
194 res = -res; <<
195 } 158 }
196 return res; 159 return res;
197 } 160 }
198 161
199 inline G4double G4Pow::log10Z(G4int Z) const { << 162 inline G4double G4Pow::log10Z(G4int Z)
200 <<
201 inline G4double G4Pow::log10A(G4double A) cons <<
202 <<
203 inline G4double G4Pow::expA(G4double A) const <<
204 { 163 {
205 G4double res; << 164 return lz[Z]/lz[10];
206 G4double a = (0.0 <= A) ? A : -A; << 165 }
207 166
208 if(a <= maxAexp) << 167 inline G4double G4Pow::log10A(G4double A)
209 { << 168 {
210 G4int i = G4int(2 * a + 0.5); << 169 return logA(A)/lz[10];
211 G4double x = a - i * 0.5; <<
212 res = fexp[i] * (1.0 + x * (1.0 + 0 <<
213 } <<
214 else <<
215 { <<
216 res = G4Exp(a); <<
217 } <<
218 if(0.0 > A) <<
219 { <<
220 res = 1.0 / res; <<
221 } <<
222 return res; <<
223 } 170 }
224 171
225 inline G4double G4Pow::powZ(G4int Z, G4double << 172 inline G4double G4Pow::powZ(G4int Z, G4double y)
226 { 173 {
227 return expA(y * lz[Z]); << 174 return std::exp(y*lz[Z]);
228 } 175 }
229 176
230 inline G4double G4Pow::powA(G4double A, G4doub << 177 inline G4double G4Pow::powA(G4double A, G4double y)
231 { 178 {
232 return (0.0 == A ? 0.0 : expA(y * logX(A))); << 179 return std::exp(y*logA(A));
233 } 180 }
234 181
235 inline G4double G4Pow::factorial(G4int Z) cons << 182 inline G4double G4Pow::factorial(G4int Z)
>> 183 {
>> 184 return fact[Z];
>> 185 }
236 186
237 inline G4double G4Pow::logfactorial(G4int Z) c << 187 // -------------------------------------------------------------------
238 188
239 #endif 189 #endif
>> 190
240 191