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

   1 /* specfunc/sinint.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_trig.h>
  26 #include <gsl/gsl_sf_expint.h>
  27
  28 #include "error.h"
  29
  30 #include "chebyshev.h"
  31 #include "cheb_eval.c"
  32
  33 /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  34
  35 /* based on SLATEC r9sifg.f, W. Fullerton */
  36
  37 /*
  38  series for f1   on the interval  2.00000e-02 to  6.25000e-02
  39                                         with weighted error   2.82e-17
  40                                          log weighted error  16.55
  41                                significant figures required  15.36
  42                                     decimal places required  17.20
  43 */
  44 static double f1_data[20] = {
  45    -0.1191081969051363610,
  46    -0.0247823144996236248,
  47     0.0011910281453357821,
  48    -0.0000927027714388562,
  49     0.0000093373141568271,
  50    -0.0000011058287820557,
  51     0.0000001464772071460,
  52    -0.0000000210694496288,
  53     0.0000000032293492367,
  54    -0.0000000005206529618,
  55     0.0000000000874878885,
  56    -0.0000000000152176187,
  57     0.0000000000027257192,
  58    -0.0000000000005007053,
  59     0.0000000000000940241,
  60    -0.0000000000000180014,
  61     0.0000000000000035063,
  62    -0.0000000000000006935,
  63     0.0000000000000001391,
  64    -0.0000000000000000282
  65 };
  66 static cheb_series f1_cs = {
  67   f1_data,
  68   19,
  69   -1, 1,
  70   10
  71 };
  72
  73 /*
  74
  75  series for f2   on the interval  0.00000e+00 to  2.00000e-02
  76                                         with weighted error   4.32e-17
  77                                          log weighted error  16.36
  78                                significant figures required  14.75
  79                                     decimal places required  17.10
  80 */
  81 static double f2_data[29] = {
  82    -0.0348409253897013234,
  83    -0.0166842205677959686,
  84     0.0006752901241237738,
  85    -0.0000535066622544701,
  86     0.0000062693421779007,
  87    -0.0000009526638801991,
  88     0.0000001745629224251,
  89    -0.0000000368795403065,
  90     0.0000000087202677705,
  91    -0.0000000022601970392,
  92     0.0000000006324624977,
  93    -0.0000000001888911889,
  94     0.0000000000596774674,
  95    -0.0000000000198044313,
  96     0.0000000000068641396,
  97    -0.0000000000024731020,
  98     0.0000000000009226360,
  99    -0.0000000000003552364,
 100     0.0000000000001407606,
 101    -0.0000000000000572623,
 102     0.0000000000000238654,
 103    -0.0000000000000101714,
 104     0.0000000000000044259,
 105    -0.0000000000000019634,
 106     0.0000000000000008868,
 107    -0.0000000000000004074,
 108     0.0000000000000001901,
 109    -0.0000000000000000900,
 110     0.0000000000000000432
 111 };
 112 static cheb_series f2_cs = {
 113   f2_data,
 114   28,
 115   -1, 1,
 116   14
 117 };
 118
 119 /*
 120
 121  series for g1   on the interval  2.00000e-02 to  6.25000e-02
 122                                         with weighted error   5.48e-17
 123                                          log weighted error  16.26
 124                                significant figures required  15.47
 125                                     decimal places required  16.92
 126 */
 127 static double g1_data[21] = {
 128    -0.3040578798253495954,
 129    -0.0566890984597120588,
 130     0.0039046158173275644,
 131    -0.0003746075959202261,
 132     0.0000435431556559844,
 133    -0.0000057417294453025,
 134     0.0000008282552104503,
 135    -0.0000001278245892595,
 136     0.0000000207978352949,
 137    -0.0000000035313205922,
 138     0.0000000006210824236,
 139    -0.0000000001125215474,
 140     0.0000000000209088918,
 141    -0.0000000000039715832,
 142     0.0000000000007690431,
 143    -0.0000000000001514697,
 144     0.0000000000000302892,
 145    -0.0000000000000061400,
 146     0.0000000000000012601,
 147    -0.0000000000000002615,
 148     0.0000000000000000548
 149 };
 150 static cheb_series g1_cs = {
 151   g1_data,
 152   20,
 153   -1, 1,
 154   13
 155 };
 156
 157 /*
 158
 159  series for g2   on the interval  0.00000e+00 to  2.00000e-02
 160                                         with weighted error   5.01e-17
 161                                          log weighted error  16.30
 162                                significant figures required  15.12
 163                                     decimal places required  17.07
 164 */
 165 static double g2_data[34] = {
 166    -0.0967329367532432218,
 167    -0.0452077907957459871,
 168     0.0028190005352706523,
 169    -0.0002899167740759160,
 170     0.0000407444664601121,
 171    -0.0000071056382192354,
 172     0.0000014534723163019,
 173    -0.0000003364116512503,
 174     0.0000000859774367886,
 175    -0.0000000238437656302,
 176     0.0000000070831906340,
 177    -0.0000000022318068154,
 178     0.0000000007401087359,
 179    -0.0000000002567171162,
 180     0.0000000000926707021,
 181    -0.0000000000346693311,
 182     0.0000000000133950573,
 183    -0.0000000000053290754,
 184     0.0000000000021775312,
 185    -0.0000000000009118621,
 186     0.0000000000003905864,
 187    -0.0000000000001708459,
 188     0.0000000000000762015,
 189    -0.0000000000000346151,
 190     0.0000000000000159996,
 191    -0.0000000000000075213,
 192     0.0000000000000035970,
 193    -0.0000000000000017530,
 194     0.0000000000000008738,
 195    -0.0000000000000004487,
 196     0.0000000000000002397,
 197    -0.0000000000000001347,
 198     0.0000000000000000801,
 199    -0.0000000000000000501
 200 };
 201 static cheb_series g2_cs = {
 202   g2_data,
 203   33,
 204   -1, 1,
 205   20
 206 };
 207
 208
 209 /* x >= 4.0 */
 210 static void fg_asymp(const double x, gsl_sf_result * f, gsl_sf_result * g)
 211 {
 212   /*
 213       xbig = sqrt (1.0/r1mach(3))
 214       xmaxf = exp (amin1(-alog(r1mach(1)), alog(r1mach(2))) - 0.01)
 215       xmaxg = 1.0/sqrt(r1mach(1))
 216       xbnd = sqrt(50.0)
 217   */
 218   const double xbig  = 1.0/GSL_SQRT_DBL_EPSILON;
 219   const double xmaxf = 1.0/GSL_DBL_MIN;
 220   const double xmaxg = 1.0/GSL_SQRT_DBL_MIN;
 221   const double xbnd  = 7.07106781187;
 222
 223   const double x2 = x*x;
 224
 225   if(x <= xbnd) {
 226     gsl_sf_result result_c1;
 227     gsl_sf_result result_c2;
 228     cheb_eval_e(&f1_cs, (1.0/x2-0.04125)/0.02125, &result_c1);
 229     cheb_eval_e(&g1_cs, (1.0/x2-0.04125)/0.02125, &result_c2);
 230     f->val = (1.0 + result_c1.val)/x;
 231     g->val = (1.0 + result_c2.val)/x2;
 232     f->err = result_c1.err/x  + 2.0 * GSL_DBL_EPSILON * fabs(f->val);
 233     g->err = result_c2.err/x2 + 2.0 * GSL_DBL_EPSILON * fabs(g->val);
 234   }
 235   else if(x <= xbig) {
 236     gsl_sf_result result_c1;
 237     gsl_sf_result result_c2;
 238     cheb_eval_e(&f2_cs, 100.0/x2-1.0, &result_c1);
 239     cheb_eval_e(&g2_cs, 100.0/x2-1.0, &result_c2);
 240     f->val = (1.0 + result_c1.val)/x;
 241     g->val = (1.0 + result_c2.val)/x2;
 242     f->err = result_c1.err/x  + 2.0 * GSL_DBL_EPSILON * fabs(f->val);
 243     g->err = result_c2.err/x2 + 2.0 * GSL_DBL_EPSILON * fabs(g->val);
 244   }
 245   else {
 246     f->val = (x < xmaxf ? 1.0/x  : 0.0);
 247     g->val = (x < xmaxg ? 1.0/x2 : 0.0);
 248     f->err = 2.0 * GSL_DBL_EPSILON * fabs(f->val);
 249     g->err = 2.0 * GSL_DBL_EPSILON * fabs(g->val);
 250   }
 251
 252   return;
 253 }
 254
 255
 256 /* based on SLATEC si.f, W. Fullerton
 257
 258  series for si   on the interval  0.00000e+00 to  1.60000e+01
 259                                         with weighted error   1.22e-17
 260                                          log weighted error  16.91
 261                                significant figures required  16.37
 262                                     decimal places required  17.45
 263 */
 264
 265 static double si_data[12] = {
 266   -0.1315646598184841929,
 267   -0.2776578526973601892,
 268    0.0354414054866659180,
 269   -0.0025631631447933978,
 270    0.0001162365390497009,
 271   -0.0000035904327241606,
 272    0.0000000802342123706,
 273   -0.0000000013562997693,
 274    0.0000000000179440722,
 275   -0.0000000000001908387,
 276    0.0000000000000016670,
 277   -0.0000000000000000122
 278 };
 279
 280 static cheb_series si_cs = {
 281   si_data,
 282   11,
 283   -1, 1,
 284   9
 285 };
 286
 287 /*
 288  series for ci   on the interval  0.00000e+00 to  1.60000e+01
 289                                         with weighted error   1.94e-18
 290                                          log weighted error  17.71
 291                                significant figures required  17.74
 292                                     decimal places required  18.27
 293 */
 294 static double ci_data[13] = {
 295    -0.34004281856055363156,
 296    -1.03302166401177456807,
 297     0.19388222659917082877,
 298    -0.01918260436019865894,
 299     0.00110789252584784967,
 300    -0.00004157234558247209,
 301     0.00000109278524300229,
 302    -0.00000002123285954183,
 303     0.00000000031733482164,
 304    -0.00000000000376141548,
 305     0.00000000000003622653,
 306    -0.00000000000000028912,
 307     0.00000000000000000194
 308 };
 309 static cheb_series ci_cs = {
 310   ci_data,
 311   12,
 312   -1, 1,
 313   9
 314 };
 315
 316
 317 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 318
 319 int gsl_sf_Si_e(const double x, gsl_sf_result * result)
 320 {
 321   double ax = fabs(x);
 322
 323   /* CHECK_POINTER(result) */
 324
 325   if(ax < GSL_SQRT_DBL_EPSILON) {
 326     result->val = x;
 327     result->err = 0.0;
 328     return GSL_SUCCESS;
 329   }
 330   else if(ax <= 4.0) {
 331     gsl_sf_result result_c;
 332     cheb_eval_e(&si_cs, (x*x-8.0)*0.125, &result_c);
 333     result->val  =  x * (0.75 + result_c.val);
 334     result->err  = ax * result_c.err;
 335     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 336     return GSL_SUCCESS;
 337   }
 338   else {
 339     /* Note there is no loss of precision
 340      * here bcause of the leading constant.
 341      */
 342     gsl_sf_result f;
 343     gsl_sf_result g;
 344     fg_asymp(ax, &f, &g);
 345     result->val  = 0.5 * M_PI - f.val*cos(ax) - g.val*sin(ax);
 346     result->err  = f.err + g.err;
 347     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 348     if(x < 0.0) result->val = -result->val;
 349     return GSL_SUCCESS;
 350   }
 351 }
 352
 353
 354 int gsl_sf_Ci_e(const double x, gsl_sf_result * result)
 355 {
 356   /* CHECK_POINTER(result) */
 357
 358   if(x <= 0.0) {
 359     DOMAIN_ERROR(result);
 360   }
 361   else if(x <= 4.0) {
 362     const double lx = log(x);
 363     const double y  = (x*x-8.0)*0.125;
 364     gsl_sf_result result_c;
 365     cheb_eval_e(&ci_cs, y, &result_c);
 366     result->val  = lx - 0.5 + result_c.val;
 367     result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(lx) + 0.5) + result_c.err;
 368     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 369     return GSL_SUCCESS;
 370   }
 371   else {
 372     gsl_sf_result sin_result;
 373     gsl_sf_result cos_result;
 374     int stat_sin = gsl_sf_sin_e(x, &sin_result);
 375     int stat_cos = gsl_sf_cos_e(x, &cos_result);
 376     gsl_sf_result f;
 377     gsl_sf_result g;
 378     fg_asymp(x, &f, &g);
 379     result->val  = f.val*sin_result.val - g.val*cos_result.val;
 380     result->err  = fabs(f.err*sin_result.val);
 381     result->err += fabs(g.err*cos_result.val);
 382     result->err += fabs(f.val*sin_result.err);
 383     result->err += fabs(g.val*cos_result.err);
 384     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 385     return GSL_ERROR_SELECT_2(stat_sin, stat_cos);
 386   }
 387 }
 388
 389
 390 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 391
 392 #include "eval.h"
 393
 394 double gsl_sf_Si(const double x)
 395 {
 396   EVAL_RESULT(gsl_sf_Si_e(x, &result));
 397 }
 398
 399 double gsl_sf_Ci(const double x)
 400 {
 401   EVAL_RESULT(gsl_sf_Ci_e(x, &result));
 402 }