Skip to content

Commit e3ee922

Browse files
perazzperazz
committed
gamma: generalize internal calculation to "next" more accurate precision
1 parent 2bdc50e commit e3ee922

File tree

1 file changed

+58
-72
lines changed

1 file changed

+58
-72
lines changed

‎src/stdlib_specialfunctions_gamma.fypp‎

Lines changed: 58 additions & 72 deletions
Original file line numberDiff line numberDiff line change
@@ -1,28 +1,24 @@
11
#:include "common.fypp"
2-
#:set R_KINDS_TYPES = [KT for KT in REAL_KINDS_TYPES if KT[0] in ["sp","dp"]]
3-
#:set C_KINDS_TYPES = [KT for KT in CMPLX_KINDS_TYPES if KT[0] in ["sp","dp"]]
4-
#:set CI_KINDS_TYPES = INT_KINDS_TYPES + C_KINDS_TYPES
2+
#:set CI_KINDS_TYPES = INT_KINDS_TYPES + CMPLX_KINDS_TYPES
3+
#:set IDX_CMPLX_KINDS_TYPES = [(i, CMPLX_KINDS[i], CMPLX_TYPES[i], CMPLX_INIT[i]) for i in range(len(CMPLX_KINDS))]
4+
#:set IDX_REAL_KINDS_TYPES = [(i, REAL_KINDS[i], REAL_TYPES[i], REAL_INIT[i]) for i in range(len(REAL_KINDS))]
55
module stdlib_specialfunctions_gamma
6-
use iso_fortran_env, only : qp => real128
76
use ieee_arithmetic, only: ieee_value, ieee_quiet_nan
8-
use stdlib_kinds, only : sp, dp, int8, int16, int32, int64
7+
use stdlib_kinds, only : sp, dp, xdp, qp, int8, int16, int32, int64
98
use stdlib_error, only : error_stop
109

1110
implicit none
1211
private
1312

14-
integer(int8), parameter :: max_fact_int8 = 6_int8
13+
integer(int8), parameter :: max_fact_int8 = 6_int8
1514
integer(int16), parameter :: max_fact_int16 = 8_int16
1615
integer(int32), parameter :: max_fact_int32 = 13_int32
1716
integer(int64), parameter :: max_fact_int64 = 21_int64
1817

19-
#:for k1, t1 in R_KINDS_TYPES
18+
#:for k1, t1 in REAL_KINDS_TYPES
2019
${t1}$, parameter :: tol_${k1}$ = epsilon(1.0_${k1}$)
2120
#:endfor
22-
real(qp), parameter :: tol_qp = epsilon(1.0_qp)
23-
24-
25-
21+
2622
public :: gamma, log_gamma, log_factorial
2723
public :: lower_incomplete_gamma, log_lower_incomplete_gamma
2824
public :: upper_incomplete_gamma, log_upper_incomplete_gamma
@@ -33,7 +29,7 @@ module stdlib_specialfunctions_gamma
3329
interface gamma
3430
!! Gamma function for integer and complex numbers
3531
!!
36-
#:for k1, t1 in CI_KINDS_TYPES
32+
#:for k1, t1 in CI_KINDS_TYPES[:-1]
3733
module procedure gamma_${t1[0]}$${k1}$
3834
#:endfor
3935
end interface gamma
@@ -43,7 +39,7 @@ module stdlib_specialfunctions_gamma
4339
interface log_gamma
4440
!! Logarithm of gamma function
4541
!!
46-
#:for k1, t1 in CI_KINDS_TYPES
42+
#:for k1, t1 in CI_KINDS_TYPES[:-1]
4743
module procedure l_gamma_${t1[0]}$${k1}$
4844
#:endfor
4945
end interface log_gamma
@@ -64,12 +60,12 @@ module stdlib_specialfunctions_gamma
6460
!! Lower incomplete gamma function
6561
!!
6662
#:for k1, t1 in INT_KINDS_TYPES
67-
#:for k2, t2 in R_KINDS_TYPES
63+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
6864
module procedure ingamma_low_${t1[0]}$${k1}$${k2}$
6965
#:endfor
7066
#:endfor
7167

72-
#:for k1, t1 in R_KINDS_TYPES
68+
#:for k1, t1 in REAL_KINDS_TYPES[:-1]
7369
module procedure ingamma_low_${t1[0]}$${k1}$
7470
#:endfor
7571
end interface lower_incomplete_gamma
@@ -80,12 +76,12 @@ module stdlib_specialfunctions_gamma
8076
!! Logarithm of lower incomplete gamma function
8177
!!
8278
#:for k1, t1 in INT_KINDS_TYPES
83-
#:for k2, t2 in R_KINDS_TYPES
79+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
8480
module procedure l_ingamma_low_${t1[0]}$${k1}$${k2}$
8581
#:endfor
8682
#:endfor
8783

88-
#:for k1, t1 in R_KINDS_TYPES
84+
#:for k1, t1 in REAL_KINDS_TYPES[:-1]
8985
module procedure l_ingamma_low_${t1[0]}$${k1}$
9086
#:endfor
9187
end interface log_lower_incomplete_gamma
@@ -96,12 +92,12 @@ module stdlib_specialfunctions_gamma
9692
!! Upper incomplete gamma function
9793
!!
9894
#:for k1, t1 in INT_KINDS_TYPES
99-
#:for k2, t2 in R_KINDS_TYPES
95+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
10096
module procedure ingamma_up_${t1[0]}$${k1}$${k2}$
10197
#:endfor
10298
#:endfor
10399

104-
#:for k1, t1 in R_KINDS_TYPES
100+
#:for k1, t1 in REAL_KINDS_TYPES[:-1]
105101
module procedure ingamma_up_${t1[0]}$${k1}$
106102
#:endfor
107103
end interface upper_incomplete_gamma
@@ -112,12 +108,12 @@ module stdlib_specialfunctions_gamma
112108
!! Logarithm of upper incomplete gamma function
113109
!!
114110
#:for k1, t1 in INT_KINDS_TYPES
115-
#:for k2, t2 in R_KINDS_TYPES
111+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
116112
module procedure l_ingamma_up_${t1[0]}$${k1}$${k2}$
117113
#:endfor
118114
#:endfor
119115

120-
#:for k1, t1 in R_KINDS_TYPES
116+
#:for k1, t1 in REAL_KINDS_TYPES[:-1]
121117
module procedure l_ingamma_up_${t1[0]}$${k1}$
122118
#:endfor
123119
end interface log_upper_incomplete_gamma
@@ -128,12 +124,12 @@ module stdlib_specialfunctions_gamma
128124
!! Regularized (normalized) lower incomplete gamma function, P
129125
!!
130126
#:for k1, t1 in INT_KINDS_TYPES
131-
#:for k2, t2 in R_KINDS_TYPES
127+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
132128
module procedure regamma_p_${t1[0]}$${k1}$${k2}$
133129
#:endfor
134130
#:endfor
135131

136-
#:for k1, t1 in R_KINDS_TYPES
132+
#:for k1, t1 in REAL_KINDS_TYPES[:-1]
137133
module procedure regamma_p_${t1[0]}$${k1}$
138134
#:endfor
139135
end interface regularized_gamma_p
@@ -144,12 +140,12 @@ module stdlib_specialfunctions_gamma
144140
!! Regularized (normalized) upper incomplete gamma function, Q
145141
!!
146142
#:for k1, t1 in INT_KINDS_TYPES
147-
#:for k2, t2 in R_KINDS_TYPES
143+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
148144
module procedure regamma_q_${t1[0]}$${k1}$${k2}$
149145
#:endfor
150146
#:endfor
151147

152-
#:for k1, t1 in R_KINDS_TYPES
148+
#:for k1, t1 in REAL_KINDS_TYPES[:-1]
153149
module procedure regamma_q_${t1[0]}$${k1}$
154150
#:endfor
155151
end interface regularized_gamma_q
@@ -160,12 +156,12 @@ module stdlib_specialfunctions_gamma
160156
! Incomplete gamma G function.
161157
! Internal use only
162158
!
163-
#:for k1, t1 in R_KINDS_TYPES
159+
#:for k1, t1 in REAL_KINDS_TYPES[:-1]
164160
module procedure gpx_${t1[0]}$${k1}$ !for real p and x
165161
#:endfor
166162

167163
#:for k1, t1 in INT_KINDS_TYPES
168-
#:for k2, t2 in R_KINDS_TYPES
164+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
169165
module procedure gpx_${t1[0]}$${k1}$${k2}$ !for integer p and real x
170166
#:endfor
171167
#:endfor
@@ -178,7 +174,7 @@ module stdlib_specialfunctions_gamma
178174
! Internal use only
179175
!
180176
#:for k1, t1 in INT_KINDS_TYPES
181-
#:for k2, t2 in R_KINDS_TYPES
177+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
182178
module procedure l_gamma_${t1[0]}$${k1}$${k2}$
183179
#:endfor
184180
#:endfor
@@ -219,14 +215,12 @@ contains
219215

220216

221217

222-
#:for k1, t1 in C_KINDS_TYPES
223-
#:if k1 == "sp"
224-
#:set k2 = "dp"
225-
#:elif k1 == "dp"
226-
#:set k2 = "qp"
227-
#:endif
228-
#:set t2 = "real({})".format(k2)
229-
218+
#! Because the KIND lists are sorted by increasing accuracy,
219+
#! gamma will use the next available more accurate KIND for the
220+
#! internal more accurate solver.
221+
#:for i, k1, t1, i1 in IDX_CMPLX_KINDS_TYPES[:-1]
222+
#:set k2 = CMPLX_KINDS[i + 1]
223+
#:set t2 = "real({})".format(k2)
230224
impure elemental function gamma_${t1[0]}$${k1}$(z) result(res)
231225
${t1}$, intent(in) :: z
232226
${t1}$ :: res
@@ -255,8 +249,8 @@ contains
255249
-2.71994908488607704e-9_${k2}$]
256250
! parameters from above referenced source.
257251

258-
#:elif k1 == "dp"
259-
#! for double precision input, using quadruple precision for calculation
252+
#:else
253+
#! for double or extended precision input, using quadruple precision for calculation
260254

261255
integer, parameter :: n = 24
262256
${t2}$, parameter :: r = 25.617904_${k2}$
@@ -290,8 +284,6 @@ contains
290284

291285
#:endif
292286

293-
294-
295287
if(abs(z % im) < tol_${k1}$) then
296288

297289
res = cmplx(gamma(z % re), kind = ${k1}$)
@@ -333,9 +325,6 @@ contains
333325

334326
#:endfor
335327

336-
337-
338-
339328
#:for k1, t1 in INT_KINDS_TYPES
340329
impure elemental function l_gamma_${t1[0]}$${k1}$(z) result(res)
341330
!
@@ -374,7 +363,7 @@ contains
374363

375364

376365
#:for k1, t1 in INT_KINDS_TYPES
377-
#:for k2, t2 in R_KINDS_TYPES
366+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
378367

379368
impure elemental function l_gamma_${t1[0]}$${k1}$${k2}$(z, x) result(res)
380369
!
@@ -415,13 +404,12 @@ contains
415404

416405

417406

418-
#:for k1, t1 in C_KINDS_TYPES
419-
#:if k1 == "sp"
420-
#:set k2 = "dp"
421-
#:elif k1 == "dp"
422-
#:set k2 = "qp"
423-
#:endif
424-
#:set t2 = "real({})".format(k2)
407+
#! Because the KIND lists are sorted by increasing accuracy,
408+
#! gamma will use the next available more accurate KIND for the
409+
#! internal more accurate solver.
410+
#:for i, k1, t1, i1 in IDX_CMPLX_KINDS_TYPES[:-1]
411+
#:set k2 = CMPLX_KINDS[i + 1]
412+
#:set t2 = "real({})".format(k2)
425413
impure elemental function l_gamma_${t1[0]}$${k1}$(z) result (res)
426414
!
427415
! log_gamma function for any complex number, excluding negative whole number
@@ -557,14 +545,12 @@ contains
557545

558546

559547

560-
#:for k1, t1 in R_KINDS_TYPES
561-
#:if k1 == "sp"
562-
#:set k2 = "dp"
563-
#:elif k1 == "dp"
564-
#:set k2 = "qp"
565-
#:endif
566-
#:set t2 = "real({})".format(k2)
567-
548+
#! Because the KIND lists are sorted by increasing accuracy,
549+
#! gamma will use the next available more accurate KIND for the
550+
#! internal more accurate solver.
551+
#:for i, k1, t1, i1 in IDX_REAL_KINDS_TYPES[:-1]
552+
#:set k2 = REAL_KINDS[i + 1]
553+
#:set t2 = REAL_TYPES[i + 1]
568554
impure elemental function gpx_${t1[0]}$${k1}$(p, x) result(res)
569555
!
570556
! Approximation of incomplete gamma G function with real argument p.
@@ -685,7 +671,7 @@ contains
685671

686672

687673
#:for k1, t1 in INT_KINDS_TYPES
688-
#:for k2, t2 in R_KINDS_TYPES
674+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
689675
impure elemental function gpx_${t1[0]}$${k1}$${k2}$(p, x) result(res)
690676
!
691677
! Approximation of incomplete gamma G function with integer argument p.
@@ -824,7 +810,7 @@ contains
824810

825811

826812

827-
#:for k1, t1 in R_KINDS_TYPES
813+
#:for k1, t1 in REAL_KINDS_TYPES[:-1]
828814
impure elemental function ingamma_low_${t1[0]}$${k1}$(p, x) result(res)
829815
!
830816
! Approximation of lower incomplete gamma function with real p.
@@ -861,7 +847,7 @@ contains
861847

862848

863849
#:for k1, t1 in INT_KINDS_TYPES
864-
#:for k2, t2 in R_KINDS_TYPES
850+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
865851
impure elemental function ingamma_low_${t1[0]}$${k1}$${k2}$(p, x) &
866852
result(res)
867853
!
@@ -901,7 +887,7 @@ contains
901887

902888

903889

904-
#:for k1, t1 in R_KINDS_TYPES
890+
#:for k1, t1 in REAL_KINDS_TYPES[:-1]
905891
impure elemental function l_ingamma_low_${t1[0]}$${k1}$(p, x) result(res)
906892

907893
${t1}$, intent(in) :: p, x
@@ -938,7 +924,7 @@ contains
938924

939925

940926
#:for k1, t1 in INT_KINDS_TYPES
941-
#:for k2, t2 in R_KINDS_TYPES
927+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
942928
impure elemental function l_ingamma_low_${t1[0]}$${k1}$${k2}$(p, x) &
943929
result(res)
944930

@@ -970,7 +956,7 @@ contains
970956

971957

972958

973-
#:for k1, t1 in R_KINDS_TYPES
959+
#:for k1, t1 in REAL_KINDS_TYPES[:-1]
974960
impure elemental function ingamma_up_${t1[0]}$${k1}$(p, x) result(res)
975961
!
976962
! Approximation of upper incomplete gamma function with real p.
@@ -1008,7 +994,7 @@ contains
1008994

1009995

1010996
#:for k1, t1 in INT_KINDS_TYPES
1011-
#:for k2, t2 in R_KINDS_TYPES
997+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
1012998
impure elemental function ingamma_up_${t1[0]}$${k1}$${k2}$(p, x) &
1013999
result(res)
10141000
!
@@ -1050,7 +1036,7 @@ contains
10501036

10511037

10521038

1053-
#:for k1, t1 in R_KINDS_TYPES
1039+
#:for k1, t1 in REAL_KINDS_TYPES[:-1]
10541040
impure elemental function l_ingamma_up_${t1[0]}$${k1}$(p, x) result(res)
10551041

10561042
${t1}$, intent(in) :: p, x
@@ -1088,7 +1074,7 @@ contains
10881074

10891075

10901076
#:for k1, t1 in INT_KINDS_TYPES
1091-
#:for k2, t2 in R_KINDS_TYPES
1077+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
10921078
impure elemental function l_ingamma_up_${t1[0]}$${k1}$${k2}$(p, x) &
10931079
result(res)
10941080

@@ -1129,7 +1115,7 @@ contains
11291115

11301116

11311117

1132-
#:for k1, t1 in R_KINDS_TYPES
1118+
#:for k1, t1 in REAL_KINDS_TYPES[:-1]
11331119
impure elemental function regamma_p_${t1[0]}$${k1}$(p, x) result(res)
11341120
!
11351121
! Approximation of regularized incomplete gamma function P(p,x) for real p
@@ -1164,7 +1150,7 @@ contains
11641150

11651151

11661152
#:for k1, t1 in INT_KINDS_TYPES
1167-
#:for k2, t2 in R_KINDS_TYPES
1153+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
11681154
impure elemental function regamma_p_${t1[0]}$${k1}$${k2}$(p, x) result(res)
11691155
!
11701156
! Approximation of regularized incomplete gamma function P(p,x) for integer p
@@ -1200,7 +1186,7 @@ contains
12001186

12011187

12021188

1203-
#:for k1, t1 in R_KINDS_TYPES
1189+
#:for k1, t1 in REAL_KINDS_TYPES[:-1]
12041190
impure elemental function regamma_q_${t1[0]}$${k1}$(p, x) result(res)
12051191
!
12061192
! Approximation of regularized incomplete gamma function Q(p,x) for real p
@@ -1235,7 +1221,7 @@ contains
12351221

12361222

12371223
#:for k1, t1 in INT_KINDS_TYPES
1238-
#:for k2, t2 in R_KINDS_TYPES
1224+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
12391225
impure elemental function regamma_q_${t1[0]}$${k1}$${k2}$(p, x) result(res)
12401226
!
12411227
! Approximation of regularized incomplet gamma function Q(p,x) for integer p

0 commit comments

Comments
 (0)