Hypergeometric 2F1() specfunc: cases x=1

[Richard Mathar <mathar at strw dot leidenuniv dot nl>]

It would be nice if the implementation of the Gaussian hypergeometric
functions would catch the special case of unit argument and compute it
via the standard formula with the Gamma functions (Abramowitz & Stegun 15.1.20).
Currently, all x>0.995 are rejects as being too close to the limit where
the branch cut prohibits the standard series to converge. So my proposal is
to implement in gsl_sf_hyperg_2F1_e() in specfunc/hyperg_2F1.c just
at the start of the function a trap for this case:

....
  const int c_neg_integer = ( c < 0.0  &&  fabs(c - rintc) < locEPS );

  result->val = 0.0;
  result->err = 0.0;

/* start proposed handling of x=1.0 RJM */
   if ( fabs(x - 1.0) < locEPS && c-a-b > 0. && c != 0. && ! c_neg_integer )
   {
        gsl_sf_result lnGamc, lnGamcab, lnGamca, lnGamcb  ;
        int stat = gsl_sf_lngamma_e(c, & lnGamc);
        if ( stat != GSL_SUCCESS )
        {
                DOMAIN_ERROR(result);
        }
        stat = gsl_sf_lngamma_e(c-a-b, & lnGamcab);
        if ( stat != GSL_SUCCESS )
        {
                DOMAIN_ERROR(result);
        }
        stat = gsl_sf_lngamma_e(c-a, & lnGamca);
        if ( stat != GSL_SUCCESS )
        {
                DOMAIN_ERROR(result);
        }
        stat = gsl_sf_lngamma_e(c-b, & lnGamcb);
        if ( stat != GSL_SUCCESS )
        {
                DOMAIN_ERROR(result);
        }
        return gsl_sf_exp_err_e(lnGamc.val+lnGamcab.val-lnGamca.val
                -lnGamcb.val,lnGamc.err+lnGamcab.err+lnGamca.err+lnGamcb.err,result) ;
   }
/*end proposed handling of x=1.0 RJM */

  if(x < -1.0 || 1.0 <= x) {
....

Richard Mathar, http://www.strw.leidenuniv.nl/~mathar