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

   1 /* specfunc/laguerre.c
   2  *
   3  * Copyright (C) 2007 Brian Gough
   4  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU General Public License as published by
   8  * the Free Software Foundation; either version 3 of the License, or (at
   9  * your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful, but
  12  * WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  19  */
  20
  21 /* Author:  G. Jungman */
  22
  23 #include <config.h>
  24 #include <gsl/gsl_math.h>
  25 #include <gsl/gsl_errno.h>
  26 #include <gsl/gsl_sf_exp.h>
  27 #include <gsl/gsl_sf_gamma.h>
  28 #include <gsl/gsl_sf_laguerre.h>
  29
  30 #include "error.h"
  31
  32 /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  33
  34
  35 /* based on the large 2b-4a asymptotic for 1F1
  36  * [Abramowitz+Stegun, 13.5.21]
  37  * L^a_n(x) = (a+1)_n / n! 1F1(-n,a+1,x)
  38  *
  39  * The second term (ser_term2) is from Slater,"The Confluent
  40  * Hypergeometric Function" p.73.  I think there may be an error in
  41  * the first term of the expression given there, comparing with AS
  42  * 13.5.21 (cf sin(a\pi+\Theta) vs sin(a\pi) + sin(\Theta)) - but the
  43  * second term appears correct.
  44  *
  45  */
  46 static
  47 int
  48 laguerre_large_n(const int n, const double alpha, const double x,
  49                  gsl_sf_result * result)
  50 {
  51   const double a = -n;
  52   const double b = alpha + 1.0;
  53   const double eta    = 2.0*b - 4.0*a;
  54   const double cos2th = x/eta;
  55   const double sin2th = 1.0 - cos2th;
  56   const double eps = asin(sqrt(cos2th));  /* theta = pi/2 - eps */
  57   const double pre_h  = 0.25*M_PI*M_PI*eta*eta*cos2th*sin2th;
  58   gsl_sf_result lg_b;
  59   gsl_sf_result lnfact;
  60   int stat_lg = gsl_sf_lngamma_e(b+n, &lg_b);
  61   int stat_lf = gsl_sf_lnfact_e(n, &lnfact);
  62   double pre_term1 = 0.5*(1.0-b)*log(0.25*x*eta);
  63   double pre_term2 = 0.25*log(pre_h);
  64   double lnpre_val = lg_b.val - lnfact.val + 0.5*x + pre_term1 - pre_term2;
  65   double lnpre_err = lg_b.err + lnfact.err + GSL_DBL_EPSILON * (fabs(pre_term1)+fabs(pre_term2));
  66
  67   double phi1 = 0.25*eta*(2*eps + sin(2.0*eps));
  68   double ser_term1 = -sin(phi1);
  69
  70   double A1 = (1.0/12.0)*(5.0/(4.0*sin2th)+(3.0*b*b-6.0*b+2.0)*sin2th - 1.0);
  71   double ser_term2 = -A1 * cos(phi1)/(0.25*eta*sin(2.0*eps));
  72
  73   double ser_val = ser_term1 + ser_term2;
  74   double ser_err = ser_term2*ser_term2 + GSL_DBL_EPSILON * (fabs(ser_term1) + fabs(ser_term2));
  75   int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err, ser_val, ser_err, result);
  76   result->err += 2.0 * GSL_SQRT_DBL_EPSILON * fabs(result->val);
  77   return GSL_ERROR_SELECT_3(stat_e, stat_lf, stat_lg);
  78 }
  79
  80
  81 /* Evaluate polynomial based on confluent hypergeometric representation.
  82  *
  83  * L^a_n(x) = (a+1)_n / n! 1F1(-n,a+1,x)
  84  *
  85  * assumes n > 0 and a != negative integer greater than -n
  86  */
  87 static
  88 int
  89 laguerre_n_cp(const int n, const double a, const double x, gsl_sf_result * result)
  90 {
  91   gsl_sf_result lnfact;
  92   gsl_sf_result lg1;
  93   gsl_sf_result lg2;
  94   double s1, s2;
  95   int stat_f = gsl_sf_lnfact_e(n, &lnfact);
  96   int stat_g1 = gsl_sf_lngamma_sgn_e(a+1.0+n, &lg1, &s1);
  97   int stat_g2 = gsl_sf_lngamma_sgn_e(a+1.0, &lg2, &s2);
  98   double poly_1F1_val = 1.0;
  99   double poly_1F1_err = 0.0;
 100   int stat_e;
 101   int k;
 102
 103   double lnpre_val = (lg1.val - lg2.val) - lnfact.val;
 104   double lnpre_err = lg1.err + lg2.err + lnfact.err + 2.0 * GSL_DBL_EPSILON * fabs(lnpre_val);
 105
 106   for(k=n-1; k>=0; k--) {
 107     double t = (-n+k)/(a+1.0+k) * (x/(k+1));
 108     double r = t + 1.0/poly_1F1_val;
 109     if(r > 0.9*GSL_DBL_MAX/poly_1F1_val) {
 110       /* internal error only, don't call the error handler */
 111       INTERNAL_OVERFLOW_ERROR(result);
 112     }
 113     else {
 114       /* Collect the Horner terms. */
 115       poly_1F1_val  = 1.0 + t * poly_1F1_val;
 116       poly_1F1_err += GSL_DBL_EPSILON + fabs(t) * poly_1F1_err;
 117     }
 118   }
 119
 120   stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
 121                                     poly_1F1_val, poly_1F1_err,
 122                                     result);
 123
 124   return GSL_ERROR_SELECT_4(stat_e, stat_f, stat_g1, stat_g2);
 125 }
 126
 127
 128 /* Evaluate the polynomial based on the confluent hypergeometric
 129  * function in a safe way, with no restriction on the arguments.
 130  *
 131  * assumes x != 0
 132  */
 133 static
 134 int
 135 laguerre_n_poly_safe(const int n, const double a, const double x, gsl_sf_result * result)
 136 {
 137   const double b = a + 1.0;
 138   const double mx = -x;
 139   const double tc_sgn = (x < 0.0 ? 1.0 : (GSL_IS_ODD(n) ? -1.0 : 1.0));
 140   gsl_sf_result tc;
 141   int stat_tc = gsl_sf_taylorcoeff_e(n, fabs(x), &tc);
 142
 143   if(stat_tc == GSL_SUCCESS) {
 144     double term = tc.val * tc_sgn;
 145     double sum_val = term;
 146     double sum_err = tc.err;
 147     int k;
 148     for(k=n-1; k>=0; k--) {
 149       term *= ((b+k)/(n-k))*(k+1.0)/mx;
 150       sum_val += term;
 151       sum_err += 4.0 * GSL_DBL_EPSILON * fabs(term);
 152     }
 153     result->val = sum_val;
 154     result->err = sum_err + 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 155     return GSL_SUCCESS;
 156   }
 157   else if(stat_tc == GSL_EOVRFLW) {
 158     result->val = 0.0; /* FIXME: should be Inf */
 159     result->err = 0.0;
 160     return stat_tc;
 161   }
 162   else {
 163     result->val = 0.0;
 164     result->err = 0.0;
 165     return stat_tc;
 166   }
 167 }
 168
 169
 170
 171 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*/
 172
 173 int
 174 gsl_sf_laguerre_1_e(const double a, const double x, gsl_sf_result * result)
 175 {
 176   /* CHECK_POINTER(result) */
 177
 178   {
 179     result->val = 1.0 + a - x;
 180     result->err = 2.0 * GSL_DBL_EPSILON * (1.0 + fabs(a) + fabs(x));
 181     return GSL_SUCCESS;
 182   }
 183 }
 184
 185 int
 186 gsl_sf_laguerre_2_e(const double a, const double x, gsl_sf_result * result)
 187 {
 188   /* CHECK_POINTER(result) */
 189
 190   if(a == -2.0) {
 191     result->val = 0.5*x*x;
 192     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 193     return GSL_SUCCESS;
 194   }
 195   else {
 196     double c0 = 0.5 * (2.0+a)*(1.0+a);
 197     double c1 = -(2.0+a);
 198     double c2 = -0.5/(2.0+a);
 199     result->val  = c0 + c1*x*(1.0 + c2*x);
 200     result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(c0) + 2.0 * fabs(c1*x) * (1.0 + 2.0 * fabs(c2*x)));
 201     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 202     return GSL_SUCCESS;
 203   }
 204 }
 205
 206 int
 207 gsl_sf_laguerre_3_e(const double a, const double x, gsl_sf_result * result)
 208 {
 209   /* CHECK_POINTER(result) */
 210
 211   if(a == -2.0) {
 212     double x2_6  = x*x/6.0;
 213     result->val  = x2_6 * (3.0 - x);
 214     result->err  = x2_6 * (3.0 + fabs(x)) * 2.0 * GSL_DBL_EPSILON;
 215     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 216     return GSL_SUCCESS;
 217   }
 218   else if(a == -3.0) {
 219     result->val = -x*x/6.0;
 220     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 221     return GSL_SUCCESS;
 222   }
 223   else {
 224     double c0 = (3.0+a)*(2.0+a)*(1.0+a) / 6.0;
 225     double c1 = -c0 * 3.0 / (1.0+a);
 226     double c2 = -1.0/(2.0+a);
 227     double c3 = -1.0/(3.0*(3.0+a));
 228     result->val  = c0 + c1*x*(1.0 + c2*x*(1.0 + c3*x));
 229     result->err  = 1.0 + 2.0 * fabs(c3*x);
 230     result->err  = 1.0 + 2.0 * fabs(c2*x) * result->err;
 231     result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(c0) + 2.0 * fabs(c1*x) * result->err);
 232     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 233     return GSL_SUCCESS;
 234   }
 235 }
 236
 237
 238 int gsl_sf_laguerre_n_e(const int n, const double a, const double x,
 239                            gsl_sf_result * result)
 240 {
 241   /* CHECK_POINTER(result) */
 242
 243   if(n < 0) {
 244     DOMAIN_ERROR(result);
 245   }
 246   else if(n == 0) {
 247     result->val = 1.0;
 248     result->err = 0.0;
 249     return GSL_SUCCESS;
 250   }
 251   else if(n == 1) {
 252     result->val = 1.0 + a - x;
 253     result->err = 2.0 * GSL_DBL_EPSILON * (1.0 + fabs(a) + fabs(x));
 254     return GSL_SUCCESS;
 255   }
 256   else if(x == 0.0) {
 257     double product = a + 1.0;
 258     int k;
 259     for(k=2; k<=n; k++) {
 260       product *= (a + k)/k;
 261     }
 262     result->val = product;
 263     result->err = 2.0 * (n + 1.0) * GSL_DBL_EPSILON * fabs(product) + GSL_DBL_EPSILON;
 264     return GSL_SUCCESS;
 265   }
 266   else if(x < 0.0 && a > -1.0) {
 267     /* In this case all the terms in the polynomial
 268      * are of the same sign. Note that this also
 269      * catches overflows correctly.
 270      */
 271     return laguerre_n_cp(n, a, x, result);
 272   }
 273   else if(n < 5 || (x > 0.0 && a < -n-1)) {
 274     /* Either the polynomial will not lose too much accuracy
 275      * or all the terms are negative. In any case,
 276      * the error estimate here is good. We try both
 277      * explicit summation methods, as they have different
 278      * characteristics. One may underflow/overflow while the
 279      * other does not.
 280      */
 281     if(laguerre_n_cp(n, a, x, result) == GSL_SUCCESS)
 282       return GSL_SUCCESS;
 283     else
 284       return laguerre_n_poly_safe(n, a, x, result);
 285   }
 286   else if(n > 1.0e+07 && x > 0.0 && a > -1.0 && x < 2.0*(a+1.0)+4.0*n) {
 287     return laguerre_large_n(n, a, x, result);
 288   }
 289   else if(a >= 0.0 || (x > 0.0 && a < -n-1)) {
 290     gsl_sf_result lg2;
 291     int stat_lg2 = gsl_sf_laguerre_2_e(a, x, &lg2);
 292     double Lkm1 = 1.0 + a - x;
 293     double Lk   = lg2.val;
 294     double Lkp1;
 295     int k;
 296
 297     for(k=2; k<n; k++) {
 298       Lkp1 = (-(k+a)*Lkm1 + (2.0*k+a+1.0-x)*Lk)/(k+1.0);
 299       Lkm1 = Lk;
 300       Lk   = Lkp1;
 301     }
 302     result->val = Lk;
 303     result->err = (fabs(lg2.err/lg2.val) + GSL_DBL_EPSILON) * n * fabs(Lk);
 304     return stat_lg2;
 305   }
 306   else {
 307     /* Despair... or magic? */
 308     return laguerre_n_poly_safe(n, a, x, result);
 309   }
 310 }
 311
 312
 313 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 314
 315 #include "eval.h"
 316
 317 double gsl_sf_laguerre_1(double a, double x)
 318 {
 319   EVAL_RESULT(gsl_sf_laguerre_1_e(a, x, &result));
 320 }
 321
 322 double gsl_sf_laguerre_2(double a, double x)
 323 {
 324   EVAL_RESULT(gsl_sf_laguerre_2_e(a, x, &result));
 325 }
 326
 327 double gsl_sf_laguerre_3(double a, double x)
 328 {
 329   EVAL_RESULT(gsl_sf_laguerre_3_e(a, x, &result));
 330 }
 331
 332 double gsl_sf_laguerre_n(int n, double a, double x)
 333 {
 334   EVAL_RESULT(gsl_sf_laguerre_n_e(n, a, x, &result));
 335 }