Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!! !!!!
4 : !!!! ParaMonte: Parallel Monte Carlo and Machine Learning Library. !!!!
5 : !!!! !!!!
6 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab !!!!
7 : !!!! !!!!
8 : !!!! This file is part of the ParaMonte library. !!!!
9 : !!!! !!!!
10 : !!!! LICENSE !!!!
11 : !!!! !!!!
12 : !!!! https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md !!!!
13 : !!!! !!!!
14 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16 :
17 : !> \brief
18 : !> This include file contains implementations of the procedures in module [pm_mathGamma](@ref pm_mathGamma).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \FatemehBagheri, Sunday 11:23 PM, September 19, 2021, Dallas, TX
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : !%%%%%%%%%%%%%%%%%%
28 : #if getGammaInc_ENABLED
29 : !%%%%%%%%%%%%%%%%%%
30 :
31 : integer(IK) :: info
32 : #if Low_ENABLED
33 : character(*, SK), parameter :: PROCEDURE_NAME = SK_"@getGammaIncLow()"
34 3045 : call setGammaIncLow(gammaIncLow, x, log_gamma(kappa), kappa, info, tol)
35 : #elif Upp_ENABLED
36 : character(*, SK), parameter :: PROCEDURE_NAME = SK_"@getGammaIncUpp()"
37 3034 : call setGammaIncUpp(gammaIncUpp, x, log_gamma(kappa), kappa, info, tol)
38 : #else
39 : #error "Unrecognized interface."
40 : #endif
41 : if (info < 0_IK) error stop SK_"@file::"//__FILE__//SK_"@line::"//getStr(__LINE__)//&! LCOV_EXCL_LINE
42 : MODULE_NAME//PROCEDURE_NAME//SK_": The Incomplete Gamma function failed to converge." ! LCOV_EXCL_LINE
43 :
44 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
45 : #elif setGammaInc_ENABLED && Def_ENABLED && Low_ENABLED
46 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
47 :
48 94394 : if (x < kappa + 1._RKC) then
49 74735 : call setGammaIncLowSeries(gammaIncLow, x, logGammaKappa, kappa, info, tol)
50 : else
51 19659 : call setGammaIncUppContFrac(gammaIncLow, x, logGammaKappa, kappa, info, tol)
52 19659 : gammaIncLow = 1._RKC - gammaIncLow
53 : end if
54 :
55 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
56 : #elif setGammaInc_ENABLED && Def_ENABLED && Upp_ENABLED
57 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
58 :
59 4223 : if (x < kappa + 1._RKC) then
60 1653 : call setGammaIncLowSeries(gammaIncUpp, x, logGammaKappa, kappa, info, tol)
61 1653 : gammaIncUpp = 1._RKC - gammaIncUpp
62 : else
63 2570 : call setGammaIncUppContFrac(gammaIncUpp, x, logGammaKappa, kappa, info, tol)
64 : end if
65 :
66 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
67 : #elif (Series_ENABLED && Low_ENABLED) || (ContFrac_ENABLED && Upp_ENABLED)
68 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
69 :
70 : integer(IK) , parameter :: ITMAX = 300_IK !< The maximum allowed number of iterations.
71 : real(RKC) , parameter :: EPS = 10 * epsilon(x) !< The default relative accuracy.
72 : real(RKC) :: tol_def, delta
73 : #if Series_ENABLED && Low_ENABLED
74 : character(*, SK), parameter :: PROCEDURE_NAME = SK_"@setGammaIncLowSeries()"
75 : real(RKC) :: kappaIncremented, summ
76 : #elif ContFrac_ENABLED && Upp_ENABLED
77 : character(*, SK), parameter :: PROCEDURE_NAME = SK_"@setGammaIncUppContFrac()"
78 : real(RKC) , parameter :: FPMIN_DEF = tiny(x) / EPS !< A number near the smallest representable floating-point number.
79 : real(RKC) :: an, b, c, d, h
80 : real(RKC) :: fpmin
81 : #else
82 : #error "Unrecognized interface."
83 : #endif
84 : ! Compute the Gamma function using the fastest method decided at runtime.
85 100623 : CHECK_ASSERTION(__LINE__, 0._RKC <= x, PROCEDURE_NAME//SK_": The input `x` must be positive. x = "//getStr(x)) ! fpp
86 100623 : CHECK_ASSERTION(__LINE__, 0._RKC < getOption(EPS, tol), PROCEDURE_NAME//SK_": The condition `0. < tol` must hold. tol = "//getStr(getOption(EPS, tol))) ! fpp
87 100623 : CHECK_ASSERTION(__LINE__, kappa > 0._RKC, PROCEDURE_NAME//SK_": The input `kappa` must be positive. kappa = "//getStr(kappa)) ! fpp
88 402492 : CHECK_ASSERTION(__LINE__, abs(log_gamma(kappa) - logGammaKappa) <= 100 * epsilon(kappa), \
89 : PROCEDURE_NAME//SK_": The input `logGammaKappa` must equal `log_gamma(kappa)`. kappa, log_gamma(kappa), logGammaKappa = "//\
90 : getStr([kappa, log_gamma(kappa), logGammaKappa])) ! fpp
91 : #if Series_ENABLED && Low_ENABLED
92 77391 : if (x > 0._RKC) then
93 77370 : if (present(tol)) then
94 6 : tol_def = tol
95 : else
96 29484 : tol_def = EPS
97 : end if
98 29490 : kappaIncremented = kappa
99 77370 : summ = 1._RKC / kappa
100 29490 : delta = summ
101 932066 : do info = 1, ITMAX
102 932066 : kappaIncremented = kappaIncremented + 1._RKC
103 932066 : delta = delta * x / kappaIncremented
104 932066 : summ = summ + delta
105 932066 : if (abs(summ) * tol_def < abs(delta)) cycle
106 77370 : gammaIncLow = summ * exp(kappa * log(x) - x - logGammaKappa)
107 932066 : return
108 : end do
109 : info = -info ! LCOV_EXCL_LINE
110 : else
111 21 : info = 0_IK
112 21 : gammaIncLow = 0._RKC
113 : end if
114 : #elif ContFrac_ENABLED && Upp_ENABLED
115 23232 : if (x > 0._RKC) then
116 23231 : if (present(tol)) then
117 0 : tol_def = tol
118 0 : fpmin = tiny(x) / tol_def
119 : else
120 2913 : tol_def = EPS
121 2913 : fpmin = FPMIN_DEF
122 : end if
123 23231 : b = x + 1._RKC - kappa
124 23231 : c = 1._RKC / fpmin
125 23231 : d = 1._RKC / b
126 2913 : h = d
127 667361 : do info = 1, ITMAX
128 667358 : an = -info * (info - kappa)
129 667358 : b = b + 2._RKC
130 667358 : d = an * d + b
131 667358 : if (abs(d) < fpmin) d = fpmin
132 667358 : c = b + an / c
133 667358 : if (abs(c) < fpmin) c = fpmin
134 667358 : d = 1._RKC / d
135 667358 : delta = d * c
136 667358 : h = h * delta
137 667358 : if (tol_def < abs(delta - 1._RKC)) cycle
138 23228 : gammaIncUpp = exp(kappa * log(x) - x - logGammaKappa) * h
139 667361 : return
140 : end do
141 : info = -info ! LCOV_EXCL_LINE
142 : else
143 1 : info = 0_IK
144 1 : gammaIncUpp = 1._RKC
145 : end if
146 : #endif
147 : #else
148 : !%%%%%%%%%%%%%%%%%%%%%%%%
149 : #error "Unrecognized interface."
150 : !%%%%%%%%%%%%%%%%%%%%%%%%
151 : #endif
|
