Skip to content

Commit 220b79f

Browse files
perazzperazz
committed
generalize l_gamma, l_factorial kinds
1 parent b4f0fcb commit 220b79f

File tree

2 files changed

+54
-74
lines changed

2 files changed

+54
-74
lines changed

‎src/stdlib_specialfunctions_gamma.fypp

+8-7
Original file line numberDiff line numberDiff line change
@@ -331,7 +331,7 @@ contains
331331
! Logarithm of gamma function for integer input
332332
!
333333
${t1}$, intent(in) :: z
334-
real :: res
334+
real(dp) :: res
335335
${t1}$ :: i
336336
${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$, two = 2_${k1}$
337337

@@ -350,7 +350,7 @@ contains
350350

351351
do i = one, z - one
352352

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

355355
end do
356356

@@ -512,14 +512,15 @@ contains
512512

513513

514514
#:for k1, t1 in INT_KINDS_TYPES
515+
#:set k2, t2 = REAL_KINDS[-2], REAL_TYPES[-2]
515516
impure elemental function l_factorial_${t1[0]}$${k1}$(n) result(res)
516517
!
517518
! Log(n!)
518519
!
519520
${t1}$, intent(in) :: n
520-
real(dp) :: res
521+
${t2}$ :: res
521522
${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$, two = 2_${k1}$
522-
real(dp), parameter :: zero_dp = 0.0_dp
523+
${t2}$, parameter :: zero_${k2}$ = 0.0_${k2}$, one_${k2}$ = 1.0_${k2}$
523524

524525
if(n < zero) call error_stop("Error(l_factorial): Logarithm of" &
525526
//" factorial function argument must be non-negative")
@@ -528,15 +529,15 @@ contains
528529

529530
case (zero)
530531

531-
res = zero_dp
532+
res = zero_${k2}$
532533

533534
case (one)
534535

535-
res = zero_dp
536+
res = zero_${k2}$
536537

537538
case (two : )
538539

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

541542
end select
542543
end function l_factorial_${t1[0]}$${k1}$

‎test/specialfunctions/test_specialfunctions_gamma.fypp

+46-67
Original file line numberDiff line numberDiff line change
@@ -80,42 +80,32 @@ contains
8080

8181

8282
#:for k1, t1 in INT_KINDS_TYPES
83-
83+
#:set k2, t2 = REAL_KINDS[-2], REAL_TYPES[-2]
8484
subroutine test_logfact_${t1[0]}$${k1}$(error)
85-
type(error_type), allocatable, intent(out) :: error
86-
integer, parameter :: n = 6
85+
type(error_type), allocatable, intent(out) :: error
8786
integer :: i
88-
89-
#:if k1 == "int8"
90-
91-
${t1}$, parameter :: x(n) = [0_${k1}$, 1_${k1}$, 2_${k1}$, 4_${k1}$, &
92-
5_${k1}$, 100_${k1}$]
93-
real(dp), parameter :: ans(n) = [0.0_dp, 0.0_dp, 0.69314718055994_dp, 3.17805383034794_dp, &
94-
4.78749174278204_dp, 3.637393755555e2_dp]
95-
96-
#:elif k1 == "int16"
97-
98-
${t1}$, parameter :: x(n) = [0_${k1}$, 1_${k1}$, 2_${k1}$, 4_${k1}$, &
99-
7_${k1}$, 500_${k1}$]
100-
real(dp), parameter :: ans(n) = [0.0_dp, 0.0_dp, 0.69314718055994_dp, 3.17805383034794_dp, &
101-
8.52516136106541_dp, 2.611330458460e3_dp]
102-
#:elif k1 == "int32"
103-
104-
${t1}$, parameter :: x(n) = [0_${k1}$, 1_${k1}$, 2_${k1}$, 4_${k1}$, &
105-
12_${k1}$, 7000_${k1}$]
106-
real(dp), parameter :: ans(n) = [0.0_dp, 0.0_dp, 0.69314718055994_dp, 3.17805383034794_dp, &
107-
1.998721449566e1_dp, 5.498100377941e4_dp]
108-
#:elif k1 == "int64"
109-
110-
${t1}$, parameter :: x(n) = [0_${k1}$, 1_${k1}$, 2_${k1}$, 4_${k1}$, &
111-
20_${k1}$, 90000_${k1}$]
112-
real(dp), parameter :: ans(n) = [0.0_dp, 0.0_dp, 0.69314718055994_dp, 3.17805383034794_dp, &
113-
4.233561646075e1_dp, 9.366874681600e5_dp]
114-
#: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)
115105

116106
do i = 1, n
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}$
@@ -201,57 +191,46 @@ contains
201191

202192
subroutine test_loggamma_${t1[0]}$${k1}$(error)
203193
type(error_type), allocatable, intent(out) :: error
204-
integer, parameter :: n = 4
205194
integer :: i
206195
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)
207212

208-
#:if k1 == "int8"
209-
210-
${t1}$, parameter :: x(n) = [1_${k1}$, 2_${k1}$, 10_${k1}$, 47_${k1}$]
211-
real(sp), parameter :: ans(n) = [0.0, 0.0, 1.28018274e1, 1.32952575e2]
212-
213-
#:elif k1 == "int16"
214-
215-
${t1}$, parameter :: x(n) = [1_${k1}$, 2_${k1}$, 111_${k1}$, 541_${k1}$]
216-
real(sp), parameter :: ans(n) = [0.0, 0.0, 4.10322777e2, 2.86151221e3]
217-
218-
#:elif k1 == "int32"
219-
220-
${t1}$, parameter :: x(n) = [1_${k1}$, 2_${k1}$, 2021_${k1}$, &
221-
42031_${k1}$]
222-
real(sp), parameter :: ans(n) = [0.0, 0.0, 1.33586470e4, 4.05433461e5]
223-
224-
#: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.)
225217

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

230220
#:elif t1[0] == "c"
231221

232-
${t1}$, parameter :: x(n) = [(0.25_${k1}$, 0.25_${k1}$), &
222+
${t1}$, parameter :: x(*) = [(0.25_${k1}$, 0.25_${k1}$), &
233223
(0.5_${k1}$, -0.5_${k1}$), &
234224
(1.0_${k1}$, 1.0_${k1}$), &
235225
(-1.254e1_${k1}$, -9.87_${k1}$)]
236226

237-
${t1}$, parameter :: ans(n) = &
227+
${t1}$, parameter :: ans(*) = &
238228
[(0.90447450949333889_${k1}$, -0.83887024394321282_${k1}$),&
239229
(0.11238724280962311_${k1}$, 0.75072920212205074_${k1}$), &
240230
(-0.65092319930185634_${k1}$, -0.30164032046753320_${k1}$),&
241231
(-4.7091788015763380e1_${k1}$, 1.4804627819235690e1_${k1}$)]
242-
#:endif
243-
244-
245-
#:if t1[0] == "i"
246-
247-
do i = 1, n
248-
249-
call check(error, log_gamma(x(i)), ans(i), "Integer kind ${k1}$ " &
250-
//"failed", thr = tol_sp, rel = .true.)
251-
252-
end do
253-
254-
#:elif t1[0] == "c"
232+
233+
integer, parameter :: n = size(x)
255234

256235
do i = 1, n
257236

0 commit comments

Comments
 (0)