16 SUBROUTINE mathlib_gausjad(fun,aa,bb,eeps,result)
28 DOUBLE PRECISION fun,aa,bb,eeps,result
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)
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/
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/
66 xplus= xmidle+range*xx(i)
67 xminu= xmidle-range*xx(i)
71 sum8 =sum8 +(fplus+fminu)*wg(i)/2d0
73 sum16=sum16 +(fplus+fminu)*wg(i)/2d0
76 calk8 = calk8 + sum8 *(x2-x1)
77 calk16= calk16+ sum16*(x2-x1)
79 erabs = abs(calk16-calk8)
81 IF(calk16 .NE. 0d0) erela= erabs/abs(calk16)
84 IF(eeps .GT. 0d0)
THEN
85 IF(erabs .LT. eps) goto 800
87 IF(erela .LT. eps) goto 800
91 WRITE(*,*)
'++++ Mathlib_GausJad: required precision to high!'
92 WRITE(*,*)
'++++ Mathlib_GausJad: eeps,erela,erabs=',eeps,erela,erabs
98 DOUBLE PRECISION FUNCTION mathlib_gauss(f,a,b,eeps)
110 DOUBLE PRECISION f,a,b,eeps
112 DOUBLE PRECISION c1,c2,bb,s8,s16,y,aa,const,delta,eps,u
115 DOUBLE PRECISION w(12),x(12)
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/
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/
134 IF(abs(y) .LE. delta)
RETURN
142 1 s8=s8+w(i)*(f(c1+u)+f(c1-u))
145 3 s16=s16+w(i)*(f(c1+u)+f(c1-u))
148 IF(eeps .LT. 0d0)
THEN
149 IF(abs(s16-s8) .GT. eps*abs(s16)) goto 4
151 IF(abs(s16-s8) .GT. eps) goto 4
153 mathlib_gauss=mathlib_gauss+s16
157 IF(abs(y) .GT. delta) goto 2
161 7
FORMAT(1x,36hgaus ... too high accuracy required)
165 DOUBLE PRECISION FUNCTION mathlib_dilogy(x)
175 DOUBLE PRECISION a,b,s,t,y,z
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
185 z=z-0.5d0*dlog(dabs(x))**2
191 3 mathlib_dilogy=1.644934066848226d0
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
228 DOUBLE PRECISION FUNCTION mathlib_dpgamm(z)
234 DOUBLE PRECISION z,z1,x,x1,x2,d1,d2,s1,s2,s3,pi,c(20),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
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
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
276 IF (x .GT. 19.d 0) goto 13
279 IF (x1 .GE. 19.d 0) goto 12
286 s1 = (x1 - 5.d-1) * dlog(x1) - x1 + const
289 IF (dabs(s2 - s1) .LT. 1.d-28) goto 21
296 IF (d1 .LT. 1.d-15) goto 31
297 30 x2 = dlog(pi/dsin(d1)) - s3
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
306 50
IF(s3 .GT. 1.74d2) goto 6666
307 mathlib_dpgamm = dexp(s3)
313 60 mathlib_dpgamm = 0.d 0
314 IF(dabs(z) .LT. 1.d-77)
RETURN
315 mathlib_dpgamm = 1.d 0/z
317 2000
FORMAT (/////, 2x, 32hdpgamm ..... argument too large., /////)
318 2001
FORMAT (/////, 2x, 32hdpgamm ..... argument is a pole., /////)
323 SUBROUTINE mathlib_gaus8(fun,aa,bb,result)
328 DOUBLE PRECISION fun,aa,bb,result
330 DOUBLE PRECISION a,b,sum8,xmidle,range,xplus,xminu
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/
342 xplus= xmidle+range*xx(i)
343 xminu= xmidle-range*xx(i)
344 sum8 =sum8 +(fun(xplus)+fun(xminu))*wg(i)/2d0
349 SUBROUTINE mathlib_gaus16(fun,aa,bb,result)
354 DOUBLE PRECISION fun,aa,bb,result
356 DOUBLE PRECISION a,b,sum16,xmidle,range,xplus,xminu
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/
372 xplus= xmidle+range*xx(i)
373 xminu= xmidle-range*xx(i)
374 sum16 =sum16 +(fun(xplus)+fun(xminu))*wg(i)/2d0