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

   1 /* specfunc/fermi_dirac.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_exp.h>
  26 #include <gsl/gsl_sf_gamma.h>
  27 #include <gsl/gsl_sf_hyperg.h>
  28 #include <gsl/gsl_sf_pow_int.h>
  29 #include <gsl/gsl_sf_zeta.h>
  30 #include <gsl/gsl_sf_fermi_dirac.h>
  31
  32 #include "error.h"
  33
  34 #include "chebyshev.h"
  35 #include "cheb_eval.c"
  36
  37 #define locEPS  (1000.0*GSL_DBL_EPSILON)
  38
  39
  40 /* Chebyshev fit for F_{1}(t);  -1 < t < 1, -1 < x < 1
  41  */
  42 static double fd_1_a_data[22] = {
  43   1.8949340668482264365,
  44   0.7237719066890052793,
  45   0.1250000000000000000,
  46   0.0101065196435973942,
  47   0.0,
  48  -0.0000600615242174119,
  49   0.0,
  50   6.816528764623e-7,
  51   0.0,
  52  -9.5895779195e-9,
  53   0.0,
  54   1.515104135e-10,
  55   0.0,
  56  -2.5785616e-12,
  57   0.0,
  58   4.62270e-14,
  59   0.0,
  60  -8.612e-16,
  61   0.0,
  62   1.65e-17,
  63   0.0,
  64  -3.e-19
  65 };
  66 static cheb_series fd_1_a_cs = {
  67   fd_1_a_data,
  68   21,
  69   -1, 1,
  70   12
  71 };
  72
  73
  74 /* Chebyshev fit for F_{1}(3/2(t+1) + 1);  -1 < t < 1, 1 < x < 4
  75  */
  76 static double fd_1_b_data[22] = {
  77   10.409136795234611872,
  78   3.899445098225161947,
  79   0.513510935510521222,
  80   0.010618736770218426,
  81  -0.001584468020659694,
  82   0.000146139297161640,
  83  -1.408095734499e-6,
  84  -2.177993899484e-6,
  85   3.91423660640e-7,
  86  -2.3860262660e-8,
  87  -4.138309573e-9,
  88   1.283965236e-9,
  89  -1.39695990e-10,
  90  -4.907743e-12,
  91   4.399878e-12,
  92  -7.17291e-13,
  93   2.4320e-14,
  94   1.4230e-14,
  95  -3.446e-15,
  96   2.93e-16,
  97   3.7e-17,
  98  -1.6e-17
  99 };
 100 static cheb_series fd_1_b_cs = {
 101   fd_1_b_data,
 102   21,
 103   -1, 1,
 104   11
 105 };
 106
 107
 108 /* Chebyshev fit for F_{1}(3(t+1) + 4);  -1 < t < 1, 4 < x < 10
 109  */
 110 static double fd_1_c_data[23] = {
 111   56.78099449124299762,
 112   21.00718468237668011,
 113   2.24592457063193457,
 114   0.00173793640425994,
 115  -0.00058716468739423,
 116   0.00016306958492437,
 117  -0.00003817425583020,
 118   7.64527252009e-6,
 119  -1.31348500162e-6,
 120   1.9000646056e-7,
 121  -2.141328223e-8,
 122   1.23906372e-9,
 123   2.1848049e-10,
 124  -1.0134282e-10,
 125   2.484728e-11,
 126  -4.73067e-12,
 127   7.3555e-13,
 128  -8.740e-14,
 129   4.85e-15,
 130   1.23e-15,
 131  -5.6e-16,
 132   1.4e-16,
 133  -3.e-17
 134 };
 135 static cheb_series fd_1_c_cs = {
 136   fd_1_c_data,
 137   22,
 138   -1, 1,
 139   13
 140 };
 141
 142
 143 /* Chebyshev fit for F_{1}(x) / x^2
 144  * 10 < x < 30
 145  * -1 < t < 1
 146  * t = 1/10 (x-10) - 1 = x/10 - 2
 147  * x = 10(t+2)
 148  */
 149 static double fd_1_d_data[30] = {
 150   1.0126626021151374442,
 151  -0.0063312525536433793,
 152   0.0024837319237084326,
 153  -0.0008764333697726109,
 154   0.0002913344438921266,
 155  -0.0000931877907705692,
 156   0.0000290151342040275,
 157  -8.8548707259955e-6,
 158   2.6603474114517e-6,
 159  -7.891415690452e-7,
 160   2.315730237195e-7,
 161  -6.73179452963e-8,
 162   1.94048035606e-8,
 163  -5.5507129189e-9,
 164   1.5766090896e-9,
 165  -4.449310875e-10,
 166   1.248292745e-10,
 167  -3.48392894e-11,
 168   9.6791550e-12,
 169  -2.6786240e-12,
 170   7.388852e-13,
 171  -2.032828e-13,
 172   5.58115e-14,
 173  -1.52987e-14,
 174   4.1886e-15,
 175  -1.1458e-15,
 176   3.132e-16,
 177  -8.56e-17,
 178   2.33e-17,
 179  -5.9e-18
 180 };
 181 static cheb_series fd_1_d_cs = {
 182   fd_1_d_data,
 183   29,
 184   -1, 1,
 185   14
 186 };
 187
 188
 189 /* Chebyshev fit for F_{1}(x) / x^2
 190  * 30 < x < Inf
 191  * -1 < t < 1
 192  * t = 60/x - 1
 193  * x = 60/(t+1)
 194  */
 195 static double fd_1_e_data[10] = {
 196   1.0013707783890401683,
 197   0.0009138522593601060,
 198   0.0002284630648400133,
 199  -1.57e-17,
 200  -1.27e-17,
 201  -9.7e-18,
 202  -6.9e-18,
 203  -4.6e-18,
 204  -2.9e-18,
 205  -1.7e-18
 206 };
 207 static cheb_series fd_1_e_cs = {
 208   fd_1_e_data,
 209   9,
 210   -1, 1,
 211   4
 212 };
 213
 214
 215 /* Chebyshev fit for F_{2}(t);  -1 < t < 1, -1 < x < 1
 216  */
 217 static double fd_2_a_data[21] = {
 218   2.1573661917148458336,
 219   0.8849670334241132182,
 220   0.1784163467613519713,
 221   0.0208333333333333333,
 222   0.0012708226459768508,
 223   0.0,
 224  -5.0619314244895e-6,
 225   0.0,
 226   4.32026533989e-8,
 227   0.0,
 228  -4.870544166e-10,
 229   0.0,
 230   6.4203740e-12,
 231   0.0,
 232  -9.37424e-14,
 233   0.0,
 234   1.4715e-15,
 235   0.0,
 236  -2.44e-17,
 237   0.0,
 238   4.e-19
 239 };
 240 static cheb_series fd_2_a_cs = {
 241   fd_2_a_data,
 242   20,
 243   -1, 1,
 244   12
 245 };
 246
 247
 248 /* Chebyshev fit for F_{2}(3/2(t+1) + 1);  -1 < t < 1, 1 < x < 4
 249  */
 250 static double fd_2_b_data[22] = {
 251   16.508258811798623599,
 252   7.421719394793067988,
 253   1.458309885545603821,
 254   0.128773850882795229,
 255   0.001963612026198147,
 256  -0.000237458988738779,
 257   0.000018539661382641,
 258  -1.92805649479e-7,
 259  -2.01950028452e-7,
 260   3.2963497518e-8,
 261  -1.885817092e-9,
 262  -2.72632744e-10,
 263   8.0554561e-11,
 264  -8.313223e-12,
 265  -2.24489e-13,
 266   2.18778e-13,
 267  -3.4290e-14,
 268   1.225e-15,
 269   5.81e-16,
 270  -1.37e-16,
 271   1.2e-17,
 272   1.e-18
 273 };
 274 static cheb_series fd_2_b_cs = {
 275   fd_2_b_data,
 276   21,
 277   -1, 1,
 278   12
 279 };
 280
 281
 282 /* Chebyshev fit for F_{1}(3(t+1) + 4);  -1 < t < 1, 4 < x < 10
 283  */
 284 static double fd_2_c_data[20] = {
 285   168.87129776686440711,
 286   81.80260488091659458,
 287   15.75408505947931513,
 288   1.12325586765966440,
 289   0.00059057505725084,
 290  -0.00016469712946921,
 291   0.00003885607810107,
 292  -7.89873660613e-6,
 293   1.39786238616e-6,
 294  -2.1534528656e-7,
 295   2.831510953e-8,
 296  -2.94978583e-9,
 297   1.6755082e-10,
 298   2.234229e-11,
 299  -1.035130e-11,
 300   2.41117e-12,
 301  -4.3531e-13,
 302   6.447e-14,
 303  -7.39e-15,
 304   4.3e-16
 305 };
 306 static cheb_series fd_2_c_cs = {
 307   fd_2_c_data,
 308   19,
 309   -1, 1,
 310   12
 311 };
 312
 313
 314 /* Chebyshev fit for F_{1}(x) / x^3
 315  * 10 < x < 30
 316  * -1 < t < 1
 317  * t = 1/10 (x-10) - 1 = x/10 - 2
 318  * x = 10(t+2)
 319  */
 320 static double fd_2_d_data[30] = {
 321   0.3459960518965277589,
 322  -0.00633136397691958024,
 323   0.00248382959047594408,
 324  -0.00087651191884005114,
 325   0.00029139255351719932,
 326  -0.00009322746111846199,
 327   0.00002904021914564786,
 328  -8.86962264810663e-6,
 329   2.66844972574613e-6,
 330  -7.9331564996004e-7,
 331   2.3359868615516e-7,
 332  -6.824790880436e-8,
 333   1.981036528154e-8,
 334  -5.71940426300e-9,
 335   1.64379426579e-9,
 336  -4.7064937566e-10,
 337   1.3432614122e-10,
 338  -3.823400534e-11,
 339   1.085771994e-11,
 340  -3.07727465e-12,
 341   8.7064848e-13,
 342  -2.4595431e-13,
 343   6.938531e-14,
 344  -1.954939e-14,
 345   5.50162e-15,
 346  -1.54657e-15,
 347   4.3429e-16,
 348  -1.2178e-16,
 349   3.394e-17,
 350  -8.81e-18
 351 };
 352 static cheb_series fd_2_d_cs = {
 353   fd_2_d_data,
 354   29,
 355   -1, 1,
 356   14
 357 };
 358
 359
 360 /* Chebyshev fit for F_{2}(x) / x^3
 361  * 30 < x < Inf
 362  * -1 < t < 1
 363  * t = 60/x - 1
 364  * x = 60/(t+1)
 365  */
 366 static double fd_2_e_data[4] = {
 367   0.3347041117223735227,
 368   0.00091385225936012645,
 369   0.00022846306484003205,
 370   5.2e-19
 371 };
 372 static cheb_series fd_2_e_cs = {
 373   fd_2_e_data,
 374   3,
 375   -1, 1,
 376   3
 377 };
 378
 379
 380 /* Chebyshev fit for F_{-1/2}(t);  -1 < t < 1, -1 < x < 1
 381  */
 382 static double fd_mhalf_a_data[20] = {
 383   1.2663290042859741974,
 384   0.3697876251911153071,
 385   0.0278131011214405055,
 386  -0.0033332848565672007,
 387  -0.0004438108265412038,
 388   0.0000616495177243839,
 389   8.7589611449897e-6,
 390  -1.2622936986172e-6,
 391  -1.837464037221e-7,
 392   2.69495091400e-8,
 393   3.9760866257e-9,
 394  -5.894468795e-10,
 395  -8.77321638e-11,
 396   1.31016571e-11,
 397   1.9621619e-12,
 398  -2.945887e-13,
 399  -4.43234e-14,
 400   6.6816e-15,
 401   1.0084e-15,
 402  -1.561e-16
 403 };
 404 static cheb_series fd_mhalf_a_cs = {
 405   fd_mhalf_a_data,
 406   19,
 407   -1, 1,
 408   12
 409 };
 410
 411
 412 /* Chebyshev fit for F_{-1/2}(3/2(t+1) + 1);  -1 < t < 1, 1 < x < 4
 413  */
 414 static double fd_mhalf_b_data[20] = {
 415   3.270796131942071484,
 416   0.5809004935853417887,
 417  -0.0299313438794694987,
 418  -0.0013287935412612198,
 419   0.0009910221228704198,
 420  -0.0001690954939688554,
 421   6.5955849946915e-6,
 422   3.5953966033618e-6,
 423  -9.430672023181e-7,
 424   8.75773958291e-8,
 425   1.06247652607e-8,
 426  -4.9587006215e-9,
 427   7.160432795e-10,
 428   4.5072219e-12,
 429  -2.3695425e-11,
 430   4.9122208e-12,
 431  -2.905277e-13,
 432  -9.59291e-14,
 433   3.00028e-14,
 434  -3.4970e-15
 435 };
 436 static cheb_series fd_mhalf_b_cs = {
 437   fd_mhalf_b_data,
 438   19,
 439   -1, 1,
 440   12
 441 };
 442
 443
 444 /* Chebyshev fit for F_{-1/2}(3(t+1) + 4);  -1 < t < 1, 4 < x < 10
 445  */
 446 static double fd_mhalf_c_data[25] = {
 447   5.828283273430595507,
 448   0.677521118293264655,
 449  -0.043946248736481554,
 450   0.005825595781828244,
 451  -0.000864858907380668,
 452   0.000110017890076539,
 453  -6.973305225404e-6,
 454  -1.716267414672e-6,
 455   8.59811582041e-7,
 456  -2.33066786976e-7,
 457   4.8503191159e-8,
 458  -8.130620247e-9,
 459   1.021068250e-9,
 460  -5.3188423e-11,
 461  -1.9430559e-11,
 462   8.750506e-12,
 463  -2.324897e-12,
 464   4.83102e-13,
 465  -8.1207e-14,
 466   1.0132e-14,
 467  -4.64e-16,
 468  -2.24e-16,
 469   9.7e-17,
 470  -2.6e-17,
 471   5.e-18
 472 };
 473 static cheb_series fd_mhalf_c_cs = {
 474   fd_mhalf_c_data,
 475   24,
 476   -1, 1,
 477   13
 478 };
 479
 480
 481 /* Chebyshev fit for F_{-1/2}(x) / x^(1/2)
 482  * 10 < x < 30
 483  * -1 < t < 1
 484  * t = 1/10 (x-10) - 1 = x/10 - 2
 485  */
 486 static double fd_mhalf_d_data[30] = {
 487   2.2530744202862438709,
 488   0.0018745152720114692,
 489  -0.0007550198497498903,
 490   0.0002759818676644382,
 491  -0.0000959406283465913,
 492   0.0000324056855537065,
 493  -0.0000107462396145761,
 494   3.5126865219224e-6,
 495  -1.1313072730092e-6,
 496   3.577454162766e-7,
 497  -1.104926666238e-7,
 498   3.31304165692e-8,
 499  -9.5837381008e-9,
 500   2.6575790141e-9,
 501  -7.015201447e-10,
 502   1.747111336e-10,
 503  -4.04909605e-11,
 504   8.5104999e-12,
 505  -1.5261885e-12,
 506   1.876851e-13,
 507   1.00574e-14,
 508  -1.82002e-14,
 509   8.6634e-15,
 510  -3.2058e-15,
 511   1.0572e-15,
 512  -3.259e-16,
 513   9.60e-17,
 514  -2.74e-17,
 515   7.6e-18,
 516  -1.9e-18
 517 };
 518 static cheb_series fd_mhalf_d_cs = {
 519   fd_mhalf_d_data,
 520   29,
 521   -1, 1,
 522   15
 523 };
 524
 525
 526 /* Chebyshev fit for F_{1/2}(t);  -1 < t < 1, -1 < x < 1
 527  */
 528 static double fd_half_a_data[23] = {
 529   1.7177138871306189157,
 530   0.6192579515822668460,
 531   0.0932802275119206269,
 532   0.0047094853246636182,
 533  -0.0004243667967864481,
 534  -0.0000452569787686193,
 535   5.2426509519168e-6,
 536   6.387648249080e-7,
 537  -8.05777004848e-8,
 538  -1.04290272415e-8,
 539   1.3769478010e-9,
 540   1.847190359e-10,
 541  -2.51061890e-11,
 542  -3.4497818e-12,
 543   4.784373e-13,
 544   6.68828e-14,
 545  -9.4147e-15,
 546  -1.3333e-15,
 547   1.898e-16,
 548   2.72e-17,
 549  -3.9e-18,
 550  -6.e-19,
 551   1.e-19
 552 };
 553 static cheb_series fd_half_a_cs = {
 554   fd_half_a_data,
 555   22,
 556   -1, 1,
 557   11
 558 };
 559
 560
 561 /* Chebyshev fit for F_{1/2}(3/2(t+1) + 1);  -1 < t < 1, 1 < x < 4
 562  */
 563 static double fd_half_b_data[20] = {
 564   7.651013792074984027,
 565   2.475545606866155737,
 566   0.218335982672476128,
 567  -0.007730591500584980,
 568  -0.000217443383867318,
 569   0.000147663980681359,
 570  -0.000021586361321527,
 571   8.07712735394e-7,
 572   3.28858050706e-7,
 573  -7.9474330632e-8,
 574   6.940207234e-9,
 575   6.75594681e-10,
 576  -3.10200490e-10,
 577   4.2677233e-11,
 578  -2.1696e-14,
 579  -1.170245e-12,
 580   2.34757e-13,
 581  -1.4139e-14,
 582  -3.864e-15,
 583   1.202e-15
 584 };
 585 static cheb_series fd_half_b_cs = {
 586   fd_half_b_data,
 587   19,
 588   -1, 1,
 589   12
 590 };
 591
 592
 593 /* Chebyshev fit for F_{1/2}(3(t+1) + 4);  -1 < t < 1, 4 < x < 10
 594  */
 595 static double fd_half_c_data[23] = {
 596   29.584339348839816528,
 597   8.808344283250615592,
 598   0.503771641883577308,
 599  -0.021540694914550443,
 600   0.002143341709406890,
 601  -0.000257365680646579,
 602   0.000027933539372803,
 603  -1.678525030167e-6,
 604  -2.78100117693e-7,
 605   1.35218065147e-7,
 606  -3.3740425009e-8,
 607   6.474834942e-9,
 608  -1.009678978e-9,
 609   1.20057555e-10,
 610  -6.636314e-12,
 611  -1.710566e-12,
 612   7.75069e-13,
 613  -1.97973e-13,
 614   3.9414e-14,
 615  -6.374e-15,
 616   7.77e-16,
 617  -4.0e-17,
 618  -1.4e-17
 619 };
 620 static cheb_series fd_half_c_cs = {
 621   fd_half_c_data,
 622   22,
 623   -1, 1,
 624   13
 625 };
 626
 627
 628 /* Chebyshev fit for F_{1/2}(x) / x^(3/2)
 629  * 10 < x < 30
 630  * -1 < t < 1
 631  * t = 1/10 (x-10) - 1 = x/10 - 2
 632  */
 633 static double fd_half_d_data[30] = {
 634   1.5116909434145508537,
 635  -0.0036043405371630468,
 636   0.0014207743256393359,
 637  -0.0005045399052400260,
 638   0.0001690758006957347,
 639  -0.0000546305872688307,
 640   0.0000172223228484571,
 641  -5.3352603788706e-6,
 642   1.6315287543662e-6,
 643  -4.939021084898e-7,
 644   1.482515450316e-7,
 645  -4.41552276226e-8,
 646   1.30503160961e-8,
 647  -3.8262599802e-9,
 648   1.1123226976e-9,
 649  -3.204765534e-10,
 650   9.14870489e-11,
 651  -2.58778946e-11,
 652   7.2550731e-12,
 653  -2.0172226e-12,
 654   5.566891e-13,
 655  -1.526247e-13,
 656   4.16121e-14,
 657  -1.12933e-14,
 658   3.0537e-15,
 659  -8.234e-16,
 660   2.215e-16,
 661  -5.95e-17,
 662   1.59e-17,
 663  -4.0e-18
 664 };
 665 static cheb_series fd_half_d_cs = {
 666   fd_half_d_data,
 667   29,
 668   -1, 1,
 669   15
 670 };
 671
 672
 673
 674 /* Chebyshev fit for F_{3/2}(t);  -1 < t < 1, -1 < x < 1
 675  */
 676 static double fd_3half_a_data[20] = {
 677   2.0404775940601704976,
 678   0.8122168298093491444,
 679   0.1536371165644008069,
 680   0.0156174323847845125,
 681   0.0005943427879290297,
 682  -0.0000429609447738365,
 683  -3.8246452994606e-6,
 684   3.802306180287e-7,
 685   4.05746157593e-8,
 686  -4.5530360159e-9,
 687  -5.306873139e-10,
 688   6.37297268e-11,
 689   7.8403674e-12,
 690  -9.840241e-13,
 691  -1.255952e-13,
 692   1.62617e-14,
 693   2.1318e-15,
 694  -2.825e-16,
 695  -3.78e-17,
 696   5.1e-18
 697 };
 698 static cheb_series fd_3half_a_cs = {
 699   fd_3half_a_data,
 700   19,
 701   -1, 1,
 702   11
 703 };
 704
 705
 706 /* Chebyshev fit for F_{3/2}(3/2(t+1) + 1);  -1 < t < 1, 1 < x < 4
 707  */
 708 static double fd_3half_b_data[22] = {
 709   13.403206654624176674,
 710   5.574508357051880924,
 711   0.931228574387527769,
 712   0.054638356514085862,
 713  -0.001477172902737439,
 714  -0.000029378553381869,
 715   0.000018357033493246,
 716  -2.348059218454e-6,
 717   8.3173787440e-8,
 718   2.6826486956e-8,
 719  -6.011244398e-9,
 720   4.94345981e-10,
 721   3.9557340e-11,
 722  -1.7894930e-11,
 723   2.348972e-12,
 724  -1.2823e-14,
 725  -5.4192e-14,
 726   1.0527e-14,
 727  -6.39e-16,
 728  -1.47e-16,
 729   4.5e-17,
 730  -5.e-18
 731 };
 732 static cheb_series fd_3half_b_cs = {
 733   fd_3half_b_data,
 734   21,
 735   -1, 1,
 736   12
 737 };
 738
 739
 740 /* Chebyshev fit for F_{3/2}(3(t+1) + 4);  -1 < t < 1, 4 < x < 10
 741  */
 742 static double fd_3half_c_data[21] = {
 743   101.03685253378877642,
 744   43.62085156043435883,
 745   6.62241373362387453,
 746   0.25081415008708521,
 747  -0.00798124846271395,
 748   0.00063462245101023,
 749  -0.00006392178890410,
 750   6.04535131939e-6,
 751  -3.4007683037e-7,
 752  -4.072661545e-8,
 753   1.931148453e-8,
 754  -4.46328355e-9,
 755   7.9434717e-10,
 756  -1.1573569e-10,
 757   1.304658e-11,
 758  -7.4114e-13,
 759  -1.4181e-13,
 760   6.491e-14,
 761  -1.597e-14,
 762   3.05e-15,
 763  -4.8e-16
 764 };
 765 static cheb_series fd_3half_c_cs = {
 766   fd_3half_c_data,
 767   20,
 768   -1, 1,
 769   12
 770 };
 771
 772
 773 /* Chebyshev fit for F_{3/2}(x) / x^(5/2)
 774  * 10 < x < 30
 775  * -1 < t < 1
 776  * t = 1/10 (x-10) - 1 = x/10 - 2
 777  */
 778 static double fd_3half_d_data[25] = {
 779   0.6160645215171852381,
 780  -0.0071239478492671463,
 781   0.0027906866139659846,
 782  -0.0009829521424317718,
 783   0.0003260229808519545,
 784  -0.0001040160912910890,
 785   0.0000322931223232439,
 786  -9.8243506588102e-6,
 787   2.9420132351277e-6,
 788  -8.699154670418e-7,
 789   2.545460071999e-7,
 790  -7.38305056331e-8,
 791   2.12545670310e-8,
 792  -6.0796532462e-9,
 793   1.7294556741e-9,
 794  -4.896540687e-10,
 795   1.380786037e-10,
 796  -3.88057305e-11,
 797   1.08753212e-11,
 798  -3.0407308e-12,
 799   8.485626e-13,
 800  -2.364275e-13,
 801   6.57636e-14,
 802  -1.81807e-14,
 803   4.6884e-15
 804 };
 805 static cheb_series fd_3half_d_cs = {
 806   fd_3half_d_data,
 807   24,
 808   -1, 1,
 809   16
 810 };
 811
 812
 813 /* Goano's modification of the Levin-u implementation.
 814  * This is a simplification of the original WHIZ algorithm.
 815  * See [Fessler et al., ACM Toms 9, 346 (1983)].
 816  */
 817 static
 818 int
 819 fd_whiz(const double term, const int iterm,
 820         double * qnum, double * qden,
 821         double * result, double * s)
 822 {
 823   if(iterm == 0) *s = 0.0;
 824
 825   *s += term;
 826
 827   qden[iterm] = 1.0/(term*(iterm+1.0)*(iterm+1.0));
 828   qnum[iterm] = *s * qden[iterm];
 829
 830   if(iterm > 0) {
 831     double factor = 1.0;
 832     double ratio  = iterm/(iterm+1.0);
 833     int j;
 834     for(j=iterm-1; j>=0; j--) {
 835       double c = factor * (j+1.0) / (iterm+1.0);
 836       factor *= ratio;
 837       qden[j] = qden[j+1] - c * qden[j];
 838       qnum[j] = qnum[j+1] - c * qnum[j];
 839     }
 840   }
 841
 842   *result = qnum[0] / qden[0];
 843   return GSL_SUCCESS;
 844 }
 845
 846
 847 /* Handle case of integer j <= -2.
 848  */
 849 static
 850 int
 851 fd_nint(const int j, const double x, gsl_sf_result * result)
 852 {
 853 /*    const int nsize = 100 + 1; */
 854   enum {
 855     nsize = 100+1
 856   };
 857   double qcoeff[nsize];
 858
 859   if(j >= -1) {
 860     result->val = 0.0;
 861     result->err = 0.0;
 862     GSL_ERROR ("error", GSL_ESANITY);
 863   }
 864   else if(j < -(nsize)) {
 865     result->val = 0.0;
 866     result->err = 0.0;
 867     GSL_ERROR ("error", GSL_EUNIMPL);
 868   }
 869   else {
 870     double a, p, f;
 871     int i, k;
 872     int n = -(j+1);
 873
 874     qcoeff[1] = 1.0;
 875
 876     for(k=2; k<=n; k++) {
 877       qcoeff[k] = -qcoeff[k-1];
 878       for(i=k-1; i>=2; i--) {
 879         qcoeff[i] = i*qcoeff[i] - (k-(i-1))*qcoeff[i-1];
 880       }
 881     }
 882
 883     if(x >= 0.0) {
 884       a = exp(-x);
 885       f = qcoeff[1];
 886       for(i=2; i<=n; i++) {
 887         f = f*a + qcoeff[i];
 888       }
 889     }
 890     else {
 891       a = exp(x);
 892       f = qcoeff[n];
 893       for(i=n-1; i>=1; i--) {
 894         f = f*a + qcoeff[i];
 895       }
 896     }
 897
 898     p = gsl_sf_pow_int(1.0+a, j);
 899     result->val = f*a*p;
 900     result->err = 3.0 * GSL_DBL_EPSILON * fabs(f*a*p);
 901     return GSL_SUCCESS;
 902   }
 903 }
 904
 905
 906 /* x < 0
 907  */
 908 static
 909 int
 910 fd_neg(const double j, const double x, gsl_sf_result * result)
 911 {
 912   enum {
 913     itmax = 100,
 914     qsize = 100+1
 915   };
 916 /*    const int itmax = 100; */
 917 /*    const int qsize = 100 + 1; */
 918   double qnum[qsize], qden[qsize];
 919
 920   if(x < GSL_LOG_DBL_MIN) {
 921     result->val = 0.0;
 922     result->err = 0.0;
 923     return GSL_SUCCESS;
 924   }
 925   else if(x < -1.0 && x < -fabs(j+1.0)) {
 926     /* Simple series implementation. Avoid the
 927      * complexity and extra work of the series
 928      * acceleration method below.
 929      */
 930     double ex   = exp(x);
 931     double term = ex;
 932     double sum  = term;
 933     int n;
 934     for(n=2; n<100; n++) {
 935       double rat = (n-1.0)/n;
 936       double p   = pow(rat, j+1.0);
 937       term *= -ex * p;
 938       sum  += term;
 939       if(fabs(term/sum) < GSL_DBL_EPSILON) break;
 940     }
 941     result->val = sum;
 942     result->err = 2.0 * GSL_DBL_EPSILON * fabs(sum);
 943     return GSL_SUCCESS;
 944   }
 945   else {
 946     double s = 0.0;
 947     double xn = x;
 948     double ex  = -exp(x);
 949     double enx = -ex;
 950     double f = 0.0;
 951     double f_previous;
 952     int jterm;
 953     for(jterm=0; jterm<=itmax; jterm++) {
 954       double p = pow(jterm+1.0, j+1.0);
 955       double term = enx/p;
 956       f_previous = f;
 957       fd_whiz(term, jterm, qnum, qden, &f, &s);
 958       xn += x;
 959       if(fabs(f-f_previous) < fabs(f)*2.0*GSL_DBL_EPSILON || xn < GSL_LOG_DBL_MIN) break;
 960       enx *= ex;
 961     }
 962
 963     result->val  = f;
 964     result->err  = fabs(f-f_previous);
 965     result->err += 2.0 * GSL_DBL_EPSILON * fabs(f);
 966
 967     if(jterm == itmax)
 968       GSL_ERROR ("error", GSL_EMAXITER);
 969     else
 970       return GSL_SUCCESS;
 971   }
 972 }
 973
 974
 975 /* asymptotic expansion
 976  * j + 2.0 > 0.0
 977  */
 978 static
 979 int
 980 fd_asymp(const double j, const double x, gsl_sf_result * result)
 981 {
 982   const int j_integer = ( fabs(j - floor(j+0.5)) < 100.0*GSL_DBL_EPSILON );
 983   const int itmax = 200;
 984   gsl_sf_result lg;
 985   int stat_lg = gsl_sf_lngamma_e(j + 2.0, &lg);
 986   double seqn_val = 0.5;
 987   double seqn_err = 0.0;
 988   double xm2  = (1.0/x)/x;
 989   double xgam = 1.0;
 990   double add = GSL_DBL_MAX;
 991   double cos_term;
 992   double ln_x;
 993   double ex_term_1;
 994   double ex_term_2;
 995   gsl_sf_result fneg;
 996   gsl_sf_result ex_arg;
 997   gsl_sf_result ex;
 998   int stat_fneg;
 999   int stat_e;
1000   int n;
1001   for(n=1; n<=itmax; n++) {
1002     double add_previous = add;
1003     gsl_sf_result eta;
1004     gsl_sf_eta_int_e(2*n, &eta);
1005     xgam = xgam * xm2 * (j + 1.0 - (2*n-2)) * (j + 1.0 - (2*n-1));
1006     add  = eta.val * xgam;
1007     if(!j_integer && fabs(add) > fabs(add_previous)) break;
1008     if(fabs(add/seqn_val) < GSL_DBL_EPSILON) break;
1009     seqn_val += add;
1010     seqn_err += 2.0 * GSL_DBL_EPSILON * fabs(add);
1011   }
1012   seqn_err += fabs(add);
1013
1014   stat_fneg = fd_neg(j, -x, &fneg);
1015   ln_x = log(x);
1016   ex_term_1 = (j+1.0)*ln_x;
1017   ex_term_2 = lg.val;
1018   ex_arg.val = ex_term_1 - ex_term_2; /*(j+1.0)*ln_x - lg.val; */
1019   ex_arg.err = GSL_DBL_EPSILON*(fabs(ex_term_1) + fabs(ex_term_2)) + lg.err;
1020   stat_e    = gsl_sf_exp_err_e(ex_arg.val, ex_arg.err, &ex);
1021   cos_term  = cos(j*M_PI);
1022   result->val  = cos_term * fneg.val + 2.0 * seqn_val * ex.val;
1023   result->err  = fabs(2.0 * ex.err * seqn_val);
1024   result->err += fabs(2.0 * ex.val * seqn_err);
1025   result->err += fabs(cos_term) * fneg.err;
1026   result->err += 4.0 * GSL_DBL_EPSILON * fabs(result->val);
1027   return GSL_ERROR_SELECT_3(stat_e, stat_fneg, stat_lg);
1028 }
1029
1030
1031 /* Series evaluation for small x, generic j.
1032  * [Goano (8)]
1033  */
1034 #if 0
1035 static
1036 int
1037 fd_series(const double j, const double x, double * result)
1038 {
1039   const int nmax = 1000;
1040   int n;
1041   double sum = 0.0;
1042   double prev;
1043   double pow_factor = 1.0;
1044   double eta_factor;
1045   gsl_sf_eta_e(j + 1.0, &eta_factor);
1046   prev = pow_factor * eta_factor;
1047   sum += prev;
1048   for(n=1; n<nmax; n++) {
1049     double term;
1050     gsl_sf_eta_e(j+1.0-n, &eta_factor);
1051     pow_factor *= x/n;
1052     term = pow_factor * eta_factor;
1053     sum += term;
1054     if(fabs(term/sum) < GSL_DBL_EPSILON && fabs(prev/sum) < GSL_DBL_EPSILON) break;
1055     prev = term;
1056   }
1057
1058   *result = sum;
1059   return GSL_SUCCESS;
1060 }
1061 #endif /* 0 */
1062
1063
1064 /* Series evaluation for small x > 0, integer j > 0; x < Pi.
1065  * [Goano (8)]
1066  */
1067 static
1068 int
1069 fd_series_int(const int j, const double x, gsl_sf_result * result)
1070 {
1071   int n;
1072   double sum = 0.0;
1073   double del;
1074   double pow_factor = 1.0;
1075   gsl_sf_result eta_factor;
1076   gsl_sf_eta_int_e(j + 1, &eta_factor);
1077   del  = pow_factor * eta_factor.val;
1078   sum += del;
1079
1080   /* Sum terms where the argument
1081    * of eta() is positive.
1082    */
1083   for(n=1; n<=j+2; n++) {
1084     gsl_sf_eta_int_e(j+1-n, &eta_factor);
1085     pow_factor *= x/n;
1086     del  = pow_factor * eta_factor.val;
1087     sum += del;
1088     if(fabs(del/sum) < GSL_DBL_EPSILON) break;
1089   }
1090
1091   /* Now sum the terms where eta() is negative.
1092    * The argument of eta() must be odd as well,
1093    * so it is convenient to transform the series
1094    * as follows:
1095    *
1096    * Sum[ eta(j+1-n) x^n / n!, {n,j+4,Infinity}]
1097    * = x^j / j! Sum[ eta(1-2m) x^(2m) j! / (2m+j)! , {m,2,Infinity}]
1098    *
1099    * We do not need to do this sum if j is large enough.
1100    */
1101   if(j < 32) {
1102     int m;
1103     gsl_sf_result jfact;
1104     double sum2;
1105     double pre2;
1106
1107     gsl_sf_fact_e((unsigned int)j, &jfact);
1108     pre2 = gsl_sf_pow_int(x, j) / jfact.val;
1109
1110     gsl_sf_eta_int_e(-3, &eta_factor);
1111     pow_factor = x*x*x*x / ((j+4)*(j+3)*(j+2)*(j+1));
1112     sum2 = eta_factor.val * pow_factor;
1113
1114     for(m=3; m<24; m++) {
1115       gsl_sf_eta_int_e(1-2*m, &eta_factor);
1116       pow_factor *= x*x / ((j+2*m)*(j+2*m-1));
1117       sum2 += eta_factor.val * pow_factor;
1118     }
1119
1120     sum += pre2 * sum2;
1121   }
1122
1123   result->val = sum;
1124   result->err = 2.0 * GSL_DBL_EPSILON * fabs(sum);
1125
1126   return GSL_SUCCESS;
1127 }
1128
1129
1130 /* series of hypergeometric functions for integer j > 0, x > 0
1131  * [Goano (7)]
1132  */
1133 static
1134 int
1135 fd_UMseries_int(const int j, const double x, gsl_sf_result * result)
1136 {
1137   const int nmax = 2000;
1138   double pre;
1139   double lnpre_val;
1140   double lnpre_err;
1141   double sum_even_val = 1.0;
1142   double sum_even_err = 0.0;
1143   double sum_odd_val  = 0.0;
1144   double sum_odd_err  = 0.0;
1145   int stat_sum;
1146   int stat_e;
1147   int stat_h = GSL_SUCCESS;
1148   int n;
1149
1150   if(x < 500.0 && j < 80) {
1151     double p = gsl_sf_pow_int(x, j+1);
1152     gsl_sf_result g;
1153     gsl_sf_fact_e(j+1, &g); /* Gamma(j+2) */
1154     lnpre_val = 0.0;
1155     lnpre_err = 0.0;
1156     pre   = p/g.val;
1157   }
1158   else {
1159     double lnx = log(x);
1160     gsl_sf_result lg;
1161     gsl_sf_lngamma_e(j + 2.0, &lg);
1162     lnpre_val = (j+1.0)*lnx - lg.val;
1163     lnpre_err = 2.0 * GSL_DBL_EPSILON * fabs((j+1.0)*lnx) + lg.err;
1164     pre = 1.0;
1165   }
1166
1167   /* Add up the odd terms of the sum.
1168    */
1169   for(n=1; n<nmax; n+=2) {
1170     double del_val;
1171     double del_err;
1172     gsl_sf_result U;
1173     gsl_sf_result M;
1174     int stat_h_U = gsl_sf_hyperg_U_int_e(1, j+2, n*x, &U);
1175     int stat_h_F = gsl_sf_hyperg_1F1_int_e(1, j+2, -n*x, &M);
1176     stat_h = GSL_ERROR_SELECT_3(stat_h, stat_h_U, stat_h_F);
1177     del_val = ((j+1.0)*U.val - M.val);
1178     del_err = (fabs(j+1.0)*U.err + M.err);
1179     sum_odd_val += del_val;
1180     sum_odd_err += del_err;
1181     if(fabs(del_val/sum_odd_val) < GSL_DBL_EPSILON) break;
1182   }
1183
1184   /* Add up the even terms of the sum.
1185    */
1186   for(n=2; n<nmax; n+=2) {
1187     double del_val;
1188     double del_err;
1189     gsl_sf_result U;
1190     gsl_sf_result M;
1191     int stat_h_U = gsl_sf_hyperg_U_int_e(1, j+2, n*x, &U);
1192     int stat_h_F = gsl_sf_hyperg_1F1_int_e(1, j+2, -n*x, &M);
1193     stat_h = GSL_ERROR_SELECT_3(stat_h, stat_h_U, stat_h_F);
1194     del_val = ((j+1.0)*U.val - M.val);
1195     del_err = (fabs(j+1.0)*U.err + M.err);
1196     sum_even_val -= del_val;
1197     sum_even_err += del_err;
1198     if(fabs(del_val/sum_even_val) < GSL_DBL_EPSILON) break;
1199   }
1200
1201   stat_sum = ( n >= nmax ? GSL_EMAXITER : GSL_SUCCESS );
1202   stat_e   = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
1203                                       pre*(sum_even_val + sum_odd_val),
1204                                       pre*(sum_even_err + sum_odd_err),
1205                                       result);
1206   result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
1207
1208   return GSL_ERROR_SELECT_3(stat_e, stat_h, stat_sum);
1209 }
1210
1211
1212 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
1213
1214 /* [Goano (4)] */
1215 int gsl_sf_fermi_dirac_m1_e(const double x, gsl_sf_result * result)
1216 {
1217   if(x < GSL_LOG_DBL_MIN) {
1218     UNDERFLOW_ERROR(result);
1219   }
1220   else if(x < 0.0) {
1221     const double ex = exp(x);
1222     result->val = ex/(1.0+ex);
1223     result->err = 2.0 * (fabs(x) + 1.0) * GSL_DBL_EPSILON * fabs(result->val);
1224     return GSL_SUCCESS;
1225   }
1226   else {
1227     double ex = exp(-x);
1228     result->val = 1.0/(1.0 + ex);
1229     result->err = 2.0 * GSL_DBL_EPSILON * (x + 1.0) * ex;
1230     return GSL_SUCCESS;
1231   }
1232 }
1233
1234
1235 /* [Goano (3)] */
1236 int gsl_sf_fermi_dirac_0_e(const double x, gsl_sf_result * result)
1237 {
1238   if(x < GSL_LOG_DBL_MIN) {
1239     UNDERFLOW_ERROR(result);
1240   }
1241   else if(x < -5.0) {
1242     double ex  = exp(x);
1243     double ser = 1.0 - ex*(0.5 - ex*(1.0/3.0 - ex*(1.0/4.0 - ex*(1.0/5.0 - ex/6.0))));
1244     result->val = ex * ser;
1245     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
1246     return GSL_SUCCESS;
1247   }
1248   else if(x < 10.0) {
1249     result->val = log(1.0 + exp(x));
1250     result->err = fabs(x * GSL_DBL_EPSILON);
1251     return GSL_SUCCESS;
1252   }
1253   else {
1254     double ex = exp(-x);
1255     result->val = x + ex * (1.0 - 0.5*ex + ex*ex/3.0 - ex*ex*ex/4.0);
1256     result->err = (x + ex) * GSL_DBL_EPSILON;
1257     return GSL_SUCCESS;
1258   }
1259 }
1260
1261
1262 int gsl_sf_fermi_dirac_1_e(const double x, gsl_sf_result * result)
1263 {
1264   if(x < GSL_LOG_DBL_MIN) {
1265     UNDERFLOW_ERROR(result);
1266   }
1267   else if(x < -1.0) {
1268     /* series [Goano (6)]
1269      */
1270     double ex   = exp(x);
1271     double term = ex;
1272     double sum  = term;
1273     int n;
1274     for(n=2; n<100 ; n++) {
1275       double rat = (n-1.0)/n;
1276       term *= -ex * rat * rat;
1277       sum  += term;
1278       if(fabs(term/sum) < GSL_DBL_EPSILON) break;
1279     }
1280     result->val = sum;
1281     result->err = 2.0 * fabs(sum) * GSL_DBL_EPSILON;
1282     return GSL_SUCCESS;
1283   }
1284   else if(x < 1.0) {
1285     return cheb_eval_e(&fd_1_a_cs, x, result);
1286   }
1287   else if(x < 4.0) {
1288     double t = 2.0/3.0*(x-1.0) - 1.0;
1289     return cheb_eval_e(&fd_1_b_cs, t, result);
1290   }
1291   else if(x < 10.0) {
1292     double t = 1.0/3.0*(x-4.0) - 1.0;
1293     return cheb_eval_e(&fd_1_c_cs, t, result);
1294   }
1295   else if(x < 30.0) {
1296     double t = 0.1*x - 2.0;
1297     gsl_sf_result c;
1298     cheb_eval_e(&fd_1_d_cs, t, &c);
1299     result->val  = c.val * x*x;
1300     result->err  = c.err * x*x + GSL_DBL_EPSILON * fabs(result->val);
1301     return GSL_SUCCESS;
1302   }
1303   else if(x < 1.0/GSL_SQRT_DBL_EPSILON) {
1304     double t = 60.0/x - 1.0;
1305     gsl_sf_result c;
1306     cheb_eval_e(&fd_1_e_cs, t, &c);
1307     result->val  = c.val * x*x;
1308     result->err  = c.err * x*x + GSL_DBL_EPSILON * fabs(result->val);
1309     return GSL_SUCCESS;
1310   }
1311   else if(x < GSL_SQRT_DBL_MAX) {
1312     result->val = 0.5 * x*x;
1313     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
1314     return GSL_SUCCESS;
1315   }
1316   else {
1317     OVERFLOW_ERROR(result);
1318   }
1319 }
1320
1321
1322 int gsl_sf_fermi_dirac_2_e(const double x, gsl_sf_result * result)
1323 {
1324   if(x < GSL_LOG_DBL_MIN) {
1325     UNDERFLOW_ERROR(result);
1326   }
1327   else if(x < -1.0) {
1328     /* series [Goano (6)]
1329      */
1330     double ex   = exp(x);
1331     double term = ex;
1332     double sum  = term;
1333     int n;
1334     for(n=2; n<100 ; n++) {
1335       double rat = (n-1.0)/n;
1336       term *= -ex * rat * rat * rat;
1337       sum  += term;
1338       if(fabs(term/sum) < GSL_DBL_EPSILON) break;
1339     }
1340     result->val = sum;
1341     result->err = 2.0 * GSL_DBL_EPSILON * fabs(sum);
1342     return GSL_SUCCESS;
1343   }
1344   else if(x < 1.0) {
1345     return cheb_eval_e(&fd_2_a_cs, x, result);
1346   }
1347   else if(x < 4.0) {
1348     double t = 2.0/3.0*(x-1.0) - 1.0;
1349     return cheb_eval_e(&fd_2_b_cs, t, result);
1350   }
1351   else if(x < 10.0) {
1352     double t = 1.0/3.0*(x-4.0) - 1.0;
1353     return cheb_eval_e(&fd_2_c_cs, t, result);
1354   }
1355   else if(x < 30.0) {
1356     double t = 0.1*x - 2.0;
1357     gsl_sf_result c;
1358     cheb_eval_e(&fd_2_d_cs, t, &c);
1359     result->val  = c.val * x*x*x;
1360     result->err  = c.err * x*x*x + 3.0 * GSL_DBL_EPSILON * fabs(result->val);
1361     return GSL_SUCCESS;
1362   }
1363   else if(x < 1.0/GSL_ROOT3_DBL_EPSILON) {
1364     double t = 60.0/x - 1.0;
1365     gsl_sf_result c;
1366     cheb_eval_e(&fd_2_e_cs, t, &c);
1367     result->val = c.val * x*x*x;
1368     result->err = c.err * x*x*x + 3.0 * GSL_DBL_EPSILON * fabs(result->val);
1369     return GSL_SUCCESS;
1370   }
1371   else if(x < GSL_ROOT3_DBL_MAX) {
1372     result->val = 1.0/6.0 * x*x*x;
1373     result->err = 3.0 * GSL_DBL_EPSILON * fabs(result->val);
1374     return GSL_SUCCESS;
1375   }
1376   else {
1377     OVERFLOW_ERROR(result);
1378   }
1379 }
1380
1381
1382 int gsl_sf_fermi_dirac_int_e(const int j, const double x, gsl_sf_result * result)
1383 {
1384   if(j < -1) {
1385     return fd_nint(j, x, result);
1386   }
1387   else if (j == -1) {
1388     return gsl_sf_fermi_dirac_m1_e(x, result);
1389   }
1390   else if(j == 0) {
1391     return gsl_sf_fermi_dirac_0_e(x, result);
1392   }
1393   else if(j == 1) {
1394     return gsl_sf_fermi_dirac_1_e(x, result);
1395   }
1396   else if(j == 2) {
1397     return gsl_sf_fermi_dirac_2_e(x, result);
1398   }
1399   else if(x < 0.0) {
1400     return fd_neg(j, x, result);
1401   }
1402   else if(x == 0.0) {
1403     return gsl_sf_eta_int_e(j+1, result);
1404   }
1405   else if(x < 1.5) {
1406     return fd_series_int(j, x, result);
1407   }
1408   else {
1409     gsl_sf_result fasymp;
1410     int stat_asymp = fd_asymp(j, x, &fasymp);
1411
1412     if(stat_asymp == GSL_SUCCESS) {
1413       result->val  = fasymp.val;
1414       result->err  = fasymp.err;
1415       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
1416       return stat_asymp;
1417     }
1418     else {
1419       return fd_UMseries_int(j, x, result);
1420     }
1421   }
1422 }
1423
1424
1425 int gsl_sf_fermi_dirac_mhalf_e(const double x, gsl_sf_result * result)
1426 {
1427   if(x < GSL_LOG_DBL_MIN) {
1428     UNDERFLOW_ERROR(result);
1429   }
1430   else if(x < -1.0) {
1431     /* series [Goano (6)]
1432      */
1433     double ex   = exp(x);
1434     double term = ex;
1435     double sum  = term;
1436     int n;
1437     for(n=2; n<200 ; n++) {
1438       double rat = (n-1.0)/n;
1439       term *= -ex * sqrt(rat);
1440       sum  += term;
1441       if(fabs(term/sum) < GSL_DBL_EPSILON) break;
1442     }
1443     result->val = sum;
1444     result->err = 2.0 * fabs(sum) * GSL_DBL_EPSILON;
1445     return GSL_SUCCESS;
1446   }
1447   else if(x < 1.0) {
1448     return cheb_eval_e(&fd_mhalf_a_cs, x, result);
1449   }
1450   else if(x < 4.0) {
1451     double t = 2.0/3.0*(x-1.0) - 1.0;
1452     return cheb_eval_e(&fd_mhalf_b_cs, t, result);
1453   }
1454   else if(x < 10.0) {
1455     double t = 1.0/3.0*(x-4.0) - 1.0;
1456     return cheb_eval_e(&fd_mhalf_c_cs, t, result);
1457   }
1458   else if(x < 30.0) {
1459     double rtx = sqrt(x);
1460     double t = 0.1*x - 2.0;
1461     gsl_sf_result c;
1462     cheb_eval_e(&fd_mhalf_d_cs, t, &c);
1463     result->val  = c.val * rtx;
1464     result->err  = c.err * rtx + 0.5 * GSL_DBL_EPSILON * fabs(result->val);
1465     return GSL_SUCCESS;
1466   }
1467   else {
1468     return fd_asymp(-0.5, x, result);
1469   }
1470 }
1471
1472
1473 int gsl_sf_fermi_dirac_half_e(const double x, gsl_sf_result * result)
1474 {
1475   if(x < GSL_LOG_DBL_MIN) {
1476     UNDERFLOW_ERROR(result);
1477   }
1478   else if(x < -1.0) {
1479     /* series [Goano (6)]
1480      */
1481     double ex   = exp(x);
1482     double term = ex;
1483     double sum  = term;
1484     int n;
1485     for(n=2; n<100 ; n++) {
1486       double rat = (n-1.0)/n;
1487       term *= -ex * rat * sqrt(rat);
1488       sum  += term;
1489       if(fabs(term/sum) < GSL_DBL_EPSILON) break;
1490     }
1491     result->val = sum;
1492     result->err = 2.0 * fabs(sum) * GSL_DBL_EPSILON;
1493     return GSL_SUCCESS;
1494   }
1495   else if(x < 1.0) {
1496     return cheb_eval_e(&fd_half_a_cs, x, result);
1497   }
1498   else if(x < 4.0) {
1499     double t = 2.0/3.0*(x-1.0) - 1.0;
1500     return cheb_eval_e(&fd_half_b_cs, t, result);
1501   }
1502   else if(x < 10.0) {
1503     double t = 1.0/3.0*(x-4.0) - 1.0;
1504     return cheb_eval_e(&fd_half_c_cs, t, result);
1505   }
1506   else if(x < 30.0) {
1507     double x32 = x*sqrt(x);
1508     double t = 0.1*x - 2.0;
1509     gsl_sf_result c;
1510     cheb_eval_e(&fd_half_d_cs, t, &c);
1511     result->val = c.val * x32;
1512     result->err = c.err * x32 + 1.5 * GSL_DBL_EPSILON * fabs(result->val);
1513     return GSL_SUCCESS;
1514   }
1515   else {
1516     return fd_asymp(0.5, x, result);
1517   }
1518 }
1519
1520
1521 int gsl_sf_fermi_dirac_3half_e(const double x, gsl_sf_result * result)
1522 {
1523   if(x < GSL_LOG_DBL_MIN) {
1524     UNDERFLOW_ERROR(result);
1525   }
1526   else if(x < -1.0) {
1527     /* series [Goano (6)]
1528      */
1529     double ex   = exp(x);
1530     double term = ex;
1531     double sum  = term;
1532     int n;
1533     for(n=2; n<100 ; n++) {
1534       double rat = (n-1.0)/n;
1535       term *= -ex * rat * rat * sqrt(rat);
1536       sum  += term;
1537       if(fabs(term/sum) < GSL_DBL_EPSILON) break;
1538     }
1539     result->val = sum;
1540     result->err = 2.0 * fabs(sum) * GSL_DBL_EPSILON;
1541     return GSL_SUCCESS;
1542   }
1543   else if(x < 1.0) {
1544     return cheb_eval_e(&fd_3half_a_cs, x, result);
1545   }
1546   else if(x < 4.0) {
1547     double t = 2.0/3.0*(x-1.0) - 1.0;
1548     return cheb_eval_e(&fd_3half_b_cs, t, result);
1549   }
1550   else if(x < 10.0) {
1551     double t = 1.0/3.0*(x-4.0) - 1.0;
1552     return cheb_eval_e(&fd_3half_c_cs, t, result);
1553   }
1554   else if(x < 30.0) {
1555     double x52 = x*x*sqrt(x);
1556     double t = 0.1*x - 2.0;
1557     gsl_sf_result c;
1558     cheb_eval_e(&fd_3half_d_cs, t, &c);
1559     result->val = c.val * x52;
1560     result->err = c.err * x52 + 2.5 * GSL_DBL_EPSILON * fabs(result->val);
1561     return GSL_SUCCESS;
1562   }
1563   else {
1564     return fd_asymp(1.5, x, result);
1565   }
1566 }
1567
1568 /* [Goano p. 222] */
1569 int gsl_sf_fermi_dirac_inc_0_e(const double x, const double b, gsl_sf_result * result)
1570 {
1571   if(b < 0.0) {
1572     DOMAIN_ERROR(result);
1573   }
1574   else {
1575     double arg = b - x;
1576     gsl_sf_result f0;
1577     int status = gsl_sf_fermi_dirac_0_e(arg, &f0);
1578     result->val = f0.val - arg;
1579     result->err = f0.err + GSL_DBL_EPSILON * (fabs(x) + fabs(b));
1580     return status;
1581   }
1582 }
1583
1584
1585
1586 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
1587
1588 #include "eval.h"
1589
1590 double gsl_sf_fermi_dirac_m1(const double x)
1591 {
1592   EVAL_RESULT(gsl_sf_fermi_dirac_m1_e(x, &result));
1593 }
1594
1595 double gsl_sf_fermi_dirac_0(const double x)
1596 {
1597   EVAL_RESULT(gsl_sf_fermi_dirac_0_e(x, &result));
1598 }
1599
1600 double gsl_sf_fermi_dirac_1(const double x)
1601 {
1602   EVAL_RESULT(gsl_sf_fermi_dirac_1_e(x, &result));
1603 }
1604
1605 double gsl_sf_fermi_dirac_2(const double x)
1606 {
1607   EVAL_RESULT(gsl_sf_fermi_dirac_2_e(x, &result));
1608 }
1609
1610 double gsl_sf_fermi_dirac_int(const int j, const double x)
1611 {
1612   EVAL_RESULT(gsl_sf_fermi_dirac_int_e(j, x, &result));
1613 }
1614
1615 double gsl_sf_fermi_dirac_mhalf(const double x)
1616 {
1617   EVAL_RESULT(gsl_sf_fermi_dirac_mhalf_e(x, &result));
1618 }
1619
1620 double gsl_sf_fermi_dirac_half(const double x)
1621 {
1622   EVAL_RESULT(gsl_sf_fermi_dirac_half_e(x, &result));
1623 }
1624
1625 double gsl_sf_fermi_dirac_3half(const double x)
1626 {
1627   EVAL_RESULT(gsl_sf_fermi_dirac_3half_e(x, &result));
1628 }
1629
1630 double gsl_sf_fermi_dirac_inc_0(const double x, const double b)
1631 {
1632   EVAL_RESULT(gsl_sf_fermi_dirac_inc_0_e(x, b, &result));
1633 }