Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!!
4 : !!!! MIT License
5 : !!!!
6 : !!!! ParaMonte: plain powerful parallel Monte Carlo library.
7 : !!!!
8 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab
9 : !!!!
10 : !!!! This file is part of the ParaMonte library.
11 : !!!!
12 : !!!! Permission is hereby granted, free of charge, to any person obtaining a
13 : !!!! copy of this software and associated documentation files (the "Software"),
14 : !!!! to deal in the Software without restriction, including without limitation
15 : !!!! the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 : !!!! and/or sell copies of the Software, and to permit persons to whom the
17 : !!!! Software is furnished to do so, subject to the following conditions:
18 : !!!!
19 : !!!! The above copyright notice and this permission notice shall be
20 : !!!! included in all copies or substantial portions of the Software.
21 : !!!!
22 : !!!! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
23 : !!!! EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
24 : !!!! MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
25 : !!!! IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
26 : !!!! DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
27 : !!!! OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
28 : !!!! OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
29 : !!!!
30 : !!!! ACKNOWLEDGMENT
31 : !!!!
32 : !!!! ParaMonte is an honor-ware and its currency is acknowledgment and citations.
33 : !!!! As per the ParaMonte library license agreement terms, if you use any parts of
34 : !!!! this library for any purposes, kindly acknowledge the use of ParaMonte in your
35 : !!!! work (education/research/industry/development/...) by citing the ParaMonte
36 : !!!! library as described on this page:
37 : !!!!
38 : !!!! https://github.com/cdslaborg/paramonte/blob/main/ACKNOWLEDGMENT.md
39 : !!!!
40 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
41 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
42 :
43 : !> \brief This module contains mathematical procedures.
44 : !> \author Amir Shahmoradi
45 :
46 : module Math_mod
47 :
48 : implicit none
49 :
50 : character(*), parameter :: MODULE_NAME = "@Math_mod"
51 :
52 : interface getCumSum
53 : module procedure :: getCumSum_IK, getCumSum_RK
54 : end interface getCumSum
55 :
56 : interface getCumSumReverse
57 : module procedure :: getCumSumReverse_IK, getCumSumReverse_RK
58 : end interface getCumSumReverse
59 :
60 : interface getLogSumExp
61 : module procedure :: getLogSumExp_RK, getLogSumExp_CK
62 : end interface getLogSumExp
63 :
64 : interface getLogSubExp
65 : module procedure :: getLogSubExp_RK
66 : end interface getLogSubExp
67 :
68 : interface getLogEggBox
69 : module procedure :: getLogEggBoxSD_RK, getLogEggBoxMD_RK
70 : module procedure :: getLogEggBoxSD_CK, getLogEggBoxMD_CK
71 : end interface getLogEggBox
72 :
73 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
74 :
75 : contains
76 :
77 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
78 :
79 : !> \brief
80 : !> Return the distance squared between the two input points.
81 : !>
82 : !> @param[in] nd : The size of the input vectors.
83 : !> @param[in] Point1 : The first point.
84 : !> @param[in] Point2 : The second point.
85 : !>
86 : !> \return
87 : !> `distanceSq` : The distance squared.
88 3 : pure function getDistanceSq(nd,Point1,Point2) result(distanceSq)
89 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
90 : !DEC$ ATTRIBUTES DLLEXPORT :: getDistanceSq
91 : #endif
92 : use Constants_mod, only: IK, RK
93 : integer(IK) , intent(in) :: nd
94 : real(RK) , intent(in) :: Point1(nd),Point2(nd)
95 : real(RK) :: distanceSq
96 : integer(IK) :: i
97 3 : distanceSq = 0._RK
98 18 : do i = 1, nd
99 18 : distanceSq = distanceSq + (Point2(i)-Point1(i))**2
100 : end do
101 3 : end function getDistanceSq
102 :
103 :
104 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
105 :
106 : !> \brief
107 : !> Return the correlation coefficient (`-1 < corCoef < 1`) corresponding to the input Fisher z-transformation.
108 : !>
109 : !> @param[in] fisherTrans : The Fisher transformation.
110 : !>
111 : !> \return
112 : !> `corCoef` : The correlation coefficient corresponding to the input Fisher transformation.
113 3 : pure elemental function getCorCeofFromFisherTrans(fisherTrans) result(corCoef)
114 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
115 : !DEC$ ATTRIBUTES DLLEXPORT :: getCorCeofFromFisherTrans
116 : #endif
117 3 : use Constants_mod, only: RK
118 : real(RK), intent(in) :: fisherTrans
119 : real(RK) :: corCoef !, twiceFisherTrans
120 3 : corCoef = tanh(fisherTrans)
121 : !twiceFisherTrans = 2._RK * corCoef
122 : !fisherTrans = (exp(twiceFisherTrans)-1._RK) / (exp(twiceFisherTrans)+1._RK)
123 6 : end function getCorCeofFromFisherTrans
124 :
125 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
126 :
127 : !> \brief
128 : !> Return Fisher z-transformation of an input correlation coefficient (`-1 < corCoef < 1`).
129 : !>
130 : !> @param[in] corCoef : The input correlation coefficient.
131 : !>
132 : !> \return
133 : !> `fisherTrans` : The output Fisher transformation.
134 3 : pure elemental function getFisherTransFromCorCoef(corCoef) result(fisherTrans)
135 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
136 : !DEC$ ATTRIBUTES DLLEXPORT :: getFisherTransFromCorCoef
137 : #endif
138 3 : use Constants_mod, only: RK
139 : real(RK), intent(in) :: corCoef
140 : real(RK) :: fisherTrans
141 3 : fisherTrans = atanh(corCoef)
142 : !fisherTrans = 0.5_RK * log(1._RK+corCoef) / (1._RK-corCoef)
143 6 : end function getFisherTransFromCorCoef
144 :
145 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
146 :
147 : !> \brief
148 : !> Return the cumulative sum of the input integer array.
149 : !>
150 : !> @param[in] vecLen : The length of the input vector.
151 : !> @param[in] Vec : The input vector.
152 : !>
153 : !> \return
154 : !> `CumSum` : An integer array of length `vecLen` representing the cumulative sum.
155 : !>
156 : !> \remark
157 : !> The first element of `CumSum` is the same as the first element of `Vec`.
158 2916 : pure function getCumSum_IK(vecLen,Vec) result(CumSum)
159 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
160 : !DEC$ ATTRIBUTES DLLEXPORT :: getCumSum_IK
161 : #endif
162 3 : use Constants_mod, only: IK
163 : integer(IK), intent(in) :: vecLen
164 : integer(IK), intent(in) :: Vec(vecLen)
165 : integer(IK) :: CumSum(vecLen)
166 : integer(IK) :: i
167 1458 : CumSum(1) = Vec(1)
168 351944 : do i = 2, vecLen
169 351944 : CumSum(i) = CumSum(i-1) + Vec(i)
170 : end do
171 1458 : end function getCumSum_IK
172 :
173 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
174 :
175 : !> \brief
176 : !> Return the cumulative sum of the input real array.
177 : !> @param[in] vecLen : The length of the input vector.
178 : !> @param[in] Vec : The input vector.
179 : !>
180 : !> \return
181 : !> `CumSum` : A real array of length `vecLen` representing the cumulative sum.
182 : !>
183 : !> \remark
184 : !> The first element of `CumSum` is the same as the first element of `Vec`.
185 1175220 : pure function getCumSum_RK(vecLen,Vec) result(CumSum)
186 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
187 : !DEC$ ATTRIBUTES DLLEXPORT :: getCumSum_RK
188 : #endif
189 1458 : use Constants_mod, only: IK, RK
190 : integer(IK), intent(in) :: vecLen
191 : real(RK) , intent(in) :: Vec(vecLen)
192 : real(RK) :: CumSum(vecLen)
193 : integer(IK) :: i
194 306 : CumSum(1) = Vec(1)
195 1174300 : do i = 2, vecLen
196 1174300 : CumSum(i) = CumSum(i-1) + Vec(i)
197 : end do
198 306 : end function getCumSum_RK
199 :
200 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
201 :
202 : !> \brief
203 : !> Return the cumulative sum of the input integer array, in the backward direction.
204 : !> @param[in] vecLen : The length of the input vector.
205 : !> @param[in] Vec : The input vector.
206 : !>
207 : !> \return
208 : !> `CumSumReverse` : An integer array of length `vecLen` representing the cumulative sum.
209 : !>
210 : !> \remark
211 : !> The last element of `CumSumReverse` is the same as the last element of `Vec`.
212 6 : pure function getCumSumReverse_IK(vecLen,Vec) result(CumSumReverse)
213 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
214 : !DEC$ ATTRIBUTES DLLEXPORT :: getCumSumReverse_IK
215 : #endif
216 306 : use Constants_mod, only: IK
217 : integer(IK), intent(in) :: vecLen
218 : integer(IK), intent(in) :: Vec(vecLen)
219 : integer(IK) :: CumSumReverse(vecLen)
220 : integer(IK) :: i, indx
221 3 : CumSumReverse(1) = Vec(vecLen)
222 30 : do i = vecLen-1,1,-1
223 27 : indx = vecLen - i
224 30 : CumSumReverse(indx+1) = CumSumReverse(indx) + Vec(i)
225 : end do
226 3 : end function getCumSumReverse_IK
227 :
228 : !> \brief
229 : !> Return the cumulative sum of the input real array, in the backward direction.
230 : !>
231 : !> @param[in] vecLen : The length of the input vector.
232 : !> @param[in] Vec : The input vector.
233 : !>
234 : !> \return
235 : !> `CumSumReverse` : A real array of length `vecLen` representing the cumulative sum.
236 : !>
237 : !> \remark
238 : !> The last element of `CumSumReverse` is the same as the last element of `Vec`.
239 39 : pure function getCumSumReverse_RK(vecLen,Vec) result(CumSumReverse)
240 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
241 : !DEC$ ATTRIBUTES DLLEXPORT :: getCumSumReverse_RK
242 : #endif
243 3 : use Constants_mod, only: IK, RK
244 : integer(IK), intent(in) :: vecLen
245 : real(RK) , intent(in) :: Vec(vecLen)
246 : real(RK) :: CumSumReverse(vecLen)
247 : integer(IK) :: i, indx
248 3 : CumSumReverse(1) = Vec(vecLen)
249 30 : do i = vecLen-1,1,-1
250 27 : indx = vecLen - i
251 30 : CumSumReverse(indx+1) = CumSumReverse(indx) + Vec(i)
252 : end do
253 3 : end function getCumSumReverse_RK
254 :
255 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
256 :
257 : !> \brief
258 : !> Return `log( exp(logValueLarger) - exp(logValueSamller) )` robustly (without overflow or underflow).
259 : !>
260 : !> @param[in] logValueLarger : A real number.
261 : !> @param[in] logValueSamller : A real number.
262 : !>
263 : !> \return
264 : !> `logSubExp` : A real number.
265 : !>
266 : !> \warning
267 : !> The onus is on the user to ensure `logValueLarger > logValueSamller`.
268 : !>
269 : !> \remark
270 : !> This function is very useful for situations where `exp(logValueLarger)` is likely to cause overflow.
271 8261 : pure function getLogSubExp_RK(logValueLarger,logValueSamller) result(logSubExp)
272 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
273 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogSubExp_RK
274 : #endif
275 3 : use Constants_mod, only: RK, LOGTINY_RK
276 : real(RK) , intent(in) :: logValueLarger, logValueSamller
277 8261 : real(RK) :: logSubExp, diff
278 8261 : logSubExp = logValueLarger
279 8261 : diff = logValueSamller - logValueLarger
280 8260 : if (diff>LOGTINY_RK) logSubExp = logSubExp + log( 1._RK - exp(diff) )
281 8261 : end function getLogSubExp_RK
282 :
283 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
284 :
285 : !> \brief
286 : !> Return the logarithm of the sum of the exponential of the input real vector robustly (without overflow or underflow).
287 : !>
288 : !> @param[in] lenLogValue : The length of the input vector.
289 : !> @param[in] LogValue : The input vector of log-values whose log-sum-exp must be computed.
290 : !> @param[in] maxLogValue : The maximum of the input `LogValue` argument (**optional**).
291 : !>
292 : !> \return
293 : !> `logSumExp` : A real number.
294 6 : pure function getLogSumExp_RK(lenLogValue, LogValue, maxLogValue) result(logSumExp)
295 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
296 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogSumExp_RK
297 : #endif
298 8261 : use Constants_mod, only: IK, RK, LOGTINY_RK
299 : integer(IK) , intent(in) :: lenLogValue
300 : real(RK) , intent(in) :: LogValue(lenLogValue)
301 : real(RK) , intent(in), optional :: maxLogValue
302 : real(RK) :: logSumExp
303 6 : real(RK) :: maxLogValueDefault
304 24 : real(RK) :: LogValueCopy(lenLogValue)
305 : integer(IK) :: i
306 6 : if (present(maxLogValue)) then
307 3 : maxLogValueDefault = maxLogValue
308 : else
309 15 : maxLogValueDefault = maxval(LogValue)
310 : end if
311 24 : LogValueCopy = LogValue - maxLogValueDefault
312 6 : do concurrent(i=1:lenLogValue)
313 24 : if (LogValueCopy(i)<LOGTINY_RK) then
314 0 : LogValueCopy(i) = 0._RK
315 : else
316 18 : LogValueCopy(i) = exp(LogValueCopy(i))
317 : end if
318 : end do
319 24 : logSumExp = maxLogValueDefault + log(sum(LogValueCopy))
320 6 : end function getLogSumExp_RK
321 :
322 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
323 :
324 : !> \brief
325 : !> Return the logarithm of the sum of the exponential of the input complex vector robustly (without overflow or underflow).
326 : !>
327 : !> @param[in] lenLogValue : The length of the input vector.
328 : !> @param[in] LogValue : The input vector of log-values whose log-sum-exp must be computed.
329 : !> @param[in] maxLogValue : The maximum of the real component of the input `LogValue` argument (**optional**).
330 : !>
331 : !> \return
332 : !> `logSumExp` : A complex number.
333 6 : pure function getLogSumExp_CK(lenLogValue, LogValue, maxLogValue) result(logSumExp)
334 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
335 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogSumExp_CK
336 : #endif
337 6 : use Constants_mod, only: IK, RK, LOGTINY_RK
338 : integer(IK) , intent(in) :: lenLogValue
339 : complex(RK) , intent(in) :: LogValue(lenLogValue)
340 : complex(RK) , intent(in), optional :: maxLogValue
341 : complex(RK) :: logSumExp
342 24 : complex(RK) :: LogValueCopy(lenLogValue)
343 : complex(RK) :: maxLogValueDefault
344 : integer(IK) :: i
345 6 : if (present(maxLogValue)) then
346 3 : maxLogValueDefault = maxLogValue
347 : else
348 15 : maxLogValueDefault = maxval(real(LogValue))
349 : end if
350 24 : LogValueCopy = LogValue - maxLogValueDefault
351 6 : do concurrent(i=1:lenLogValue)
352 24 : if (real(LogValueCopy(i))<LOGTINY_RK) then
353 0 : LogValueCopy(i) = 0._RK
354 : else
355 18 : LogValueCopy(i) = exp(LogValueCopy(i))
356 : end if
357 : end do
358 24 : logSumExp = maxLogValueDefault + log(sum(LogValueCopy))
359 6 : end function getLogSumExp_CK
360 :
361 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
362 :
363 : !> \brief
364 : !> Return the logarithm of the egg-box probability density function in one dimension.
365 3 : pure function getLogEggBoxSD_RK(constant,exponent,coef,point) result(logEggBox)
366 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
367 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogEggBoxSD_RK
368 : #endif
369 6 : use Constants_mod, only: IK, RK
370 : implicit none
371 : real(RK), intent(in) :: constant
372 : real(RK), intent(in) :: exponent
373 : real(RK), intent(in) :: coef
374 : real(RK), intent(in) :: point
375 : real(RK) :: logEggBox
376 3 : logEggBox = exponent * ( constant + cos(coef*point) )
377 6 : end function getLogEggBoxSD_RK
378 :
379 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
380 :
381 : !> \brief
382 : !> Return the logarithm of the egg-box probability density function in one dimension, as a complex number.
383 3 : pure function getLogEggBoxSD_CK(constant,exponent,coef,point) result(logEggBox)
384 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
385 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogEggBoxSD_CK
386 : #endif
387 3 : use Constants_mod, only: IK, RK
388 : implicit none
389 : complex(RK), intent(in) :: constant
390 : complex(RK), intent(in) :: exponent
391 : complex(RK), intent(in) :: coef
392 : complex(RK), intent(in) :: Point
393 : complex(RK) :: logEggBox
394 3 : logEggBox = exponent * ( constant + cos(coef*point) )
395 6 : end function getLogEggBoxSD_CK
396 :
397 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
398 :
399 : !> \brief
400 : !> Return the logarithm of the egg-box probability density function in multiple dimensions, as a real number.
401 3 : pure function getLogEggBoxMD_RK(nd,constant,exponent,Coef,Point) result(logEggBox)
402 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
403 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogEggBoxMD_RK
404 : #endif
405 3 : use Constants_mod, only: IK, RK
406 : implicit none
407 : integer(IK), intent(in) :: nd
408 : real(RK), intent(in) :: constant
409 : real(RK), intent(in) :: exponent
410 : real(RK), intent(in) :: Coef(nd)
411 : real(RK), intent(in) :: Point(nd)
412 : real(RK) :: logEggBox
413 : integer(IK) :: i
414 3 : logEggBox = 1._RK
415 12 : do i = 1, nd
416 12 : logEggBox = logEggBox * cos(Coef(i)*Point(i))
417 : end do
418 3 : logEggBox = exponent * ( constant + logEggBox )
419 3 : end function getLogEggBoxMD_RK
420 :
421 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
422 :
423 : !> \brief
424 : !> Return the logarithm of the egg-box probability density function in multiple dimensions, as a complex number.
425 : !> \remark
426 : !> The multidimensional EggBox function follows an equation of the form,
427 : !> \f{equation}
428 : !> f({\mathbf{\theta}}) = \bigg[ \exp( \text{constant} + \prod_{i=1)^{\text{nd}} ~ \cos(\text{Coef}_i^2)
429 : !> \f}
430 3 : pure function getLogEggBoxMD_CK(nd,constant,exponent,Coef,Point) result(logEggBox)
431 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
432 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogEggBoxMD_CK
433 : #endif
434 3 : use Constants_mod, only: IK, RK
435 : implicit none
436 : integer(IK), intent(in) :: nd
437 : complex(RK), intent(in) :: constant
438 : complex(RK), intent(in) :: exponent
439 : complex(RK), intent(in) :: Coef(nd)
440 : complex(RK), intent(in) :: Point(nd)
441 : complex(RK) :: logEggBox
442 : integer(IK) :: i
443 3 : logEggBox = 1._RK
444 12 : do i = 1, nd
445 12 : logEggBox = logEggBox * cos(Coef(i)*Point(i))
446 : end do
447 3 : logEggBox = exponent * ( constant + logEggBox )
448 3 : end function getLogEggBoxMD_CK
449 :
450 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
451 :
452 : !> \brief
453 : !> Return the log(factorial) for a whole integer input. This is basically `logGamma( positiveInteger + 1 )`.
454 : !>
455 : !> \param[in] positiveInteger : The input positive integer whose log factorial must be computed.
456 : !>
457 : !> \return
458 : !> `logFactorial` : The natural logarithm of the factorial.
459 : !>
460 : !> \remark
461 : !> This function is mostly useful for large input integers.
462 123 : pure function getLogFactorial(positiveInteger) result(logFactorial)
463 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
464 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogFactorial
465 : #endif
466 : use Constants_mod, only: IK, RK ! LCOV_EXCL_LINE
467 : implicit none
468 : integer(IK), intent(in) :: positiveInteger
469 : integer(IK) :: i
470 : real(RK) :: logFactorial
471 123 : logFactorial = 0._RK
472 174 : do i = 2, positiveInteger
473 174 : logFactorial = logFactorial + log(real(i,kind=RK))
474 : end do
475 123 : end function getLogFactorial
476 :
477 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
478 :
479 : !> \brief
480 : !> Return the factorial for a whole integer input. This is basically `Gamma( intNum + 1 )`.
481 : !>
482 : !> \param[in] positiveInteger : The input positive integer whose log factorial must be computed.
483 : !>
484 : !> \return
485 : !> `factorial` : The factorial.
486 3 : pure function getFactorial(positiveInteger) result(factorial)
487 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
488 : !DEC$ ATTRIBUTES DLLEXPORT :: getFactorial
489 : #endif
490 123 : use Constants_mod, only: IK, RK
491 : implicit none
492 : integer(IK), intent(in) :: positiveInteger
493 : integer(IK) :: i
494 : real(RK) :: factorial
495 3 : factorial = 1._RK
496 30 : do i = 2, positiveInteger
497 30 : factorial = factorial * i
498 : end do
499 3 : end function getFactorial
500 :
501 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
502 :
503 : !> \brief
504 : !> Return the Gamma function for a half-integer input as real of kind `RK`.
505 : !>
506 : !> \param[in] positiveHalfInteger : The input half integer as a real number
507 : !>
508 : !> \remark
509 : !> The equation for half-integer Gamma-function is given as,
510 : !> \f{equation}
511 : !> \Gamma \left( \frac{n}{2} \right) = \sqrt \pi \frac{ (n-2)!! }{ 2^\frac{n-1}{2} } ~,
512 : !> \f}
513 : !>
514 : !> \return
515 : !> `gammaHalfInt` : The Gamma function for a half integer input.
516 3 : pure function getGammaHalfInt(positiveHalfInteger) result(gammaHalfInt)
517 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
518 : !DEC$ ATTRIBUTES DLLEXPORT :: getGammaHalfInt
519 : #endif
520 3 : use Constants_mod, only: IK, RK, SQRTPI
521 : implicit none
522 : real(RK), intent(in) :: positiveHalfInteger
523 : real(RK) :: gammaHalfInt
524 : integer(IK) :: i,k
525 3 : gammaHalfInt = SQRTPI
526 3 : k = nint(positiveHalfInteger-0.5_RK,kind=IK) ! positiveHalfInteger = k + 1/2
527 3 : do i = k+1, 2*k
528 3 : gammaHalfInt = gammaHalfInt * 0.25_RK * i
529 : end do
530 3 : end function getGammaHalfInt
531 :
532 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
533 :
534 : !> \brief
535 : !> Return the natural logarithm of the Gamma function for a half-integer input as real of kind `RK`.
536 : !>
537 : !> \param[in] positiveHalfInteger : The input half integer as a real number
538 : !>
539 : !> \remark
540 : !> The equation for half-integer Gamma-function is given as,
541 : !> \f{equation}
542 : !> \Gamma \left( \frac{n}{2} \right) = \sqrt \pi \frac{ (n-2)!! }{ 2^\frac{n-1}{2} } ~,
543 : !> \f}
544 : !>
545 : !> \return
546 : !> `gammaHalfInt` : The Gamma function for a half integer input.
547 3 : pure function getLogGammaHalfInt(positiveHalfInteger) result(logGammaHalfInt)
548 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
549 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogGammaHalfInt
550 : #endif
551 3 : use Constants_mod, only: IK, RK, SQRTPI
552 : implicit none
553 : real(RK), intent(in) :: positiveHalfInteger
554 : real(RK), parameter :: COEF = log(0.25_RK)
555 : real(RK), parameter :: LOG_SQRTPI = log(SQRTPI)
556 : real(RK) :: logGammaHalfInt
557 : integer(IK) :: i, k
558 3 : k = nint(positiveHalfInteger-0.5_RK,kind=IK) ! positiveHalfInteger = k + 1/2
559 3 : logGammaHalfInt = LOG_SQRTPI
560 3 : do i = k+1, 2*k
561 3 : logGammaHalfInt = logGammaHalfInt + COEF + log(real(i,kind=RK))
562 : end do
563 3 : end function getLogGammaHalfInt
564 :
565 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
566 :
567 : ! see for example: http://math.stackexchange.com/questions/606184/volume-of-n-dimensional-ellipsoid
568 : ! see for example: https://en.wikipedia.org/wiki/Volume_of_an_n-ball
569 : ! see for example: https://en.wikipedia.org/wiki/Particular_values_of_the_Gamma_function
570 : ! getEllVolCoef = PI^(nd/2) / gamma(nd/2+1) where n is just a positive integer
571 : !
572 : !> \brief
573 : !> Return the coefficient of the volume of an `nd`-dimensional ellipsoid.
574 : !>
575 : !> param[in] nd : The number of dimensions.
576 : !>
577 : !> \return
578 : !> `ellVolCoef` : The coefficient of the volume of an `nd`-dimensional ellipsoid.
579 6 : pure function getEllVolCoef(nd) result(ellVolCoef)
580 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
581 : !DEC$ ATTRIBUTES DLLEXPORT :: getEllVolCoef
582 : #endif
583 3 : use Constants_mod, only: IK, RK, PI
584 : implicit none
585 : integer(IK), intent(in) :: nd
586 : integer(IK) :: i,k
587 : real(RK), parameter :: FOUR_PI = PI * 4._RK
588 : real(RK) :: ellVolCoef
589 6 : if (mod(nd,2_IK)==0_IK) then ! nd is even
590 :
591 3 : ellVolCoef = PI
592 15 : do i = 2_IK, nd / 2_IK
593 15 : ellVolCoef = ellVolCoef * PI / i ! nd = 2k ; ellVolCoef = PI^(k) / Factorial(k)
594 : end do
595 :
596 : else ! nd is an odd integer
597 :
598 : ! nd = 2k-1 ; gamma(nd/2 + 1) = gamma(k + 1/2) ; gamma(k+1/2) = sqrt(PI) * (2k)! / (4^k * k!)
599 3 : k = (nd + 1_IK) / 2_IK
600 :
601 : ! This is to avoid an extra unnecessary division of ellVolCoef by PI
602 3 : ellVolCoef = 4._RK / (k + 1_IK)
603 :
604 18 : do i = k+2_IK, 2_IK*k
605 : ! ellVolCoef = PI^(k-1/2) / gamma(k+1/2) = PI^(k+1) * 4^k * k! / (2k)!
606 18 : ellVolCoef = ellVolCoef * FOUR_PI / i
607 : end do
608 :
609 : end if
610 6 : end function getEllVolCoef
611 :
612 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
613 :
614 : ! see for example: http://math.stackexchange.com/questions/606184/volume-of-n-dimensional-ellipsoid
615 : ! see for example: https://en.wikipedia.org/wiki/Volume_of_an_n-ball
616 : ! see for example: https://en.wikipedia.org/wiki/Particular_values_of_the_Gamma_function
617 : ! getEllVolCoef = PI^(nd/2) / gamma(nd/2+1) where n is just a positive integer
618 : !
619 : !> \brief
620 : !> Return the coefficient of the volume of an `nd`-dimensional ellipsoid.
621 : !>
622 : !> param[in] nd : The number of dimensions.
623 : !>
624 : !> \return
625 : !> `logEllVolCoef` : The natural logarithm of the volume of an `nd`-dimensional unit-ball.
626 6 : pure function getLogEllVolCoef(nd) result(logEllVolCoef)
627 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
628 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogEllVolCoef
629 : #endif
630 6 : use Constants_mod, only: IK, RK, PI
631 : implicit none
632 : integer(IK), intent(in) :: nd
633 : integer(IK) :: i,k
634 : real(RK), parameter :: FOUR_PI = PI * 4._RK
635 : real(RK) :: logEllVolCoef
636 6 : if (mod(nd,2_IK)==0_IK) then ! nd is even
637 :
638 3 : logEllVolCoef = PI
639 15 : do i = 2_IK, nd / 2_IK
640 15 : logEllVolCoef = logEllVolCoef * PI / i ! nd = 2k ; ellVolCoef = PI^(k) / Factorial(k)
641 : end do
642 :
643 : else ! nd is an odd integer
644 :
645 : ! nd = 2k-1 ; gamma(nd/2 + 1) = gamma(k + 1/2) ; gamma(k+1/2) = sqrt(PI) * (2k)! / (4^k * k!)
646 3 : k = (nd + 1_IK) / 2_IK
647 :
648 : ! This is to avoid an extra unnecessary division of ellVolCoef by PI
649 3 : logEllVolCoef = 4._RK / (k + 1_IK)
650 :
651 18 : do i = k+2_IK, 2_IK*k
652 : ! ellVolCoef = PI^(k-1/2) / gamma(k+1/2) = PI^(k+1) * 4^k * k! / (2k)!
653 18 : logEllVolCoef = logEllVolCoef * FOUR_PI / i
654 : end do
655 :
656 : end if
657 6 : logEllVolCoef = log(logEllVolCoef)
658 6 : end function getLogEllVolCoef
659 :
660 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
661 :
662 : !> \brief
663 : !> Return the logarithm of the volume of an `nd`-dimensional ball of unit-radius.
664 : !>
665 : !> param[in] nd : The number of dimensions.
666 : !>
667 : !> \return
668 : !> `logVolUnitBall` : The logarithm of the volume of an `nd`-dimensional ball of unit-radius.
669 : !>
670 : !> \brief
671 : !> This routine is an exact replacement for [getLogEllVolCoef](@ref getlogellvolcoef).
672 : !> However, it is, on average, 10%-15% faster than [getLogEllVolCoef](@ref getlogellvolcoef) based
673 : !> on benchmarks with Intel ifort 19.4 with all optimization flags on.
674 : !
675 : ! see for example: http://math.stackexchange.com/questions/606184/volume-of-n-dimensional-ellipsoid
676 : ! see for example: https://en.wikipedia.org/wiki/Volume_of_an_n-ball
677 : ! see for example: https://en.wikipedia.org/wiki/Particular_values_of_the_Gamma_function
678 : ! getEllVolCoef = PI^(nd/2) / gamma(nd/2+1) where n is just a positive integer
679 123 : pure function getLogVolUnitBall(nd) result(logVolUnitBall)
680 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
681 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogVolUnitBall
682 : #endif
683 6 : use Constants_mod, only: IK, RK, PI
684 : implicit none
685 : integer(IK) , intent(in) :: nd
686 : real(RK) , parameter :: LOG_PI = log(PI)
687 : integer(IK) :: ndHalfInteger
688 : real(RK) :: logVolUnitBall
689 123 : real(RK) :: ndHalfReal
690 123 : if (mod(nd,2_IK)==0_IK) then ! nd is even
691 120 : ndHalfInteger = nd / 2_IK
692 120 : logVolUnitBall = ndHalfInteger*LOG_PI - getLogFactorial(ndHalfInteger) ! nd = 2k ; logVolUnitBall = PI^(k) / Factorial(k)
693 : else ! nd is an odd integer
694 3 : ndHalfReal = 0.5_RK * real(nd,kind=RK)
695 3 : logVolUnitBall = ndHalfReal * LOG_PI - log_gamma(ndHalfReal+1._RK)
696 : end if
697 123 : end function getLogVolUnitBall
698 :
699 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
700 :
701 : !> \brief
702 : !> Return the logarithm of the volume of an `nd`-dimensional hyper-ellipsoid.
703 : !>
704 : !> param[in] nd : The number of dimensions.
705 : !> param[in] logSqrtDetCovMat : The logarithm of the square root of the determinant of the representative covariance matrix of the ellipsoid.
706 : !>
707 : !> \return
708 : !> `logVolEllipsoid` : The logarithm of the volume of an `nd`-dimensional hyper-ellipsoid.
709 : !
710 : ! see for example: https://math.stackexchange.com/questions/332391/volume-of-hyperellipsoid/332434
711 : ! see for example: https://math.stackexchange.com/questions/2854930/volume-of-an-n-dimensional-ellipsoid
712 6 : pure function getLogVolEllipsoid(nd,logSqrtDetCovMat) result(logVolEllipsoid)
713 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
714 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogVolEllipsoid
715 : #endif
716 123 : use Constants_mod, only: IK, RK
717 : implicit none
718 : integer(IK) , intent(in) :: nd
719 : real(RK) , intent(in) :: logSqrtDetCovMat
720 : real(RK) :: logVolEllipsoid
721 6 : logVolEllipsoid = getLogVolUnitBall(nd) + logSqrtDetCovMat
722 12 : end function getLogVolEllipsoid
723 :
724 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
725 :
726 : !> \brief
727 : !> Return the volume of a list of `nEllipsoid` `nd`-dimensional hyper-ellipsoids.
728 : !>
729 : !> @param[in] nd : The number of dimensions of the ellipsoids.
730 : !> @param[in] nEllipsoid : The number of ellipsoids.
731 : !> @param[in] LogSqrtDetCovMat : The vector of natural logarithms of the square roots of the determinants of the covariance matrices of the ellipsoids.
732 : !>
733 : !> \return
734 : !> `logVolEllipsoid` : The logarithm of the volume of an `nd`-dimensional hyper-ellipsoid.
735 : !>
736 : !> \remark
737 : !> There is really no reason to use this function for HPC solutions, as it can be likely done more efficiently.
738 15 : pure function getLogVolEllipsoids(nd,nEllipsoid,LogSqrtDetCovMat) result(LogVolEllipsoids)
739 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
740 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogVolEllipsoids
741 : #endif
742 6 : use Constants_mod, only: IK, RK
743 : implicit none
744 : integer(IK) :: i
745 : integer(IK) , intent(in) :: nd
746 : integer(IK) , intent(in) :: nEllipsoid
747 : real(RK) , intent(in) :: LogSqrtDetCovMat(nEllipsoid)
748 : real(RK) :: LogVolEllipsoids(nEllipsoid)
749 3 : real(RK) :: logVolUnitBall
750 3 : logVolUnitBall = getLogVolUnitBall(nd)
751 9 : do i = 1, nEllipsoid
752 9 : LogVolEllipsoids(i) = logVolUnitBall + LogSqrtDetCovMat(i)
753 : end do
754 3 : end function getLogVolEllipsoids
755 :
756 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
757 :
758 : !> \brief
759 : !> Return the lower incomplete gamma function for the specified `exponent` and upper limit.
760 : !> `tolerance` represents the relative accuracy.
761 : !>
762 : !> \warning
763 : !> Do not set tolerance to a number larger than 1, or else the code will crash spectacularly.
764 : !>
765 : !> \remark
766 : !> `logGammaExponent = log_gamma(exponent)`
767 : !>
768 : !> \warning
769 : !> On input, both `exponent` and `upperLim` must be positive, otherwise, a `lowerGamma = -infinity` will be returned to signal error.
770 : ! This algorithm borrows from the `gammp` implementation of the Numerical Recipes by Press et al 1992.
771 39 : pure function getLowerGamma(exponent,logGammaExponent,upperLim,tolerance) result(lowerGamma)
772 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
773 : !DEC$ ATTRIBUTES DLLEXPORT :: getLowerGamma
774 : #endif
775 3 : use Constants_mod, only: IK, RK, HUGE_RK
776 : implicit none
777 : real(RK), intent(in) :: exponent, logGammaExponent, upperLim
778 : real(RK), intent(in), optional :: tolerance
779 : real(RK) :: lowerGamma
780 39 : if ( upperLim < 0._RK .or. exponent <= 0._RK ) then
781 3 : lowerGamma = -HUGE_RK
782 3 : return
783 36 : elseif ( upperLim < exponent + 1._RK ) then
784 24 : lowerGamma = getGammaSeries(exponent,log_gamma(exponent),upperLim,tolerance)
785 : else
786 12 : lowerGamma = 1._RK - getGammaContFrac(exponent,logGammaExponent,upperLim,tolerance)
787 : end if
788 39 : end function getLowerGamma
789 :
790 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
791 :
792 : !> \brief
793 : !> Return the lower incomplete gamma function for the specified exponent and upper limit.
794 : !>
795 : !> `tolerance` represents the relative accuracy.
796 : !>
797 : !> \warning
798 : !> Do not set `tolerance` to a number larger than 1, or else the code will crash spectacularly.
799 : !>
800 : !> \warning
801 : !> On input, both `exponent` and `lowerLim` must be positive, otherwise, a `upperGamma = -infinity` will be returned to signal error.
802 : ! This algorithm borrows from the gammq implementation of the Numerical Recipes by Press et al 1992.
803 39 : pure function getUpperGamma(exponent,logGammaExponent,lowerLim,tolerance) result(upperGamma)
804 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
805 : !DEC$ ATTRIBUTES DLLEXPORT :: getUpperGamma
806 : #endif
807 39 : use Constants_mod, only: IK, RK, HUGE_RK
808 : implicit none
809 : real(RK), intent(in) :: exponent, logGammaExponent, lowerLim
810 : real(RK), intent(in), optional :: tolerance
811 : real(RK) :: upperGamma
812 39 : if ( lowerLim < 0._RK .or. exponent <= 0._RK ) then
813 3 : upperGamma = -HUGE_RK
814 3 : return
815 36 : elseif (lowerLim < exponent + 1._RK) then
816 24 : upperGamma = 1._RK - getGammaSeries(exponent,logGammaExponent,lowerLim,tolerance)
817 : else
818 12 : upperGamma = getGammaContFrac(exponent,logGammaExponent,lowerLim,tolerance)
819 : end if
820 39 : end function getUpperGamma
821 :
822 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
823 :
824 : !> \brief
825 : !> Return the lower incomplete gamma function `P(exponent, upperLim)` evaluated by its series representation as `gammaSeries`.
826 : !>
827 : !> \param[in] exponent : The exponent in lower incomplete gamma function `P(exponent, upperLim)`.
828 : !> \param[in] logGammaExponent : The input `log( gamma(exponent) )`.
829 : !> \param[in] upperLim : The upper limit in lower incomplete gamma function `P(exponent, upperLim)`.
830 : !> \param[in] tolerance : The input relative accuracy.
831 : !>
832 : !> \warning
833 : !> Do not set `tolerance` to a number larger than 1, or else the code will crash spectacularly.
834 : !>
835 : !> \warning
836 : !> If the algorithm fails to converge, `gammaSeries` will be set to negative infinity on output to signal error.
837 : ! This algorithm borrows from the gser implementation of the Numerical Recipes by Press et al 1992.
838 51 : pure function getGammaSeries(exponent,logGammaExponent,upperLim,tolerance) result(gammaSeries)
839 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
840 : !DEC$ ATTRIBUTES DLLEXPORT :: getGammaSeries
841 : #endif
842 39 : use Constants_mod, only: IK, RK, HUGE_RK
843 : implicit none
844 : real(RK), intent(in) :: exponent, logGammaExponent, upperLim
845 : real(RK), intent(in), optional :: tolerance
846 : real(RK) :: gammaSeries
847 : integer(IK), parameter :: ITMAX = 100
848 : real(RK), parameter :: EPS_DEFAULT = epsilon(upperLim)
849 51 : real(RK) :: ap,del,summ,eps
850 : integer(IK) :: iter
851 51 : if (present(tolerance)) then
852 39 : eps = tolerance
853 : else
854 12 : eps = EPS_DEFAULT
855 : end if
856 51 : if (upperLim == 0._RK) then
857 3 : gammaSeries = 0._RK
858 3 : return
859 : end if
860 48 : ap = exponent
861 48 : summ = 1.0_RK / exponent
862 48 : del = summ
863 606 : do iter = 1, ITMAX
864 606 : ap = ap + 1.0_RK
865 606 : del = del*upperLim / ap
866 606 : summ = summ + del
867 606 : if (abs(del) < abs(summ)*eps) exit
868 : end do
869 48 : if (iter>ITMAX) then
870 0 : gammaSeries = -HUGE_RK
871 : else
872 48 : gammaSeries = summ * exp( exponent*log(upperLim) - upperLim - logGammaExponent )
873 : end if
874 51 : end function getGammaSeries
875 :
876 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
877 :
878 : !> \brief
879 : !> Return the incomplete gamma function `Q(exponent, lowerLim)` evaluated by its continued fraction representation as `gammaContFrac`.
880 : !>
881 : !> @param[in] logGammaExponent : This is the `log( gamma(exponent) )`.
882 : !> @param[in] tolerance : Represents the relative accuracy (**optional**).
883 : !> @param[in] exponent : The exponent.
884 : !> @param[in] lowerLim : The lower limit.
885 : !>
886 : !> \warning
887 : !> If the algorithm fails to converge, `gammaContFrac` will be set to negative infinity on output to signal error.
888 : ! This algorithm borrows from the gcf implementation of the Numerical Recipes by Press et al 1992.
889 27 : pure function getGammaContFrac(exponent,logGammaExponent,lowerLim,tolerance) result(gammaContFrac)
890 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
891 : !DEC$ ATTRIBUTES DLLEXPORT :: getGammaContFrac
892 : #endif
893 51 : use Constants_mod, only: IK, RK, HUGE_RK
894 : implicit none
895 : real(RK), intent(in) :: exponent, logGammaExponent, lowerLim
896 : real(RK), optional, intent(in) :: tolerance
897 : real(RK) :: gammaContFrac
898 : integer(IK), parameter :: ITMAX = 100_IK !< The maximum allowed number of iterations.
899 : real(RK), parameter :: EPS_DEFAULT = epsilon(lowerLim) !< The default relative accuracy.
900 : real(RK), parameter :: FPMIN_DEFAULT = tiny(lowerLim) / EPS_DEFAULT !< A number near the smallest representable floating-point number.
901 27 : real(RK) :: eps, fpmin
902 27 : real(RK) :: an,b,c,d,del,h
903 : integer(IK) :: iter
904 27 : if (lowerLim == 0._RK) then
905 3 : gammaContFrac = 1._RK
906 3 : return
907 : end if
908 24 : if (present(tolerance)) then
909 18 : eps = tolerance
910 18 : fpmin = tiny(lowerLim) / eps
911 : else
912 6 : eps = EPS_DEFAULT
913 6 : fpmin = FPMIN_DEFAULT
914 : end if
915 24 : b = lowerLim + 1._RK - exponent
916 24 : c = 1._RK / fpmin
917 24 : d = 1._RK / b
918 24 : h = d
919 198 : do iter = 1, ITMAX
920 198 : an = -iter * (iter-exponent)
921 198 : b = b + 2.0_RK
922 198 : d = an*d + b
923 198 : if (abs(d) < fpmin) d = fpmin
924 198 : c = b + an/c
925 198 : if (abs(c) < fpmin) c = fpmin
926 198 : d = 1.0_RK / d
927 198 : del = d*c
928 198 : h = h*del
929 198 : if (abs(del-1._RK) <= eps) exit
930 : end do
931 : !if (iter > ITMAX) call nrerror('exponent too large, ITMAX too small in getGammaContFrac')
932 24 : if (iter>ITMAX) then
933 0 : gammaContFrac = -HUGE_RK
934 : else
935 24 : gammaContFrac = exp(exponent*log(lowerLim) - lowerLim - logGammaExponent)*h
936 : end if
937 27 : end function getGammaContFrac
938 :
939 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
940 :
941 : end module Math_mod ! LCOV_EXCL_LINE
|
