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

   1 /* specfunc/log.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_log.h>
  26
  27 #include "error.h"
  28
  29 #include "chebyshev.h"
  30 #include "cheb_eval.c"
  31
  32 /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  33
  34 /* Chebyshev expansion for log(1 + x(t))/x(t)
  35  *
  36  * x(t) = (4t-1)/(2(4-t))
  37  * t(x) = (8x+1)/(2(x+2))
  38  * -1/2 < x < 1/2
  39  * -1 < t < 1
  40  */
  41 static double lopx_data[21] = {
  42   2.16647910664395270521272590407,
  43  -0.28565398551049742084877469679,
  44   0.01517767255690553732382488171,
  45  -0.00200215904941415466274422081,
  46   0.00019211375164056698287947962,
  47  -0.00002553258886105542567601400,
  48   2.9004512660400621301999384544e-06,
  49  -3.8873813517057343800270917900e-07,
  50   4.7743678729400456026672697926e-08,
  51  -6.4501969776090319441714445454e-09,
  52   8.2751976628812389601561347296e-10,
  53  -1.1260499376492049411710290413e-10,
  54   1.4844576692270934446023686322e-11,
  55  -2.0328515972462118942821556033e-12,
  56   2.7291231220549214896095654769e-13,
  57  -3.7581977830387938294437434651e-14,
  58   5.1107345870861673561462339876e-15,
  59  -7.0722150011433276578323272272e-16,
  60   9.7089758328248469219003866867e-17,
  61  -1.3492637457521938883731579510e-17,
  62   1.8657327910677296608121390705e-18
  63 };
  64 static cheb_series lopx_cs = {
  65   lopx_data,
  66   20,
  67   -1, 1,
  68   10
  69 };
  70
  71 /* Chebyshev expansion for (log(1 + x(t)) - x(t))/x(t)^2
  72  *
  73  * x(t) = (4t-1)/(2(4-t))
  74  * t(x) = (8x+1)/(2(x+2))
  75  * -1/2 < x < 1/2
  76  * -1 < t < 1
  77  */
  78 static double lopxmx_data[20] = {
  79  -1.12100231323744103373737274541,
  80   0.19553462773379386241549597019,
  81  -0.01467470453808083971825344956,
  82   0.00166678250474365477643629067,
  83  -0.00018543356147700369785746902,
  84   0.00002280154021771635036301071,
  85  -2.8031253116633521699214134172e-06,
  86   3.5936568872522162983669541401e-07,
  87  -4.6241857041062060284381167925e-08,
  88   6.0822637459403991012451054971e-09,
  89  -8.0339824424815790302621320732e-10,
  90   1.0751718277499375044851551587e-10,
  91  -1.4445310914224613448759230882e-11,
  92   1.9573912180610336168921438426e-12,
  93  -2.6614436796793061741564104510e-13,
  94   3.6402634315269586532158344584e-14,
  95  -4.9937495922755006545809120531e-15,
  96   6.8802890218846809524646902703e-16,
  97  -9.5034129794804273611403251480e-17,
  98   1.3170135013050997157326965813e-17
  99 };
 100 static cheb_series lopxmx_cs = {
 101   lopxmx_data,
 102   19,
 103   -1, 1,
 104   9
 105 };
 106
 107
 108 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 109
 110 int
 111 gsl_sf_log_e(const double x, gsl_sf_result * result)
 112 {
 113   /* CHECK_POINTER(result) */
 114
 115   if(x <= 0.0) {
 116     DOMAIN_ERROR(result);
 117   }
 118   else {
 119     result->val = log(x);
 120     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 121     return GSL_SUCCESS;
 122   }
 123 }
 124
 125
 126 int
 127 gsl_sf_log_abs_e(const double x, gsl_sf_result * result)
 128 {
 129   /* CHECK_POINTER(result) */
 130
 131   if(x == 0.0) {
 132     DOMAIN_ERROR(result);
 133   }
 134   else {
 135     result->val = log(fabs(x));
 136     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 137     return GSL_SUCCESS;
 138   }
 139 }
 140
 141 int
 142 gsl_sf_complex_log_e(const double zr, const double zi, gsl_sf_result * lnr, gsl_sf_result * theta)
 143 {
 144   /* CHECK_POINTER(lnr) */
 145   /* CHECK_POINTER(theta) */
 146
 147   if(zr != 0.0 || zi != 0.0) {
 148     const double ax = fabs(zr);
 149     const double ay = fabs(zi);
 150     const double min = GSL_MIN(ax, ay);
 151     const double max = GSL_MAX(ax, ay);
 152     lnr->val = log(max) + 0.5 * log(1.0 + (min/max)*(min/max));
 153     lnr->err = 2.0 * GSL_DBL_EPSILON * fabs(lnr->val);
 154     theta->val = atan2(zi, zr);
 155     theta->err = GSL_DBL_EPSILON * fabs(lnr->val);
 156     return GSL_SUCCESS;
 157   }
 158   else {
 159     DOMAIN_ERROR_2(lnr, theta);
 160   }
 161 }
 162
 163
 164 int
 165 gsl_sf_log_1plusx_e(const double x, gsl_sf_result * result)
 166 {
 167   /* CHECK_POINTER(result) */
 168
 169   if(x <= -1.0) {
 170     DOMAIN_ERROR(result);
 171   }
 172   else if(fabs(x) < GSL_ROOT6_DBL_EPSILON) {
 173     const double c1 = -0.5;
 174     const double c2 =  1.0/3.0;
 175     const double c3 = -1.0/4.0;
 176     const double c4 =  1.0/5.0;
 177     const double c5 = -1.0/6.0;
 178     const double c6 =  1.0/7.0;
 179     const double c7 = -1.0/8.0;
 180     const double c8 =  1.0/9.0;
 181     const double c9 = -1.0/10.0;
 182     const double t  =  c5 + x*(c6 + x*(c7 + x*(c8 + x*c9)));
 183     result->val = x * (1.0 + x*(c1 + x*(c2 + x*(c3 + x*(c4 + x*t)))));
 184     result->err = GSL_DBL_EPSILON * fabs(result->val);
 185     return GSL_SUCCESS;
 186   }
 187   else if(fabs(x) < 0.5) {
 188     double t = 0.5*(8.0*x + 1.0)/(x+2.0);
 189     gsl_sf_result c;
 190     cheb_eval_e(&lopx_cs, t, &c);
 191     result->val = x * c.val;
 192     result->err = fabs(x * c.err);
 193     return GSL_SUCCESS;
 194   }
 195   else {
 196     result->val = log(1.0 + x);
 197     result->err = GSL_DBL_EPSILON * fabs(result->val);
 198     return GSL_SUCCESS;
 199   }
 200 }
 201
 202
 203 int
 204 gsl_sf_log_1plusx_mx_e(const double x, gsl_sf_result * result)
 205 {
 206   /* CHECK_POINTER(result) */
 207
 208   if(x <= -1.0) {
 209     DOMAIN_ERROR(result);
 210   }
 211   else if(fabs(x) < GSL_ROOT5_DBL_EPSILON) {
 212     const double c1 = -0.5;
 213     const double c2 =  1.0/3.0;
 214     const double c3 = -1.0/4.0;
 215     const double c4 =  1.0/5.0;
 216     const double c5 = -1.0/6.0;
 217     const double c6 =  1.0/7.0;
 218     const double c7 = -1.0/8.0;
 219     const double c8 =  1.0/9.0;
 220     const double c9 = -1.0/10.0;
 221     const double t  =  c5 + x*(c6 + x*(c7 + x*(c8 + x*c9)));
 222     result->val = x*x * (c1 + x*(c2 + x*(c3 + x*(c4 + x*t))));
 223     result->err = GSL_DBL_EPSILON * fabs(result->val);
 224     return GSL_SUCCESS;
 225   }
 226   else if(fabs(x) < 0.5) {
 227     double t = 0.5*(8.0*x + 1.0)/(x+2.0);
 228     gsl_sf_result c;
 229     cheb_eval_e(&lopxmx_cs, t, &c);
 230     result->val = x*x * c.val;
 231     result->err = x*x * c.err;
 232     return GSL_SUCCESS;
 233   }
 234   else {
 235     const double lterm = log(1.0 + x);
 236     result->val = lterm - x;
 237     result->err = GSL_DBL_EPSILON * (fabs(lterm) + fabs(x));
 238     return GSL_SUCCESS;
 239   }
 240 }
 241
 242
 243
 244 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 245
 246 #include "eval.h"
 247
 248 double gsl_sf_log(const double x)
 249 {
 250   EVAL_RESULT(gsl_sf_log_e(x, &result));
 251 }
 252
 253 double gsl_sf_log_abs(const double x)
 254 {
 255   EVAL_RESULT(gsl_sf_log_abs_e(x, &result));
 256 }
 257
 258 double gsl_sf_log_1plusx(const double x)
 259 {
 260   EVAL_RESULT(gsl_sf_log_1plusx_e(x, &result));
 261 }
 262
 263 double gsl_sf_log_1plusx_mx(const double x)
 264 {
 265   EVAL_RESULT(gsl_sf_log_1plusx_mx_e(x, &result));
 266 }