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