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

   1 /* specfunc/shint.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_expint.h>
  26
  27 #include "error.h"
  28
  29 #include "chebyshev.h"
  30 #include "cheb_eval.c"
  31
  32 /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  33
  34 /* based on SLATEC shi.f, W. Fullerton
  35
  36  series for shi  on the interval  0.00000e+00 to  1.40625e-01
  37                                         with weighted error   4.67e-20
  38                                          log weighted error  19.33
  39                                significant figures required  17.07
  40                                     decimal places required  19.75
  41 */
  42 static double shi_data[7] = {
  43    0.0078372685688900950695,
  44    0.0039227664934234563973,
  45    0.0000041346787887617267,
  46    0.0000000024707480372883,
  47    0.0000000000009379295591,
  48    0.0000000000000002451817,
  49    0.0000000000000000000467
  50 };
  51 static cheb_series shi_cs = {
  52   shi_data,
  53   6,
  54   -1, 1,
  55   6
  56 };
  57
  58
  59 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  60
  61 int gsl_sf_Shi_e(const double x, gsl_sf_result * result)
  62 {
  63   const double xsml = GSL_SQRT_DBL_EPSILON;  /* sqrt (r1mach(3)) */
  64   const double ax   = fabs(x);
  65
  66   if(ax < xsml) {
  67     result->val = x;
  68     result->err = 0.0;
  69     return GSL_SUCCESS;
  70   }
  71   else if(ax <= 0.375) {
  72     gsl_sf_result result_c;
  73     cheb_eval_e(&shi_cs, 128.0*x*x/9.0-1.0, &result_c);
  74     result->val  = x * (1.0 + result_c.val);
  75     result->err  = x * result_c.err;
  76     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  77     return GSL_SUCCESS;
  78   }
  79   else {
  80     gsl_sf_result result_Ei;
  81     gsl_sf_result result_E1;
  82     int status_Ei = gsl_sf_expint_Ei_e(x, &result_Ei);
  83     int status_E1 = gsl_sf_expint_E1_e(x, &result_E1);
  84     result->val  = 0.5*(result_Ei.val + result_E1.val);
  85     result->err  = 0.5*(result_Ei.err + result_E1.err);
  86     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  87     if(status_Ei == GSL_EUNDRFLW && status_E1 == GSL_EUNDRFLW) {
  88       GSL_ERROR ("underflow", GSL_EUNDRFLW);
  89     }
  90     else if(status_Ei == GSL_EOVRFLW || status_E1 == GSL_EOVRFLW) {
  91       GSL_ERROR ("overflow", GSL_EOVRFLW);
  92     }
  93     else {
  94       return GSL_SUCCESS;
  95     }
  96   }
  97 }
  98
  99
 100 int gsl_sf_Chi_e(const double x, gsl_sf_result * result)
 101 {
 102   gsl_sf_result result_Ei;
 103   gsl_sf_result result_E1;
 104   int status_Ei = gsl_sf_expint_Ei_e(x, &result_Ei);
 105   int status_E1 = gsl_sf_expint_E1_e(x, &result_E1);
 106   if(status_Ei == GSL_EDOM || status_E1 == GSL_EDOM) {
 107     DOMAIN_ERROR(result);
 108   }
 109   else if(status_Ei == GSL_EUNDRFLW && status_E1 == GSL_EUNDRFLW) {
 110     UNDERFLOW_ERROR(result);
 111   }
 112   else if(status_Ei == GSL_EOVRFLW || status_E1 == GSL_EOVRFLW) {
 113     OVERFLOW_ERROR(result);
 114   }
 115   else {
 116     result->val  = 0.5 * (result_Ei.val - result_E1.val);
 117     result->err  = 0.5 * (result_Ei.err + result_E1.err);
 118     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 119     return GSL_SUCCESS;
 120   }
 121 }
 122
 123 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 124
 125 #include "eval.h"
 126
 127 double gsl_sf_Shi(const double x)
 128 {
 129   EVAL_RESULT(gsl_sf_Shi_e(x, &result));
 130 }
 131
 132 double gsl_sf_Chi(const double x)
 133 {
 134   EVAL_RESULT(gsl_sf_Chi_e(x, &result));
 135 }