C++InterfacetoTauola
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 *//////////////////////////////////////////////////////////////////////////////