summaryrefslogtreecommitdiff
path: root/gsl-1.9/specfunc/poch.c
diff options
context:
space:
mode:
Diffstat (limited to 'gsl-1.9/specfunc/poch.c')
-rw-r--r--gsl-1.9/specfunc/poch.c435
1 files changed, 435 insertions, 0 deletions
diff --git a/gsl-1.9/specfunc/poch.c b/gsl-1.9/specfunc/poch.c
new file mode 100644
index 0000000..782900e
--- /dev/null
+++ b/gsl-1.9/specfunc/poch.c
@@ -0,0 +1,435 @@
+/* specfunc/poch.c
+ *
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+/* Author: G. Jungman */
+
+#include <config.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_sf_exp.h>
+#include <gsl/gsl_sf_log.h>
+#include <gsl/gsl_sf_pow_int.h>
+#include <gsl/gsl_sf_psi.h>
+#include <gsl/gsl_sf_gamma.h>
+
+#include "error.h"
+
+static const double bern[21] = {
+ 0.0 /* no element 0 */,
+ +0.833333333333333333333333333333333e-01,
+ -0.138888888888888888888888888888888e-02,
+ +0.330687830687830687830687830687830e-04,
+ -0.826719576719576719576719576719576e-06,
+ +0.208767569878680989792100903212014e-07,
+ -0.528419013868749318484768220217955e-09,
+ +0.133825365306846788328269809751291e-10,
+ -0.338968029632258286683019539124944e-12,
+ +0.858606205627784456413590545042562e-14,
+ -0.217486869855806187304151642386591e-15,
+ +0.550900282836022951520265260890225e-17,
+ -0.139544646858125233407076862640635e-18,
+ +0.353470703962946747169322997780379e-20,
+ -0.895351742703754685040261131811274e-22,
+ +0.226795245233768306031095073886816e-23,
+ -0.574472439520264523834847971943400e-24,
+ +0.145517247561486490186626486727132e-26,
+ -0.368599494066531017818178247990866e-28,
+ +0.933673425709504467203255515278562e-30,
+ -0.236502241570062993455963519636983e-31
+};
+
+
+/* ((a)_x - 1)/x in the "small x" region where
+ * cancellation must be controlled.
+ *
+ * Based on SLATEC DPOCH1().
+ */
+/*
+C When ABS(X) is so small that substantial cancellation will occur if
+C the straightforward formula is used, we use an expansion due
+C to Fields and discussed by Y. L. Luke, The Special Functions and Their
+C Approximations, Vol. 1, Academic Press, 1969, page 34.
+C
+C The ratio POCH(A,X) = GAMMA(A+X)/GAMMA(A) is written by Luke as
+C (A+(X-1)/2)**X * polynomial in (A+(X-1)/2)**(-2) .
+C In order to maintain significance in POCH1, we write for positive a
+C (A+(X-1)/2)**X = EXP(X*LOG(A+(X-1)/2)) = EXP(Q)
+C = 1.0 + Q*EXPREL(Q) .
+C Likewise the polynomial is written
+C POLY = 1.0 + X*POLY1(A,X) .
+C Thus,
+C POCH1(A,X) = (POCH(A,X) - 1) / X
+C = EXPREL(Q)*(Q/X + Q*POLY1(A,X)) + POLY1(A,X)
+C
+*/
+static
+int
+pochrel_smallx(const double a, const double x, gsl_sf_result * result)
+{
+ /*
+ SQTBIG = 1.0D0/SQRT(24.0D0*D1MACH(1))
+ ALNEPS = LOG(D1MACH(3))
+ */
+ const double SQTBIG = 1.0/(2.0*M_SQRT2*M_SQRT3*GSL_SQRT_DBL_MIN);
+ const double ALNEPS = GSL_LOG_DBL_EPSILON - M_LN2;
+
+ if(x == 0.0) {
+ return gsl_sf_psi_e(a, result);
+ }
+ else {
+ const double bp = ( (a < -0.5) ? 1.0-a-x : a );
+ const int incr = ( (bp < 10.0) ? 11.0-bp : 0 );
+ const double b = bp + incr;
+ double dpoch1;
+ gsl_sf_result dexprl;
+ int stat_dexprl;
+ int i;
+
+ double var = b + 0.5*(x-1.0);
+ double alnvar = log(var);
+ double q = x*alnvar;
+
+ double poly1 = 0.0;
+
+ if(var < SQTBIG) {
+ const int nterms = (int)(-0.5*ALNEPS/alnvar + 1.0);
+ const double var2 = (1.0/var)/var;
+ const double rho = 0.5 * (x + 1.0);
+ double term = var2;
+ double gbern[24];
+ int k, j;
+
+ gbern[1] = 1.0;
+ gbern[2] = -rho/12.0;
+ poly1 = gbern[2] * term;
+
+ if(nterms > 20) {
+ /* NTERMS IS TOO BIG, MAYBE D1MACH(3) IS BAD */
+ /* nterms = 20; */
+ result->val = 0.0;
+ result->err = 0.0;
+ GSL_ERROR ("error", GSL_ESANITY);
+ }
+
+ for(k=2; k<=nterms; k++) {
+ double gbk = 0.0;
+ for(j=1; j<=k; j++) {
+ gbk += bern[k-j+1]*gbern[j];
+ }
+ gbern[k+1] = -rho*gbk/k;
+
+ term *= (2*k-2-x)*(2*k-1-x)*var2;
+ poly1 += gbern[k+1]*term;
+ }
+ }
+
+ stat_dexprl = gsl_sf_expm1_e(q, &dexprl);
+ if(stat_dexprl != GSL_SUCCESS) {
+ result->val = 0.0;
+ result->err = 0.0;
+ return stat_dexprl;
+ }
+ dexprl.val = dexprl.val/q;
+ poly1 *= (x - 1.0);
+ dpoch1 = dexprl.val * (alnvar + q * poly1) + poly1;
+
+ for(i=incr-1; i >= 0; i--) {
+ /*
+ C WE HAVE DPOCH1(B,X), BUT BP IS SMALL, SO WE USE BACKWARDS RECURSION
+ C TO OBTAIN DPOCH1(BP,X).
+ */
+ double binv = 1.0/(bp+i);
+ dpoch1 = (dpoch1 - binv) / (1.0 + x*binv);
+ }
+
+ if(bp == a) {
+ result->val = dpoch1;
+ result->err = 2.0 * GSL_DBL_EPSILON * (fabs(incr) + 1.0) * fabs(result->val);
+ return GSL_SUCCESS;
+ }
+ else {
+ /*
+ C WE HAVE DPOCH1(BP,X), BUT A IS LT -0.5. WE THEREFORE USE A
+ C REFLECTION FORMULA TO OBTAIN DPOCH1(A,X).
+ */
+ double sinpxx = sin(M_PI*x)/x;
+ double sinpx2 = sin(0.5*M_PI*x);
+ double t1 = sinpxx/tan(M_PI*b);
+ double t2 = 2.0*sinpx2*(sinpx2/x);
+ double trig = t1 - t2;
+ result->val = dpoch1 * (1.0 + x*trig) + trig;
+ result->err = (fabs(dpoch1*x) + 1.0) * GSL_DBL_EPSILON * (fabs(t1) + fabs(t2));
+ result->err += 2.0 * GSL_DBL_EPSILON * (fabs(incr) + 1.0) * fabs(result->val);
+ return GSL_SUCCESS;
+ }
+ }
+}
+
+
+/* Assumes a>0 and a+x>0.
+ */
+static
+int
+lnpoch_pos(const double a, const double x, gsl_sf_result * result)
+{
+ double absx = fabs(x);
+
+ if(absx > 0.1*a || absx*log(GSL_MAX_DBL(a,2.0)) > 0.1) {
+ if(a < GSL_SF_GAMMA_XMAX && a+x < GSL_SF_GAMMA_XMAX) {
+ /* If we can do it by calculating the gamma functions
+ * directly, then that will be more accurate than
+ * doing the subtraction of the logs.
+ */
+ gsl_sf_result g1;
+ gsl_sf_result g2;
+ gsl_sf_gammainv_e(a, &g1);
+ gsl_sf_gammainv_e(a+x, &g2);
+ result->val = -log(g2.val/g1.val);
+ result->err = g1.err/fabs(g1.val) + g2.err/fabs(g2.val);
+ result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
+ return GSL_SUCCESS;
+ }
+ else {
+ /* Otherwise we must do the subtraction.
+ */
+ gsl_sf_result lg1;
+ gsl_sf_result lg2;
+ int stat_1 = gsl_sf_lngamma_e(a, &lg1);
+ int stat_2 = gsl_sf_lngamma_e(a+x, &lg2);
+ result->val = lg2.val - lg1.val;
+ result->err = lg2.err + lg1.err;
+ result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
+ return GSL_ERROR_SELECT_2(stat_1, stat_2);
+ }
+ }
+ else if(absx < 0.1*a && a > 15.0) {
+ /* Be careful about the implied subtraction.
+ * Note that both a+x and and a must be
+ * large here since a is not small
+ * and x is not relatively large.
+ * So we calculate using Stirling for Log[Gamma(z)].
+ *
+ * Log[Gamma(a+x)/Gamma(a)] = x(Log[a]-1) + (x+a-1/2)Log[1+x/a]
+ * + (1/(1+eps) - 1) / (12 a)
+ * - (1/(1+eps)^3 - 1) / (360 a^3)
+ * + (1/(1+eps)^5 - 1) / (1260 a^5)
+ * - (1/(1+eps)^7 - 1) / (1680 a^7)
+ * + ...
+ */
+ const double eps = x/a;
+ const double den = 1.0 + eps;
+ const double d3 = den*den*den;
+ const double d5 = d3*den*den;
+ const double d7 = d5*den*den;
+ const double c1 = -eps/den;
+ const double c3 = -eps*(3.0+eps*(3.0+eps))/d3;
+ const double c5 = -eps*(5.0+eps*(10.0+eps*(10.0+eps*(5.0+eps))))/d5;
+ const double c7 = -eps*(7.0+eps*(21.0+eps*(35.0+eps*(35.0+eps*(21.0+eps*(7.0+eps))))))/d7;
+ const double p8 = gsl_sf_pow_int(1.0+eps,8);
+ const double c8 = 1.0/p8 - 1.0; /* these need not */
+ const double c9 = 1.0/(p8*(1.0+eps)) - 1.0; /* be very accurate */
+ const double a4 = a*a*a*a;
+ const double a6 = a4*a*a;
+ const double ser_1 = c1 + c3/(30.0*a*a) + c5/(105.0*a4) + c7/(140.0*a6);
+ const double ser_2 = c8/(99.0*a6*a*a) - 691.0/360360.0 * c9/(a6*a4);
+ const double ser = (ser_1 + ser_2)/ (12.0*a);
+
+ double term1 = x * log(a/M_E);
+ double term2;
+ gsl_sf_result ln_1peps;
+ gsl_sf_log_1plusx_e(eps, &ln_1peps); /* log(1 + x/a) */
+ term2 = (x + a - 0.5) * ln_1peps.val;
+
+ result->val = term1 + term2 + ser;
+ result->err = GSL_DBL_EPSILON*fabs(term1);
+ result->err += fabs((x + a - 0.5)*ln_1peps.err);
+ result->err += fabs(ln_1peps.val) * GSL_DBL_EPSILON * (fabs(x) + fabs(a) + 0.5);
+ result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
+ return GSL_SUCCESS;
+ }
+ else {
+ gsl_sf_result poch_rel;
+ int stat_p = pochrel_smallx(a, x, &poch_rel);
+ double eps = x*poch_rel.val;
+ int stat_e = gsl_sf_log_1plusx_e(eps, result);
+ result->err = 2.0 * fabs(x * poch_rel.err / (1.0 + eps));
+ result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
+ return GSL_ERROR_SELECT_2(stat_e, stat_p);
+ }
+}
+
+
+/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
+
+int
+gsl_sf_lnpoch_e(const double a, const double x, gsl_sf_result * result)
+{
+ /* CHECK_POINTER(result) */
+
+ if(a <= 0.0 || a+x <= 0.0) {
+ DOMAIN_ERROR(result);
+ }
+ else if(x == 0.0) {
+ result->val = 0.0;
+ result->err = 0.0;
+ return GSL_SUCCESS;
+ }
+ else {
+ return lnpoch_pos(a, x, result);
+ }
+}
+
+
+int
+gsl_sf_lnpoch_sgn_e(const double a, const double x,
+ gsl_sf_result * result, double * sgn)
+{
+ if(a == 0.0 || a+x == 0.0) {
+ *sgn = 0.0;
+ DOMAIN_ERROR(result);
+ }
+ else if(x == 0.0) {
+ *sgn = 1.0;
+ result->val = 0.0;
+ result->err = 0.0;
+ return GSL_SUCCESS;
+ }
+ else if(a > 0.0 && a+x > 0.0) {
+ *sgn = 1.0;
+ return lnpoch_pos(a, x, result);
+ }
+ else if(a < 0.0 && a+x < 0.0) {
+ /* Reduce to positive case using reflection.
+ */
+ double sin_1 = sin(M_PI * (1.0 - a));
+ double sin_2 = sin(M_PI * (1.0 - a - x));
+ if(sin_1 == 0.0 || sin_2 == 0.0) {
+ *sgn = 0.0;
+ DOMAIN_ERROR(result);
+ }
+ else {
+ gsl_sf_result lnp_pos;
+ int stat_pp = lnpoch_pos(1.0-a, -x, &lnp_pos);
+ double lnterm = log(fabs(sin_1/sin_2));
+ result->val = lnterm - lnp_pos.val;
+ result->err = lnp_pos.err;
+ result->err += 2.0 * GSL_DBL_EPSILON * (fabs(1.0-a) + fabs(1.0-a-x)) * fabs(lnterm);
+ result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
+ *sgn = GSL_SIGN(sin_1*sin_2);
+ return stat_pp;
+ }
+ }
+ else {
+ /* Evaluate gamma ratio directly.
+ */
+ gsl_sf_result lg_apn;
+ gsl_sf_result lg_a;
+ double s_apn, s_a;
+ int stat_apn = gsl_sf_lngamma_sgn_e(a+x, &lg_apn, &s_apn);
+ int stat_a = gsl_sf_lngamma_sgn_e(a, &lg_a, &s_a);
+ if(stat_apn == GSL_SUCCESS && stat_a == GSL_SUCCESS) {
+ result->val = lg_apn.val - lg_a.val;
+ result->err = lg_apn.err + lg_a.err;
+ result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
+ *sgn = s_a * s_apn;
+ return GSL_SUCCESS;
+ }
+ else if(stat_apn == GSL_EDOM || stat_a == GSL_EDOM){
+ *sgn = 0.0;
+ DOMAIN_ERROR(result);
+ }
+ else {
+ result->val = 0.0;
+ result->err = 0.0;
+ *sgn = 0.0;
+ return GSL_FAILURE;
+ }
+ }
+}
+
+
+int
+gsl_sf_poch_e(const double a, const double x, gsl_sf_result * result)
+{
+ /* CHECK_POINTER(result) */
+
+ if(x == 0.0) {
+ result->val = 1.0;
+ result->err = 0.0;
+ return GSL_SUCCESS;
+ }
+ else {
+ gsl_sf_result lnpoch;
+ double sgn;
+ int stat_lnpoch = gsl_sf_lnpoch_sgn_e(a, x, &lnpoch, &sgn);
+ int stat_exp = gsl_sf_exp_err_e(lnpoch.val, lnpoch.err, result);
+ result->val *= sgn;
+ result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
+ return GSL_ERROR_SELECT_2(stat_exp, stat_lnpoch);
+ }
+}
+
+
+int
+gsl_sf_pochrel_e(const double a, const double x, gsl_sf_result * result)
+{
+ const double absx = fabs(x);
+ const double absa = fabs(a);
+
+ /* CHECK_POINTER(result) */
+
+ if(absx > 0.1*absa || absx*log(GSL_MAX(absa,2.0)) > 0.1) {
+ gsl_sf_result lnpoch;
+ double sgn;
+ int stat_poch = gsl_sf_lnpoch_sgn_e(a, x, &lnpoch, &sgn);
+ if(lnpoch.val > GSL_LOG_DBL_MAX) {
+ OVERFLOW_ERROR(result);
+ }
+ else {
+ const double el = exp(lnpoch.val);
+ result->val = (sgn*el - 1.0)/x;
+ result->err = fabs(result->val) * (lnpoch.err + 2.0 * GSL_DBL_EPSILON);
+ result->err += 2.0 * GSL_DBL_EPSILON * (fabs(sgn*el) + 1.0) / fabs(x);
+ return stat_poch;
+ }
+ }
+ else {
+ return pochrel_smallx(a, x, result);
+ }
+}
+
+
+/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
+
+#include "eval.h"
+
+double gsl_sf_lnpoch(const double a, const double x)
+{
+ EVAL_RESULT(gsl_sf_lnpoch_e(a, x, &result));
+}
+
+double gsl_sf_poch(const double a, const double x)
+{
+ EVAL_RESULT(gsl_sf_poch_e(a, x, &result));
+}
+
+double gsl_sf_pochrel(const double a, const double x)
+{
+ EVAL_RESULT(gsl_sf_pochrel_e(a, x, &result));
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-01 12:07:02 +0000