2 ! PR 33683 - we used to pick up the wrong gamma function
3 ! from the library on some systems.
6 integer, parameter :: n_max = 20
7 double precision, dimension(0:n_max) :: c
10 double precision :: td, xd
16 c(n) = (2*n-1)*c(n-1)*0.5d0
24 if (abs(gamma(xs)-ts)/ts > 9e-6) call abort
25 if (abs(gamma(xd)-td)/td > 5e-14) call abort
27 call tst_s(2.3, gamma(2.3))
28 call tst_s(3.7, gamma(3.7))
29 call tst_s(5.5, gamma(5.5))
30 call tst_d(4.2d0, gamma(4.2d0))
31 call tst_d(8.1d0, gamma(8.1d0))
33 subroutine tst_s(a, b)
35 if (abs(gamma(a) - b)/b > 1e-6) call abort
38 subroutine tst_d(a, b)
39 double precision :: a,b
40 if (abs(gamma(a) - b)/b > 5e-14) call abort