C++ Interface to Tauola
MathLib.f
1*//////////////////////////////////////////////////////////////////////////////////////////////////
2*// //
3*// //
4*// Pseudo-CLASS Mathlib //
5*// //
6*// Purpose: library of math utilies //
7*// //
8*// SUBROUTINE Mathlib_GausJad(fun,aa,bb,eeps,result) : Gauss integration //
9*// DOUBLE PRECISION FUNCTION Mathlib_Gauss(f,a,b,eeps) : Gauss integration //
10*// DOUBLE PRECISION FUNCTION Mathlib_dilogy(x) : Dilog function Li_2 //
11*// DOUBLE PRECISION FUNCTION Mathlib_dpgamm(z) : Euler Gamma function //
12*// //
13*//////////////////////////////////////////////////////////////////////////////////////////////////
14
15
16 SUBROUTINE mathlib_gausjad(fun,aa,bb,eeps,result)
17*//////////////////////////////////////////////////////////////////////////////
18*// //
19*// Gauss-type integration by S. Jadach, Oct. 1990, June 1997 //
20*// this is non-adaptive (!!!!) unoptimized (!!!) integration subprogram. //
21*// //
22*// Eeps>0 treated as ABSOLUTE requested error //
23*// Eeps<0 treated as RELATIVE requested error //
24*// //
25*//////////////////////////////////////////////////////////////////////////////
26 IMPLICIT NONE
27*
28 DOUBLE PRECISION fun,aa,bb,eeps,result
29 EXTERNAL fun
30*
31 DOUBLE PRECISION a,b,xplus,sum16,range,sum8,erabs,erela,fminu,xminu
32 DOUBLE PRECISION fplus,xmidle,calk8,eps,x1,x2,delta,calk16
33 INTEGER iter,ndivi,itermx,k,i
34 DOUBLE PRECISION wg(12),xx(12)
35 DATA wg
36 $/0.101228536290376d0, 0.222381034453374d0, 0.313706645877887d0,
37 $ 0.362683783378362d0, 0.027152459411754d0, 0.062253523938648d0,
38 $ 0.095158511682493d0, 0.124628971255534d0, 0.149595988816577d0,
39 $ 0.169156519395003d0, 0.182603415044924d0, 0.189450610455069d0/
40 DATA xx
41 $/0.960289856497536d0, 0.796666477413627d0, 0.525532409916329d0,
42 $ 0.183434642495650d0, 0.989400934991650d0, 0.944575023073233d0,
43 $ 0.865631202387832d0, 0.755404408355003d0, 0.617876244402644d0,
44 $ 0.458016777657227d0, 0.281603550779259d0, 0.095012509837637d0/
45 DATA itermx / 15/
46*-------------------------------------------------------------------------------
47 a = aa
48 b = bb
49 eps= abs(eeps)
50 ndivi=1
51*** iteration over subdivisions terminated by precision requirement
52 DO iter=1,itermx
53 calk8 =0d0
54 calk16 =0d0
55*** sum over delta subintegrals
56 DO k = 1,ndivi
57 delta = (b-a)/ndivi
58 x1 = a + (k-1)*delta
59 x2 = x1+ delta
60 xmidle= 0.5d0*(x2+x1)
61 range = 0.5d0*(x2-x1)
62 sum8 =0d0
63 sum16=0d0
64*** 8- and 12-point gauss integration over single delta subinterval
65 DO i=1,12
66 xplus= xmidle+range*xx(i)
67 xminu= xmidle-range*xx(i)
68 fplus=fun(xplus)
69 fminu=fun(xminu)
70 IF(i .LE. 4) THEN
71 sum8 =sum8 +(fplus+fminu)*wg(i)/2d0
72 ELSE
73 sum16=sum16 +(fplus+fminu)*wg(i)/2d0
74 ENDIF
75 ENDDO
76 calk8 = calk8 + sum8 *(x2-x1)
77 calk16= calk16+ sum16*(x2-x1)
78 ENDDO
79 erabs = abs(calk16-calk8)
80 erela = 0d0
81 IF(calk16 .NE. 0d0) erela= erabs/abs(calk16)
82*** WRITE(*,*) 'gausjad: calk8,calk16=',iter,calk8,calk16,erela
83*** precision check to terminate integration
84 IF(eeps .GT. 0d0) THEN
85 IF(erabs .LT. eps) GOTO 800
86 ELSE
87 IF(erela .LT. eps) GOTO 800
88 ENDIF
89 ndivi=ndivi*2
90 ENDDO
91 WRITE(*,*) '++++ Mathlib_GausJad: required precision to high!'
92 WRITE(*,*) '++++ Mathlib_GausJad: eeps,erela,erabs=',eeps,erela,erabs
93 800 CONTINUE
94 result = calk16
95 END
96
97
98 DOUBLE PRECISION FUNCTION mathlib_gauss(f,a,b,eeps)
99*//////////////////////////////////////////////////////////////////////////////
100*// //
101*// this is iterative integration procedure //
102*// originates probably from cern library //
103*// it subdivides inegration range until required PRECISION is reached //
104*// PRECISION is a difference from 8 and 16 point gauss itegr. result //
105*// eeps positive treated as absolute PRECISION //
106*// eeps negative treated as relative PRECISION //
107*// //
108*//////////////////////////////////////////////////////////////////////////////
109 IMPLICIT NONE
110 DOUBLE PRECISION f,a,b,eeps
111*
112 DOUBLE PRECISION c1,c2,bb,s8,s16,y,aa,const,delta,eps,u
113 INTEGER i
114*
115 DOUBLE PRECISION w(12),x(12)
116 EXTERNAL f
117 DATA const /1.0d-19/
118 DATA w
119 1/0.10122 85362 90376, 0.22238 10344 53374, 0.31370 66458 77887,
120 2 0.36268 37833 78362, 0.02715 24594 11754, 0.06225 35239 38648,
121 3 0.09515 85116 82493, 0.12462 89712 55534, 0.14959 59888 16577,
122 4 0.16915 65193 95003, 0.18260 34150 44924, 0.18945 06104 55069/
123 DATA x
124 1/0.96028 98564 97536, 0.79666 64774 13627, 0.52553 24099 16329,
125 2 0.18343 46424 95650, 0.98940 09349 91650, 0.94457 50230 73233,
126 3 0.86563 12023 87832, 0.75540 44083 55003, 0.61787 62444 02644,
127 4 0.45801 67776 57227, 0.28160 35507 79259, 0.09501 25098 37637/
128*-----------------------------------------------------------------------------
129 eps=abs(eeps)
130 delta=const*abs(a-b)
131 mathlib_gauss=0d0
132 aa=a
133 5 y=b-aa
134 IF(abs(y) .LE. delta) RETURN
135 2 bb=aa+y
136 c1=0.5d0*(aa+bb)
137 c2=c1-aa
138 s8=0d0
139 s16=0d0
140 DO 1 i=1,4
141 u=x(i)*c2
142 1 s8=s8+w(i)*(f(c1+u)+f(c1-u))
143 DO 3 i=5,12
144 u=x(i)*c2
145 3 s16=s16+w(i)*(f(c1+u)+f(c1-u))
146 s8=s8*c2
147 s16=s16*c2
148 IF(eeps .LT. 0d0) THEN
149 IF(abs(s16-s8) .GT. eps*abs(s16)) GOTO 4
150 ELSE
151 IF(abs(s16-s8) .GT. eps) GOTO 4
152 ENDIF
153 mathlib_gauss=mathlib_gauss+s16
154 aa=bb
155 GOTO 5
156 4 y=0.5d0*y
157 IF(abs(y) .GT. delta) GOTO 2
158 WRITE(*,7)
159 mathlib_gauss=0d0
160 RETURN
161 7 FORMAT(1x,36hgaus ... too high accuracy required)
162 END
163
164
165 DOUBLE PRECISION FUNCTION mathlib_dilogy(x)
166*//////////////////////////////////////////////////////////////////////////////
167*// //
168*// dilogarithm FUNCTION: dilog(x)=int( -ln(1-z)/z ) , 0 < z < x . //
169*// this is the cernlib version. //
170*// //
171*//////////////////////////////////////////////////////////////////////////////
172 IMPLICIT NONE
173 DOUBLE PRECISION x
174* locals
175 DOUBLE PRECISION a,b,s,t,y,z
176*------------------------------------------------------------------------------
177 z=-1.644934066848226d0
178 IF(x .LT. -1.d0) go to 1
179 IF(x .LE. 0.5d0) go to 2
180 IF(x .EQ. 1.d0) go to 3
181 IF(x .LE. 2.d0) go to 4
182 z=3.289868133696453d0
183 1 t=1.d0/x
184 s=-0.5d0
185 z=z-0.5d0*dlog(dabs(x))**2
186 go to 5
187 2 t=x
188 s=0.5d0
189 z=0.d0
190 go to 5
191 3 mathlib_dilogy=1.644934066848226d0
192 RETURN
193 4 t=1.d0-x
194 s=-0.5d0
195 z=1.644934066848226d0-dlog(x)*dlog(dabs(t))
196 5 y=2.666666666666667d0*t+0.666666666666667d0
197 b= 0.000000000000001d0
198 a=y*b +0.000000000000004d0
199 b=y*a-b+0.000000000000011d0
200 a=y*b-a+0.000000000000037d0
201 b=y*a-b+0.000000000000121d0
202 a=y*b-a+0.000000000000398d0
203 b=y*a-b+0.000000000001312d0
204 a=y*b-a+0.000000000004342d0
205 b=y*a-b+0.000000000014437d0
206 a=y*b-a+0.000000000048274d0
207 b=y*a-b+0.000000000162421d0
208 a=y*b-a+0.000000000550291d0
209 b=y*a-b+0.000000001879117d0
210 a=y*b-a+0.000000006474338d0
211 b=y*a-b+0.000000022536705d0
212 a=y*b-a+0.000000079387055d0
213 b=y*a-b+0.000000283575385d0
214 a=y*b-a+0.000001029904264d0
215 b=y*a-b+0.000003816329463d0
216 a=y*b-a+0.000014496300557d0
217 b=y*a-b+0.000056817822718d0
218 a=y*b-a+0.000232002196094d0
219 b=y*a-b+0.001001627496164d0
220 a=y*b-a+0.004686361959447d0
221 b=y*a-b+0.024879322924228d0
222 a=y*b-a+0.166073032927855d0
223 a=y*a-b+1.935064300869969d0
224 mathlib_dilogy = s*t*(a-b)+z
225 END
226
227
228 DOUBLE PRECISION FUNCTION mathlib_dpgamm(z)
229*//////////////////////////////////////////////////////////////////////////////
230*// //
231*// Double precision gamma function //
232*// //
233*//////////////////////////////////////////////////////////////////////////////
234 DOUBLE PRECISION z,z1,x,x1,x2,d1,d2,s1,s2,s3,pi,c(20),const
235 SAVE c,pi,const
236 DATA c( 1) / 8.3333333333333333333333333332d-02/
237 DATA c( 2) /-2.7777777777777777777777777777d-03/
238 DATA c( 3) / 7.9365079365079365079365079364d-04/
239 DATA c( 4) /-5.9523809523809523809523809523d-04/
240 DATA c( 5) / 8.4175084175084175084175084175d-04/
241 DATA c( 6) /-1.9175269175269175269175269175d-03/
242 DATA c( 7) / 6.4102564102564102564102564102d-03/
243 DATA c( 8) /-2.9550653594771241830065359477d-02/
244 DATA c( 9) / 1.7964437236883057316493849001d-01/
245 DATA c(10) /-1.3924322169059011164274322169d+00/
246 DATA c(11) / 1.3402864044168391994478951001d+01/
247 DATA c(12) /-1.5684828462600201730636513245d+02/
248 DATA c(13) / 2.1931033333333333333333333333d+03/
249 DATA c(14) /-3.6108771253724989357173265219d+04/
250 DATA c(15) / 6.9147226885131306710839525077d+05/
251 DATA c(16) /-1.5238221539407416192283364959d+07/
252 DATA c(17) / 3.8290075139141414141414141414d+08/
253 DATA c(18) /-1.0882266035784391089015149165d+10/
254 DATA c(19) / 3.4732028376500225225225225224d+11/
255 DATA c(20) /-1.2369602142269274454251710349d+13/
256 DATA pi / 3.1415926535897932384626433832d+00/
257 DATA const / 9.1893853320467274178032973641d-01/
258 IF(z .GT. 5.75d 1) GOTO 6666
259 nn = z
260 IF (z - dble(float(nn))) 3,1,3
261 1 IF (z .LE. 0.d 0) GOTO 6667
262 mathlib_dpgamm = 1.d 0
263 IF (z .LE. 2.d 0) RETURN
264 z1 = z
265 2 z1 = z1 - 1.d 0
266 mathlib_dpgamm = mathlib_dpgamm * z1
267 IF (z1 - 2.d 0) 61,61,2
268 3 IF (dabs(z) .LT. 1.d-29) GOTO 60
269 IF (z .LT. 0.d 0) GOTO 4
270 x = z
271 kk = 1
272 GOTO 10
273 4 x = 1.d 0 - z
274 kk = 2
275 10 x1 = x
276 IF (x .GT. 19.d 0) GOTO 13
277 d1 = x
278 11 x1 = x1 + 1.d 0
279 IF (x1 .GE. 19.d 0) GOTO 12
280 d1 = d1 * x1
281 GOTO 11
282 12 s3 = -dlog(d1)
283 GOTO 14
284 13 s3 = 0.d 0
285 14 d1 = x1 * x1
286 s1 = (x1 - 5.d-1) * dlog(x1) - x1 + const
287 DO 20 k=1,20
288 s2 = s1 + c(k)/x1
289 IF (dabs(s2 - s1) .LT. 1.d-28) GOTO 21
290 x1 = x1 * d1
291 20 s1 = s2
292 21 s3 = s3 + s2
293 GOTO (50,22), kk
294 22 d2 = dabs(z - nn)
295 d1 = d2 * pi
296 IF (d1 .LT. 1.d-15) GOTO 31
297 30 x2 = dlog(pi/dsin(d1)) - s3
298 GOTO 40
299 31 x2 = -dlog(d2)
300 40 mm = dabs(z)
301 IF(x2 .GT. 1.74d2) GOTO 6666
302 mathlib_dpgamm = dexp(x2)
303 IF (mm .ne. (mm/2) * 2) RETURN
304 mathlib_dpgamm = -mathlib_dpgamm
305 RETURN
306 50 IF(s3 .GT. 1.74d2) GOTO 6666
307 mathlib_dpgamm = dexp(s3)
308 RETURN
309 6666 print *, 2000
310 RETURN
311 6667 print *, 2001
312 RETURN
313 60 mathlib_dpgamm = 0.d 0
314 IF(dabs(z) .LT. 1.d-77) RETURN
315 mathlib_dpgamm = 1.d 0/z
316 61 RETURN
317 2000 FORMAT (/////, 2x, 32hdpgamm ..... argument too large., /////)
318 2001 FORMAT (/////, 2x, 32hdpgamm ..... argument is a pole., /////)
319 END
320
321
322
323 SUBROUTINE mathlib_gaus8(fun,aa,bb,result)
324*//////////////////////////////////////////////////////////////////////////////
325*// 8-point Gauss //
326*//////////////////////////////////////////////////////////////////////////////
327 IMPLICIT NONE
328 DOUBLE PRECISION fun,aa,bb,result
329 EXTERNAL fun
330 DOUBLE PRECISION a,b,sum8,xmidle,range,xplus,xminu
331 INTEGER k,i
332 DOUBLE PRECISION wg(4),xx(4)
333 DATA wg /0.101228536290376d0, 0.222381034453374d0, 0.313706645877887d0, 0.362683783378362d0/
334 DATA xx /0.960289856497536d0, 0.796666477413627d0, 0.525532409916329d0, 0.183434642495650d0/
335*-------------------------------------------------------------------------------
336 a = aa
337 b = bb
338 xmidle= 0.5d0*(a+b)
339 range = 0.5d0*(b-a)
340 sum8 =0d0
341 DO i=1,4
342 xplus= xmidle+range*xx(i)
343 xminu= xmidle-range*xx(i)
344 sum8 =sum8 +(fun(xplus)+fun(xminu))*wg(i)/2d0
345 ENDDO
346 result = sum8*(b-a)
347 END
348
349 SUBROUTINE mathlib_gaus16(fun,aa,bb,result)
350*//////////////////////////////////////////////////////////////////////////////
351*// 12-point Gauss //
352*//////////////////////////////////////////////////////////////////////////////
353 IMPLICIT NONE
354 DOUBLE PRECISION fun,aa,bb,result
355 EXTERNAL fun
356 DOUBLE PRECISION a,b,sum16,xmidle,range,xplus,xminu
357 INTEGER k,i
358 DOUBLE PRECISION wg(8),xx(8)
359 DATA wg /0.027152459411754d0, 0.062253523938648d0,
360 $ 0.095158511682493d0, 0.124628971255534d0, 0.149595988816577d0,
361 $ 0.169156519395003d0, 0.182603415044924d0, 0.189450610455069d0/
362 DATA xx /0.989400934991650d0, 0.944575023073233d0,
363 $ 0.865631202387832d0, 0.755404408355003d0, 0.617876244402644d0,
364 $ 0.458016777657227d0, 0.281603550779259d0, 0.095012509837637d0/
365*-------------------------------------------------------------------------------
366 a = aa
367 b = bb
368 xmidle= 0.5d0*(a+b)
369 range = 0.5d0*(b-a)
370 sum16 =0d0
371 DO i=1,8
372 xplus= xmidle+range*xx(i)
373 xminu= xmidle-range*xx(i)
374 sum16 =sum16 +(fun(xplus)+fun(xminu))*wg(i)/2d0
375 ENDDO
376 result = sum16*(b-a)
377 END
378
379
380
381
382*//////////////////////////////////////////////////////////////////////////////
383*// //
384*// End Pseudo-CLASS Mathlib //
385*//////////////////////////////////////////////////////////////////////////////