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

   1 /* specfunc/dilog.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 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_clausen.h>
  26 #include <gsl/gsl_sf_trig.h>
  27 #include <gsl/gsl_sf_log.h>
  28 #include <gsl/gsl_sf_dilog.h>
  29
  30
  31 /* Evaluate series for real dilog(x)
  32  * Sum[ x^k / k^2, {k,1,Infinity}]
  33  *
  34  * Converges rapidly for |x| < 1/2.
  35  */
  36 static
  37 int
  38 dilog_series_1(const double x, gsl_sf_result * result)
  39 {
  40   const int kmax = 1000;
  41   double sum  = x;
  42   double term = x;
  43   int k;
  44   for(k=2; k<kmax; k++) {
  45     const double rk = (k-1.0)/k;
  46     term *= x;
  47     term *= rk*rk;
  48     sum += term;
  49     if(fabs(term/sum) < GSL_DBL_EPSILON) break;
  50   }
  51
  52   result->val  = sum;
  53   result->err  = 2.0 * fabs(term);
  54   result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  55
  56   if(k == kmax)
  57     GSL_ERROR ("error", GSL_EMAXITER);
  58   else
  59     return GSL_SUCCESS;
  60 }
  61
  62
  63 /* Compute the associated series
  64  *
  65  *   sum_{k=1}{infty} r^k / (k^2 (k+1))
  66  *
  67  * This is a series which appears in the one-step accelerated
  68  * method, which splits out one elementary function from the
  69  * full definition of Li_2(x). See below.
  70  */
  71 static int
  72 series_2(double r, gsl_sf_result * result)
  73 {
  74   static const int kmax = 100;
  75   double rk = r;
  76   double sum = 0.5 * r;
  77   int k;
  78   for(k=2; k<10; k++)
  79   {
  80     double ds;
  81     rk *= r;
  82     ds = rk/(k*k*(k+1.0));
  83     sum += ds;
  84   }
  85   for(; k<kmax; k++)
  86   {
  87     double ds;
  88     rk *= r;
  89     ds = rk/(k*k*(k+1.0));
  90     sum += ds;
  91     if(fabs(ds/sum) < 0.5*GSL_DBL_EPSILON) break;
  92   }
  93
  94   result->val = sum;
  95   result->err = 2.0 * kmax * GSL_DBL_EPSILON * fabs(sum);
  96
  97   return GSL_SUCCESS;
  98 }
  99
 100
 101 /* Compute Li_2(x) using the accelerated series representation.
 102  *
 103  * Li_2(x) = 1 + (1-x)ln(1-x)/x + series_2(x)
 104  *
 105  * assumes: -1 < x < 1
 106  */
 107 static int
 108 dilog_series_2(double x, gsl_sf_result * result)
 109 {
 110   const int stat_s3 = series_2(x, result);
 111   double t;
 112   if(x > 0.01)
 113     t = (1.0 - x) * log(1.0-x) / x;
 114   else
 115   {
 116     static const double c3 = 1.0/3.0;
 117     static const double c4 = 1.0/4.0;
 118     static const double c5 = 1.0/5.0;
 119     static const double c6 = 1.0/6.0;
 120     static const double c7 = 1.0/7.0;
 121     static const double c8 = 1.0/8.0;
 122     const double t68 = c6 + x*(c7 + x*c8);
 123     const double t38 = c3 + x *(c4 + x *(c5 + x * t68));
 124     t = (x - 1.0) * (1.0 + x*(0.5 + x*t38));
 125   }
 126   result->val += 1.0 + t;
 127   result->err += 2.0 * GSL_DBL_EPSILON * fabs(t);
 128   return stat_s3;
 129 }
 130
 131
 132 /* Calculates Li_2(x) for real x. Assumes x >= 0.0.
 133  */
 134 static
 135 int
 136 dilog_xge0(const double x, gsl_sf_result * result)
 137 {
 138   if(x > 2.0) {
 139     gsl_sf_result ser;
 140     const int stat_ser = dilog_series_2(1.0/x, &ser);
 141     const double log_x = log(x);
 142     const double t1 = M_PI*M_PI/3.0;
 143     const double t2 = ser.val;
 144     const double t3 = 0.5*log_x*log_x;
 145     result->val  = t1 - t2 - t3;
 146     result->err  = GSL_DBL_EPSILON * fabs(log_x) + ser.err;
 147     result->err += GSL_DBL_EPSILON * (fabs(t1) + fabs(t2) + fabs(t3));
 148     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 149     return stat_ser;
 150   }
 151   else if(x > 1.01) {
 152     gsl_sf_result ser;
 153     const int stat_ser = dilog_series_2(1.0 - 1.0/x, &ser);
 154     const double log_x    = log(x);
 155     const double log_term = log_x * (log(1.0-1.0/x) + 0.5*log_x);
 156     const double t1 = M_PI*M_PI/6.0;
 157     const double t2 = ser.val;
 158     const double t3 = log_term;
 159     result->val  = t1 + t2 - t3;
 160     result->err  = GSL_DBL_EPSILON * fabs(log_x) + ser.err;
 161     result->err += GSL_DBL_EPSILON * (fabs(t1) + fabs(t2) + fabs(t3));
 162     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 163     return stat_ser;
 164   }
 165   else if(x > 1.0) {
 166     /* series around x = 1.0 */
 167     const double eps = x - 1.0;
 168     const double lne = log(eps);
 169     const double c0 = M_PI*M_PI/6.0;
 170     const double c1 =   1.0 - lne;
 171     const double c2 = -(1.0 - 2.0*lne)/4.0;
 172     const double c3 =  (1.0 - 3.0*lne)/9.0;
 173     const double c4 = -(1.0 - 4.0*lne)/16.0;
 174     const double c5 =  (1.0 - 5.0*lne)/25.0;
 175     const double c6 = -(1.0 - 6.0*lne)/36.0;
 176     const double c7 =  (1.0 - 7.0*lne)/49.0;
 177     const double c8 = -(1.0 - 8.0*lne)/64.0;
 178     result->val = c0+eps*(c1+eps*(c2+eps*(c3+eps*(c4+eps*(c5+eps*(c6+eps*(c7+eps*c8)))))));
 179     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 180     return GSL_SUCCESS;
 181   }
 182   else if(x == 1.0) {
 183     result->val = M_PI*M_PI/6.0;
 184     result->err = 2.0 * GSL_DBL_EPSILON * M_PI*M_PI/6.0;
 185     return GSL_SUCCESS;
 186   }
 187   else if(x > 0.5) {
 188     gsl_sf_result ser;
 189     const int stat_ser = dilog_series_2(1.0-x, &ser);
 190     const double log_x = log(x);
 191     const double t1 = M_PI*M_PI/6.0;
 192     const double t2 = ser.val;
 193     const double t3 = log_x*log(1.0-x);
 194     result->val  = t1 - t2 - t3;
 195     result->err  = GSL_DBL_EPSILON * fabs(log_x) + ser.err;
 196     result->err += GSL_DBL_EPSILON * (fabs(t1) + fabs(t2) + fabs(t3));
 197     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 198     return stat_ser;
 199   }
 200   else if(x > 0.25) {
 201     return dilog_series_2(x, result);
 202   }
 203   else if(x > 0.0) {
 204     return dilog_series_1(x, result);
 205   }
 206   else {
 207     /* x == 0.0 */
 208     result->val = 0.0;
 209     result->err = 0.0;
 210     return GSL_SUCCESS;
 211   }
 212 }
 213
 214
 215 /* Evaluate the series representation for Li2(z):
 216  *
 217  *   Li2(z) = Sum[ |z|^k / k^2 Exp[i k arg(z)], {k,1,Infinity}]
 218  *   |z|    = r
 219  *   arg(z) = theta
 220  *
 221  * Assumes 0 < r < 1.
 222  * It is used only for small r.
 223  */
 224 static
 225 int
 226 dilogc_series_1(
 227   const double r,
 228   const double x,
 229   const double y,
 230   gsl_sf_result * real_result,
 231   gsl_sf_result * imag_result
 232   )
 233 {
 234   const double cos_theta = x/r;
 235   const double sin_theta = y/r;
 236   const double alpha = 1.0 - cos_theta;
 237   const double beta  = sin_theta;
 238   double ck = cos_theta;
 239   double sk = sin_theta;
 240   double rk = r;
 241   double real_sum = r*ck;
 242   double imag_sum = r*sk;
 243   const int kmax = 50 + (int)(22.0/(-log(r))); /* tuned for double-precision */
 244   int k;
 245   for(k=2; k<kmax; k++) {
 246     double dr, di;
 247     double ck_tmp = ck;
 248     ck = ck - (alpha*ck + beta*sk);
 249     sk = sk - (alpha*sk - beta*ck_tmp);
 250     rk *= r;
 251     dr = rk/((double)k*k) * ck;
 252     di = rk/((double)k*k) * sk;
 253     real_sum += dr;
 254     imag_sum += di;
 255     if(fabs((dr*dr + di*di)/(real_sum*real_sum + imag_sum*imag_sum)) < GSL_DBL_EPSILON*GSL_DBL_EPSILON) break;
 256   }
 257
 258   real_result->val = real_sum;
 259   real_result->err = 2.0 * kmax * GSL_DBL_EPSILON * fabs(real_sum);
 260   imag_result->val = imag_sum;
 261   imag_result->err = 2.0 * kmax * GSL_DBL_EPSILON * fabs(imag_sum);
 262
 263   return GSL_SUCCESS;
 264 }
 265
 266
 267 /* Compute
 268  *
 269  *   sum_{k=1}{infty} z^k / (k^2 (k+1))
 270  *
 271  * This is a series which appears in the one-step accelerated
 272  * method, which splits out one elementary function from the
 273  * full definition of Li_2.
 274  */
 275 static int
 276 series_2_c(
 277   double r,
 278   double x,
 279   double y,
 280   gsl_sf_result * sum_re,
 281   gsl_sf_result * sum_im
 282   )
 283 {
 284   const double cos_theta = x/r;
 285   const double sin_theta = y/r;
 286   const double alpha = 1.0 - cos_theta;
 287   const double beta  = sin_theta;
 288   double ck = cos_theta;
 289   double sk = sin_theta;
 290   double rk = r;
 291   double real_sum = 0.5 * r*ck;
 292   double imag_sum = 0.5 * r*sk;
 293   const int kmax = 30 + (int)(18.0/(-log(r))); /* tuned for double-precision */
 294   int k;
 295   for(k=2; k<kmax; k++)
 296   {
 297     double dr, di;
 298     const double ck_tmp = ck;
 299     ck = ck - (alpha*ck + beta*sk);
 300     sk = sk - (alpha*sk - beta*ck_tmp);
 301     rk *= r;
 302     dr = rk/((double)k*k*(k+1.0)) * ck;
 303     di = rk/((double)k*k*(k+1.0)) * sk;
 304     real_sum += dr;
 305     imag_sum += di;
 306     if(fabs((dr*dr + di*di)/(real_sum*real_sum + imag_sum*imag_sum)) < GSL_DBL_EPSILON*GSL_DBL_EPSILON) break;
 307   }
 308
 309   sum_re->val = real_sum;
 310   sum_re->err = 2.0 * kmax * GSL_DBL_EPSILON * fabs(real_sum);
 311   sum_im->val = imag_sum;
 312   sum_im->err = 2.0 * kmax * GSL_DBL_EPSILON * fabs(imag_sum);
 313
 314   return GSL_SUCCESS;
 315 }
 316
 317
 318 /* Compute Li_2(z) using the one-step accelerated series.
 319  *
 320  * Li_2(z) = 1 + (1-z)ln(1-z)/z + series_2_c(z)
 321  *
 322  * z = r exp(i theta)
 323  * assumes: r < 1
 324  * assumes: r > epsilon, so that we take no special care with log(1-z)
 325  */
 326 static
 327 int
 328 dilogc_series_2(
 329   const double r,
 330   const double x,
 331   const double y,
 332   gsl_sf_result * real_dl,
 333   gsl_sf_result * imag_dl
 334   )
 335 {
 336   if(r == 0.0)
 337   {
 338     real_dl->val = 0.0;
 339     imag_dl->val = 0.0;
 340     real_dl->err = 0.0;
 341     imag_dl->err = 0.0;
 342     return GSL_SUCCESS;
 343   }
 344   else
 345   {
 346     gsl_sf_result sum_re;
 347     gsl_sf_result sum_im;
 348     const int stat_s3 = series_2_c(r, x, y, &sum_re, &sum_im);
 349
 350     /* t = ln(1-z)/z */
 351     gsl_sf_result ln_omz_r;
 352     gsl_sf_result ln_omz_theta;
 353     const int stat_log = gsl_sf_complex_log_e(1.0-x, -y, &ln_omz_r, &ln_omz_theta);
 354     const double t_x = ( ln_omz_r.val * x + ln_omz_theta.val * y)/(r*r);
 355     const double t_y = (-ln_omz_r.val * y + ln_omz_theta.val * x)/(r*r);
 356
 357     /* r = (1-z) ln(1-z)/z */
 358     const double r_x = (1.0 - x) * t_x + y * t_y;
 359     const double r_y = (1.0 - x) * t_y - y * t_x;
 360
 361     real_dl->val = sum_re.val + r_x + 1.0;
 362     imag_dl->val = sum_im.val + r_y;
 363     real_dl->err = sum_re.err + 2.0*GSL_DBL_EPSILON*(fabs(real_dl->val) + fabs(r_x));
 364     imag_dl->err = sum_im.err + 2.0*GSL_DBL_EPSILON*(fabs(imag_dl->val) + fabs(r_y));
 365     return GSL_ERROR_SELECT_2(stat_s3, stat_log);
 366   }
 367 }
 368
 369
 370 /* Evaluate a series for Li_2(z) when |z| is near 1.
 371  * This is uniformly good away from z=1.
 372  *
 373  *   Li_2(z) = Sum[ a^n/n! H_n(theta), {n, 0, Infinity}]
 374  *
 375  * where
 376  *   H_n(theta) = Sum[ e^(i m theta) m^n / m^2, {m, 1, Infinity}]
 377  *   a = ln(r)
 378  *
 379  *  H_0(t) = Gl_2(t) + i Cl_2(t)
 380  *  H_1(t) = 1/2 ln(2(1-c)) + I atan2(-s, 1-c)
 381  *  H_2(t) = -1/2 + I/2 s/(1-c)
 382  *  H_3(t) = -1/2 /(1-c)
 383  *  H_4(t) = -I/2 s/(1-c)^2
 384  *  H_5(t) = 1/2 (2 + c)/(1-c)^2
 385  *  H_6(t) = I/2 s/(1-c)^5 (8(1-c) - s^2 (3 + c))
 386  */
 387 static
 388 int
 389 dilogc_series_3(
 390   const double r,
 391   const double x,
 392   const double y,
 393   gsl_sf_result * real_result,
 394   gsl_sf_result * imag_result
 395   )
 396 {
 397   const double theta = atan2(y, x);
 398   const double cos_theta = x/r;
 399   const double sin_theta = y/r;
 400   const double a = log(r);
 401   const double omc = 1.0 - cos_theta;
 402   const double omc2 = omc*omc;
 403   double H_re[7];
 404   double H_im[7];
 405   double an, nfact;
 406   double sum_re, sum_im;
 407   gsl_sf_result Him0;
 408   int n;
 409
 410   H_re[0] = M_PI*M_PI/6.0 + 0.25*(theta*theta - 2.0*M_PI*fabs(theta));
 411   gsl_sf_clausen_e(theta, &Him0);
 412   H_im[0] = Him0.val;
 413
 414   H_re[1] = -0.5*log(2.0*omc);
 415   H_im[1] = -atan2(-sin_theta, omc);
 416
 417   H_re[2] = -0.5;
 418   H_im[2] = 0.5 * sin_theta/omc;
 419
 420   H_re[3] = -0.5/omc;
 421   H_im[3] = 0.0;
 422
 423   H_re[4] = 0.0;
 424   H_im[4] = -0.5*sin_theta/omc2;
 425
 426   H_re[5] = 0.5 * (2.0 + cos_theta)/omc2;
 427   H_im[5] = 0.0;
 428
 429   H_re[6] = 0.0;
 430   H_im[6] = 0.5 * sin_theta/(omc2*omc2*omc) * (8.0*omc - sin_theta*sin_theta*(3.0 + cos_theta));
 431
 432   sum_re = H_re[0];
 433   sum_im = H_im[0];
 434   an = 1.0;
 435   nfact = 1.0;
 436   for(n=1; n<=6; n++) {
 437     double t;
 438     an *= a;
 439     nfact *= n;
 440     t = an/nfact;
 441     sum_re += t * H_re[n];
 442     sum_im += t * H_im[n];
 443   }
 444
 445   real_result->val = sum_re;
 446   real_result->err = 2.0 * 6.0 * GSL_DBL_EPSILON * fabs(sum_re) + fabs(an/nfact);
 447   imag_result->val = sum_im;
 448   imag_result->err = 2.0 * 6.0 * GSL_DBL_EPSILON * fabs(sum_im) + Him0.err + fabs(an/nfact);
 449
 450   return GSL_SUCCESS;
 451 }
 452
 453
 454 /* Calculate complex dilogarithm Li_2(z) in the fundamental region,
 455  * which we take to be the intersection of the unit disk with the
 456  * half-space x < MAGIC_SPLIT_VALUE. It turns out that 0.732 is a
 457  * nice choice for MAGIC_SPLIT_VALUE since then points mapped out
 458  * of the x > MAGIC_SPLIT_VALUE region and into another part of the
 459  * unit disk are bounded in radius by MAGIC_SPLIT_VALUE itself.
 460  *
 461  * If |z| < 0.98 we use a direct series summation. Otherwise z is very
 462  * near the unit circle, and the series_2 expansion is used; see above.
 463  * Because the fundamental region is bounded away from z = 1, this
 464  * works well.
 465  */
 466 static
 467 int
 468 dilogc_fundamental(double r, double x, double y, gsl_sf_result * real_dl, gsl_sf_result * imag_dl)
 469 {
 470   if(r > 0.98)
 471     return dilogc_series_3(r, x, y, real_dl, imag_dl);
 472   else if(r > 0.25)
 473     return dilogc_series_2(r, x, y, real_dl, imag_dl);
 474   else
 475     return dilogc_series_1(r, x, y, real_dl, imag_dl);
 476 }
 477
 478
 479 /* Compute Li_2(z) for z in the unit disk, |z| < 1. If z is outside
 480  * the fundamental region, which means that it is too close to z = 1,
 481  * then it is reflected into the fundamental region using the identity
 482  *
 483  *   Li2(z) = -Li2(1-z) + zeta(2) - ln(z) ln(1-z).
 484  */
 485 static
 486 int
 487 dilogc_unitdisk(double x, double y, gsl_sf_result * real_dl, gsl_sf_result * imag_dl)
 488 {
 489   static const double MAGIC_SPLIT_VALUE = 0.732;
 490   static const double zeta2 = M_PI*M_PI/6.0;
 491   const double r = hypot(x, y);
 492
 493   if(x > MAGIC_SPLIT_VALUE)
 494   {
 495     /* Reflect away from z = 1 if we are too close. The magic value
 496      * insures that the reflected value of the radius satisfies the
 497      * related inequality r_tmp < MAGIC_SPLIT_VALUE.
 498      */
 499     const double x_tmp = 1.0 - x;
 500     const double y_tmp =     - y;
 501     const double r_tmp = hypot(x_tmp, y_tmp);
 502     /* const double cos_theta_tmp = x_tmp/r_tmp; */
 503     /* const double sin_theta_tmp = y_tmp/r_tmp; */
 504
 505     gsl_sf_result result_re_tmp;
 506     gsl_sf_result result_im_tmp;
 507
 508     const int stat_dilog = dilogc_fundamental(r_tmp, x_tmp, y_tmp, &result_re_tmp, &result_im_tmp);
 509
 510     const double lnz    =  log(r);               /*  log(|z|)   */
 511     const double lnomz  =  log(r_tmp);           /*  log(|1-z|) */
 512     const double argz   =  atan2(y, x);          /*  arg(z) assuming principal branch */
 513     const double argomz =  atan2(y_tmp, x_tmp);  /*  arg(1-z)   */
 514     real_dl->val  = -result_re_tmp.val + zeta2 - lnz*lnomz + argz*argomz;
 515     real_dl->err  =  result_re_tmp.err;
 516     real_dl->err +=  2.0 * GSL_DBL_EPSILON * (zeta2 + fabs(lnz*lnomz) + fabs(argz*argomz));
 517     imag_dl->val  = -result_im_tmp.val - argz*lnomz - argomz*lnz;
 518     imag_dl->err  =  result_im_tmp.err;
 519     imag_dl->err +=  2.0 * GSL_DBL_EPSILON * (fabs(argz*lnomz) + fabs(argomz*lnz));
 520
 521     return stat_dilog;
 522   }
 523   else
 524   {
 525     return dilogc_fundamental(r, x, y, real_dl, imag_dl);
 526   }
 527 }
 528
 529
 530
 531 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 532
 533
 534 int
 535 gsl_sf_dilog_e(const double x, gsl_sf_result * result)
 536 {
 537   if(x >= 0.0) {
 538     return dilog_xge0(x, result);
 539   }
 540   else {
 541     gsl_sf_result d1, d2;
 542     int stat_d1 = dilog_xge0( -x, &d1);
 543     int stat_d2 = dilog_xge0(x*x, &d2);
 544     result->val  = -d1.val + 0.5 * d2.val;
 545     result->err  =  d1.err + 0.5 * d2.err;
 546     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 547     return GSL_ERROR_SELECT_2(stat_d1, stat_d2);
 548   }
 549 }
 550
 551
 552 int
 553 gsl_sf_complex_dilog_xy_e(
 554   const double x,
 555   const double y,
 556   gsl_sf_result * real_dl,
 557   gsl_sf_result * imag_dl
 558   )
 559 {
 560   const double zeta2 = M_PI*M_PI/6.0;
 561   const double r2 = x*x + y*y;
 562
 563   if(y == 0.0)
 564   {
 565     if(x >= 1.0)
 566     {
 567       imag_dl->val = -M_PI * log(x);
 568       imag_dl->err = 2.0 * GSL_DBL_EPSILON * fabs(imag_dl->val);
 569     }
 570     else
 571     {
 572       imag_dl->val = 0.0;
 573       imag_dl->err = 0.0;
 574     }
 575     return gsl_sf_dilog_e(x, real_dl);
 576   }
 577   else if(fabs(r2 - 1.0) < GSL_DBL_EPSILON)
 578   {
 579     /* Lewin A.2.4.1 and A.2.4.2 */
 580
 581     const double theta = atan2(y, x);
 582     const double term1 = theta*theta/4.0;
 583     const double term2 = M_PI*fabs(theta)/2.0;
 584     real_dl->val = zeta2 + term1 - term2;
 585     real_dl->err = 2.0 * GSL_DBL_EPSILON * (zeta2 + term1 + term2);
 586     return gsl_sf_clausen_e(theta, imag_dl);
 587   }
 588   else if(r2 < 1.0)
 589   {
 590     return dilogc_unitdisk(x, y, real_dl, imag_dl);
 591   }
 592   else
 593   {
 594     /* Reduce argument to unit disk. */
 595     const double r = sqrt(r2);
 596     const double x_tmp =  x/r2;
 597     const double y_tmp = -y/r2;
 598     /* const double r_tmp = 1.0/r; */
 599     gsl_sf_result result_re_tmp, result_im_tmp;
 600
 601     const int stat_dilog =
 602       dilogc_unitdisk(x_tmp, y_tmp, &result_re_tmp, &result_im_tmp);
 603
 604     /* Unwind the inversion.
 605      *
 606      *  Li_2(z) + Li_2(1/z) = -zeta(2) - 1/2 ln(-z)^2
 607      */
 608     const double theta = atan2(y, x);
 609     const double theta_abs = fabs(theta);
 610     const double theta_sgn = ( theta < 0.0 ? -1.0 : 1.0 );
 611     const double ln_minusz_re = log(r);
 612     const double ln_minusz_im = theta_sgn * (theta_abs - M_PI);
 613     const double lmz2_re = ln_minusz_re*ln_minusz_re - ln_minusz_im*ln_minusz_im;
 614     const double lmz2_im = 2.0*ln_minusz_re*ln_minusz_im;
 615     real_dl->val = -result_re_tmp.val - 0.5 * lmz2_re - zeta2;
 616     real_dl->err =  result_re_tmp.err + 2.0*GSL_DBL_EPSILON*(0.5 * fabs(lmz2_re) + zeta2);
 617     imag_dl->val = -result_im_tmp.val - 0.5 * lmz2_im;
 618     imag_dl->err =  result_im_tmp.err + 2.0*GSL_DBL_EPSILON*fabs(lmz2_im);
 619     return stat_dilog;
 620   }
 621 }
 622
 623
 624 int
 625 gsl_sf_complex_dilog_e(
 626   const double r,
 627   const double theta,
 628   gsl_sf_result * real_dl,
 629   gsl_sf_result * imag_dl
 630   )
 631 {
 632   const double cos_theta = cos(theta);
 633   const double sin_theta = sin(theta);
 634   const double x = r * cos_theta;
 635   const double y = r * sin_theta;
 636   return gsl_sf_complex_dilog_xy_e(x, y, real_dl, imag_dl);
 637 }
 638
 639
 640 int
 641 gsl_sf_complex_spence_xy_e(
 642   const double x,
 643   const double y,
 644   gsl_sf_result * real_sp,
 645   gsl_sf_result * imag_sp
 646   )
 647 {
 648   const double oms_x = 1.0 - x;
 649   const double oms_y =     - y;
 650   return gsl_sf_complex_dilog_xy_e(oms_x, oms_y, real_sp, imag_sp);
 651 }
 652
 653
 654
 655 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 656
 657 #include "eval.h"
 658
 659 double gsl_sf_dilog(const double x)
 660 {
 661   EVAL_RESULT(gsl_sf_dilog_e(x, &result));
 662 }