Skip to content
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.

Commit 1ba499c

Browse files
perazzperazz
authoredApr 8, 2025··
specialfunctions: generalize gamma precision (#969)
2 parents def840d + f9d999f commit 1ba499c

File tree

4 files changed

+147
-177
lines changed

4 files changed

+147
-177
lines changed
 

‎example/specialfunctions_gamma/example_gamma.f90

Lines changed: 2 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,21 +1,19 @@
11
program example_gamma
2-
use stdlib_kinds, only: dp, int64
2+
use stdlib_kinds, only: sp, dp, int64
33
use stdlib_specialfunctions_gamma, only: gamma
44
implicit none
55

66
integer :: i
77
integer(int64) :: n
88
real :: x
99
real(dp) :: y
10-
complex :: z
11-
complex(dp) :: z1
10+
complex(sp) :: z
1211

1312
i = 10
1413
n = 15_int64
1514
x = 2.5
1615
y = 4.3_dp
1716
z = (2.3, 0.6)
18-
z1 = (-4.2_dp, 3.1_dp)
1917

2018
print *, gamma(i) !integer gives exact result
2119
! 362880
@@ -32,6 +30,4 @@ program example_gamma
3230
print *, gamma(z)
3331
! (0.988054395, 0.383354813)
3432

35-
print *, gamma(z1)
36-
! (-2.78916032990983999E-005, 9.83164600163221218E-006)
3733
end program example_gamma

‎example/specialfunctions_gamma/example_log_gamma.f90

Lines changed: 3 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1,19 +1,18 @@
11
program example_log_gamma
2-
use stdlib_kinds, only: dp
2+
use stdlib_kinds, only: sp, dp
33
use stdlib_specialfunctions_gamma, only: log_gamma
44
implicit none
55

66
integer :: i
77
real :: x
88
real(dp) :: y
9-
complex :: z
10-
complex(dp) :: z1
9+
complex(sp) :: z
1110

1211
i = 10
1312
x = 8.76
1413
y = x
1514
z = (5.345, -3.467)
16-
z1 = z
15+
1716
print *, log_gamma(i) !default single precision output
1817
!12.8018274
1918

@@ -29,7 +28,4 @@ program example_log_gamma
2928

3029
!(2.56165648, -5.73382425)
3130

32-
print *, log_gamma(z1)
33-
34-
!(2.5616575105114614, -5.7338247782852498)
3531
end program example_log_gamma

‎src/stdlib_specialfunctions_gamma.fypp

Lines changed: 68 additions & 81 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] if k1 == "sp" else CMPLX_KINDS[-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,16 +325,13 @@ 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
!
342331
! Logarithm of gamma function for integer input
343332
!
344333
${t1}$, intent(in) :: z
345-
real :: res
334+
real(dp) :: res
346335
${t1}$ :: i
347336
${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$, two = 2_${k1}$
348337

@@ -361,7 +350,7 @@ contains
361350

362351
do i = one, z - one
363352

364-
res = res + log(real(i))
353+
res = res + log(real(i,dp))
365354

366355
end do
367356

@@ -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] if k1 == "sp" else CMPLX_KINDS[-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
@@ -436,7 +424,7 @@ contains
436424
real(${k1}$) :: d
437425
integer :: m, i
438426
complex(${k2}$) :: zr, zr2, sum, s
439-
real(${k1}$), parameter :: z_limit = 10_${k1}$, zero_k1 = 0.0_${k1}$
427+
real(${k1}$), parameter :: z_limit = 10.0_${k1}$, zero_k1 = 0.0_${k1}$
440428
integer, parameter :: n = 20
441429
${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$, &
442430
pi = acos(-one), ln2pi = log(2 * pi)
@@ -524,14 +512,15 @@ contains
524512

525513

526514
#:for k1, t1 in INT_KINDS_TYPES
515+
#:set k2, t2 = REAL_KINDS[-2], REAL_TYPES[-2]
527516
impure elemental function l_factorial_${t1[0]}$${k1}$(n) result(res)
528517
!
529518
! Log(n!)
530519
!
531520
${t1}$, intent(in) :: n
532-
real(dp) :: res
521+
${t2}$ :: res
533522
${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$, two = 2_${k1}$
534-
real(dp), parameter :: zero_dp = 0.0_dp
523+
${t2}$, parameter :: zero_${k2}$ = 0.0_${k2}$, one_${k2}$ = 1.0_${k2}$
535524

536525
if(n < zero) call error_stop("Error(l_factorial): Logarithm of" &
537526
//" factorial function argument must be non-negative")
@@ -540,15 +529,15 @@ contains
540529

541530
case (zero)
542531

543-
res = zero_dp
532+
res = zero_${k2}$
544533

545534
case (one)
546535

547-
res = zero_dp
536+
res = zero_${k2}$
548537

549538
case (two : )
550539

551-
res = l_gamma(n + 1, 1.0_dp)
540+
res = l_gamma(n + 1, one_${k2}$)
552541

553542
end select
554543
end function l_factorial_${t1[0]}$${k1}$
@@ -557,14 +546,12 @@ contains
557546

558547

559548

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-
549+
#! Because the KIND lists are sorted by increasing accuracy,
550+
#! gamma will use the next available more accurate KIND for the
551+
#! internal more accurate solver.
552+
#:for i, k1, t1, i1 in IDX_REAL_KINDS_TYPES[:-1]
553+
#:set k2 = REAL_KINDS[i + 1] if k1 == "sp" else REAL_KINDS[-1]
554+
#:set t2 = REAL_TYPES[i + 1]
568555
impure elemental function gpx_${t1[0]}$${k1}$(p, x) result(res)
569556
!
570557
! Approximation of incomplete gamma G function with real argument p.
@@ -685,7 +672,7 @@ contains
685672

686673

687674
#:for k1, t1 in INT_KINDS_TYPES
688-
#:for k2, t2 in R_KINDS_TYPES
675+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
689676
impure elemental function gpx_${t1[0]}$${k1}$${k2}$(p, x) result(res)
690677
!
691678
! Approximation of incomplete gamma G function with integer argument p.
@@ -732,7 +719,7 @@ contains
732719

733720
if(mod(n, 2) == 0) then
734721

735-
a = (1 - p - n / 2) * x
722+
a = (one - p - n / 2) * x
736723

737724
else
738725

@@ -824,7 +811,7 @@ contains
824811

825812

826813

827-
#:for k1, t1 in R_KINDS_TYPES
814+
#:for k1, t1 in REAL_KINDS_TYPES[:-1]
828815
impure elemental function ingamma_low_${t1[0]}$${k1}$(p, x) result(res)
829816
!
830817
! Approximation of lower incomplete gamma function with real p.
@@ -861,7 +848,7 @@ contains
861848

862849

863850
#:for k1, t1 in INT_KINDS_TYPES
864-
#:for k2, t2 in R_KINDS_TYPES
851+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
865852
impure elemental function ingamma_low_${t1[0]}$${k1}$${k2}$(p, x) &
866853
result(res)
867854
!
@@ -901,7 +888,7 @@ contains
901888

902889

903890

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

907894
${t1}$, intent(in) :: p, x
@@ -938,7 +925,7 @@ contains
938925

939926

940927
#:for k1, t1 in INT_KINDS_TYPES
941-
#:for k2, t2 in R_KINDS_TYPES
928+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
942929
impure elemental function l_ingamma_low_${t1[0]}$${k1}$${k2}$(p, x) &
943930
result(res)
944931

@@ -970,7 +957,7 @@ contains
970957

971958

972959

973-
#:for k1, t1 in R_KINDS_TYPES
960+
#:for k1, t1 in REAL_KINDS_TYPES[:-1]
974961
impure elemental function ingamma_up_${t1[0]}$${k1}$(p, x) result(res)
975962
!
976963
! Approximation of upper incomplete gamma function with real p.
@@ -1008,7 +995,7 @@ contains
1008995

1009996

1010997
#:for k1, t1 in INT_KINDS_TYPES
1011-
#:for k2, t2 in R_KINDS_TYPES
998+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
1012999
impure elemental function ingamma_up_${t1[0]}$${k1}$${k2}$(p, x) &
10131000
result(res)
10141001
!
@@ -1050,7 +1037,7 @@ contains
10501037

10511038

10521039

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

10561043
${t1}$, intent(in) :: p, x
@@ -1088,7 +1075,7 @@ contains
10881075

10891076

10901077
#:for k1, t1 in INT_KINDS_TYPES
1091-
#:for k2, t2 in R_KINDS_TYPES
1078+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
10921079
impure elemental function l_ingamma_up_${t1[0]}$${k1}$${k2}$(p, x) &
10931080
result(res)
10941081

@@ -1129,7 +1116,7 @@ contains
11291116

11301117

11311118

1132-
#:for k1, t1 in R_KINDS_TYPES
1119+
#:for k1, t1 in REAL_KINDS_TYPES[:-1]
11331120
impure elemental function regamma_p_${t1[0]}$${k1}$(p, x) result(res)
11341121
!
11351122
! Approximation of regularized incomplete gamma function P(p,x) for real p
@@ -1164,7 +1151,7 @@ contains
11641151

11651152

11661153
#:for k1, t1 in INT_KINDS_TYPES
1167-
#:for k2, t2 in R_KINDS_TYPES
1154+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
11681155
impure elemental function regamma_p_${t1[0]}$${k1}$${k2}$(p, x) result(res)
11691156
!
11701157
! Approximation of regularized incomplete gamma function P(p,x) for integer p
@@ -1200,7 +1187,7 @@ contains
12001187

12011188

12021189

1203-
#:for k1, t1 in R_KINDS_TYPES
1190+
#:for k1, t1 in REAL_KINDS_TYPES[:-1]
12041191
impure elemental function regamma_q_${t1[0]}$${k1}$(p, x) result(res)
12051192
!
12061193
! Approximation of regularized incomplete gamma function Q(p,x) for real p
@@ -1235,7 +1222,7 @@ contains
12351222

12361223

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

‎test/specialfunctions/test_specialfunctions_gamma.fypp

Lines changed: 74 additions & 83 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,11 @@
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 test_specialfunctions_gamma
66
use testdrive, only : new_unittest, unittest_type, error_type, check
7-
use stdlib_kinds, only: sp, dp, int8, int16, int32, int64
7+
use stdlib_kinds, only: sp, dp, xdp, qp, int8, int16, int32, int64
8+
use stdlib_error, only: state_type, STDLIB_VALUE_ERROR
89
use stdlib_specialfunctions_gamma, only: gamma, log_gamma, log_factorial, &
910
lower_incomplete_gamma, &
1011
upper_incomplete_gamma, &
@@ -18,8 +19,8 @@ module test_specialfunctions_gamma
1819

1920
public :: collect_specialfunctions_gamma
2021

21-
#:for k1, t1 in R_KINDS_TYPES
22-
${t1}$, parameter :: tol_${k1}$ = 1000 * epsilon(1.0_${k1}$)
22+
#:for k1, t1 in REAL_KINDS_TYPES
23+
${t1}$, parameter :: tol_${k1}$ = sqrt(epsilon(1.0_${k1}$))
2324
#:endfor
2425

2526
contains
@@ -35,15 +36,15 @@ contains
3536
test_logfact_${t1[0]}$${k1}$) &
3637
#:endfor
3738

38-
#:for k1, t1 in CI_KINDS_TYPES
39+
#:for k1, t1 in CI_KINDS_TYPES[:-1]
3940
, new_unittest("gamma_${t1[0]}$${k1}$", &
4041
test_gamma_${t1[0]}$${k1}$) &
4142
, new_unittest("log_gamma_${t1[0]}$${k1}$", &
4243
test_loggamma_${t1[0]}$${k1}$) &
4344
#:endfor
4445

4546
#:for k1, t1 in INT_KINDS_TYPES
46-
#:for k2, t2 in R_KINDS_TYPES
47+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
4748
, new_unittest("lower_incomplete_gamma_${t1[0]}$${k1}$${k2}$", &
4849
test_lincgamma_${t1[0]}$${k1}$${k2}$) &
4950
, new_unittest("log_lower_incomplete_gamma_${t1[0]}$${k1}$${k2}$", &
@@ -59,7 +60,7 @@ contains
5960
#:endfor
6061
#:endfor
6162

62-
#:for k1, t1 in R_KINDS_TYPES
63+
#:for k1, t1 in REAL_KINDS_TYPES[:-1]
6364
, new_unittest("lower_incomplete_gamma_${t1[0]}$${k1}$", &
6465
test_lincgamma_${t1[0]}$${k1}$) &
6566
, new_unittest("log_lower_incomplete_gamma_${t1[0]}$${k1}$", &
@@ -79,43 +80,32 @@ contains
7980

8081

8182
#:for k1, t1 in INT_KINDS_TYPES
82-
83+
#:set k2, t2 = REAL_KINDS[-2], REAL_TYPES[-2]
8384
subroutine test_logfact_${t1[0]}$${k1}$(error)
84-
type(error_type), allocatable, intent(out) :: error
85-
integer, parameter :: n = 6
85+
type(error_type), allocatable, intent(out) :: error
8686
integer :: i
87-
88-
#:if k1 == "int8"
89-
90-
${t1}$, parameter :: x(n) = [0_${k1}$, 1_${k1}$, 2_${k1}$, 4_${k1}$, &
91-
5_${k1}$, 100_${k1}$]
92-
real(dp), parameter :: ans(n) = [0.0_dp, 0.0_dp, 0.69314718055994_dp, 3.17805383034794_dp, &
93-
4.78749174278204_dp, 3.637393755555e2_dp]
94-
95-
#:elif k1 == "int16"
96-
97-
${t1}$, parameter :: x(n) = [0_${k1}$, 1_${k1}$, 2_${k1}$, 4_${k1}$, &
98-
7_${k1}$, 500_${k1}$]
99-
real(dp), parameter :: ans(n) = [0.0_dp, 0.0_dp, 0.69314718055994_dp, 3.17805383034794_dp, &
100-
8.52516136106541_dp, 2.611330458460e3_dp]
101-
#:elif k1 == "int32"
102-
103-
${t1}$, parameter :: x(n) = [0_${k1}$, 1_${k1}$, 2_${k1}$, 4_${k1}$, &
104-
12_${k1}$, 7000_${k1}$]
105-
real(dp), parameter :: ans(n) = [0.0_dp, 0.0_dp, 0.69314718055994_dp, 3.17805383034794_dp, &
106-
1.998721449566e1_dp, 5.498100377941e4_dp]
107-
#:elif k1 == "int64"
108-
109-
${t1}$, parameter :: x(n) = [0_${k1}$, 1_${k1}$, 2_${k1}$, 4_${k1}$, &
110-
20_${k1}$, 90000_${k1}$]
111-
real(dp), parameter :: ans(n) = [0.0_dp, 0.0_dp, 0.69314718055994_dp, 3.17805383034794_dp, &
112-
4.233561646075e1_dp, 9.366874681600e5_dp]
113-
#:endif
87+
88+
integer, parameter :: xtest(*) = [0,1,2,4,5,7,12,20,100,500,7000,90000]
89+
${t2}$, parameter :: res(*) = [0.0_${k2}$, &
90+
0.0_${k2}$, &
91+
0.69314718055994_${k2}$, &
92+
3.17805383034794_${k2}$, &
93+
4.78749174278204_${k2}$, &
94+
8.52516136106541_${k2}$, &
95+
1.998721449566e1_${k2}$, &
96+
4.233561646075e1_${k2}$, &
97+
3.637393755555e2_${k2}$, &
98+
2.611330458460e3_${k2}$, &
99+
5.498100377941e4_${k2}$, &
100+
9.366874681600e5_${k2}$]
101+
102+
${t1}$, parameter :: x(*) = pack(xtest, xtest<huge(0_${k1}$))
103+
${t2}$, parameter :: ans(*) = pack(res , xtest<huge(0_${k1}$))
104+
integer, parameter :: n = size(x)
114105

115106
do i = 1, n
116-
117107
call check(error, log_factorial(x(i)), ans(i), "Integer kind " &
118-
//"${k1}$ failed", thr = tol_dp, rel = .true.)
108+
//"${k1}$ failed", thr = tol_${k2}$, rel = .true.)
119109

120110
end do
121111
end subroutine test_logfact_${t1[0]}$${k1}$
@@ -124,12 +114,13 @@ contains
124114

125115

126116

127-
#:for k1, t1 in CI_KINDS_TYPES
117+
#:for k1, t1 in CI_KINDS_TYPES[:-1]
128118

129119
subroutine test_gamma_${t1[0]}$${k1}$(error)
130120
type(error_type), allocatable, intent(out) :: error
131121
integer, parameter :: n = 5
132122
integer :: i
123+
type(state_type) :: err
133124

134125
#:if k1 == "int8"
135126

@@ -182,74 +173,74 @@ contains
182173
#:elif t1[0] == "c"
183174

184175
do i = 1, n
176+
177+
err = state_type(STDLIB_VALUE_ERROR,'Complex ${k1}$ failed: x=',x(i), &
178+
' gamma=',gamma(x(i)), &
179+
'expected=',ans(i), &
180+
' tol=',tol_${k1}$)
185181

186-
call check(error, gamma(x(i)), ans(i), "Complex kind ${k1}$ failed",&
182+
call check(error, gamma(x(i)), ans(i), err%print(),&
187183
thr = tol_${k1}$, rel = .true.)
188184

189185
end do
190-
191186
#:endif
187+
192188
end subroutine test_gamma_${t1[0]}$${k1}$
193189

194190

195191

196192
subroutine test_loggamma_${t1[0]}$${k1}$(error)
197193
type(error_type), allocatable, intent(out) :: error
198-
integer, parameter :: n = 4
199194
integer :: i
195+
type(state_type) :: err
196+
197+
#:if t1[0] == "i"
198+
199+
integer, parameter :: xtest(*) = [1,2,10,47,111,541,2021,42031]
200+
real(dp), parameter :: res(*) = [0.0_dp, &
201+
0.0_dp, &
202+
1.28018274e1_dp, &
203+
1.32952575e2_dp, &
204+
4.10322777e2_dp, &
205+
2.86151221e3_dp, &
206+
1.33586470e4_dp, &
207+
4.05433461e5_dp]
208+
209+
integer, parameter :: x(*) = pack(xtest,xtest<huge(0_${k1}$))
210+
real(dp), parameter :: ans(*) = pack(res,xtest<huge(0_${k1}$))
211+
integer, parameter :: n = size(x)
200212

201-
#:if k1 == "int8"
202-
203-
${t1}$, parameter :: x(n) = [1_${k1}$, 2_${k1}$, 10_${k1}$, 47_${k1}$]
204-
real(sp), parameter :: ans(n) = [0.0, 0.0, 1.28018274e1, 1.32952575e2]
205-
206-
#:elif k1 == "int16"
207-
208-
${t1}$, parameter :: x(n) = [1_${k1}$, 2_${k1}$, 111_${k1}$, 541_${k1}$]
209-
real(sp), parameter :: ans(n) = [0.0, 0.0, 4.10322777e2, 2.86151221e3]
210-
211-
#:elif k1 == "int32"
212-
213-
${t1}$, parameter :: x(n) = [1_${k1}$, 2_${k1}$, 2021_${k1}$, &
214-
42031_${k1}$]
215-
real(sp), parameter :: ans(n) = [0.0, 0.0, 1.33586470e4, 4.05433461e5]
216-
217-
#:elif k1 == "int64"
213+
do i = 1, n
214+
print *, 'log ',log_gamma(x(i)),' ans=',ans(i),' tol=',tol_dp
215+
call check(error, log_gamma(x(i)), ans(i), "Integer kind ${k1}$ " &
216+
//"failed", thr = tol_dp, rel = .true.)
218217

219-
${t1}$, parameter :: x(n) = [1_${k1}$, 2_${k1}$, 2021_${k1}$, &
220-
42031_${k1}$]
221-
real(sp), parameter :: ans(n) = [0.0, 0.0, 1.33586470e4, 4.05433461e5]
218+
end do
222219

223220
#:elif t1[0] == "c"
224221

225-
${t1}$, parameter :: x(n) = [(0.25_${k1}$, 0.25_${k1}$), &
222+
${t1}$, parameter :: x(*) = [(0.25_${k1}$, 0.25_${k1}$), &
226223
(0.5_${k1}$, -0.5_${k1}$), &
227224
(1.0_${k1}$, 1.0_${k1}$), &
228225
(-1.254e1_${k1}$, -9.87_${k1}$)]
229226

230-
${t1}$, parameter :: ans(n) = &
227+
${t1}$, parameter :: ans(*) = &
231228
[(0.90447450949333889_${k1}$, -0.83887024394321282_${k1}$),&
232229
(0.11238724280962311_${k1}$, 0.75072920212205074_${k1}$), &
233230
(-0.65092319930185634_${k1}$, -0.30164032046753320_${k1}$),&
234231
(-4.7091788015763380e1_${k1}$, 1.4804627819235690e1_${k1}$)]
235-
#:endif
236-
237-
238-
#:if t1[0] == "i"
232+
233+
integer, parameter :: n = size(x)
239234

240235
do i = 1, n
241236

242-
call check(error, log_gamma(x(i)), ans(i), "Integer kind ${k1}$ " &
243-
//"failed", thr = tol_sp, rel = .true.)
244-
245-
end do
237+
err = state_type(STDLIB_VALUE_ERROR,'Complex ${k1}$ failed: x=',x(i), &
238+
' log(gamma)=',log_gamma(x(i)), &
239+
'expected=',ans(i), &
240+
' tol=',tol_${k1}$)
246241

247-
#:elif t1[0] == "c"
248-
249-
do i = 1, n
250-
251-
call check(error, log_gamma(x(i)), ans(i), "Complex kind ${k1}$ " &
252-
//"failed", thr = tol_${k1}$, rel = .true.)
242+
call check(error, log_gamma(x(i)), ans(i), err%print(),&
243+
thr = tol_${k1}$, rel = .true.)
253244

254245
end do
255246

@@ -262,7 +253,7 @@ contains
262253

263254

264255
#:for k1, t1 in INT_KINDS_TYPES
265-
#:for k2, t2 in R_KINDS_TYPES
256+
#:for k2, t2 in REAL_KINDS_TYPES[:-1]
266257

267258
subroutine test_lincgamma_${t1[0]}$${k1}$${k2}$(error)
268259
type(error_type), allocatable, intent(out) :: error
@@ -411,7 +402,7 @@ contains
411402

412403

413404

414-
#:for k1, t1 in R_KINDS_TYPES
405+
#:for k1, t1 in REAL_KINDS_TYPES[:-1]
415406

416407
subroutine test_lincgamma_${t1[0]}$${k1}$(error)
417408
type(error_type), allocatable, intent(out) :: error

0 commit comments

Comments
 (0)
Please sign in to comment.