chguit Subroutine

   subroutine chguit(a,b,x,Hu,Id)

      real(wp),intent(in) :: a !! Parameter ( `a > 0` )
      real(wp),intent(in) :: b !! Parameter
      real(wp),intent(in) :: x !! Argument ( `x > 0` )
      real(wp),intent(out) :: Hu !! `U(a,b,z)`
      integer,intent(out) :: Id !! Estimated number of significant digits

      real(wp) :: a1 , b1 , c , d , f1 , f2 , g , ga ,&
                  hu0 , hu1 , hu2 , s, t1 , t2 , t3 , t4
      integer :: j , k , m

      real(wp),dimension(30),parameter :: t = &
            [.259597723012478e-01_wp , .778093339495366e-01_wp , &
             .129449135396945e+00_wp , .180739964873425e+00_wp , &
             .231543551376029e+00_wp , .281722937423262e+00_wp , &
             .331142848268448e+00_wp , .379670056576798e+00_wp , &
             .427173741583078e+00_wp , .473525841761707e+00_wp , &
             .518601400058570e+00_wp , .562278900753945e+00_wp , &
             .604440597048510e+00_wp , .644972828489477e+00_wp , &
             .683766327381356e+00_wp , .720716513355730e+00_wp , &
             .755723775306586e+00_wp , .788693739932264e+00_wp , &
             .819537526162146e+00_wp , .848171984785930e+00_wp , &
             .874519922646898e+00_wp , .898510310810046e+00_wp , &
             .920078476177628e+00_wp , .939166276116423e+00_wp , &
             .955722255839996e+00_wp , .969701788765053e+00_wp , &
             .981067201752598e+00_wp , .989787895222222e+00_wp , &
             .995840525118838e+00_wp , .999210123227436e+00_wp ]

      real(wp),dimension(30),parameter :: w = &
         [ .519078776312206e-01_wp , .517679431749102e-01_wp , &
           .514884515009810e-01_wp , .510701560698557e-01_wp , &
           .505141845325094e-01_wp , .498220356905502e-01_wp , &
           .489955754557568e-01_wp , .480370318199712e-01_wp , &
           .469489888489122e-01_wp , .457343797161145e-01_wp , &
           .443964787957872e-01_wp , .429388928359356e-01_wp , &
           .413655512355848e-01_wp , .396806954523808e-01_wp , &
           .378888675692434e-01_wp , .359948980510845e-01_wp , &
           .340038927249464e-01_wp , .319212190192963e-01_wp , &
           .297524915007890e-01_wp , .275035567499248e-01_wp , &
           .251804776215213e-01_wp , .227895169439978e-01_wp , &
           .203371207294572e-01_wp , .178299010142074e-01_wp , &
           .152746185967848e-01_wp , .126781664768159e-01_wp , &
           .100475571822880e-01_wp , .738993116334531e-02_wp , &
           .471272992695363e-02_wp , .202681196887362e-02_wp ]

      Id = 9
      ! DLMF 13.4.4, integration up to C=12/X
      a1 = a - 1.0_wp
      b1 = b - a - 1.0_wp
      c = 12.0_wp/x
      hu0 = 0.0_wp
      do m = 10 , 100 , 5
         hu1 = 0.0_wp
         g = 0.5_wp*c/m
         d = g
         do j = 1 , m
            s = 0.0_wp
            do k = 1 , 30
               t1 = d + g*t(k)
               t2 = d - g*t(k)
               f1 = exp(-x*t1)*t1**a1*(1.0_wp+t1)**b1
               f2 = exp(-x*t2)*t2**a1*(1.0_wp+t2)**b1
               s = s + w(k)*(f1+f2)
            enddo
            hu1 = hu1 + s*g
            d = d + 2.0_wp*g
         enddo
         if ( abs(1.0_wp-hu0/hu1)<1.0e-9_wp ) exit
         hu0 = hu1
      enddo
      call gamma2(a,ga)
      hu1 = hu1/ga
      ! DLMF 13.4.4 with substitution t=C/(1-u)
      ! integration u from 0 to 1, i.e. t from C=12/X to infinity
      do m = 2 , 10 , 2
         hu2 = 0.0_wp
         g = 0.5_wp/m
         d = g
         do j = 1 , m
            s = 0.0_wp
            do k = 1 , 30
               t1 = d + g*t(k)
               t2 = d - g*t(k)
               t3 = c/(1.0_wp-t1)
               t4 = c/(1.0_wp-t2)
               f1 = t3*t3/c*exp(-x*t3)*t3**a1*(1.0_wp+t3)**b1
               f2 = t4*t4/c*exp(-x*t4)*t4**a1*(1.0_wp+t4)**b1
               s = s + w(k)*(f1+f2)
            enddo
            hu2 = hu2 + s*g
            d = d + 2.0_wp*g
         enddo
         if ( abs(1.0_wp-hu0/hu2)<1.0e-9_wp ) exit
         hu0 = hu2
      enddo
      call gamma2(a,ga)
      hu2 = hu2/ga
      Hu = hu1 + hu2

   end subroutine chguit