gam0 Subroutine

   subroutine gam0(x,Ga)

      real(wp),intent(in) :: x !! Argument of `Г(x) ( |x| ≤ 1 )`
      real(wp),intent(out) :: Ga !! `Г(x)`

      real(wp) :: gr
      integer :: k

      real(wp),dimension(25),parameter :: g = [ 1.0_wp , &
                                                gamma , &
                                               -0.6558780715202538_wp , &
                                               -0.420026350340952e-1_wp , &
                                                0.1665386113822915_wp , &
                                               -0.421977345555443e-1_wp , &
                                               -0.96219715278770e-2_wp , &
                                                0.72189432466630e-2_wp , &
                                               -0.11651675918591e-2_wp , &
                                               -0.2152416741149e-3_wp , &
                                                0.1280502823882e-3_wp , &
                                               -0.201348547807e-4_wp , &
                                               -0.12504934821e-5_wp , &
                                                0.11330272320e-5_wp , &
                                               -0.2056338417e-6_wp , &
                                                0.61160950e-8_wp , &
                                                0.50020075e-8_wp , &
                                               -0.11812746e-8_wp , &
                                                0.1043427e-9_wp ,  &
                                                0.77823e-11_wp , &
                                               -0.36968e-11_wp , &
                                                0.51e-12_wp , &
                                               -0.206e-13_wp , &
                                               -0.54e-14_wp ,  &
                                                0.14e-14_wp]
      gr = g(25) ! JW: typo in the original code
      do k = 24 , 1 , -1
         gr = gr*x + g(k)
      enddo
      Ga = 1.0_wp/(gr*x)

   end subroutine gam0