htswanalysis/MACS/lib/gsl/gsl-1.11/specfunc/poch.c

   1 /* specfunc/poch.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Author:  G. Jungman */
  21
  22 #include <config.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_sf_exp.h>
  26 #include <gsl/gsl_sf_log.h>
  27 #include <gsl/gsl_sf_pow_int.h>
  28 #include <gsl/gsl_sf_psi.h>
  29 #include <gsl/gsl_sf_gamma.h>
  30
  31 #include "error.h"
  32
  33 static const double bern[21] = {
  34    0.0   /* no element 0 */,
  35   +0.833333333333333333333333333333333e-01,
  36   -0.138888888888888888888888888888888e-02,
  37   +0.330687830687830687830687830687830e-04,
  38   -0.826719576719576719576719576719576e-06,
  39   +0.208767569878680989792100903212014e-07,
  40   -0.528419013868749318484768220217955e-09,
  41   +0.133825365306846788328269809751291e-10,
  42   -0.338968029632258286683019539124944e-12,
  43   +0.858606205627784456413590545042562e-14,
  44   -0.217486869855806187304151642386591e-15,
  45   +0.550900282836022951520265260890225e-17,
  46   -0.139544646858125233407076862640635e-18,
  47   +0.353470703962946747169322997780379e-20,
  48   -0.895351742703754685040261131811274e-22,
  49   +0.226795245233768306031095073886816e-23,
  50   -0.574472439520264523834847971943400e-24,
  51   +0.145517247561486490186626486727132e-26,
  52   -0.368599494066531017818178247990866e-28,
  53   +0.933673425709504467203255515278562e-30,
  54   -0.236502241570062993455963519636983e-31
  55 };
  56
  57
  58 /* ((a)_x - 1)/x in the "small x" region where
  59  * cancellation must be controlled.
  60  *
  61  * Based on SLATEC DPOCH1().
  62  */
  63 /*
  64 C When ABS(X) is so small that substantial cancellation will occur if
  65 C the straightforward formula is used, we use an expansion due
  66 C to Fields and discussed by Y. L. Luke, The Special Functions and Their
  67 C Approximations, Vol. 1, Academic Press, 1969, page 34.
  68 C
  69 C The ratio POCH(A,X) = GAMMA(A+X)/GAMMA(A) is written by Luke as
  70 C        (A+(X-1)/2)**X * polynomial in (A+(X-1)/2)**(-2) .
  71 C In order to maintain significance in POCH1, we write for positive a
  72 C        (A+(X-1)/2)**X = EXP(X*LOG(A+(X-1)/2)) = EXP(Q)
  73 C                       = 1.0 + Q*EXPREL(Q) .
  74 C Likewise the polynomial is written
  75 C        POLY = 1.0 + X*POLY1(A,X) .
  76 C Thus,
  77 C        POCH1(A,X) = (POCH(A,X) - 1) / X
  78 C                   = EXPREL(Q)*(Q/X + Q*POLY1(A,X)) + POLY1(A,X)
  79 C
  80 */
  81 static
  82 int
  83 pochrel_smallx(const double a, const double x, gsl_sf_result * result)
  84 {
  85   /*
  86    SQTBIG = 1.0D0/SQRT(24.0D0*D1MACH(1))
  87    ALNEPS = LOG(D1MACH(3))
  88    */
  89   const double SQTBIG = 1.0/(2.0*M_SQRT2*M_SQRT3*GSL_SQRT_DBL_MIN);
  90   const double ALNEPS = GSL_LOG_DBL_EPSILON - M_LN2;
  91
  92   if(x == 0.0) {
  93     return gsl_sf_psi_e(a, result);
  94   }
  95   else {
  96     const double bp   = (  (a < -0.5) ? 1.0-a-x : a );
  97     const int    incr = ( (bp < 10.0) ? 11.0-bp : 0 );
  98     const double b    = bp + incr;
  99     double dpoch1;
 100     gsl_sf_result dexprl;
 101     int stat_dexprl;
 102     int i;
 103
 104     double var    = b + 0.5*(x-1.0);
 105     double alnvar = log(var);
 106     double q = x*alnvar;
 107
 108     double poly1 = 0.0;
 109
 110     if(var < SQTBIG) {
 111       const int nterms = (int)(-0.5*ALNEPS/alnvar + 1.0);
 112       const double var2 = (1.0/var)/var;
 113       const double rho  = 0.5 * (x + 1.0);
 114       double term = var2;
 115       double gbern[24];
 116       int k, j;
 117
 118       gbern[1] = 1.0;
 119       gbern[2] = -rho/12.0;
 120       poly1 = gbern[2] * term;
 121
 122       if(nterms > 20) {
 123         /* NTERMS IS TOO BIG, MAYBE D1MACH(3) IS BAD */
 124         /* nterms = 20; */
 125         result->val = 0.0;
 126         result->err = 0.0;
 127         GSL_ERROR ("error", GSL_ESANITY);
 128       }
 129
 130       for(k=2; k<=nterms; k++) {
 131         double gbk = 0.0;
 132         for(j=1; j<=k; j++) {
 133           gbk += bern[k-j+1]*gbern[j];
 134         }
 135         gbern[k+1] = -rho*gbk/k;
 136
 137         term  *= (2*k-2-x)*(2*k-1-x)*var2;
 138         poly1 += gbern[k+1]*term;
 139       }
 140     }
 141
 142     stat_dexprl = gsl_sf_expm1_e(q, &dexprl);
 143     if(stat_dexprl != GSL_SUCCESS) {
 144       result->val = 0.0;
 145       result->err = 0.0;
 146       return stat_dexprl;
 147     }
 148     dexprl.val = dexprl.val/q;
 149     poly1 *= (x - 1.0);
 150     dpoch1 = dexprl.val * (alnvar + q * poly1) + poly1;
 151
 152     for(i=incr-1; i >= 0; i--) {
 153       /*
 154        C WE HAVE DPOCH1(B,X), BUT BP IS SMALL, SO WE USE BACKWARDS RECURSION
 155        C TO OBTAIN DPOCH1(BP,X).
 156        */
 157       double binv = 1.0/(bp+i);
 158       dpoch1 = (dpoch1 - binv) / (1.0 + x*binv);
 159     }
 160
 161     if(bp == a) {
 162       result->val = dpoch1;
 163       result->err = 2.0 * GSL_DBL_EPSILON * (fabs(incr) + 1.0) * fabs(result->val);
 164       return GSL_SUCCESS;
 165     }
 166     else {
 167       /*
 168        C WE HAVE DPOCH1(BP,X), BUT A IS LT -0.5.  WE THEREFORE USE A
 169        C REFLECTION FORMULA TO OBTAIN DPOCH1(A,X).
 170        */
 171       double sinpxx = sin(M_PI*x)/x;
 172       double sinpx2 = sin(0.5*M_PI*x);
 173       double t1 = sinpxx/tan(M_PI*b);
 174       double t2 = 2.0*sinpx2*(sinpx2/x);
 175       double trig  = t1 - t2;
 176       result->val  = dpoch1 * (1.0 + x*trig) + trig;
 177       result->err  = (fabs(dpoch1*x) + 1.0) * GSL_DBL_EPSILON * (fabs(t1) + fabs(t2));
 178       result->err += 2.0 * GSL_DBL_EPSILON * (fabs(incr) + 1.0) * fabs(result->val);
 179       return GSL_SUCCESS;
 180     }
 181   }
 182 }
 183
 184
 185 /* Assumes a>0 and a+x>0.
 186  */
 187 static
 188 int
 189 lnpoch_pos(const double a, const double x, gsl_sf_result * result)
 190 {
 191   double absx = fabs(x);
 192
 193   if(absx > 0.1*a || absx*log(GSL_MAX_DBL(a,2.0)) > 0.1) {
 194     if(a < GSL_SF_GAMMA_XMAX && a+x < GSL_SF_GAMMA_XMAX) {
 195       /* If we can do it by calculating the gamma functions
 196        * directly, then that will be more accurate than
 197        * doing the subtraction of the logs.
 198        */
 199       gsl_sf_result g1;
 200       gsl_sf_result g2;
 201       gsl_sf_gammainv_e(a,   &g1);
 202       gsl_sf_gammainv_e(a+x, &g2);
 203       result->val  = -log(g2.val/g1.val);
 204       result->err  = g1.err/fabs(g1.val) + g2.err/fabs(g2.val);
 205       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 206       return GSL_SUCCESS;
 207     }
 208     else {
 209       /* Otherwise we must do the subtraction.
 210        */
 211       gsl_sf_result lg1;
 212       gsl_sf_result lg2;
 213       int stat_1 = gsl_sf_lngamma_e(a,   &lg1);
 214       int stat_2 = gsl_sf_lngamma_e(a+x, &lg2);
 215       result->val  = lg2.val - lg1.val;
 216       result->err  = lg2.err + lg1.err;
 217       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 218       return GSL_ERROR_SELECT_2(stat_1, stat_2);
 219     }
 220   }
 221   else if(absx < 0.1*a && a > 15.0) {
 222     /* Be careful about the implied subtraction.
 223      * Note that both a+x and and a must be
 224      * large here since a is not small
 225      * and x is not relatively large.
 226      * So we calculate using Stirling for Log[Gamma(z)].
 227      *
 228      *   Log[Gamma(a+x)/Gamma(a)] = x(Log[a]-1) + (x+a-1/2)Log[1+x/a]
 229      *                              + (1/(1+eps)   - 1) / (12 a)
 230      *                              - (1/(1+eps)^3 - 1) / (360 a^3)
 231      *                              + (1/(1+eps)^5 - 1) / (1260 a^5)
 232      *                              - (1/(1+eps)^7 - 1) / (1680 a^7)
 233      *                              + ...
 234      */
 235     const double eps = x/a;
 236     const double den = 1.0 + eps;
 237     const double d3 = den*den*den;
 238     const double d5 = d3*den*den;
 239     const double d7 = d5*den*den;
 240     const double c1 = -eps/den;
 241     const double c3 = -eps*(3.0+eps*(3.0+eps))/d3;
 242     const double c5 = -eps*(5.0+eps*(10.0+eps*(10.0+eps*(5.0+eps))))/d5;
 243     const double c7 = -eps*(7.0+eps*(21.0+eps*(35.0+eps*(35.0+eps*(21.0+eps*(7.0+eps))))))/d7;
 244     const double p8 = gsl_sf_pow_int(1.0+eps,8);
 245     const double c8 = 1.0/p8             - 1.0;  /* these need not   */
 246     const double c9 = 1.0/(p8*(1.0+eps)) - 1.0;  /* be very accurate */
 247     const double a4 = a*a*a*a;
 248     const double a6 = a4*a*a;
 249     const double ser_1 = c1 + c3/(30.0*a*a) + c5/(105.0*a4) + c7/(140.0*a6);
 250     const double ser_2 = c8/(99.0*a6*a*a) - 691.0/360360.0 * c9/(a6*a4);
 251     const double ser = (ser_1 + ser_2)/ (12.0*a);
 252
 253     double term1 = x * log(a/M_E);
 254     double term2;
 255     gsl_sf_result ln_1peps;
 256     gsl_sf_log_1plusx_e(eps, &ln_1peps);  /* log(1 + x/a) */
 257     term2 = (x + a - 0.5) * ln_1peps.val;
 258
 259     result->val  = term1 + term2 + ser;
 260     result->err  = GSL_DBL_EPSILON*fabs(term1);
 261     result->err += fabs((x + a - 0.5)*ln_1peps.err);
 262     result->err += fabs(ln_1peps.val) * GSL_DBL_EPSILON * (fabs(x) + fabs(a) + 0.5);
 263     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 264     return GSL_SUCCESS;
 265   }
 266   else {
 267     gsl_sf_result poch_rel;
 268     int stat_p = pochrel_smallx(a, x, &poch_rel);
 269     double eps = x*poch_rel.val;
 270     int stat_e = gsl_sf_log_1plusx_e(eps, result);
 271     result->err  = 2.0 * fabs(x * poch_rel.err / (1.0 + eps));
 272     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 273     return GSL_ERROR_SELECT_2(stat_e, stat_p);
 274   }
 275 }
 276
 277
 278 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 279
 280 int
 281 gsl_sf_lnpoch_e(const double a, const double x, gsl_sf_result * result)
 282 {
 283   /* CHECK_POINTER(result) */
 284
 285   if(a <= 0.0 || a+x <= 0.0) {
 286     DOMAIN_ERROR(result);
 287   }
 288   else if(x == 0.0) {
 289     result->val = 0.0;
 290     result->err = 0.0;
 291     return GSL_SUCCESS;
 292   }
 293   else {
 294     return lnpoch_pos(a, x, result);
 295   }
 296 }
 297
 298
 299 int
 300 gsl_sf_lnpoch_sgn_e(const double a, const double x,
 301                        gsl_sf_result * result, double * sgn)
 302 {
 303   if(a == 0.0 || a+x == 0.0) {
 304     *sgn = 0.0;
 305     DOMAIN_ERROR(result);
 306   }
 307   else if(x == 0.0) {
 308     *sgn = 1.0;
 309     result->val = 0.0;
 310     result->err = 0.0;
 311     return GSL_SUCCESS;
 312   }
 313   else if(a > 0.0 && a+x > 0.0) {
 314     *sgn = 1.0;
 315     return lnpoch_pos(a, x, result);
 316   }
 317   else if(a < 0.0 && a+x < 0.0) {
 318     /* Reduce to positive case using reflection.
 319      */
 320     double sin_1 = sin(M_PI * (1.0 - a));
 321     double sin_2 = sin(M_PI * (1.0 - a - x));
 322     if(sin_1 == 0.0 || sin_2 == 0.0) {
 323       *sgn = 0.0;
 324       DOMAIN_ERROR(result);
 325     }
 326     else {
 327       gsl_sf_result lnp_pos;
 328       int stat_pp   = lnpoch_pos(1.0-a, -x, &lnp_pos);
 329       double lnterm = log(fabs(sin_1/sin_2));
 330       result->val  = lnterm - lnp_pos.val;
 331       result->err  = lnp_pos.err;
 332       result->err += 2.0 * GSL_DBL_EPSILON * (fabs(1.0-a) + fabs(1.0-a-x)) * fabs(lnterm);
 333       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 334       *sgn = GSL_SIGN(sin_1*sin_2);
 335       return stat_pp;
 336     }
 337   }
 338   else {
 339     /* Evaluate gamma ratio directly.
 340      */
 341     gsl_sf_result lg_apn;
 342     gsl_sf_result lg_a;
 343     double s_apn, s_a;
 344     int stat_apn = gsl_sf_lngamma_sgn_e(a+x, &lg_apn, &s_apn);
 345     int stat_a   = gsl_sf_lngamma_sgn_e(a,   &lg_a,   &s_a);
 346     if(stat_apn == GSL_SUCCESS && stat_a == GSL_SUCCESS) {
 347       result->val  = lg_apn.val - lg_a.val;
 348       result->err  = lg_apn.err + lg_a.err;
 349       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 350       *sgn = s_a * s_apn;
 351       return GSL_SUCCESS;
 352     }
 353     else if(stat_apn == GSL_EDOM || stat_a == GSL_EDOM){
 354       *sgn = 0.0;
 355       DOMAIN_ERROR(result);
 356     }
 357     else {
 358       result->val = 0.0;
 359       result->err = 0.0;
 360       *sgn = 0.0;
 361       return GSL_FAILURE;
 362     }
 363   }
 364 }
 365
 366
 367 int
 368 gsl_sf_poch_e(const double a, const double x, gsl_sf_result * result)
 369 {
 370   /* CHECK_POINTER(result) */
 371
 372   if(x == 0.0) {
 373     result->val = 1.0;
 374     result->err = 0.0;
 375     return GSL_SUCCESS;
 376   }
 377   else {
 378     gsl_sf_result lnpoch;
 379     double sgn;
 380     int stat_lnpoch = gsl_sf_lnpoch_sgn_e(a, x, &lnpoch, &sgn);
 381     int stat_exp    = gsl_sf_exp_err_e(lnpoch.val, lnpoch.err, result);
 382     result->val *= sgn;
 383     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 384     return GSL_ERROR_SELECT_2(stat_exp, stat_lnpoch);
 385   }
 386 }
 387
 388
 389 int
 390 gsl_sf_pochrel_e(const double a, const double x, gsl_sf_result * result)
 391 {
 392   const double absx = fabs(x);
 393   const double absa = fabs(a);
 394
 395   /* CHECK_POINTER(result) */
 396
 397   if(absx > 0.1*absa || absx*log(GSL_MAX(absa,2.0)) > 0.1) {
 398     gsl_sf_result lnpoch;
 399     double sgn;
 400     int stat_poch = gsl_sf_lnpoch_sgn_e(a, x, &lnpoch, &sgn);
 401     if(lnpoch.val > GSL_LOG_DBL_MAX) {
 402       OVERFLOW_ERROR(result);
 403     }
 404     else {
 405       const double el = exp(lnpoch.val);
 406       result->val  = (sgn*el - 1.0)/x;
 407       result->err  = fabs(result->val) * (lnpoch.err + 2.0 * GSL_DBL_EPSILON);
 408       result->err += 2.0 * GSL_DBL_EPSILON * (fabs(sgn*el) + 1.0) / fabs(x);
 409       return stat_poch;
 410     }
 411   }
 412   else {
 413     return pochrel_smallx(a, x, result);
 414   }
 415 }
 416
 417
 418 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 419
 420 #include "eval.h"
 421
 422 double gsl_sf_lnpoch(const double a, const double x)
 423 {
 424   EVAL_RESULT(gsl_sf_lnpoch_e(a, x, &result));
 425 }
 426
 427 double gsl_sf_poch(const double a, const double x)
 428 {
 429   EVAL_RESULT(gsl_sf_poch_e(a, x, &result));
 430 }
 431
 432 double gsl_sf_pochrel(const double a, const double x)
 433 {
 434   EVAL_RESULT(gsl_sf_pochrel_e(a, x, &result));
 435 }