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

   1 /* specfunc/psi.c
   2  *
   3  * Copyright (C) 2007 Brian Gough
   4  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004, 2005, 2006 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_zeta.h>
  29 #include <gsl/gsl_sf_psi.h>
  30 #include <gsl/gsl_complex_math.h>
  31
  32 #include <stdio.h>
  33
  34 #include "error.h"
  35
  36 #include "chebyshev.h"
  37 #include "cheb_eval.c"
  38
  39 /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  40
  41
  42 /* Chebyshev fit for f(y) = Re(Psi(1+Iy)) + M_EULER - y^2/(1+y^2) - y^2/(2(4+y^2))
  43  * 1 < y < 10
  44  *   ==>
  45  * y(x) = (9x + 11)/2,  -1 < x < 1
  46  * x(y) = (2y - 11)/9
  47  *
  48  * g(x) := f(y(x))
  49  */
  50 static double r1py_data[] = {
  51    1.59888328244976954803168395603,
  52    0.67905625353213463845115658455,
  53   -0.068485802980122530009506482524,
  54   -0.005788184183095866792008831182,
  55    0.008511258167108615980419855648,
  56   -0.004042656134699693434334556409,
  57    0.001352328406159402601778462956,
  58   -0.000311646563930660566674525382,
  59    0.000018507563785249135437219139,
  60    0.000028348705427529850296492146,
  61   -0.000019487536014574535567541960,
  62    8.0709788710834469408621587335e-06,
  63   -2.2983564321340518037060346561e-06,
  64    3.0506629599604749843855962658e-07,
  65    1.3042238632418364610774284846e-07,
  66   -1.2308657181048950589464690208e-07,
  67    5.7710855710682427240667414345e-08,
  68   -1.8275559342450963966092636354e-08,
  69    3.1020471300626589420759518930e-09,
  70    6.8989327480593812470039430640e-10,
  71   -8.7182290258923059852334818997e-10,
  72    4.4069147710243611798213548777e-10,
  73   -1.4727311099198535963467200277e-10,
  74    2.7589682523262644748825844248e-11,
  75    4.1871826756975856411554363568e-12,
  76   -6.5673460487260087541400767340e-12,
  77    3.4487900886723214020103638000e-12,
  78   -1.1807251417448690607973794078e-12,
  79    2.3798314343969589258709315574e-13,
  80    2.1663630410818831824259465821e-15
  81 };
  82 static cheb_series r1py_cs = {
  83   r1py_data,
  84   29,
  85   -1,1,
  86   18
  87 };
  88
  89
  90 /* Chebyshev fits from SLATEC code for psi(x)
  91
  92  Series for PSI        on the interval  0.         to  1.00000D+00
  93                                        with weighted error   2.03E-17
  94                                         log weighted error  16.69
  95                               significant figures required  16.39
  96                                    decimal places required  17.37
  97
  98  Series for APSI       on the interval  0.         to  2.50000D-01
  99                                        with weighted error   5.54E-17
 100                                         log weighted error  16.26
 101                               significant figures required  14.42
 102                                    decimal places required  16.86
 103
 104 */
 105
 106 static double psics_data[23] = {
 107   -.038057080835217922,
 108    .491415393029387130,
 109   -.056815747821244730,
 110    .008357821225914313,
 111   -.001333232857994342,
 112    .000220313287069308,
 113   -.000037040238178456,
 114    .000006283793654854,
 115   -.000001071263908506,
 116    .000000183128394654,
 117   -.000000031353509361,
 118    .000000005372808776,
 119   -.000000000921168141,
 120    .000000000157981265,
 121   -.000000000027098646,
 122    .000000000004648722,
 123   -.000000000000797527,
 124    .000000000000136827,
 125   -.000000000000023475,
 126    .000000000000004027,
 127   -.000000000000000691,
 128    .000000000000000118,
 129   -.000000000000000020
 130 };
 131 static cheb_series psi_cs = {
 132   psics_data,
 133   22,
 134   -1, 1,
 135   17
 136 };
 137
 138 static double apsics_data[16] = {
 139   -.0204749044678185,
 140   -.0101801271534859,
 141    .0000559718725387,
 142   -.0000012917176570,
 143    .0000000572858606,
 144   -.0000000038213539,
 145    .0000000003397434,
 146   -.0000000000374838,
 147    .0000000000048990,
 148   -.0000000000007344,
 149    .0000000000001233,
 150   -.0000000000000228,
 151    .0000000000000045,
 152   -.0000000000000009,
 153    .0000000000000002,
 154   -.0000000000000000
 155 };
 156 static cheb_series apsi_cs = {
 157   apsics_data,
 158   15,
 159   -1, 1,
 160   9
 161 };
 162
 163 #define PSI_TABLE_NMAX 100
 164 static double psi_table[PSI_TABLE_NMAX+1] = {
 165   0.0,  /* Infinity */              /* psi(0) */
 166  -M_EULER,                          /* psi(1) */
 167   0.42278433509846713939348790992,  /* ...    */
 168   0.92278433509846713939348790992,
 169   1.25611766843180047272682124325,
 170   1.50611766843180047272682124325,
 171   1.70611766843180047272682124325,
 172   1.87278433509846713939348790992,
 173   2.01564147795560999653634505277,
 174   2.14064147795560999653634505277,
 175   2.25175258906672110764745616389,
 176   2.35175258906672110764745616389,
 177   2.44266167997581201673836525479,
 178   2.52599501330914535007169858813,
 179   2.60291809023222227314862166505,
 180   2.67434666166079370172005023648,
 181   2.74101332832746036838671690315,
 182   2.80351332832746036838671690315,
 183   2.86233685773922507426906984432,
 184   2.91789241329478062982462539988,
 185   2.97052399224214905087725697883,
 186   3.02052399224214905087725697883,
 187   3.06814303986119666992487602645,
 188   3.11359758531574212447033057190,
 189   3.15707584618530734186163491973,
 190   3.1987425128519740085283015864,
 191   3.2387425128519740085283015864,
 192   3.2772040513135124700667631249,
 193   3.3142410883505495071038001619,
 194   3.3499553740648352213895144476,
 195   3.3844381326855248765619282407,
 196   3.4177714660188582098952615740,
 197   3.4500295305349872421533260902,
 198   3.4812795305349872421533260902,
 199   3.5115825608380175451836291205,
 200   3.5409943255438998981248055911,
 201   3.5695657541153284695533770196,
 202   3.5973435318931062473311547974,
 203   3.6243705589201332743581818244,
 204   3.6506863483938174848844976139,
 205   3.6763273740348431259101386396,
 206   3.7013273740348431259101386396,
 207   3.7257176179372821503003825420,
 208   3.7495271417468059598241920658,
 209   3.7727829557002943319172153216,
 210   3.7955102284275670591899425943,
 211   3.8177324506497892814121648166,
 212   3.8394715810845718901078169905,
 213   3.8607481768292527411716467777,
 214   3.8815815101625860745049801110,
 215   3.9019896734278921969539597029,
 216   3.9219896734278921969539597029,
 217   3.9415975165651470989147440166,
 218   3.9608282857959163296839747858,
 219   3.9796962103242182164764276160,
 220   3.9982147288427367349949461345,
 221   4.0163965470245549168131279527,
 222   4.0342536898816977739559850956,
 223   4.0517975495308205809735289552,
 224   4.0690389288411654085597358518,
 225   4.0859880813835382899156680552,
 226   4.1026547480502049565823347218,
 227   4.1190481906731557762544658694,
 228   4.1351772229312202923834981274,
 229   4.1510502388042361653993711433,
 230   4.1666752388042361653993711433,
 231   4.1820598541888515500147557587,
 232   4.1972113693403667015299072739,
 233   4.2121367424746950597388624977,
 234   4.2268426248276362362094507330,
 235   4.2413353784508246420065521823,
 236   4.2556210927365389277208378966,
 237   4.2697055997787924488475984600,
 238   4.2835944886676813377364873489,
 239   4.2972931188046676391063503626,
 240   4.3108066323181811526198638761,
 241   4.3241399656515144859531972094,
 242   4.3372978603883565912163551041,
 243   4.3502848733753695782293421171,
 244   4.3631053861958823987421626300,
 245   4.3757636140439836645649474401,
 246   4.3882636140439836645649474401,
 247   4.4006092930563293435772931191,
 248   4.4128044150075488557724150703,
 249   4.4248526077786331931218126607,
 250   4.4367573696833950978837174226,
 251   4.4485220755657480390601880108,
 252   4.4601499825424922251066996387,
 253   4.4716442354160554434975042364,
 254   4.4830078717796918071338678728,
 255   4.4942438268358715824147667492,
 256   4.5053549379469826935258778603,
 257   4.5163439489359936825368668713,
 258   4.5272135141533849868846929582,
 259   4.5379662023254279976373811303,
 260   4.5486045001977684231692960239,
 261   4.5591308159872421073798223397,
 262   4.5695474826539087740464890064,
 263   4.5798567610044242379640147796,
 264   4.5900608426370772991885045755,
 265   4.6001618527380874001986055856
 266 };
 267
 268
 269 #define PSI_1_TABLE_NMAX 100
 270 static double psi_1_table[PSI_1_TABLE_NMAX+1] = {
 271   0.0,  /* Infinity */              /* psi(1,0) */
 272   M_PI*M_PI/6.0,                    /* psi(1,1) */
 273   0.644934066848226436472415,       /* ...      */
 274   0.394934066848226436472415,
 275   0.2838229557371153253613041,
 276   0.2213229557371153253613041,
 277   0.1813229557371153253613041,
 278   0.1535451779593375475835263,
 279   0.1331370146940314251345467,
 280   0.1175120146940314251345467,
 281   0.1051663356816857461222010,
 282   0.0951663356816857461222010,
 283   0.0869018728717683907503002,
 284   0.0799574284273239463058557,
 285   0.0740402686640103368384001,
 286   0.0689382278476838062261552,
 287   0.0644937834032393617817108,
 288   0.0605875334032393617817108,
 289   0.0571273257907826143768665,
 290   0.0540409060376961946237801,
 291   0.0512708229352031198315363,
 292   0.0487708229352031198315363,
 293   0.0465032492390579951149830,
 294   0.0444371335365786562720078,
 295   0.0425467743683366902984728,
 296   0.0408106632572255791873617,
 297   0.0392106632572255791873617,
 298   0.0377313733163971768204978,
 299   0.0363596312039143235969038,
 300   0.0350841209998326909438426,
 301   0.0338950603577399442137594,
 302   0.0327839492466288331026483,
 303   0.0317433665203020901265817,
 304   0.03076680402030209012658168,
 305   0.02984853037475571730748159,
 306   0.02898347847164153045627052,
 307   0.02816715194102928555831133,
 308   0.02739554700275768062003973,
 309   0.02666508681283803124093089,
 310   0.02597256603721476254286995,
 311   0.02531510384129102815759710,
 312   0.02469010384129102815759710,
 313   0.02409521984367056414807896,
 314   0.02352832641963428296894063,
 315   0.02298749353699501850166102,
 316   0.02247096461137518379091722,
 317   0.02197713745088135663042339,
 318   0.02150454765882086513703965,
 319   0.02105185413233829383780923,
 320   0.02061782635456051606003145,
 321   0.02020133322669712580597065,
 322   0.01980133322669712580597065,
 323   0.01941686571420193164987683,
 324   0.01904704322899483105816086,
 325   0.01869104465298913508094477,
 326   0.01834810912486842177504628,
 327   0.01801753061247172756017024,
 328   0.01769865306145131939690494,
 329   0.01739086605006319997554452,
 330   0.01709360088954001329302371,
 331   0.01680632711763538818529605,
 332   0.01652854933985761040751827,
 333   0.01625980437882562975715546,
 334   0.01599965869724394401313881,
 335   0.01574770606433893015574400,
 336   0.01550356543933893015574400,
 337   0.01526687904880638577704578,
 338   0.01503731063741979257227076,
 339   0.01481454387422086185273411,
 340   0.01459828089844231513993134,
 341   0.01438824099085987447620523,
 342   0.01418415935820681325171544,
 343   0.01398578601958352422176106,
 344   0.01379288478501562298719316,
 345   0.01360523231738567365335942,
 346   0.01342261726990576130858221,
 347   0.01324483949212798353080444,
 348   0.01307170929822216635628920,
 349   0.01290304679189732236910755,
 350   0.01273868124291638877278934,
 351   0.01257845051066194236996928,
 352   0.01242220051066194236996928,
 353   0.01226978472038606978956995,
 354   0.01212106372098095378719041,
 355   0.01197590477193174490346273,
 356   0.01183418141592267460867815,
 357   0.01169577311142440471248438,
 358   0.01156056489076458859566448,
 359   0.01142844704164317229232189,
 360   0.01129931481023821361463594,
 361   0.01117306812421372175754719,
 362   0.01104961133409026496742374,
 363   0.01092885297157366069257770,
 364   0.01081070552355853781923177,
 365   0.01069508522063334415522437,
 366   0.01058191183901270133041676,
 367   0.01047110851491297833872701,
 368   0.01036260157046853389428257,
 369   0.01025632035036012704977199,  /* ...        */
 370   0.01015219706839427948625679,  /* psi(1,99)  */
 371   0.01005016666333357139524567   /* psi(1,100) */
 372 };
 373
 374
 375 /* digamma for x both positive and negative; we do both
 376  * cases here because of the way we use even/odd parts
 377  * of the function
 378  */
 379 static int
 380 psi_x(const double x, gsl_sf_result * result)
 381 {
 382   const double y = fabs(x);
 383
 384   if(x == 0.0 || x == -1.0 || x == -2.0) {
 385     DOMAIN_ERROR(result);
 386   }
 387   else if(y >= 2.0) {
 388     const double t = 8.0/(y*y)-1.0;
 389     gsl_sf_result result_c;
 390     cheb_eval_e(&apsi_cs, t, &result_c);
 391     if(x < 0.0) {
 392       const double s = sin(M_PI*x);
 393       const double c = cos(M_PI*x);
 394       if(fabs(s) < 2.0*GSL_SQRT_DBL_MIN) {
 395         DOMAIN_ERROR(result);
 396       }
 397       else {
 398         result->val  = log(y) - 0.5/x + result_c.val - M_PI * c/s;
 399         result->err  = M_PI*fabs(x)*GSL_DBL_EPSILON/(s*s);
 400         result->err += result_c.err;
 401         result->err += GSL_DBL_EPSILON * fabs(result->val);
 402         return GSL_SUCCESS;
 403       }
 404     }
 405     else {
 406       result->val  = log(y) - 0.5/x + result_c.val;
 407       result->err  = result_c.err;
 408       result->err += GSL_DBL_EPSILON * fabs(result->val);
 409       return GSL_SUCCESS;
 410     }
 411   }
 412   else { /* -2 < x < 2 */
 413     gsl_sf_result result_c;
 414
 415     if(x < -1.0) { /* x = -2 + v */
 416       const double v  = x + 2.0;
 417       const double t1 = 1.0/x;
 418       const double t2 = 1.0/(x+1.0);
 419       const double t3 = 1.0/v;
 420       cheb_eval_e(&psi_cs, 2.0*v-1.0, &result_c);
 421
 422       result->val  = -(t1 + t2 + t3) + result_c.val;
 423       result->err  = GSL_DBL_EPSILON * (fabs(t1) + fabs(x/(t2*t2)) + fabs(x/(t3*t3)));
 424       result->err += result_c.err;
 425       result->err += GSL_DBL_EPSILON * fabs(result->val);
 426       return GSL_SUCCESS;
 427     }
 428     else if(x < 0.0) { /* x = -1 + v */
 429       const double v  = x + 1.0;
 430       const double t1 = 1.0/x;
 431       const double t2 = 1.0/v;
 432       cheb_eval_e(&psi_cs, 2.0*v-1.0, &result_c);
 433
 434       result->val  = -(t1 + t2) + result_c.val;
 435       result->err  = GSL_DBL_EPSILON * (fabs(t1) + fabs(x/(t2*t2)));
 436       result->err += result_c.err;
 437       result->err += GSL_DBL_EPSILON * fabs(result->val);
 438       return GSL_SUCCESS;
 439     }
 440     else if(x < 1.0) { /* x = v */
 441       const double t1 = 1.0/x;
 442       cheb_eval_e(&psi_cs, 2.0*x-1.0, &result_c);
 443
 444       result->val  = -t1 + result_c.val;
 445       result->err  = GSL_DBL_EPSILON * t1;
 446       result->err += result_c.err;
 447       result->err += GSL_DBL_EPSILON * fabs(result->val);
 448       return GSL_SUCCESS;
 449     }
 450     else { /* x = 1 + v */
 451       const double v = x - 1.0;
 452       return cheb_eval_e(&psi_cs, 2.0*v-1.0, result);
 453     }
 454   }
 455 }
 456
 457
 458 /* psi(z) for large |z| in the right half-plane; [Abramowitz + Stegun, 6.3.18] */
 459 static
 460 gsl_complex
 461 psi_complex_asymp(gsl_complex z)
 462 {
 463   /* coefficients in the asymptotic expansion for large z;
 464    * let w = z^(-2) and write the expression in the form
 465    *
 466    *   ln(z) - 1/(2z) - 1/12 w (1 + c1 w + c2 w + c3 w + ... )
 467    */
 468   static const double c1 = -0.1;
 469   static const double c2 =  1.0/21.0;
 470   static const double c3 = -0.05;
 471
 472   gsl_complex zi = gsl_complex_inverse(z);
 473   gsl_complex w  = gsl_complex_mul(zi, zi);
 474   gsl_complex cs;
 475
 476   /* Horner method evaluation of term in parentheses */
 477   gsl_complex sum;
 478   sum = gsl_complex_mul_real(w, c3/c2);
 479   sum = gsl_complex_add_real(sum, 1.0);
 480   sum = gsl_complex_mul_real(sum, c2/c1);
 481   sum = gsl_complex_mul(sum, w);
 482   sum = gsl_complex_add_real(sum, 1.0);
 483   sum = gsl_complex_mul_real(sum, c1);
 484   sum = gsl_complex_mul(sum, w);
 485   sum = gsl_complex_add_real(sum, 1.0);
 486
 487   /* correction added to log(z) */
 488   cs = gsl_complex_mul(sum, w);
 489   cs = gsl_complex_mul_real(cs, -1.0/12.0);
 490   cs = gsl_complex_add(cs, gsl_complex_mul_real(zi, -0.5));
 491
 492   return gsl_complex_add(gsl_complex_log(z), cs);
 493 }
 494
 495
 496
 497 /* psi(z) for complex z in the right half-plane */
 498 static int
 499 psi_complex_rhp(
 500   gsl_complex z,
 501   gsl_sf_result * result_re,
 502   gsl_sf_result * result_im
 503   )
 504 {
 505   int n_recurse = 0;
 506   int i;
 507   gsl_complex a;
 508
 509   if(GSL_REAL(z) == 0.0 && GSL_IMAG(z) == 0.0)
 510   {
 511     result_re->val = 0.0;
 512     result_im->val = 0.0;
 513     result_re->err = 0.0;
 514     result_im->err = 0.0;
 515     return GSL_EDOM;
 516   }
 517
 518   /* compute the number of recurrences to apply */
 519   if(GSL_REAL(z) < 20.0 && fabs(GSL_IMAG(z)) < 20.0)
 520   {
 521     const double sp = sqrt(20.0 + GSL_IMAG(z));
 522     const double sn = sqrt(20.0 - GSL_IMAG(z));
 523     const double rhs = sp*sn - GSL_REAL(z);
 524     if(rhs > 0.0) n_recurse = ceil(rhs);
 525   }
 526
 527   /* compute asymptotic at the large value z + n_recurse */
 528   a = psi_complex_asymp(gsl_complex_add_real(z, n_recurse));
 529
 530   result_re->err = 2.0 * GSL_DBL_EPSILON * fabs(GSL_REAL(a));
 531   result_im->err = 2.0 * GSL_DBL_EPSILON * fabs(GSL_IMAG(a));
 532
 533   /* descend recursively, if necessary */
 534   for(i = n_recurse; i >= 1; --i)
 535   {
 536     gsl_complex zn = gsl_complex_add_real(z, i - 1.0);
 537     gsl_complex zn_inverse = gsl_complex_inverse(zn);
 538     a = gsl_complex_sub(a, zn_inverse);
 539
 540     /* accumulate the error, to catch cancellations */
 541     result_re->err += 2.0 * GSL_DBL_EPSILON * fabs(GSL_REAL(zn_inverse));
 542     result_im->err += 2.0 * GSL_DBL_EPSILON * fabs(GSL_IMAG(zn_inverse));
 543   }
 544
 545   result_re->val = GSL_REAL(a);
 546   result_im->val = GSL_IMAG(a);
 547
 548   result_re->err += 2.0 * GSL_DBL_EPSILON * fabs(result_re->val);
 549   result_im->err += 2.0 * GSL_DBL_EPSILON * fabs(result_im->val);
 550
 551   return GSL_SUCCESS;
 552 }
 553
 554
 555
 556 /* generic polygamma; assumes n >= 0 and x > 0
 557  */
 558 static int
 559 psi_n_xg0(const int n, const double x, gsl_sf_result * result)
 560 {
 561   if(n == 0) {
 562     return gsl_sf_psi_e(x, result);
 563   }
 564   else {
 565     /* Abramowitz + Stegun 6.4.10 */
 566     gsl_sf_result ln_nf;
 567     gsl_sf_result hzeta;
 568     int stat_hz = gsl_sf_hzeta_e(n+1.0, x, &hzeta);
 569     int stat_nf = gsl_sf_lnfact_e((unsigned int) n, &ln_nf);
 570     int stat_e  = gsl_sf_exp_mult_err_e(ln_nf.val, ln_nf.err,
 571                                            hzeta.val, hzeta.err,
 572                                            result);
 573     if(GSL_IS_EVEN(n)) result->val = -result->val;
 574     return GSL_ERROR_SELECT_3(stat_e, stat_nf, stat_hz);
 575   }
 576 }
 577
 578
 579
 580 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 581
 582 int gsl_sf_psi_int_e(const int n, gsl_sf_result * result)
 583 {
 584   /* CHECK_POINTER(result) */
 585
 586   if(n <= 0) {
 587     DOMAIN_ERROR(result);
 588   }
 589   else if(n <= PSI_TABLE_NMAX) {
 590     result->val = psi_table[n];
 591     result->err = GSL_DBL_EPSILON * fabs(result->val);
 592     return GSL_SUCCESS;
 593   }
 594   else {
 595     /* Abramowitz+Stegun 6.3.18 */
 596     const double c2 = -1.0/12.0;
 597     const double c3 =  1.0/120.0;
 598     const double c4 = -1.0/252.0;
 599     const double c5 =  1.0/240.0;
 600     const double ni2 = (1.0/n)*(1.0/n);
 601     const double ser = ni2 * (c2 + ni2 * (c3 + ni2 * (c4 + ni2*c5)));
 602     result->val  = log(n) - 0.5/n + ser;
 603     result->err  = GSL_DBL_EPSILON * (fabs(log(n)) + fabs(0.5/n) + fabs(ser));
 604     result->err += GSL_DBL_EPSILON * fabs(result->val);
 605     return GSL_SUCCESS;
 606   }
 607 }
 608
 609
 610 int gsl_sf_psi_e(const double x, gsl_sf_result * result)
 611 {
 612   /* CHECK_POINTER(result) */
 613   return psi_x(x, result);
 614 }
 615
 616
 617 int
 618 gsl_sf_psi_1piy_e(const double y, gsl_sf_result * result)
 619 {
 620   const double ay = fabs(y);
 621
 622   /* CHECK_POINTER(result) */
 623
 624   if(ay > 1000.0) {
 625     /* [Abramowitz+Stegun, 6.3.19] */
 626     const double yi2 = 1.0/(ay*ay);
 627     const double lny = log(ay);
 628     const double sum = yi2 * (1.0/12.0 + 1.0/120.0 * yi2 + 1.0/252.0 * yi2*yi2);
 629     result->val = lny + sum;
 630     result->err = 2.0 * GSL_DBL_EPSILON * (fabs(lny) + fabs(sum));
 631     return GSL_SUCCESS;
 632   }
 633   else if(ay > 10.0) {
 634     /* [Abramowitz+Stegun, 6.3.19] */
 635     const double yi2 = 1.0/(ay*ay);
 636     const double lny = log(ay);
 637     const double sum = yi2 * (1.0/12.0 +
 638                          yi2 * (1.0/120.0 +
 639                            yi2 * (1.0/252.0 +
 640                              yi2 * (1.0/240.0 +
 641                                yi2 * (1.0/132.0 + 691.0/32760.0 * yi2)))));
 642     result->val = lny + sum;
 643     result->err = 2.0 * GSL_DBL_EPSILON * (fabs(lny) + fabs(sum));
 644     return GSL_SUCCESS;
 645   }
 646   else if(ay > 1.0){
 647     const double y2 = ay*ay;
 648     const double x  = (2.0*ay - 11.0)/9.0;
 649     const double v  = y2*(1.0/(1.0+y2) + 0.5/(4.0+y2));
 650     gsl_sf_result result_c;
 651     cheb_eval_e(&r1py_cs, x, &result_c);
 652     result->val  = result_c.val - M_EULER + v;
 653     result->err  = result_c.err;
 654     result->err += 2.0 * GSL_DBL_EPSILON * (fabs(v) + M_EULER + fabs(result_c.val));
 655     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 656     result->err *= 5.0; /* FIXME: losing a digit somewhere... maybe at x=... ? */
 657     return GSL_SUCCESS;
 658   }
 659   else {
 660     /* [Abramowitz+Stegun, 6.3.17]
 661      *
 662      * -M_EULER + y^2 Sum[1/n 1/(n^2 + y^2), {n,1,M}]
 663      *   +     Sum[1/n^3, {n,M+1,Infinity}]
 664      *   - y^2 Sum[1/n^5, {n,M+1,Infinity}]
 665      *   + y^4 Sum[1/n^7, {n,M+1,Infinity}]
 666      *   - y^6 Sum[1/n^9, {n,M+1,Infinity}]
 667      *   + O(y^8)
 668      *
 669      * We take M=50 for at least 15 digit precision.
 670      */
 671     const int M = 50;
 672     const double y2 = y*y;
 673     const double c0 = 0.00019603999466879846570;
 674     const double c2 = 3.8426659205114376860e-08;
 675     const double c4 = 1.0041592839497643554e-11;
 676     const double c6 = 2.9516743763500191289e-15;
 677     const double p  = c0 + y2 *(-c2 + y2*(c4 - y2*c6));
 678     double sum = 0.0;
 679     double v;
 680
 681     int n;
 682     for(n=1; n<=M; n++) {
 683       sum += 1.0/(n * (n*n + y*y));
 684     }
 685
 686     v = y2 * (sum + p);
 687     result->val  = -M_EULER + v;
 688     result->err  = GSL_DBL_EPSILON * (M_EULER + fabs(v));
 689     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 690     return GSL_SUCCESS;
 691   }
 692 }
 693
 694
 695 int gsl_sf_psi_1_int_e(const int n, gsl_sf_result * result)
 696 {
 697   /* CHECK_POINTER(result) */
 698   if(n <= 0) {
 699     DOMAIN_ERROR(result);
 700   }
 701   else if(n <= PSI_1_TABLE_NMAX) {
 702     result->val = psi_1_table[n];
 703     result->err = GSL_DBL_EPSILON * result->val;
 704     return GSL_SUCCESS;
 705   }
 706   else {
 707     /* Abramowitz+Stegun 6.4.12
 708      * double-precision for n > 100
 709      */
 710     const double c0 = -1.0/30.0;
 711     const double c1 =  1.0/42.0;
 712     const double c2 = -1.0/30.0;
 713     const double ni2 = (1.0/n)*(1.0/n);
 714     const double ser =  ni2*ni2 * (c0 + ni2*(c1 + c2*ni2));
 715     result->val = (1.0 + 0.5/n + 1.0/(6.0*n*n) + ser) / n;
 716     result->err = GSL_DBL_EPSILON * result->val;
 717     return GSL_SUCCESS;
 718   }
 719 }
 720
 721
 722 int gsl_sf_psi_1_e(const double x, gsl_sf_result * result)
 723 {
 724   /* CHECK_POINTER(result) */
 725
 726   if(x == 0.0 || x == -1.0 || x == -2.0) {
 727     DOMAIN_ERROR(result);
 728   }
 729   else if(x > 0.0)
 730   {
 731     return psi_n_xg0(1, x, result);
 732   }
 733   else if(x > -5.0)
 734   {
 735     /* Abramowitz + Stegun 6.4.6 */
 736     int M = -floor(x);
 737     double fx = x + M;
 738     double sum = 0.0;
 739     int m;
 740
 741     if(fx == 0.0)
 742       DOMAIN_ERROR(result);
 743
 744     for(m = 0; m < M; ++m)
 745       sum += 1.0/((x+m)*(x+m));
 746
 747     {
 748       int stat_psi = psi_n_xg0(1, fx, result);
 749       result->val += sum;
 750       result->err += M * GSL_DBL_EPSILON * sum;
 751       return stat_psi;
 752     }
 753   }
 754   else
 755   {
 756     /* Abramowitz + Stegun 6.4.7 */
 757     const double sin_px = sin(M_PI * x);
 758     const double d = M_PI*M_PI/(sin_px*sin_px);
 759     gsl_sf_result r;
 760     int stat_psi = psi_n_xg0(1, 1.0-x, &r);
 761     result->val = d - r.val;
 762     result->err = r.err + 2.0*GSL_DBL_EPSILON*d;
 763     return stat_psi;
 764   }
 765 }
 766
 767
 768 int gsl_sf_psi_n_e(const int n, const double x, gsl_sf_result * result)
 769 {
 770   /* CHECK_POINTER(result) */
 771
 772   if(n == 0)
 773   {
 774     return gsl_sf_psi_e(x, result);
 775   }
 776   else if(n == 1)
 777   {
 778     return gsl_sf_psi_1_e(x, result);
 779   }
 780   else if(n < 0 || x <= 0.0) {
 781     DOMAIN_ERROR(result);
 782   }
 783   else {
 784     gsl_sf_result ln_nf;
 785     gsl_sf_result hzeta;
 786     int stat_hz = gsl_sf_hzeta_e(n+1.0, x, &hzeta);
 787     int stat_nf = gsl_sf_lnfact_e((unsigned int) n, &ln_nf);
 788     int stat_e  = gsl_sf_exp_mult_err_e(ln_nf.val, ln_nf.err,
 789                                            hzeta.val, hzeta.err,
 790                                            result);
 791     if(GSL_IS_EVEN(n)) result->val = -result->val;
 792     return GSL_ERROR_SELECT_3(stat_e, stat_nf, stat_hz);
 793   }
 794 }
 795
 796
 797 int
 798 gsl_sf_complex_psi_e(
 799   const double x,
 800   const double y,
 801   gsl_sf_result * result_re,
 802   gsl_sf_result * result_im
 803   )
 804 {
 805   if(x >= 0.0)
 806   {
 807     gsl_complex z = gsl_complex_rect(x, y);
 808     return psi_complex_rhp(z, result_re, result_im);
 809   }
 810   else
 811   {
 812     /* reflection formula [Abramowitz+Stegun, 6.3.7] */
 813     gsl_complex z = gsl_complex_rect(x, y);
 814     gsl_complex omz = gsl_complex_rect(1.0 - x, -y);
 815     gsl_complex zpi = gsl_complex_mul_real(z, M_PI);
 816     gsl_complex cotzpi = gsl_complex_cot(zpi);
 817     int ret_val = psi_complex_rhp(omz, result_re, result_im);
 818
 819     if(GSL_IS_REAL(GSL_REAL(cotzpi)) && GSL_IS_REAL(GSL_IMAG(cotzpi)))
 820     {
 821       result_re->val -= M_PI * GSL_REAL(cotzpi);
 822       result_im->val -= M_PI * GSL_IMAG(cotzpi);
 823       return ret_val;
 824     }
 825     else
 826     {
 827       GSL_ERROR("singularity", GSL_EDOM);
 828     }
 829   }
 830 }
 831
 832
 833
 834 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 835
 836 #include "eval.h"
 837
 838 double gsl_sf_psi_int(const int n)
 839 {
 840   EVAL_RESULT(gsl_sf_psi_int_e(n, &result));
 841 }
 842
 843 double gsl_sf_psi(const double x)
 844 {
 845   EVAL_RESULT(gsl_sf_psi_e(x, &result));
 846 }
 847
 848 double gsl_sf_psi_1piy(const double x)
 849 {
 850   EVAL_RESULT(gsl_sf_psi_1piy_e(x, &result));
 851 }
 852
 853 double gsl_sf_psi_1_int(const int n)
 854 {
 855   EVAL_RESULT(gsl_sf_psi_1_int_e(n, &result));
 856 }
 857
 858 double gsl_sf_psi_1(const double x)
 859 {
 860   EVAL_RESULT(gsl_sf_psi_1_e(x, &result));
 861 }
 862
 863 double gsl_sf_psi_n(const int n, const double x)
 864 {
 865   EVAL_RESULT(gsl_sf_psi_n_e(n, x, &result));
 866 }