Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : MODULE beta_gamma_psi
10 : ! not tested in the case where dp would stand for single precision
11 : ! Routines for the beta function are untested
12 :
13 : USE kinds, ONLY: dp
14 : #include "./base/base_uses.f90"
15 :
16 : IMPLICIT NONE
17 :
18 : PRIVATE
19 :
20 : PUBLIC :: psi
21 :
22 : CONTAINS
23 :
24 : ! **************************************************************************************************
25 : !> \brief ...
26 : !> \param i ...
27 : !> \return ...
28 : ! **************************************************************************************************
29 7564 : FUNCTION ipmpar(i) RESULT(fn_val)
30 : !-----------------------------------------------------------------------
31 :
32 : ! IPMPAR PROVIDES THE INTEGER MACHINE CONSTANTS FOR THE COMPUTER
33 : ! THAT IS USED. IT IS ASSUMED THAT THE ARGUMENT I IS AN INTEGER
34 : ! HAVING ONE OF THE VALUES 1-10. IPMPAR(I) HAS THE VALUE ...
35 :
36 : ! INTEGERS.
37 :
38 : ! ASSUME INTEGERS ARE REPRESENTED IN THE N-DIGIT, BASE-A FORM
39 :
40 : ! SIGN ( X(N-1)*A**(N-1) + ... + X(1)*A + X(0) )
41 :
42 : ! WHERE 0 .LE. X(I) .LT. A FOR I=0,...,N-1.
43 :
44 : ! IPMPAR(1) = A, THE BASE (radix).
45 :
46 : ! IPMPAR(2) = N, THE NUMBER OF BASE-A DIGITS (digits).
47 :
48 : ! IPMPAR(3) = A**N - 1, THE LARGEST MAGNITUDE (huge).
49 :
50 : ! FLOATING-POINT NUMBERS.
51 :
52 : ! IT IS ASSUMED THAT THE SINGLE AND DOUBLE PRECISION FLOATING
53 : ! POINT ARITHMETICS HAVE THE SAME BASE, SAY B, AND THAT THE
54 : ! NONZERO NUMBERS ARE REPRESENTED IN THE FORM
55 :
56 : ! SIGN (B**E) * (X(1)/B + ... + X(M)/B**M)
57 :
58 : ! WHERE X(I) = 0,1,...,B-1 FOR I=1,...,M,
59 : ! X(1) .GE. 1, AND EMIN .LE. E .LE. EMAX.
60 :
61 : ! IPMPAR(4) = B, THE BASE.
62 :
63 : ! SINGLE-PRECISION
64 :
65 : ! IPMPAR(5) = M, THE NUMBER OF BASE-B DIGITS.
66 :
67 : ! IPMPAR(6) = EMIN, THE SMALLEST EXPONENT E.
68 :
69 : ! IPMPAR(7) = EMAX, THE LARGEST EXPONENT E.
70 :
71 : ! DOUBLE-PRECISION
72 :
73 : ! IPMPAR(8) = M, THE NUMBER OF BASE-B DIGITS.
74 :
75 : ! IPMPAR(9) = EMIN, THE SMALLEST EXPONENT E.
76 :
77 : ! IPMPAR(10) = EMAX, THE LARGEST EXPONENT E.
78 :
79 : !-----------------------------------------------------------------------
80 :
81 : INTEGER, INTENT(IN) :: i
82 : INTEGER :: fn_val
83 :
84 7564 : SELECT CASE (i)
85 : CASE (1)
86 0 : fn_val = RADIX(i)
87 : CASE (2)
88 0 : fn_val = DIGITS(i)
89 : CASE (3)
90 7564 : fn_val = HUGE(i)
91 : CASE (4)
92 0 : fn_val = RADIX(1.0)
93 : CASE (5)
94 0 : fn_val = DIGITS(1.0)
95 : CASE (6)
96 0 : fn_val = MINEXPONENT(1.0)
97 : CASE (7)
98 0 : fn_val = MAXEXPONENT(1.0)
99 : CASE (8)
100 0 : fn_val = DIGITS(1.0e0_dp)
101 : CASE (9)
102 0 : fn_val = MINEXPONENT(1.0e0_dp)
103 : CASE (10)
104 0 : fn_val = MAXEXPONENT(1.0e0_dp)
105 : CASE DEFAULT
106 7564 : CPABORT("unknown case")
107 : END SELECT
108 :
109 7564 : END FUNCTION ipmpar
110 :
111 : ! **************************************************************************************************
112 : !> \brief ...
113 : !> \param i ...
114 : !> \return ...
115 : ! **************************************************************************************************
116 7564 : FUNCTION dpmpar(i) RESULT(fn_val)
117 : !-----------------------------------------------------------------------
118 :
119 : ! DPMPAR PROVIDES THE DOUBLE PRECISION MACHINE CONSTANTS FOR
120 : ! THE COMPUTER BEING USED. IT IS ASSUMED THAT THE ARGUMENT
121 : ! I IS AN INTEGER HAVING ONE OF THE VALUES 1, 2, OR 3. IF THE
122 : ! DOUBLE PRECISION ARITHMETIC BEING USED HAS M BASE B DIGITS AND
123 : ! ITS SMALLEST AND LARGEST EXPONENTS ARE EMIN AND EMAX, THEN
124 :
125 : ! DPMPAR(1) = B**(1 - M), THE MACHINE PRECISION,
126 :
127 : ! DPMPAR(2) = B**(EMIN - 1), THE SMALLEST MAGNITUDE,
128 :
129 : ! DPMPAR(3) = B**EMAX*(1 - B**(-M)), THE LARGEST MAGNITUDE.
130 : !-----------------------------------------------------------------------
131 :
132 : INTEGER, INTENT(IN) :: i
133 : REAL(dp) :: fn_val
134 :
135 : REAL(dp), PARAMETER :: one = 1._dp
136 :
137 : ! Local variable
138 :
139 15128 : SELECT CASE (i)
140 : CASE (1)
141 7564 : fn_val = EPSILON(one)
142 : CASE (2)
143 0 : fn_val = TINY(one)
144 : CASE (3)
145 7564 : fn_val = HUGE(one)
146 : END SELECT
147 :
148 7564 : END FUNCTION dpmpar
149 :
150 : ! **************************************************************************************************
151 : !> \brief ...
152 : !> \param l ...
153 : !> \return ...
154 : ! **************************************************************************************************
155 0 : FUNCTION dxparg(l) RESULT(fn_val)
156 : !--------------------------------------------------------------------
157 : ! IF L = 0 THEN DXPARG(L) = THE LARGEST POSITIVE W FOR WHICH
158 : ! DEXP(W) CAN BE COMPUTED.
159 : !
160 : ! IF L IS NONZERO THEN DXPARG(L) = THE LARGEST NEGATIVE W FOR
161 : ! WHICH THE COMPUTED VALUE OF DEXP(W) IS NONZERO.
162 : !
163 : ! NOTE... ONLY AN APPROXIMATE VALUE FOR DXPARG(L) IS NEEDED.
164 : !--------------------------------------------------------------------
165 : INTEGER, INTENT(IN) :: l
166 : REAL(dp) :: fn_val
167 :
168 : REAL(dp), PARAMETER :: one = 1._dp
169 :
170 : ! Local variable
171 :
172 0 : IF (l == 0) THEN
173 : fn_val = LOG(HUGE(one))
174 : ELSE
175 0 : fn_val = LOG(TINY(one))
176 : END IF
177 :
178 0 : END FUNCTION dxparg
179 :
180 : ! **************************************************************************************************
181 : !> \brief ...
182 : !> \param a ...
183 : !> \return ...
184 : ! **************************************************************************************************
185 0 : FUNCTION alnrel(a) RESULT(fn_val)
186 : !-----------------------------------------------------------------------
187 : ! EVALUATION OF THE FUNCTION LN(1 + A)
188 : !-----------------------------------------------------------------------
189 : REAL(dp), INTENT(IN) :: a
190 : REAL(dp) :: fn_val
191 :
192 : REAL(dp), PARAMETER :: half = 0.5e0_dp, one = 1.e0_dp, p1 = -.129418923021993e+01_dp, &
193 : p2 = .405303492862024e+00_dp, p3 = -.178874546012214e-01_dp, &
194 : q1 = -.162752256355323e+01_dp, q2 = .747811014037616e+00_dp, &
195 : q3 = -.845104217945565e-01_dp, two = 2.e0_dp, zero = 0.e0_dp
196 :
197 : REAL(dp) :: t, t2, w, x
198 :
199 : !--------------------------
200 :
201 0 : IF (ABS(a) <= 0.375e0_dp) THEN
202 0 : t = a/(a + two)
203 0 : t2 = t*t
204 0 : w = (((p3*t2 + p2)*t2 + p1)*t2 + one)/(((q3*t2 + q2)*t2 + q1)*t2 + one)
205 0 : fn_val = two*t*w
206 : ELSE
207 0 : x = one + a
208 0 : IF (a < zero) x = (a + half) + half
209 0 : fn_val = LOG(x)
210 : END IF
211 :
212 0 : END FUNCTION alnrel
213 :
214 : ! **************************************************************************************************
215 : !> \brief ...
216 : !> \param a ...
217 : !> \return ...
218 : ! **************************************************************************************************
219 0 : FUNCTION gam1(a) RESULT(fn_val)
220 : !-----------------------------------------------------------------------
221 : ! COMPUTATION OF 1/GAMMA(A+1) - 1 FOR -0.5 <= A <= 1.5
222 : !-----------------------------------------------------------------------
223 : REAL(dp), INTENT(IN) :: a
224 : REAL(dp) :: fn_val
225 :
226 : REAL(dp), PARAMETER :: p(7) = (/.577215664901533e+00_dp, -.409078193005776e+00_dp, &
227 : -.230975380857675e+00_dp, .597275330452234e-01_dp, .766968181649490e-02_dp, &
228 : -.514889771323592e-02_dp, .589597428611429e-03_dp/), q(5) = (/.100000000000000e+01_dp, &
229 : .427569613095214e+00_dp, .158451672430138e+00_dp, .261132021441447e-01_dp, &
230 : .423244297896961e-02_dp/), r(9) = (/-.422784335098468e+00_dp, -.771330383816272e+00_dp, &
231 : -.244757765222226e+00_dp, .118378989872749e+00_dp, .930357293360349e-03_dp, &
232 : -.118290993445146e-01_dp, .223047661158249e-02_dp, .266505979058923e-03_dp, &
233 : -.132674909766242e-03_dp/)
234 : REAL(dp), PARAMETER :: s1 = .273076135303957e+00_dp, s2 = .559398236957378e-01_dp
235 :
236 : REAL(dp) :: bot, d, t, top, w
237 :
238 0 : t = a
239 0 : d = a - 0.5e0_dp
240 0 : IF (d > 0.0e0_dp) t = d - 0.5e0_dp
241 :
242 0 : IF (t > 0.e0_dp) THEN
243 0 : top = (((((p(7)*t + p(6))*t + p(5))*t + p(4))*t + p(3))*t + p(2))*t + p(1)
244 0 : bot = (((q(5)*t + q(4))*t + q(3))*t + q(2))*t + 1.0e0_dp
245 0 : w = top/bot
246 0 : IF (d > 0.0e0_dp) THEN
247 0 : fn_val = (t/a)*((w - 0.5e0_dp) - 0.5e0_dp)
248 : ELSE
249 0 : fn_val = a*w
250 : END IF
251 0 : ELSE IF (t < 0.e0_dp) THEN
252 : top = (((((((r(9)*t + r(8))*t + r(7))*t + r(6))*t + r(5))*t &
253 0 : + r(4))*t + r(3))*t + r(2))*t + r(1)
254 0 : bot = (s2*t + s1)*t + 1.0e0_dp
255 0 : w = top/bot
256 0 : IF (d > 0.0e0_dp) THEN
257 0 : fn_val = t*w/a
258 : ELSE
259 0 : fn_val = a*((w + 0.5e0_dp) + 0.5e0_dp)
260 : END IF
261 : ELSE
262 : fn_val = 0.0e0_dp
263 : END IF
264 :
265 0 : END FUNCTION gam1
266 :
267 : ! **************************************************************************************************
268 : !> \brief ...
269 : !> \param a ...
270 : !> \param b ...
271 : !> \return ...
272 : ! **************************************************************************************************
273 0 : FUNCTION algdiv(a, b) RESULT(fn_val)
274 : !-----------------------------------------------------------------------
275 :
276 : ! COMPUTATION OF LN(GAMMA(B)/GAMMA(A+B)) WHEN B >= 8
277 :
278 : ! --------
279 :
280 : ! IN THIS ALGORITHM, DEL(X) IS THE FUNCTION DEFINED BY
281 : ! LN(GAMMA(X)) = (X - 0.5)*LN(X) - X + 0.5*LN(2*PI) + DEL(X).
282 :
283 : !-----------------------------------------------------------------------
284 : REAL(dp), INTENT(IN) :: a, b
285 : REAL(dp) :: fn_val
286 :
287 : REAL(dp), PARAMETER :: c0 = .833333333333333e-01_dp, c1 = -.277777777760991e-02_dp, &
288 : c2 = .793650666825390e-03_dp, c3 = -.595202931351870e-03_dp, &
289 : c4 = .837308034031215e-03_dp, c5 = -.165322962780713e-02_dp
290 :
291 : REAL(dp) :: c, d, h, s11, s3, s5, s7, s9, t, u, v, &
292 : w, x, x2
293 :
294 0 : IF (a > b) THEN
295 0 : h = b/a
296 0 : c = 1.0e0_dp/(1.0e0_dp + h)
297 0 : x = h/(1.0e0_dp + h)
298 0 : d = a + (b - 0.5e0_dp)
299 : ELSE
300 0 : h = a/b
301 0 : c = h/(1.0e0_dp + h)
302 0 : x = 1.0e0_dp/(1.0e0_dp + h)
303 0 : d = b + (a - 0.5e0_dp)
304 : END IF
305 :
306 : ! SET SN = (1 - X**N)/(1 - X)
307 :
308 0 : x2 = x*x
309 0 : s3 = 1.0e0_dp + (x + x2)
310 0 : s5 = 1.0e0_dp + (x + x2*s3)
311 0 : s7 = 1.0e0_dp + (x + x2*s5)
312 0 : s9 = 1.0e0_dp + (x + x2*s7)
313 0 : s11 = 1.0e0_dp + (x + x2*s9)
314 :
315 : ! SET W = DEL(B) - DEL(A + B)
316 :
317 0 : t = (1.0e0_dp/b)**2
318 0 : w = ((((c5*s11*t + c4*s9)*t + c3*s7)*t + c2*s5)*t + c1*s3)*t + c0
319 0 : w = w*(c/b)
320 :
321 : ! COMBINE THE RESULTS
322 :
323 0 : u = d*alnrel(a/b)
324 0 : v = a*(LOG(b) - 1.0e0_dp)
325 0 : IF (u > v) THEN
326 0 : fn_val = (w - v) - u
327 : ELSE
328 0 : fn_val = (w - u) - v
329 : END IF
330 :
331 0 : END FUNCTION algdiv
332 :
333 : ! **************************************************************************************************
334 : !> \brief ...
335 : !> \param x ...
336 : !> \return ...
337 : ! **************************************************************************************************
338 0 : FUNCTION rexp(x) RESULT(fn_val)
339 : !-----------------------------------------------------------------------
340 : ! EVALUATION OF THE FUNCTION EXP(X) - 1
341 : !-----------------------------------------------------------------------
342 : REAL(dp), INTENT(IN) :: x
343 : REAL(dp) :: fn_val
344 :
345 : REAL(dp), PARAMETER :: p1 = .914041914819518e-09_dp, p2 = .238082361044469e-01_dp, &
346 : q1 = -.499999999085958e+00_dp, q2 = .107141568980644e+00_dp, &
347 : q3 = -.119041179760821e-01_dp, q4 = .595130811860248e-03_dp
348 :
349 : REAL(dp) :: e
350 :
351 0 : IF (ABS(x) < 0.15e0_dp) THEN
352 0 : fn_val = x*(((p2*x + p1)*x + 1.0e0_dp)/((((q4*x + q3)*x + q2)*x + q1)*x + 1.0e0_dp))
353 0 : RETURN
354 : END IF
355 :
356 0 : IF (x < 0.0e0_dp) THEN
357 0 : IF (x > -37.0e0_dp) THEN
358 0 : fn_val = (EXP(x) - 0.5e0_dp) - 0.5e0_dp
359 : ELSE
360 : fn_val = -1.0e0_dp
361 : END IF
362 : ELSE
363 0 : e = EXP(x)
364 0 : fn_val = e*(0.5e0_dp + (0.5e0_dp - 1.0e0_dp/e))
365 : END IF
366 :
367 : END FUNCTION rexp
368 :
369 : ! **************************************************************************************************
370 : !> \brief ...
371 : !> \param a ...
372 : !> \param b ...
373 : !> \param x ...
374 : !> \param y ...
375 : !> \param w ...
376 : !> \param eps ...
377 : !> \param ierr ...
378 : ! **************************************************************************************************
379 0 : SUBROUTINE bgrat(a, b, x, y, w, eps, ierr)
380 : !-----------------------------------------------------------------------
381 : ! ASYMPTOTIC EXPANSION FOR IX(A,B) WHEN A IS LARGER THAN B.
382 : ! THE RESULT OF THE EXPANSION IS ADDED TO W. IT IS ASSUMED
383 : ! THAT A <= 15 AND B <= 1. EPS IS THE TOLERANCE USED.
384 : ! IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
385 : !-----------------------------------------------------------------------
386 : REAL(dp), INTENT(IN) :: a, b, x, y
387 : REAL(dp), INTENT(INOUT) :: w
388 : REAL(dp), INTENT(IN) :: eps
389 : INTEGER, INTENT(OUT) :: ierr
390 :
391 : REAL(dp), PARAMETER :: half = 0.5e0_dp, one = 1.e0_dp, &
392 : quarter = 0.25e0_dp, zero = 0.e0_dp
393 :
394 : INTEGER :: i, n
395 : REAL(dp) :: bm1, bp2n, c(30), cn, coef, d(30), dj, &
396 : j, l, lnx, n2, nu, q, r, s, sum, t, &
397 : t2, tol, u, v, z
398 :
399 0 : bm1 = (b - half) - half
400 0 : nu = a + half*bm1
401 0 : IF (y > 0.375e0_dp) THEN
402 0 : lnx = LOG(x)
403 : ELSE
404 0 : lnx = alnrel(-y)
405 : END IF
406 0 : z = -nu*lnx
407 0 : IF (b*z == zero) THEN
408 : ! THE EXPANSION CANNOT BE COMPUTED
409 0 : ierr = 1
410 0 : RETURN
411 : END IF
412 :
413 : ! COMPUTATION OF THE EXPANSION
414 : ! SET R = EXP(-Z)*Z**B/GAMMA(B)
415 :
416 0 : r = b*(one + gam1(b))*EXP(b*LOG(z))
417 0 : r = r*EXP(a*lnx)*EXP(half*bm1*lnx)
418 0 : u = algdiv(b, a) + b*LOG(nu)
419 0 : u = r*EXP(-u)
420 0 : IF (u == zero) THEN
421 : ! THE EXPANSION CANNOT BE COMPUTED
422 0 : ierr = 1
423 0 : RETURN
424 : END IF
425 0 : CALL grat1(b, z, r, q=q, eps=eps)
426 :
427 0 : tol = 15.0e0_dp*eps
428 0 : v = quarter*(one/nu)**2
429 0 : t2 = quarter*lnx*lnx
430 0 : l = w/u
431 0 : j = q/r
432 0 : sum = j
433 0 : t = one
434 0 : cn = one
435 0 : n2 = zero
436 0 : DO n = 1, 30
437 0 : bp2n = b + n2
438 0 : j = (bp2n*(bp2n + one)*j + (z + bp2n + one)*t)*v
439 0 : n2 = n2 + 2.0e0_dp
440 0 : t = t*t2
441 0 : cn = cn/(n2*(n2 + one))
442 0 : c(n) = cn
443 0 : s = zero
444 0 : IF (.NOT. (n == 1)) THEN
445 0 : coef = b - n
446 0 : DO i = 1, n - 1
447 0 : s = s + coef*c(i)*d(n - i)
448 0 : coef = coef + b
449 : END DO
450 : END IF
451 0 : d(n) = bm1*cn + s/n
452 0 : dj = d(n)*j
453 0 : sum = sum + dj
454 0 : IF (sum <= zero) THEN
455 : ! THE EXPANSION CANNOT BE COMPUTED
456 0 : ierr = 1
457 0 : RETURN
458 : END IF
459 0 : IF (ABS(dj) <= tol*(sum + l)) EXIT
460 : END DO
461 :
462 : ! ADD THE RESULTS TO W
463 :
464 0 : ierr = 0
465 0 : w = w + u*sum
466 :
467 : END SUBROUTINE bgrat
468 :
469 : ! **************************************************************************************************
470 : !> \brief ...
471 : !> \param a ...
472 : !> \param x ...
473 : !> \param r ...
474 : !> \param p ...
475 : !> \param q ...
476 : !> \param eps ...
477 : ! **************************************************************************************************
478 0 : SUBROUTINE grat1(a, x, r, p, q, eps)
479 : !-----------------------------------------------------------------------
480 : ! EVALUATION OF P(A,X) AND Q(A,X) WHERE A <= 1 AND
481 : ! THE INPUT ARGUMENT R HAS THE VALUE E**(-X)*X**A/GAMMA(A)
482 : !-----------------------------------------------------------------------
483 : REAL(dp), INTENT(IN) :: a, x, r
484 : REAL(dp), INTENT(OUT), OPTIONAL :: p, q
485 : REAL(dp), INTENT(IN) :: eps
486 :
487 : REAL(dp), PARAMETER :: half = 0.5e0_dp, one = 1.e0_dp, &
488 : quarter = 0.25e0_dp, three = 3.e0_dp, &
489 : two = 2.e0_dp, zero = 0.e0_dp
490 :
491 : REAL(dp) :: a2n, a2nm1, an, b2n, b2nm1, c, g, h, j, &
492 : l, pp, qq, sum, t, tol, w, z
493 :
494 0 : IF (a*x == zero) THEN
495 0 : IF (x <= a) THEN
496 : pp = zero
497 : qq = one
498 : ELSE
499 0 : pp = one
500 0 : qq = zero
501 : END IF
502 0 : IF (PRESENT(p)) p = pp
503 0 : IF (PRESENT(q)) q = qq
504 0 : RETURN
505 : END IF
506 0 : IF (a == half) THEN
507 0 : IF (x < quarter) THEN
508 0 : pp = ERF(SQRT(x))
509 0 : qq = half + (half - pp)
510 : ELSE
511 0 : qq = ERFC(SQRT(x))
512 0 : pp = half + (half - qq)
513 : END IF
514 0 : IF (PRESENT(p)) p = pp
515 0 : IF (PRESENT(q)) q = qq
516 0 : RETURN
517 : END IF
518 0 : IF (x < 1.1e0_dp) THEN
519 : ! TAYLOR SERIES FOR P(A,X)/X**A
520 :
521 0 : an = three
522 0 : c = x
523 0 : sum = x/(a + three)
524 0 : tol = three*eps/(a + one)
525 0 : an = an + one
526 0 : c = -c*(x/an)
527 0 : t = c/(a + an)
528 0 : sum = sum + t
529 0 : DO WHILE (ABS(t) > tol)
530 0 : an = an + one
531 0 : c = -c*(x/an)
532 0 : t = c/(a + an)
533 0 : sum = sum + t
534 : END DO
535 0 : j = a*x*((sum/6.0e0_dp - half/(a + two))*x + one/(a + one))
536 :
537 0 : z = a*LOG(x)
538 0 : h = gam1(a)
539 0 : g = one + h
540 0 : IF ((x < quarter .AND. z > -.13394e0_dp) .OR. a < x/2.59e0_dp) THEN
541 0 : l = rexp(z)
542 0 : qq = ((half + (half + l))*j - l)*g - h
543 0 : IF (qq <= zero) THEN
544 : pp = one
545 : qq = zero
546 : ELSE
547 0 : pp = half + (half - qq)
548 : END IF
549 : ELSE
550 0 : w = EXP(z)
551 0 : pp = w*g*(half + (half - j))
552 0 : qq = half + (half - pp)
553 : END IF
554 : ELSE
555 : ! CONTINUED FRACTION EXPANSION
556 :
557 0 : tol = 8.0e0_dp*eps
558 0 : a2nm1 = one
559 0 : a2n = one
560 0 : b2nm1 = x
561 0 : b2n = x + (one - a)
562 0 : c = one
563 : DO
564 0 : a2nm1 = x*a2n + c*a2nm1
565 0 : b2nm1 = x*b2n + c*b2nm1
566 0 : c = c + one
567 0 : a2n = a2nm1 + (c - a)*a2n
568 0 : b2n = b2nm1 + (c - a)*b2n
569 0 : a2nm1 = a2nm1/b2n
570 0 : b2nm1 = b2nm1/b2n
571 0 : a2n = a2n/b2n
572 0 : b2n = one
573 0 : IF (ABS(a2n - a2nm1/b2nm1) < tol*a2n) EXIT
574 : END DO
575 :
576 0 : qq = r*a2n
577 0 : pp = half + (half - qq)
578 : END IF
579 :
580 0 : IF (PRESENT(p)) p = pp
581 0 : IF (PRESENT(q)) q = qq
582 :
583 : END SUBROUTINE grat1
584 :
585 : ! **************************************************************************************************
586 : !> \brief ...
587 : !> \param mu ...
588 : !> \param x ...
589 : !> \return ...
590 : ! **************************************************************************************************
591 0 : FUNCTION esum(mu, x) RESULT(fn_val)
592 : !-----------------------------------------------------------------------
593 : ! EVALUATION OF EXP(MU + X)
594 : !-----------------------------------------------------------------------
595 : INTEGER, INTENT(IN) :: mu
596 : REAL(dp), INTENT(IN) :: x
597 : REAL(dp) :: fn_val
598 :
599 : REAL(dp) :: w
600 :
601 0 : IF (x > 0.0e0_dp) THEN
602 0 : IF (mu > 0 .OR. mu + x < 0.0_dp) THEN
603 0 : w = mu
604 0 : fn_val = EXP(w)*EXP(x)
605 : ELSE
606 0 : w = mu + x
607 0 : fn_val = EXP(w)
608 : END IF
609 : ELSE
610 0 : IF (mu < 0 .OR. mu + x < 0.0_dp) THEN
611 0 : w = mu
612 0 : fn_val = EXP(w)*EXP(x)
613 : ELSE
614 0 : w = mu + x
615 0 : fn_val = EXP(w)
616 : END IF
617 : END IF
618 :
619 0 : END FUNCTION esum
620 :
621 : ! **************************************************************************************************
622 : !> \brief ...
623 : !> \param x ...
624 : !> \return ...
625 : ! **************************************************************************************************
626 0 : FUNCTION rlog1(x) RESULT(fn_val)
627 : !-----------------------------------------------------------------------
628 : ! EVALUATION OF THE FUNCTION X - LN(1 + X)
629 : !-----------------------------------------------------------------------
630 : ! A = RLOG (0.7)
631 : ! B = RLOG (4/3)
632 : !------------------------
633 : REAL(dp), INTENT(IN) :: x
634 : REAL(dp) :: fn_val
635 :
636 : REAL(dp), PARAMETER :: a = .566749439387324e-01_dp, b = .456512608815524e-01_dp, &
637 : p0 = .333333333333333e+00_dp, p1 = -.224696413112536e+00_dp, &
638 : p2 = .620886815375787e-02_dp, q1 = -.127408923933623e+01_dp, q2 = .354508718369557e+00_dp
639 :
640 : REAL(dp) :: r, t, u, up2, w, w1
641 :
642 0 : IF (x < -0.39e0_dp .OR. x > 0.57e0_dp) THEN
643 0 : w = (x + 0.5e0_dp) + 0.5e0_dp
644 0 : fn_val = x - LOG(w)
645 0 : RETURN
646 : END IF
647 :
648 : ! ARGUMENT REDUCTION
649 :
650 0 : IF (x < -0.18e0_dp) THEN
651 0 : u = (x + 0.3e0_dp)/0.7e0_dp
652 0 : up2 = u + 2.0e0_dp
653 0 : w1 = a - u*0.3e0_dp
654 0 : ELSEIF (x > 0.18e0_dp) THEN
655 0 : t = 0.75e0_dp*x
656 0 : u = t - 0.25e0_dp
657 0 : up2 = t + 1.75e0_dp
658 0 : w1 = b + u/3.0e0_dp
659 : ELSE
660 0 : u = x
661 0 : up2 = u + 2.0e0_dp
662 0 : w1 = 0.0e0_dp
663 : END IF
664 :
665 : ! SERIES EXPANSION
666 :
667 0 : r = u/up2
668 0 : t = r*r
669 0 : w = ((p2*t + p1)*t + p0)/((q2*t + q1)*t + 1.0e0_dp)
670 0 : fn_val = r*(u - 2.0e0_dp*t*w) + w1
671 :
672 0 : END FUNCTION rlog1
673 :
674 : ! **************************************************************************************************
675 : !> \brief ...
676 : !> \param a ...
677 : !> \return ...
678 : ! **************************************************************************************************
679 0 : FUNCTION gamln1(a) RESULT(fn_val)
680 : !-----------------------------------------------------------------------
681 : ! EVALUATION OF LN(GAMMA(1 + A)) FOR -0.2 .LE. A .LE. 1.25
682 : !-----------------------------------------------------------------------
683 : REAL(dp), INTENT(IN) :: a
684 : REAL(dp) :: fn_val
685 :
686 : REAL(dp), PARAMETER :: p0 = .577215664901533e+00_dp, p1 = .844203922187225e+00_dp, &
687 : p2 = -.168860593646662e+00_dp, p3 = -.780427615533591e+00_dp, &
688 : p4 = -.402055799310489e+00_dp, p5 = -.673562214325671e-01_dp, &
689 : p6 = -.271935708322958e-02_dp, q1 = .288743195473681e+01_dp, &
690 : q2 = .312755088914843e+01_dp, q3 = .156875193295039e+01_dp, q4 = .361951990101499e+00_dp, &
691 : q5 = .325038868253937e-01_dp, q6 = .667465618796164e-03_dp, r0 = .422784335098467e+00_dp, &
692 : r1 = .848044614534529e+00_dp, r2 = .565221050691933e+00_dp, r3 = .156513060486551e+00_dp, &
693 : r4 = .170502484022650e-01_dp, r5 = .497958207639485e-03_dp
694 : REAL(dp), PARAMETER :: s1 = .124313399877507e+01_dp, s2 = .548042109832463e+00_dp, &
695 : s3 = .101552187439830e+00_dp, s4 = .713309612391000e-02_dp, s5 = .116165475989616e-03_dp
696 :
697 : REAL(dp) :: w, x
698 :
699 0 : IF (a < 0.6e0_dp) THEN
700 : w = ((((((p6*a + p5)*a + p4)*a + p3)*a + p2)*a + p1)*a + p0)/ &
701 0 : ((((((q6*a + q5)*a + q4)*a + q3)*a + q2)*a + q1)*a + 1.0e0_dp)
702 0 : fn_val = -a*w
703 : ELSE
704 0 : x = (a - 0.5e0_dp) - 0.5e0_dp
705 : w = (((((r5*x + r4)*x + r3)*x + r2)*x + r1)*x + r0)/ &
706 0 : (((((s5*x + s4)*x + s3)*x + s2)*x + s1)*x + 1.0e0_dp)
707 0 : fn_val = x*w
708 : END IF
709 :
710 0 : END FUNCTION gamln1
711 :
712 : ! **************************************************************************************************
713 : !> \brief ...
714 : !> \param xx ...
715 : !> \return ...
716 : ! **************************************************************************************************
717 7564 : FUNCTION psi(xx) RESULT(fn_val)
718 : !---------------------------------------------------------------------
719 :
720 : ! EVALUATION OF THE DIGAMMA FUNCTION
721 :
722 : ! -----------
723 :
724 : ! PSI(XX) IS ASSIGNED THE VALUE 0 WHEN THE DIGAMMA FUNCTION CANNOT
725 : ! BE COMPUTED.
726 :
727 : ! THE MAIN COMPUTATION INVOLVES EVALUATION OF RATIONAL CHEBYSHEV
728 : ! APPROXIMATIONS PUBLISHED IN MATH. COMP. 27, 123-127(1973) BY
729 : ! CODY, STRECOK AND THACHER.
730 :
731 : !---------------------------------------------------------------------
732 : ! PSI WAS WRITTEN AT ARGONNE NATIONAL LABORATORY FOR THE FUNPACK
733 : ! PACKAGE OF SPECIAL FUNCTION SUBROUTINES. PSI WAS MODIFIED BY
734 : ! A.H. MORRIS (NSWC).
735 : !---------------------------------------------------------------------
736 : REAL(dp), INTENT(IN) :: xx
737 : REAL(dp) :: fn_val
738 :
739 : REAL(dp), PARAMETER :: dx0 = 1.461632144968362341262659542325721325e0_dp, p1(7) = (/ &
740 : .895385022981970e-02_dp, .477762828042627e+01_dp, .142441585084029e+03_dp, &
741 : .118645200713425e+04_dp, .363351846806499e+04_dp, .413810161269013e+04_dp, &
742 : .130560269827897e+04_dp/), p2(4) = (/-.212940445131011e+01_dp, -.701677227766759e+01_dp, &
743 : -.448616543918019e+01_dp, -.648157123766197e+00_dp/), piov4 = .785398163397448e0_dp, q1(6)&
744 : = (/.448452573429826e+02_dp, .520752771467162e+03_dp, .221000799247830e+04_dp, &
745 : .364127349079381e+04_dp, .190831076596300e+04_dp, .691091682714533e-05_dp/)
746 : REAL(dp), PARAMETER :: q2(4) = (/.322703493791143e+02_dp, .892920700481861e+02_dp, &
747 : .546117738103215e+02_dp, .777788548522962e+01_dp/)
748 :
749 : INTEGER :: i, m, n, nq
750 : REAL(dp) :: aug, den, sgn, upper, w, x, xmax1, xmx0, &
751 : xsmall, z
752 :
753 : !---------------------------------------------------------------------
754 : ! PIOV4 = PI/4
755 : ! DX0 = ZERO OF PSI TO EXTENDED PRECISION
756 : !---------------------------------------------------------------------
757 : !---------------------------------------------------------------------
758 : ! COEFFICIENTS FOR RATIONAL APPROXIMATION OF
759 : ! PSI(X) / (X - X0), 0.5 <= X <= 3.0
760 : !---------------------------------------------------------------------
761 : !---------------------------------------------------------------------
762 : ! COEFFICIENTS FOR RATIONAL APPROXIMATION OF
763 : ! PSI(X) - LN(X) + 1 / (2*X), X > 3.0
764 : !---------------------------------------------------------------------
765 : !---------------------------------------------------------------------
766 : ! MACHINE DEPENDENT CONSTANTS ...
767 : ! XMAX1 = THE SMALLEST POSITIVE FLOATING POINT CONSTANT
768 : ! WITH ENTIRELY INTEGER REPRESENTATION. ALSO USED
769 : ! AS NEGATIVE OF LOWER BOUND ON ACCEPTABLE NEGATIVE
770 : ! ARGUMENTS AND AS THE POSITIVE ARGUMENT BEYOND WHICH
771 : ! PSI MAY BE REPRESENTED AS ALOG(X).
772 : ! XSMALL = ABSOLUTE ARGUMENT BELOW WHICH PI*COTAN(PI*X)
773 : ! MAY BE REPRESENTED BY 1/X.
774 : !---------------------------------------------------------------------
775 :
776 7564 : xmax1 = ipmpar(3)
777 7564 : xmax1 = MIN(xmax1, 1.0e0_dp/dpmpar(1))
778 7564 : xsmall = 1.e-9_dp
779 : !---------------------------------------------------------------------
780 7564 : x = xx
781 7564 : aug = 0.0e0_dp
782 7564 : IF (x < 0.5e0_dp) THEN
783 : !---------------------------------------------------------------------
784 : ! X .LT. 0.5, USE REFLECTION FORMULA
785 : ! PSI(1-X) = PSI(X) + PI * COTAN(PI*X)
786 : !---------------------------------------------------------------------
787 0 : IF (ABS(x) <= xsmall) THEN
788 0 : IF (x == 0.0e0_dp) THEN
789 : ! ERROR RETURN
790 7564 : fn_val = 0.0e0_dp
791 : RETURN
792 : END IF
793 : !---------------------------------------------------------------------
794 : ! 0 .LT. ABS(X) .LE. XSMALL. USE 1/X AS A SUBSTITUTE
795 : ! FOR PI*COTAN(PI*X)
796 : !---------------------------------------------------------------------
797 0 : aug = -1.0e0_dp/x
798 0 : x = 1.0e0_dp - x
799 : ELSE
800 : !---------------------------------------------------------------------
801 : ! REDUCTION OF ARGUMENT FOR COTAN
802 : !---------------------------------------------------------------------
803 0 : w = -x
804 0 : sgn = piov4
805 0 : IF (w <= 0.0e0_dp) THEN
806 0 : w = -w
807 0 : sgn = -sgn
808 : END IF
809 : !---------------------------------------------------------------------
810 : ! MAKE AN ERROR EXIT IF X .LE. -XMAX1
811 : !---------------------------------------------------------------------
812 0 : IF (w >= xmax1) THEN
813 : ! ERROR RETURN
814 7564 : fn_val = 0.0e0_dp
815 : RETURN
816 : END IF
817 0 : nq = INT(w)
818 0 : w = w - nq
819 0 : nq = INT(w*4.0e0_dp)
820 0 : w = 4.0e0_dp*(w - nq*.25e0_dp)
821 : !---------------------------------------------------------------------
822 : ! W IS NOW RELATED TO THE FRACTIONAL PART OF 4.0 * X.
823 : ! ADJUST ARGUMENT TO CORRESPOND TO VALUES IN FIRST
824 : ! QUADRANT AND DETERMINE SIGN
825 : !---------------------------------------------------------------------
826 0 : n = nq/2
827 0 : IF ((n + n) /= nq) w = 1.0e0_dp - w
828 0 : z = piov4*w
829 0 : m = n/2
830 0 : IF ((m + m) /= n) sgn = -sgn
831 : !---------------------------------------------------------------------
832 : ! DETERMINE FINAL VALUE FOR -PI*COTAN(PI*X)
833 : !---------------------------------------------------------------------
834 0 : n = (nq + 1)/2
835 0 : m = n/2
836 0 : m = m + m
837 0 : IF (m /= n) THEN
838 0 : aug = sgn*((SIN(z)/COS(z))*4.0e0_dp)
839 : ELSE
840 : !---------------------------------------------------------------------
841 : ! CHECK FOR SINGULARITY
842 : !---------------------------------------------------------------------
843 0 : IF (z == 0.0e0_dp) THEN
844 : ! ERROR RETURN
845 7564 : fn_val = 0.0e0_dp
846 : RETURN
847 : END IF
848 : !---------------------------------------------------------------------
849 : ! USE COS/SIN AS A SUBSTITUTE FOR COTAN, AND
850 : ! SIN/COS AS A SUBSTITUTE FOR TAN
851 : !---------------------------------------------------------------------
852 0 : aug = sgn*((COS(z)/SIN(z))*4.0e0_dp)
853 : END IF
854 0 : x = 1.0e0_dp - x
855 : END IF
856 : END IF
857 7564 : IF (x <= 3.0e0_dp) THEN
858 : !---------------------------------------------------------------------
859 : ! 0.5 .LE. X .LE. 3.0
860 : !---------------------------------------------------------------------
861 0 : den = x
862 0 : upper = p1(1)*x
863 :
864 0 : DO i = 1, 5
865 0 : den = (den + q1(i))*x
866 0 : upper = (upper + p1(i + 1))*x
867 : END DO
868 :
869 0 : den = (upper + p1(7))/(den + q1(6))
870 0 : xmx0 = x - dx0
871 0 : fn_val = den*xmx0 + aug
872 0 : RETURN
873 : END IF
874 : !---------------------------------------------------------------------
875 : ! IF X .GE. XMAX1, PSI = LN(X)
876 : !---------------------------------------------------------------------
877 7564 : IF (x < xmax1) THEN
878 : !---------------------------------------------------------------------
879 : ! 3.0 .LT. X .LT. XMAX1
880 : !---------------------------------------------------------------------
881 7564 : w = 1.0e0_dp/(x*x)
882 7564 : den = w
883 7564 : upper = p2(1)*w
884 :
885 30256 : DO i = 1, 3
886 22692 : den = (den + q2(i))*w
887 30256 : upper = (upper + p2(i + 1))*w
888 : END DO
889 :
890 7564 : aug = upper/(den + q2(4)) - 0.5e0_dp/x + aug
891 : END IF
892 7564 : fn_val = aug + LOG(x)
893 :
894 7564 : END FUNCTION psi
895 :
896 : ! **************************************************************************************************
897 : !> \brief ...
898 : !> \param a0 ...
899 : !> \param b0 ...
900 : !> \return ...
901 : ! **************************************************************************************************
902 0 : FUNCTION betaln(a0, b0) RESULT(fn_val)
903 : !-----------------------------------------------------------------------
904 : ! EVALUATION OF THE LOGARITHM OF THE BETA FUNCTION
905 : !-----------------------------------------------------------------------
906 : ! E = 0.5*LN(2*PI)
907 : !--------------------------
908 : REAL(dp), INTENT(IN) :: a0, b0
909 : REAL(dp) :: fn_val
910 :
911 : REAL(dp), PARAMETER :: e = .918938533204673e0_dp
912 :
913 : INTEGER :: i, n
914 : REAL(dp) :: a, b, c, h, u, v, w, z
915 :
916 : !--------------------------
917 :
918 0 : a = MIN(a0, b0)
919 0 : b = MAX(a0, b0)
920 : !-----------------------------------------------------------------------
921 : ! PROCEDURE WHEN A .GE. 8
922 : !-----------------------------------------------------------------------
923 0 : IF (a >= 8.0e0_dp) THEN
924 0 : w = bcorr(a, b)
925 0 : h = a/b
926 0 : c = h/(1.0e0_dp + h)
927 0 : u = -(a - 0.5e0_dp)*LOG(c)
928 0 : v = b*alnrel(h)
929 0 : IF (u > v) THEN
930 0 : fn_val = (((-0.5e0_dp*LOG(b) + e) + w) - v) - u
931 : ELSE
932 0 : fn_val = (((-0.5e0_dp*LOG(b) + e) + w) - u) - v
933 : END IF
934 : !-----------------------------------------------------------------------
935 : ! PROCEDURE WHEN A .LT. 1
936 : !-----------------------------------------------------------------------
937 0 : ELSEIF (a < 1.0e0_dp) THEN
938 0 : IF (b < 8.0e0_dp) THEN
939 0 : fn_val = LOG_GAMMA(a) + (LOG_GAMMA(b) - LOG_GAMMA(a + b))
940 : ELSE
941 0 : fn_val = LOG_GAMMA(a) + algdiv(a, b)
942 : END IF
943 : !-----------------------------------------------------------------------
944 : ! PROCEDURE WHEN 1 .LE. A .LT. 8
945 : !-----------------------------------------------------------------------
946 0 : ELSEIF (a <= 2.0e0_dp) THEN
947 0 : IF (b <= 2.0e0_dp) THEN
948 0 : fn_val = LOG_GAMMA(a) + LOG_GAMMA(b) - gsumln(a, b)
949 0 : RETURN
950 : END IF
951 0 : w = 0.0e0_dp
952 0 : IF (b < 8.0e0_dp) THEN
953 : ! REDUCTION OF B WHEN B .LT. 8
954 :
955 0 : n = INT(b - 1.0e0_dp)
956 0 : z = 1.0e0_dp
957 0 : DO i = 1, n
958 0 : b = b - 1.0e0_dp
959 0 : z = z*(b/(a + b))
960 : END DO
961 0 : fn_val = w + LOG(z) + (LOG_GAMMA(a) + (LOG_GAMMA(b) - gsumln(a, b)))
962 0 : RETURN
963 : END IF
964 0 : fn_val = LOG_GAMMA(a) + algdiv(a, b)
965 : ! REDUCTION OF A WHEN B .LE. 1000
966 0 : ELSEIF (b <= 1000.0e0_dp) THEN
967 0 : n = INT(a - 1.0e0_dp)
968 0 : w = 1.0e0_dp
969 0 : DO i = 1, n
970 0 : a = a - 1.0e0_dp
971 0 : h = a/b
972 0 : w = w*(h/(1.0e0_dp + h))
973 : END DO
974 0 : w = LOG(w)
975 0 : IF (b >= 8.0e0_dp) THEN
976 0 : fn_val = w + LOG_GAMMA(a) + algdiv(a, b)
977 0 : RETURN
978 : END IF
979 :
980 : ! REDUCTION OF B WHEN B .LT. 8
981 :
982 0 : n = INT(b - 1.0e0_dp)
983 0 : z = 1.0e0_dp
984 0 : DO i = 1, n
985 0 : b = b - 1.0e0_dp
986 0 : z = z*(b/(a + b))
987 : END DO
988 0 : fn_val = w + LOG(z) + (LOG_GAMMA(a) + (LOG_GAMMA(b) - gsumln(a, b)))
989 : ELSE
990 : ! REDUCTION OF A WHEN B .GT. 1000
991 :
992 0 : n = INT(a - 1.0e0_dp)
993 0 : w = 1.0e0_dp
994 0 : DO i = 1, n
995 0 : a = a - 1.0e0_dp
996 0 : w = w*(a/(1.0e0_dp + a/b))
997 : END DO
998 0 : fn_val = (LOG(w) - n*LOG(b)) + (LOG_GAMMA(a) + algdiv(a, b))
999 : END IF
1000 :
1001 : END FUNCTION betaln
1002 :
1003 : ! **************************************************************************************************
1004 : !> \brief ...
1005 : !> \param a ...
1006 : !> \param b ...
1007 : !> \return ...
1008 : ! **************************************************************************************************
1009 0 : FUNCTION gsumln(a, b) RESULT(fn_val)
1010 : !-----------------------------------------------------------------------
1011 : ! EVALUATION OF THE FUNCTION LN(GAMMA(A + B))
1012 : ! FOR 1 .LE. A .LE. 2 AND 1 .LE. B .LE. 2
1013 : !-----------------------------------------------------------------------
1014 : REAL(dp), INTENT(IN) :: a, b
1015 : REAL(dp) :: fn_val
1016 :
1017 : REAL(dp) :: x
1018 :
1019 0 : x = a + b - 2.e0_dp
1020 0 : IF (x <= 0.25e0_dp) THEN
1021 0 : fn_val = gamln1(1.0e0_dp + x)
1022 0 : ELSEIF (x <= 1.25e0_dp) THEN
1023 0 : fn_val = gamln1(x) + alnrel(x)
1024 : ELSE
1025 0 : fn_val = gamln1(x - 1.0e0_dp) + LOG(x*(1.0e0_dp + x))
1026 : END IF
1027 0 : END FUNCTION gsumln
1028 :
1029 : ! **************************************************************************************************
1030 : !> \brief ...
1031 : !> \param a0 ...
1032 : !> \param b0 ...
1033 : !> \return ...
1034 : ! **************************************************************************************************
1035 0 : FUNCTION bcorr(a0, b0) RESULT(fn_val)
1036 : !-----------------------------------------------------------------------
1037 :
1038 : ! EVALUATION OF DEL(A0) + DEL(B0) - DEL(A0 + B0) WHERE
1039 : ! LN(GAMMA(A)) = (A - 0.5)*LN(A) - A + 0.5*LN(2*PI) + DEL(A).
1040 : ! IT IS ASSUMED THAT A0 .GE. 8 AND B0 .GE. 8.
1041 :
1042 : !-----------------------------------------------------------------------
1043 : REAL(dp), INTENT(IN) :: a0, b0
1044 : REAL(dp) :: fn_val
1045 :
1046 : REAL(dp), PARAMETER :: c0 = .833333333333333e-01_dp, c1 = -.277777777760991e-02_dp, &
1047 : c2 = .793650666825390e-03_dp, c3 = -.595202931351870e-03_dp, &
1048 : c4 = .837308034031215e-03_dp, c5 = -.165322962780713e-02_dp
1049 :
1050 : REAL(dp) :: a, b, c, h, s11, s3, s5, s7, s9, t, w, &
1051 : x, x2
1052 :
1053 0 : a = MIN(a0, b0)
1054 0 : b = MAX(a0, b0)
1055 :
1056 0 : h = a/b
1057 0 : c = h/(1.0e0_dp + h)
1058 0 : x = 1.0e0_dp/(1.0e0_dp + h)
1059 0 : x2 = x*x
1060 :
1061 : ! SET SN = (1 - X**N)/(1 - X)
1062 :
1063 0 : s3 = 1.0e0_dp + (x + x2)
1064 0 : s5 = 1.0e0_dp + (x + x2*s3)
1065 0 : s7 = 1.0e0_dp + (x + x2*s5)
1066 0 : s9 = 1.0e0_dp + (x + x2*s7)
1067 0 : s11 = 1.0e0_dp + (x + x2*s9)
1068 :
1069 : ! SET W = DEL(B) - DEL(A + B)
1070 :
1071 0 : t = (1.0e0_dp/b)**2
1072 0 : w = ((((c5*s11*t + c4*s9)*t + c3*s7)*t + c2*s5)*t + c1*s3)*t + c0
1073 0 : w = w*(c/b)
1074 :
1075 : ! COMPUTE DEL(A) + W
1076 :
1077 0 : t = (1.0e0_dp/a)**2
1078 0 : fn_val = (((((c5*t + c4)*t + c3)*t + c2)*t + c1)*t + c0)/a + w
1079 : RETURN
1080 : END FUNCTION bcorr
1081 :
1082 : ! **************************************************************************************************
1083 : !> \brief ...
1084 : !> \param a ...
1085 : !> \param b ...
1086 : !> \param x ...
1087 : !> \param y ...
1088 : !> \param w ...
1089 : !> \param w1 ...
1090 : !> \param ierr ...
1091 : ! **************************************************************************************************
1092 0 : SUBROUTINE bratio(a, b, x, y, w, w1, ierr)
1093 : !-----------------------------------------------------------------------
1094 :
1095 : ! EVALUATION OF THE INCOMPLETE BETA FUNCTION IX(A,B)
1096 :
1097 : ! --------------------
1098 :
1099 : ! IT IS ASSUMED THAT A AND B ARE NONNEGATIVE, AND THAT X <= 1
1100 : ! AND Y = 1 - X. BRATIO ASSIGNS W AND W1 THE VALUES
1101 :
1102 : ! W = IX(A,B)
1103 : ! W1 = 1 - IX(A,B)
1104 :
1105 : ! IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
1106 : ! IF NO INPUT ERRORS ARE DETECTED THEN IERR IS SET TO 0 AND
1107 : ! W AND W1 ARE COMPUTED. OTHERWISE, IF AN ERROR IS DETECTED,
1108 : ! THEN W AND W1 ARE ASSIGNED THE VALUE 0 AND IERR IS SET TO
1109 : ! ONE OF THE FOLLOWING VALUES ...
1110 :
1111 : ! IERR = 1 IF A OR B IS NEGATIVE
1112 : ! IERR = 2 IF A = B = 0
1113 : ! IERR = 3 IF X .LT. 0 OR X .GT. 1
1114 : ! IERR = 4 IF Y .LT. 0 OR Y .GT. 1
1115 : ! IERR = 5 IF X + Y .NE. 1
1116 : ! IERR = 6 IF X = A = 0
1117 : ! IERR = 7 IF Y = B = 0
1118 :
1119 : !--------------------
1120 : ! WRITTEN BY ALFRED H. MORRIS, JR.
1121 : ! NAVAL SURFACE WARFARE CENTER
1122 : ! DAHLGREN, VIRGINIA
1123 : ! REVISED ... APRIL 1993
1124 : !-----------------------------------------------------------------------
1125 : REAL(dp), INTENT(IN) :: a, b, x, y
1126 : REAL(dp), INTENT(OUT) :: w, w1
1127 : INTEGER, INTENT(OUT) :: ierr
1128 :
1129 : INTEGER :: ierr1, ind, n
1130 : REAL(dp) :: a0, b0, eps, lambda, t, x0, y0, z
1131 :
1132 : !-----------------------------------------------------------------------
1133 : ! ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE SMALLEST
1134 : ! FLOATING POINT NUMBER FOR WHICH 1.0 + EPS .GT. 1.0
1135 :
1136 0 : eps = dpmpar(1)
1137 :
1138 0 : w = 0.0e0_dp
1139 0 : w1 = 0.0e0_dp
1140 0 : IF (a < 0.0e0_dp .OR. b < 0.0e0_dp) THEN
1141 0 : ierr = 1
1142 0 : RETURN
1143 : END IF
1144 0 : IF (a == 0.0e0_dp .AND. b == 0.0e0_dp) THEN
1145 0 : ierr = 2
1146 0 : RETURN
1147 : END IF
1148 0 : IF (x < 0.0e0_dp .OR. x > 1.0e0_dp) THEN
1149 0 : ierr = 3
1150 0 : RETURN
1151 : END IF
1152 0 : IF (y < 0.0e0_dp .OR. y > 1.0e0_dp) THEN
1153 0 : ierr = 4
1154 0 : RETURN
1155 : END IF
1156 0 : z = ((x + y) - 0.5e0_dp) - 0.5e0_dp
1157 0 : IF (ABS(z) > 3.0e0_dp*eps) THEN
1158 0 : ierr = 5
1159 0 : RETURN
1160 : END IF
1161 :
1162 0 : ierr = 0
1163 :
1164 0 : IF (x == 0.0e0_dp .OR. a == 0.0e0_dp) THEN
1165 0 : IF (x == 0.0e0_dp .AND. a == 0.0e0_dp) THEN
1166 0 : ierr = 6
1167 : ELSE
1168 : w = 0.0e0_dp
1169 0 : w1 = 1.0e0_dp
1170 : END IF
1171 0 : RETURN
1172 : END IF
1173 :
1174 0 : IF (y == 0.0e0_dp .OR. b == 0.0e0_dp) THEN
1175 0 : IF (y == 0.0e0_dp .AND. b == 0.0e0_dp) THEN
1176 0 : ierr = 7
1177 : ELSE
1178 0 : w = 1.0e0_dp
1179 : w1 = 0.0e0_dp
1180 : END IF
1181 0 : RETURN
1182 : END IF
1183 :
1184 0 : eps = MAX(eps, 1.e-15_dp)
1185 0 : IF (MAX(a, b) < 1.e-3_dp*eps) THEN
1186 : ! PROCEDURE FOR A AND B .LT. 1.E-3*EPS
1187 0 : w = b/(a + b)
1188 0 : w1 = a/(a + b)
1189 0 : RETURN
1190 : END IF
1191 :
1192 0 : ind = 0
1193 0 : a0 = a
1194 0 : b0 = b
1195 0 : x0 = x
1196 0 : y0 = y
1197 0 : IF (MIN(a0, b0) > 1.0e0_dp) GO TO 30
1198 :
1199 : ! PROCEDURE FOR A0 .LE. 1 OR B0 .LE. 1
1200 :
1201 0 : IF (x > 0.5e0_dp) THEN
1202 0 : ind = 1
1203 0 : a0 = b
1204 0 : b0 = a
1205 0 : x0 = y
1206 0 : y0 = x
1207 : END IF
1208 :
1209 0 : IF (b0 < MIN(eps, eps*a0)) GO TO 80
1210 0 : IF (a0 < MIN(eps, eps*b0) .AND. b0*x0 <= 1.0e0_dp) GO TO 90
1211 0 : IF (MAX(a0, b0) > 1.0e0_dp) GO TO 20
1212 0 : IF (a0 >= MIN(0.2e0_dp, b0)) GO TO 100
1213 0 : IF (x0**a0 <= 0.9e0_dp) GO TO 100
1214 0 : IF (x0 >= 0.3e0_dp) GO TO 110
1215 0 : n = 20
1216 0 : GO TO 130
1217 :
1218 0 : 20 IF (b0 <= 1.0e0_dp) GO TO 100
1219 0 : IF (x0 >= 0.3e0_dp) GO TO 110
1220 0 : IF (x0 < 0.1e0_dp .AND. (x0*b0)**a0 <= 0.7e0_dp) GO TO 100
1221 0 : IF (b0 > 15.0e0_dp) GO TO 131
1222 0 : n = 20
1223 0 : GO TO 130
1224 :
1225 : ! PROCEDURE FOR A0 .GT. 1 AND B0 .GT. 1
1226 :
1227 0 : 30 IF (a > b) THEN
1228 0 : lambda = (a + b)*y - b
1229 : ELSE
1230 0 : lambda = a - (a + b)*x
1231 : END IF
1232 :
1233 0 : IF (lambda < 0.0e0_dp) THEN
1234 0 : ind = 1
1235 0 : a0 = b
1236 0 : b0 = a
1237 0 : x0 = y
1238 0 : y0 = x
1239 0 : lambda = ABS(lambda)
1240 : END IF
1241 :
1242 0 : IF (b0 < 40.0e0_dp .AND. b0*x0 <= 0.7e0_dp) GO TO 100
1243 0 : IF (b0 < 40.0e0_dp) GO TO 140
1244 0 : IF (a0 > b0) GO TO 50
1245 0 : IF (a0 <= 100.0e0_dp) GO TO 120
1246 0 : IF (lambda > 0.03e0_dp*a0) GO TO 120
1247 0 : GO TO 180
1248 0 : 50 IF (b0 <= 100.0e0_dp) GO TO 120
1249 0 : IF (lambda > 0.03e0_dp*b0) GO TO 120
1250 0 : GO TO 180
1251 :
1252 : ! EVALUATION OF THE APPROPRIATE ALGORITHM
1253 :
1254 0 : 80 w = fpser(a0, b0, x0, eps)
1255 0 : w1 = 0.5e0_dp + (0.5e0_dp - w)
1256 0 : GO TO 220
1257 :
1258 0 : 90 w1 = apser(a0, b0, x0, eps)
1259 0 : w = 0.5e0_dp + (0.5e0_dp - w1)
1260 0 : GO TO 220
1261 :
1262 0 : 100 w = bpser(a0, b0, x0, eps)
1263 0 : w1 = 0.5e0_dp + (0.5e0_dp - w)
1264 0 : GO TO 220
1265 :
1266 0 : 110 w1 = bpser(b0, a0, y0, eps)
1267 0 : w = 0.5e0_dp + (0.5e0_dp - w1)
1268 0 : GO TO 220
1269 :
1270 0 : 120 w = bfrac(a0, b0, x0, y0, lambda, 15.0e0_dp*eps)
1271 0 : w1 = 0.5e0_dp + (0.5e0_dp - w)
1272 0 : GO TO 220
1273 :
1274 0 : 130 w1 = bup(b0, a0, y0, x0, n, eps)
1275 0 : b0 = b0 + n
1276 0 : 131 CALL bgrat(b0, a0, y0, x0, w1, eps, ierr1)
1277 0 : IF (ierr1 > 0) CPABORT("Error in BGRAT")
1278 0 : w = 0.5e0_dp + (0.5e0_dp - w1)
1279 0 : GO TO 220
1280 :
1281 0 : 140 n = INT(b0)
1282 0 : b0 = b0 - n
1283 0 : IF (b0 /= 0.0e0_dp) GO TO 141
1284 0 : n = n - 1
1285 0 : b0 = 1.0e0_dp
1286 0 : 141 w = bup(b0, a0, y0, x0, n, eps)
1287 0 : IF (x0 > 0.7e0_dp) GO TO 150
1288 0 : w = w + bpser(a0, b0, x0, eps)
1289 0 : w1 = 0.5e0_dp + (0.5e0_dp - w)
1290 0 : GO TO 220
1291 :
1292 0 : 150 IF (a0 > 15.0e0_dp) GO TO 151
1293 0 : n = 20
1294 0 : w = w + bup(a0, b0, x0, y0, n, eps)
1295 0 : a0 = a0 + n
1296 0 : 151 CALL bgrat(a0, b0, x0, y0, w, eps, ierr1)
1297 0 : w1 = 0.5e0_dp + (0.5e0_dp - w)
1298 0 : GO TO 220
1299 :
1300 0 : 180 w = basym(a0, b0, lambda, 100.0e0_dp*eps)
1301 0 : w1 = 0.5e0_dp + (0.5e0_dp - w)
1302 0 : GO TO 220
1303 :
1304 : ! TERMINATION OF THE PROCEDURE
1305 :
1306 0 : 220 IF (ind == 0) RETURN
1307 0 : t = w
1308 0 : w = w1
1309 0 : w1 = t
1310 :
1311 : END SUBROUTINE bratio
1312 :
1313 : ! **************************************************************************************************
1314 : !> \brief ...
1315 : !> \param a ...
1316 : !> \param b ...
1317 : !> \param x ...
1318 : !> \param eps ...
1319 : !> \return ...
1320 : ! **************************************************************************************************
1321 0 : FUNCTION fpser(a, b, x, eps) RESULT(fn_val)
1322 : !-----------------------------------------------------------------------
1323 :
1324 : ! EVALUATION OF I (A,B)
1325 : ! X
1326 :
1327 : ! FOR B .LT. MIN(EPS,EPS*A) AND X .LE. 0.5.
1328 :
1329 : !-----------------------------------------------------------------------
1330 : REAL(dp), INTENT(IN) :: a, b, x, eps
1331 : REAL(dp) :: fn_val
1332 :
1333 : REAL(dp) :: an, c, s, t, tol
1334 :
1335 : ! SET FPSER = X**A
1336 :
1337 0 : fn_val = 1.0e0_dp
1338 0 : IF (a > 1.e-3_dp*eps) THEN
1339 0 : fn_val = 0.0e0_dp
1340 0 : t = a*LOG(x)
1341 0 : IF (t < dxparg(1)) RETURN
1342 0 : fn_val = EXP(t)
1343 : END IF
1344 :
1345 : ! NOTE THAT 1/B(A,B) = B
1346 :
1347 0 : fn_val = (b/a)*fn_val
1348 0 : tol = eps/a
1349 0 : an = a + 1.0e0_dp
1350 0 : t = x
1351 0 : s = t/an
1352 0 : an = an + 1.0e0_dp
1353 0 : t = x*t
1354 0 : c = t/an
1355 0 : s = s + c
1356 0 : DO WHILE (ABS(c) > tol)
1357 0 : an = an + 1.0e0_dp
1358 0 : t = x*t
1359 0 : c = t/an
1360 0 : s = s + c
1361 : END DO
1362 :
1363 0 : fn_val = fn_val*(1.0e0_dp + a*s)
1364 :
1365 0 : END FUNCTION fpser
1366 :
1367 : ! **************************************************************************************************
1368 : !> \brief ...
1369 : !> \param a ...
1370 : !> \param b ...
1371 : !> \param x ...
1372 : !> \param eps ...
1373 : !> \return ...
1374 : ! **************************************************************************************************
1375 0 : FUNCTION apser(a, b, x, eps) RESULT(fn_val)
1376 : !-----------------------------------------------------------------------
1377 : ! APSER YIELDS THE INCOMPLETE BETA RATIO I(SUB(1-X))(B,A) FOR
1378 : ! A .LE. MIN(EPS,EPS*B), B*X .LE. 1, AND X .LE. 0.5. USED WHEN
1379 : ! A IS VERY SMALL. USE ONLY IF ABOVE INEQUALITIES ARE SATISFIED.
1380 : !-----------------------------------------------------------------------
1381 : REAL(dp), INTENT(IN) :: a, b, x, eps
1382 : REAL(dp) :: fn_val
1383 :
1384 : REAL(dp), PARAMETER :: g = .577215664901533e0_dp
1385 :
1386 : REAL(dp) :: aj, bx, c, j, s, t, tol
1387 :
1388 0 : bx = b*x
1389 0 : t = x - bx
1390 0 : IF (b*eps > 2.e-2_dp) THEN
1391 0 : c = LOG(bx) + g + t
1392 : ELSE
1393 0 : c = LOG(x) + psi(b) + g + t
1394 : END IF
1395 :
1396 0 : tol = 5.0e0_dp*eps*ABS(c)
1397 : j = 1.0e0_dp
1398 0 : s = 0.0e0_dp
1399 0 : j = j + 1.0e0_dp
1400 0 : t = t*(x - bx/j)
1401 0 : aj = t/j
1402 0 : s = s + aj
1403 0 : DO WHILE (ABS(aj) > tol)
1404 0 : t = t*(x - bx/j)
1405 0 : aj = t/j
1406 0 : s = s + aj
1407 : END DO
1408 :
1409 0 : fn_val = -a*(c + s)
1410 :
1411 0 : END FUNCTION apser
1412 :
1413 : ! **************************************************************************************************
1414 : !> \brief ...
1415 : !> \param a ...
1416 : !> \param b ...
1417 : !> \param x ...
1418 : !> \param eps ...
1419 : !> \return ...
1420 : ! **************************************************************************************************
1421 0 : FUNCTION bpser(a, b, x, eps) RESULT(fn_val)
1422 : !-----------------------------------------------------------------------
1423 : ! POWER SERIES EXPANSION FOR EVALUATING IX(A,B) WHEN B .LE. 1
1424 : ! OR B*X .LE. 0.7. EPS IS THE TOLERANCE USED.
1425 : !-----------------------------------------------------------------------
1426 : REAL(dp), INTENT(IN) :: a, b, x, eps
1427 : REAL(dp) :: fn_val
1428 :
1429 : INTEGER :: i, m
1430 : REAL(dp) :: a0, apb, b0, c, n, sum, t, tol, u, w, z
1431 :
1432 0 : fn_val = 0.0e0_dp
1433 0 : IF (x == 0.0e0_dp) RETURN
1434 : !-----------------------------------------------------------------------
1435 : ! COMPUTE THE FACTOR X**A/(A*BETA(A,B))
1436 : !-----------------------------------------------------------------------
1437 0 : a0 = MIN(a, b)
1438 0 : b0 = MAX(a, b)
1439 0 : IF (a0 >= 1.0e0_dp) THEN
1440 0 : z = a*LOG(x) - betaln(a, b)
1441 0 : fn_val = EXP(z)/a
1442 0 : ELSEIF (b0 >= 8.0e0_dp) THEN
1443 0 : u = gamln1(a0) + algdiv(a0, b0)
1444 0 : z = a*LOG(x) - u
1445 0 : fn_val = (a0/a)*EXP(z)
1446 0 : ELSEIF (b0 > 1.0e0_dp) THEN
1447 : ! PROCEDURE FOR A0 .LT. 1 AND 1 .LT. B0 .LT. 8
1448 :
1449 0 : u = gamln1(a0)
1450 0 : m = INT(b0 - 1.0e0_dp)
1451 0 : IF (m >= 1) THEN
1452 : c = 1.0e0_dp
1453 0 : DO i = 1, m
1454 0 : b0 = b0 - 1.0e0_dp
1455 0 : c = c*(b0/(a0 + b0))
1456 : END DO
1457 0 : u = LOG(c) + u
1458 : END IF
1459 :
1460 0 : z = a*LOG(x) - u
1461 0 : b0 = b0 - 1.0e0_dp
1462 0 : apb = a0 + b0
1463 0 : IF (apb > 1.0e0_dp) THEN
1464 0 : u = a0 + b0 - 1.e0_dp
1465 0 : t = (1.0e0_dp + gam1(u))/apb
1466 : ELSE
1467 0 : t = 1.0e0_dp + gam1(apb)
1468 : END IF
1469 0 : fn_val = EXP(z)*(a0/a)*(1.0e0_dp + gam1(b0))/t
1470 : ELSE
1471 :
1472 : ! PROCEDURE FOR A0 .LT. 1 AND B0 .LE. 1
1473 :
1474 0 : fn_val = x**a
1475 0 : IF (fn_val == 0.0e0_dp) RETURN
1476 :
1477 0 : apb = a + b
1478 0 : IF (apb > 1.0e0_dp) THEN
1479 0 : u = a + b - 1.e0_dp
1480 0 : z = (1.0e0_dp + gam1(u))/apb
1481 : ELSE
1482 0 : z = 1.0e0_dp + gam1(apb)
1483 : END IF
1484 :
1485 0 : c = (1.0e0_dp + gam1(a))*(1.0e0_dp + gam1(b))/z
1486 0 : fn_val = fn_val*c*(b/apb)
1487 : END IF
1488 :
1489 : ! PROCEDURE FOR A0 .LT. 1 AND B0 .GE. 8
1490 :
1491 0 : IF (fn_val == 0.0e0_dp .OR. a <= 0.1e0_dp*eps) RETURN
1492 : !-----------------------------------------------------------------------
1493 : ! COMPUTE THE SERIES
1494 : !-----------------------------------------------------------------------
1495 0 : sum = 0.0e0_dp
1496 0 : n = 0.0e0_dp
1497 0 : c = 1.0e0_dp
1498 0 : tol = eps/a
1499 0 : n = n + 1.0e0_dp
1500 0 : c = c*(0.5e0_dp + (0.5e0_dp - b/n))*x
1501 0 : w = c/(a + n)
1502 0 : sum = sum + w
1503 0 : DO WHILE (ABS(w) > tol)
1504 0 : n = n + 1.0e0_dp
1505 0 : c = c*(0.5e0_dp + (0.5e0_dp - b/n))*x
1506 0 : w = c/(a + n)
1507 0 : sum = sum + w
1508 : END DO
1509 0 : fn_val = fn_val*(1.0e0_dp + a*sum)
1510 :
1511 0 : END FUNCTION bpser
1512 :
1513 : ! **************************************************************************************************
1514 : !> \brief ...
1515 : !> \param a ...
1516 : !> \param b ...
1517 : !> \param x ...
1518 : !> \param y ...
1519 : !> \param n ...
1520 : !> \param eps ...
1521 : !> \return ...
1522 : ! **************************************************************************************************
1523 0 : FUNCTION bup(a, b, x, y, n, eps) RESULT(fn_val)
1524 : !-----------------------------------------------------------------------
1525 : ! EVALUATION OF IX(A,B) - IX(A+N,B) WHERE N IS A POSITIVE INTEGER.
1526 : ! EPS IS THE TOLERANCE USED.
1527 : !-----------------------------------------------------------------------
1528 : REAL(dp), INTENT(IN) :: a, b, x, y
1529 : INTEGER, INTENT(IN) :: n
1530 : REAL(dp), INTENT(IN) :: eps
1531 : REAL(dp) :: fn_val
1532 :
1533 : INTEGER :: i, k, kp1, mu, nm1
1534 : REAL(dp) :: ap1, apb, d, l, r, t, w
1535 :
1536 : ! OBTAIN THE SCALING FACTOR EXP(-MU) AND
1537 : ! EXP(MU)*(X**A*Y**B/BETA(A,B))/A
1538 :
1539 0 : apb = a + b
1540 0 : ap1 = a + 1.0e0_dp
1541 0 : mu = 0
1542 0 : d = 1.0e0_dp
1543 0 : IF (.NOT. (n == 1 .OR. a < 1.0e0_dp .OR. apb < 1.1e0_dp*ap1)) THEN
1544 0 : mu = INT(ABS(dxparg(1)))
1545 0 : k = INT(dxparg(0))
1546 0 : IF (k < mu) mu = k
1547 0 : t = mu
1548 0 : d = EXP(-t)
1549 : END IF
1550 :
1551 0 : fn_val = brcmp1(mu, a, b, x, y)/a
1552 0 : IF (n == 1 .OR. fn_val == 0.0e0_dp) RETURN
1553 0 : nm1 = n - 1
1554 0 : w = d
1555 :
1556 : ! LET K BE THE INDEX OF THE MAXIMUM TERM
1557 :
1558 0 : k = 0
1559 0 : IF (b > 1.0e0_dp) THEN
1560 0 : IF (y <= 1.e-4_dp) THEN
1561 0 : k = nm1
1562 0 : DO i = 1, k
1563 0 : l = i - 1
1564 0 : d = ((apb + l)/(ap1 + l))*x*d
1565 0 : w = w + d
1566 : END DO
1567 : IF (k == nm1) THEN
1568 0 : fn_val = fn_val*w
1569 0 : RETURN
1570 : END IF
1571 : ELSE
1572 0 : r = (b - 1.0e0_dp)*x/y - a
1573 0 : IF (r >= 1.0e0_dp) THEN
1574 0 : k = nm1
1575 0 : t = nm1
1576 0 : IF (r < t) k = INT(r)
1577 :
1578 : ! ADD THE INCREASING TERMS OF THE SERIES
1579 :
1580 0 : DO i = 1, k
1581 0 : l = i - 1
1582 0 : d = ((apb + l)/(ap1 + l))*x*d
1583 0 : w = w + d
1584 : END DO
1585 0 : IF (k == nm1) THEN
1586 0 : fn_val = fn_val*w
1587 0 : RETURN
1588 : END IF
1589 : END IF
1590 : END IF
1591 : END IF
1592 :
1593 : ! ADD THE REMAINING TERMS OF THE SERIES
1594 :
1595 0 : kp1 = k + 1
1596 0 : DO i = kp1, nm1
1597 0 : l = i - 1
1598 0 : d = ((apb + l)/(ap1 + l))*x*d
1599 0 : w = w + d
1600 0 : IF (d <= eps*w) EXIT
1601 : END DO
1602 :
1603 : ! TERMINATE THE PROCEDURE
1604 :
1605 0 : fn_val = fn_val*w
1606 :
1607 0 : END FUNCTION bup
1608 :
1609 : ! **************************************************************************************************
1610 : !> \brief ...
1611 : !> \param a ...
1612 : !> \param b ...
1613 : !> \param x ...
1614 : !> \param y ...
1615 : !> \param lambda ...
1616 : !> \param eps ...
1617 : !> \return ...
1618 : ! **************************************************************************************************
1619 0 : FUNCTION bfrac(a, b, x, y, lambda, eps) RESULT(fn_val)
1620 : !-----------------------------------------------------------------------
1621 : ! CONTINUED FRACTION EXPANSION FOR IX(A,B) WHEN A,B .GT. 1.
1622 : ! IT IS ASSUMED THAT LAMBDA = (A + B)*Y - B.
1623 : !-----------------------------------------------------------------------
1624 : REAL(dp), INTENT(IN) :: a, b, x, y, lambda, eps
1625 : REAL(dp) :: fn_val
1626 :
1627 : REAL(dp) :: alpha, an, anp1, beta, bn, bnp1, c, c0, &
1628 : c1, e, n, p, r, r0, s, t, w, yp1
1629 :
1630 0 : fn_val = brcomp(a, b, x, y)
1631 0 : IF (fn_val == 0.0e0_dp) RETURN
1632 :
1633 0 : c = 1.0e0_dp + lambda
1634 0 : c0 = b/a
1635 0 : c1 = 1.0e0_dp + 1.0e0_dp/a
1636 0 : yp1 = y + 1.0e0_dp
1637 :
1638 0 : n = 0.0e0_dp
1639 0 : p = 1.0e0_dp
1640 0 : s = a + 1.0e0_dp
1641 0 : an = 0.0e0_dp
1642 0 : bn = 1.0e0_dp
1643 0 : anp1 = 1.0e0_dp
1644 0 : bnp1 = c/c1
1645 0 : r = c1/c
1646 :
1647 : ! CONTINUED FRACTION CALCULATION
1648 :
1649 0 : DO WHILE (.TRUE.)
1650 0 : n = n + 1.0e0_dp
1651 0 : t = n/a
1652 0 : w = n*(b - n)*x
1653 0 : e = a/s
1654 0 : alpha = (p*(p + c0)*e*e)*(w*x)
1655 0 : IF (alpha <= 0.0e0_dp) THEN
1656 : ! TERMINATION
1657 0 : fn_val = fn_val*r
1658 0 : RETURN
1659 : END IF
1660 0 : e = (1.0e0_dp + t)/(c1 + t + t)
1661 0 : beta = n + w/s + e*(c + n*yp1)
1662 0 : p = 1.0e0_dp + t
1663 0 : s = s + 2.0e0_dp
1664 :
1665 : ! UPDATE AN, BN, ANP1, AND BNP1
1666 :
1667 0 : t = alpha*an + beta*anp1
1668 0 : an = anp1
1669 0 : anp1 = t
1670 0 : t = alpha*bn + beta*bnp1
1671 0 : bn = bnp1
1672 0 : bnp1 = t
1673 0 : r0 = r
1674 0 : r = anp1/bnp1
1675 0 : IF (ABS(r - r0) <= eps*r) THEN
1676 : ! TERMINATION
1677 0 : fn_val = fn_val*r
1678 0 : RETURN
1679 : END IF
1680 :
1681 : ! RESCALE AN, BN, ANP1, AND BNP1
1682 :
1683 0 : an = an/bnp1
1684 0 : bn = bn/bnp1
1685 0 : anp1 = r
1686 0 : bnp1 = 1.0e0_dp
1687 : END DO
1688 :
1689 : END FUNCTION bfrac
1690 :
1691 : ! **************************************************************************************************
1692 : !> \brief ...
1693 : !> \param a ...
1694 : !> \param b ...
1695 : !> \param x ...
1696 : !> \param y ...
1697 : !> \return ...
1698 : ! **************************************************************************************************
1699 0 : FUNCTION brcomp(a, b, x, y) RESULT(fn_val)
1700 : !-----------------------------------------------------------------------
1701 : ! EVALUATION OF X**A*Y**B/BETA(A,B)
1702 : !-----------------------------------------------------------------------
1703 : REAL(dp), INTENT(IN) :: a, b, x, y
1704 : REAL(dp) :: fn_val
1705 :
1706 : REAL(dp), PARAMETER :: const = 0.398942280401433e0_dp
1707 :
1708 : INTEGER :: i, n
1709 : REAL(dp) :: a0, apb, b0, c, e, h, lambda, lnx, lny, &
1710 : t, u, v, x0, y0, z
1711 :
1712 : !-----------------
1713 : ! CONST = 1/SQRT(2*PI)
1714 : !-----------------
1715 :
1716 0 : fn_val = 0.0e0_dp
1717 0 : IF (x == 0.0e0_dp .OR. y == 0.0e0_dp) RETURN
1718 0 : a0 = MIN(a, b)
1719 0 : IF (a0 >= 8.0e0_dp) THEN
1720 : !-----------------------------------------------------------------------
1721 : ! PROCEDURE FOR A .GE. 8 AND B .GE. 8
1722 : !-----------------------------------------------------------------------
1723 0 : IF (a > b) THEN
1724 0 : h = b/a
1725 0 : x0 = 1.0e0_dp/(1.0e0_dp + h)
1726 0 : y0 = h/(1.0e0_dp + h)
1727 0 : lambda = (a + b)*y - b
1728 : ELSE
1729 0 : h = a/b
1730 0 : x0 = h/(1.0e0_dp + h)
1731 0 : y0 = 1.0e0_dp/(1.0e0_dp + h)
1732 0 : lambda = a - (a + b)*x
1733 : END IF
1734 :
1735 0 : e = -lambda/a
1736 0 : IF (ABS(e) > 0.6e0_dp) THEN
1737 0 : u = e - LOG(x/x0)
1738 : ELSE
1739 0 : u = rlog1(e)
1740 : END IF
1741 :
1742 0 : e = lambda/b
1743 0 : IF (ABS(e) > 0.6e0_dp) THEN
1744 0 : v = e - LOG(y/y0)
1745 : ELSE
1746 0 : v = rlog1(e)
1747 : END IF
1748 :
1749 0 : z = EXP(-(a*u + b*v))
1750 0 : fn_val = const*SQRT(b*x0)*z*EXP(-bcorr(a, b))
1751 0 : RETURN
1752 : END IF
1753 :
1754 0 : IF (x > 0.375e0_dp) THEN
1755 0 : IF (y > 0.375e0_dp) THEN
1756 0 : lnx = LOG(x)
1757 0 : lny = LOG(y)
1758 : ELSE
1759 0 : lnx = alnrel(-y)
1760 0 : lny = LOG(y)
1761 : END IF
1762 : ELSE
1763 0 : lnx = LOG(x)
1764 0 : lny = alnrel(-x)
1765 : END IF
1766 :
1767 0 : z = a*lnx + b*lny
1768 0 : IF (a0 >= 1.0e0_dp) THEN
1769 0 : z = z - betaln(a, b)
1770 0 : fn_val = EXP(z)
1771 0 : RETURN
1772 : END IF
1773 : !-----------------------------------------------------------------------
1774 : ! PROCEDURE FOR A .LT. 1 OR B .LT. 1
1775 : !-----------------------------------------------------------------------
1776 0 : b0 = MAX(a, b)
1777 0 : IF (b0 >= 8.0e0_dp) THEN
1778 : ! ALGORITHM FOR B0 .GE. 8
1779 0 : u = gamln1(a0) + algdiv(a0, b0)
1780 0 : fn_val = a0*EXP(z - u)
1781 : END IF
1782 0 : IF (b0 <= 1.0e0_dp) THEN
1783 :
1784 : ! ALGORITHM FOR B0 .LE. 1
1785 :
1786 0 : fn_val = EXP(z)
1787 0 : IF (fn_val == 0.0e0_dp) RETURN
1788 :
1789 0 : apb = a + b
1790 0 : IF (apb > 1.0e0_dp) THEN
1791 0 : u = a + b - 1.e0_dp
1792 0 : z = (1.0e0_dp + gam1(u))/apb
1793 : ELSE
1794 0 : z = 1.0e0_dp + gam1(apb)
1795 : END IF
1796 :
1797 0 : c = (1.0e0_dp + gam1(a))*(1.0e0_dp + gam1(b))/z
1798 0 : fn_val = fn_val*(a0*c)/(1.0e0_dp + a0/b0)
1799 0 : RETURN
1800 : END IF
1801 :
1802 : ! ALGORITHM FOR 1 .LT. B0 .LT. 8
1803 :
1804 0 : u = gamln1(a0)
1805 0 : n = INT(b0 - 1.0e0_dp)
1806 0 : IF (n >= 1) THEN
1807 : c = 1.0e0_dp
1808 0 : DO i = 1, n
1809 0 : b0 = b0 - 1.0e0_dp
1810 0 : c = c*(b0/(a0 + b0))
1811 : END DO
1812 0 : u = LOG(c) + u
1813 : END IF
1814 :
1815 0 : z = z - u
1816 0 : b0 = b0 - 1.0e0_dp
1817 0 : apb = a0 + b0
1818 0 : IF (apb > 1.0e0_dp) THEN
1819 0 : u = a0 + b0 - 1.e0_dp
1820 0 : t = (1.0e0_dp + gam1(u))/apb
1821 : ELSE
1822 0 : t = 1.0e0_dp + gam1(apb)
1823 : END IF
1824 0 : fn_val = a0*EXP(z)*(1.0e0_dp + gam1(b0))/t
1825 :
1826 0 : END FUNCTION brcomp
1827 :
1828 : ! **************************************************************************************************
1829 : !> \brief ...
1830 : !> \param mu ...
1831 : !> \param a ...
1832 : !> \param b ...
1833 : !> \param x ...
1834 : !> \param y ...
1835 : !> \return ...
1836 : ! **************************************************************************************************
1837 0 : FUNCTION brcmp1(mu, a, b, x, y) RESULT(fn_val)
1838 : !-----------------------------------------------------------------------
1839 : ! EVALUATION OF EXP(MU) * (X**A*Y**B/BETA(A,B))
1840 : !-----------------------------------------------------------------------
1841 : INTEGER, INTENT(IN) :: mu
1842 : REAL(dp), INTENT(IN) :: a, b, x, y
1843 : REAL(dp) :: fn_val
1844 :
1845 : REAL(dp), PARAMETER :: const = 0.398942280401433e0_dp
1846 :
1847 : INTEGER :: i, n
1848 : REAL(dp) :: a0, apb, b0, c, e, h, lambda, lnx, lny, &
1849 : t, u, v, x0, y0, z
1850 :
1851 : !-----------------
1852 : ! CONST = 1/SQRT(2*PI)
1853 : !-----------------
1854 :
1855 0 : a0 = MIN(a, b)
1856 0 : IF (a0 >= 8.0e0_dp) THEN
1857 : !-----------------------------------------------------------------------
1858 : ! PROCEDURE FOR A .GE. 8 AND B .GE. 8
1859 : !-----------------------------------------------------------------------
1860 : IF (a > b) THEN
1861 0 : h = b/a
1862 0 : x0 = 1.0e0_dp/(1.0e0_dp + h)
1863 0 : y0 = h/(1.0e0_dp + h)
1864 0 : lambda = (a + b)*y - b
1865 : END IF
1866 0 : h = a/b
1867 0 : x0 = h/(1.0e0_dp + h)
1868 0 : y0 = 1.0e0_dp/(1.0e0_dp + h)
1869 0 : lambda = a - (a + b)*x
1870 :
1871 0 : e = -lambda/a
1872 0 : IF (ABS(e) > 0.6e0_dp) THEN
1873 0 : u = e - LOG(x/x0)
1874 : ELSE
1875 0 : u = rlog1(e)
1876 : END IF
1877 :
1878 0 : e = lambda/b
1879 0 : IF (ABS(e) <= 0.6e0_dp) THEN
1880 0 : v = rlog1(e)
1881 : ELSE
1882 0 : v = e - LOG(y/y0)
1883 : END IF
1884 :
1885 0 : z = esum(mu, -(a*u + b*v))
1886 0 : fn_val = const*SQRT(b*x0)*z*EXP(-bcorr(a, b))
1887 : END IF
1888 :
1889 0 : IF (x > 0.375e0_dp) THEN
1890 0 : IF (y > 0.375e0_dp) THEN
1891 0 : lnx = LOG(x)
1892 0 : lny = LOG(y)
1893 : ELSE
1894 0 : lnx = alnrel(-y)
1895 0 : lny = LOG(y)
1896 : END IF
1897 : ELSE
1898 0 : lnx = LOG(x)
1899 0 : lny = alnrel(-x)
1900 : END IF
1901 0 : z = a*lnx + b*lny
1902 0 : IF (a0 >= 1.0e0_dp) THEN
1903 0 : z = z - betaln(a, b)
1904 0 : fn_val = esum(mu, z)
1905 0 : RETURN
1906 : END IF
1907 : !-----------------------------------------------------------------------
1908 : ! PROCEDURE FOR A .LT. 1 OR B .LT. 1
1909 : !-----------------------------------------------------------------------
1910 0 : b0 = MAX(a, b)
1911 0 : IF (b0 >= 8.0e0_dp) THEN
1912 : ! ALGORITHM FOR B0 .GE. 8
1913 :
1914 0 : u = gamln1(a0) + algdiv(a0, b0)
1915 0 : fn_val = a0*esum(mu, z - u)
1916 0 : RETURN
1917 : END IF
1918 0 : IF (b0 <= 1.0e0_dp) THEN
1919 :
1920 : ! ALGORITHM FOR B0 .LE. 1
1921 :
1922 0 : fn_val = esum(mu, z)
1923 0 : IF (fn_val == 0.0e0_dp) RETURN
1924 :
1925 0 : apb = a + b
1926 0 : IF (apb > 1.0e0_dp) THEN
1927 0 : u = a + b - 1.e0_dp
1928 0 : z = (1.0e0_dp + gam1(u))/apb
1929 : ELSE
1930 0 : z = 1.0e0_dp + gam1(apb)
1931 : END IF
1932 :
1933 0 : c = (1.0e0_dp + gam1(a))*(1.0e0_dp + gam1(b))/z
1934 0 : fn_val = fn_val*(a0*c)/(1.0e0_dp + a0/b0)
1935 0 : RETURN
1936 : END IF
1937 :
1938 : ! ALGORITHM FOR 1 .LT. B0 .LT. 8
1939 :
1940 0 : u = gamln1(a0)
1941 0 : n = INT(b0 - 1.0e0_dp)
1942 0 : IF (n >= 1) THEN
1943 : c = 1.0e0_dp
1944 0 : DO i = 1, n
1945 0 : b0 = b0 - 1.0e0_dp
1946 0 : c = c*(b0/(a0 + b0))
1947 : END DO
1948 0 : u = LOG(c) + u
1949 : END IF
1950 :
1951 0 : z = z - u
1952 0 : b0 = b0 - 1.0e0_dp
1953 0 : apb = a0 + b0
1954 0 : IF (apb > 1.0e0_dp) THEN
1955 0 : u = a0 + b0 - 1.e0_dp
1956 0 : t = (1.0e0_dp + gam1(u))/apb
1957 : ELSE
1958 0 : t = 1.0e0_dp + gam1(apb)
1959 : END IF
1960 0 : fn_val = a0*esum(mu, z)*(1.0e0_dp + gam1(b0))/t
1961 :
1962 0 : END FUNCTION brcmp1
1963 :
1964 : ! **************************************************************************************************
1965 : !> \brief ...
1966 : !> \param a ...
1967 : !> \param b ...
1968 : !> \param lambda ...
1969 : !> \param eps ...
1970 : !> \return ...
1971 : ! **************************************************************************************************
1972 0 : FUNCTION basym(a, b, lambda, eps) RESULT(fn_val)
1973 : !-----------------------------------------------------------------------
1974 : ! ASYMPTOTIC EXPANSION FOR IX(A,B) FOR LARGE A AND B.
1975 : ! LAMBDA = (A + B)*Y - B AND EPS IS THE TOLERANCE USED.
1976 : ! IT IS ASSUMED THAT LAMBDA IS NONNEGATIVE AND THAT
1977 : ! A AND B ARE GREATER THAN OR EQUAL TO 15.
1978 : !-----------------------------------------------------------------------
1979 : REAL(dp), INTENT(IN) :: a, b, lambda, eps
1980 : REAL(dp) :: fn_val
1981 :
1982 : INTEGER, PARAMETER :: num = 20
1983 : REAL(dp), PARAMETER :: e0 = 1.12837916709551e0_dp, &
1984 : e1 = .353553390593274e0_dp
1985 :
1986 : INTEGER :: i, im1, imj, j, m, mm1, mmj, n, np1
1987 : REAL(dp) :: a0(21), b0(21), bsum, c(21), d(21), &
1988 : dsum, f, h, h2, hn, j0, j1, r, r0, r1, &
1989 : s, sum, t, t0, t1, u, w, w0, z, z0, &
1990 : z2, zn, znm1
1991 :
1992 : !------------------------
1993 : ! ****** NUM IS THE MAXIMUM VALUE THAT N CAN TAKE IN THE DO LOOP
1994 : ! ENDING AT STATEMENT 50. IT IS REQUIRED THAT NUM BE EVEN.
1995 : ! THE ARRAYS A0, B0, C, D HAVE DIMENSION NUM + 1.
1996 : !------------------------
1997 : ! E0 = 2/SQRT(PI)
1998 : ! E1 = 2**(-3/2)
1999 : !------------------------
2000 :
2001 0 : fn_val = 0.0e0_dp
2002 0 : IF (a >= b) THEN
2003 0 : h = b/a
2004 0 : r0 = 1.0e0_dp/(1.0e0_dp + h)
2005 0 : r1 = (b - a)/a
2006 0 : w0 = 1.0e0_dp/SQRT(b*(1.0e0_dp + h))
2007 : ELSE
2008 0 : h = a/b
2009 0 : r0 = 1.0e0_dp/(1.0e0_dp + h)
2010 0 : r1 = (b - a)/b
2011 0 : w0 = 1.0e0_dp/SQRT(a*(1.0e0_dp + h))
2012 : END IF
2013 :
2014 0 : f = a*rlog1(-lambda/a) + b*rlog1(lambda/b)
2015 0 : t = EXP(-f)
2016 0 : IF (t == 0.0e0_dp) RETURN
2017 0 : z0 = SQRT(f)
2018 0 : z = 0.5e0_dp*(z0/e1)
2019 0 : z2 = f + f
2020 :
2021 0 : a0(1) = (2.0e0_dp/3.0e0_dp)*r1
2022 0 : c(1) = -0.5e0_dp*a0(1)
2023 0 : d(1) = -c(1)
2024 0 : j0 = (0.5e0_dp/e0)*ERFC_SCALED(z0)
2025 0 : j1 = e1
2026 0 : sum = j0 + d(1)*w0*j1
2027 :
2028 0 : s = 1.0e0_dp
2029 0 : h2 = h*h
2030 0 : hn = 1.0e0_dp
2031 0 : w = w0
2032 0 : znm1 = z
2033 0 : zn = z2
2034 0 : DO n = 2, num, 2
2035 0 : hn = h2*hn
2036 0 : a0(n) = 2.0e0_dp*r0*(1.0e0_dp + h*hn)/(n + 2.0e0_dp)
2037 0 : np1 = n + 1
2038 0 : s = s + hn
2039 0 : a0(np1) = 2.0e0_dp*r1*s/(n + 3.0e0_dp)
2040 :
2041 0 : DO i = n, np1
2042 0 : r = -0.5e0_dp*(i + 1.0e0_dp)
2043 0 : b0(1) = r*a0(1)
2044 0 : DO m = 2, i
2045 : bsum = 0.0e0_dp
2046 : mm1 = m - 1
2047 0 : DO j = 1, mm1
2048 0 : mmj = m - j
2049 0 : bsum = bsum + (j*r - mmj)*a0(j)*b0(mmj)
2050 : END DO
2051 0 : b0(m) = r*a0(m) + bsum/m
2052 : END DO
2053 0 : c(i) = b0(i)/(i + 1.0e0_dp)
2054 :
2055 0 : dsum = 0.0e0_dp
2056 0 : im1 = i - 1
2057 0 : DO j = 1, im1
2058 0 : imj = i - j
2059 0 : dsum = dsum + d(imj)*c(j)
2060 : END DO
2061 0 : d(i) = -(dsum + c(i))
2062 : END DO
2063 :
2064 0 : j0 = e1*znm1 + (n - 1.0e0_dp)*j0
2065 0 : j1 = e1*zn + n*j1
2066 0 : znm1 = z2*znm1
2067 0 : zn = z2*zn
2068 0 : w = w0*w
2069 0 : t0 = d(n)*w*j0
2070 0 : w = w0*w
2071 0 : t1 = d(np1)*w*j1
2072 0 : sum = sum + (t0 + t1)
2073 0 : IF ((ABS(t0) + ABS(t1)) <= eps*sum) EXIT
2074 : END DO
2075 :
2076 0 : u = EXP(-bcorr(a, b))
2077 0 : fn_val = e0*t*u*sum
2078 :
2079 0 : END FUNCTION basym
2080 :
2081 : END MODULE beta_gamma_psi
|
