From mboxrd@z Thu Jan  1 00:00:00 1970
Return-Path: <gsl-discuss-return-5211-listarch-gsl-discuss=sources.redhat.com@sourceware.org>
Received: (qmail 26477 invoked by alias); 11 Dec 2007 21:16:27 -0000
Received: (qmail 25966 invoked by uid 22791); 11 Dec 2007 21:14:46 -0000
X-Spam-Check-By: sourceware.org
Date: Tue, 11 Dec 2007 21:16:00 -0000
From: Richard Mathar <mathar@strw.leidenuniv.nl>
Message-Id: <200712112114.lBBLERNe005280@amer.strw.leidenuniv.nl>
To: gsl-discuss@sourceware.org
Subject: Hypergeometric 2F1() specfunc: cases x=1
Cc: mathar@mail.strw.leidenuniv.nl
X-LeidenObservatory-MailScanner-Information: Processed through MailScanner at Leiden Observatory
X-LeidenObservatory-MailScanner: Virusscan: Found to be clean
X-LeidenObservatory-MailScanner-SpamCheck: not spam,
	SpamAssassin (not cached, score=-2.599, required 4.5,
	autolearn=not spam, BAYES_00 -2.60)
X-LeidenObservatory-MailScanner-From: mathar@mail.strw.leidenuniv.nl
Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm
Precedence: bulk
List-Id: <gsl-discuss.sourceware.org>
List-Subscribe: <mailto:gsl-discuss-subscribe@sourceware.org>
List-Archive: <http://sourceware.org/ml/gsl-discuss/>
List-Post: <mailto:gsl-discuss@sourceware.org>
List-Help: <mailto:gsl-discuss-help@sourceware.org>, <http://sourceware.org/ml/#faqs>
Sender: gsl-discuss-owner@sourceware.org
X-SW-Source: 2007-q4/txt/msg00038.txt.bz2


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