|
Packit |
67cb25 |
/* specfunc/dilog.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Author: G. Jungman */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_clausen.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_trig.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_log.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_dilog.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Evaluate series for real dilog(x)
|
|
Packit |
67cb25 |
* Sum[ x^k / k^2, {k,1,Infinity}]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Converges rapidly for |x| < 1/2.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
dilog_series_1(const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const int kmax = 1000;
|
|
Packit |
67cb25 |
double sum = x;
|
|
Packit |
67cb25 |
double term = x;
|
|
Packit |
67cb25 |
int k;
|
|
Packit |
67cb25 |
for(k=2; k
|
|
Packit |
67cb25 |
const double rk = (k-1.0)/k;
|
|
Packit |
67cb25 |
term *= x;
|
|
Packit |
67cb25 |
term *= rk*rk;
|
|
Packit |
67cb25 |
sum += term;
|
|
Packit |
67cb25 |
if(fabs(term/sum) < GSL_DBL_EPSILON) break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = sum;
|
|
Packit |
67cb25 |
result->err = 2.0 * fabs(term);
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(k == kmax)
|
|
Packit |
67cb25 |
GSL_ERROR ("error", GSL_EMAXITER);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the associated series
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* sum_{k=1}{infty} r^k / (k^2 (k+1))
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This is a series which appears in the one-step accelerated
|
|
Packit |
67cb25 |
* method, which splits out one elementary function from the
|
|
Packit |
67cb25 |
* full definition of Li_2(x). See below.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
series_2(double r, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
static const int kmax = 100;
|
|
Packit |
67cb25 |
double rk = r;
|
|
Packit |
67cb25 |
double sum = 0.5 * r;
|
|
Packit |
67cb25 |
int k;
|
|
Packit |
67cb25 |
for(k=2; k<10; k++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double ds;
|
|
Packit |
67cb25 |
rk *= r;
|
|
Packit |
67cb25 |
ds = rk/(k*k*(k+1.0));
|
|
Packit |
67cb25 |
sum += ds;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
for(; k
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double ds;
|
|
Packit |
67cb25 |
rk *= r;
|
|
Packit |
67cb25 |
ds = rk/(k*k*(k+1.0));
|
|
Packit |
67cb25 |
sum += ds;
|
|
Packit |
67cb25 |
if(fabs(ds/sum) < 0.5*GSL_DBL_EPSILON) break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = sum;
|
|
Packit |
67cb25 |
result->err = 2.0 * kmax * GSL_DBL_EPSILON * fabs(sum);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute Li_2(x) using the accelerated series representation.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Li_2(x) = 1 + (1-x)ln(1-x)/x + series_2(x)
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* assumes: -1 < x < 1
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
dilog_series_2(double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const int stat_s3 = series_2(x, result);
|
|
Packit |
67cb25 |
double t;
|
|
Packit |
67cb25 |
if(x > 0.01)
|
|
Packit |
67cb25 |
t = (1.0 - x) * log(1.0-x) / x;
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
static const double c3 = 1.0/3.0;
|
|
Packit |
67cb25 |
static const double c4 = 1.0/4.0;
|
|
Packit |
67cb25 |
static const double c5 = 1.0/5.0;
|
|
Packit |
67cb25 |
static const double c6 = 1.0/6.0;
|
|
Packit |
67cb25 |
static const double c7 = 1.0/7.0;
|
|
Packit |
67cb25 |
static const double c8 = 1.0/8.0;
|
|
Packit |
67cb25 |
const double t68 = c6 + x*(c7 + x*c8);
|
|
Packit |
67cb25 |
const double t38 = c3 + x *(c4 + x *(c5 + x * t68));
|
|
Packit |
67cb25 |
t = (x - 1.0) * (1.0 + x*(0.5 + x*t38));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
result->val += 1.0 + t;
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(t);
|
|
Packit |
67cb25 |
return stat_s3;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Calculates Li_2(x) for real x. Assumes x >= 0.0.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
dilog_xge0(const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(x > 2.0) {
|
|
Packit |
67cb25 |
gsl_sf_result ser;
|
|
Packit |
67cb25 |
const int stat_ser = dilog_series_2(1.0/x, &ser;;
|
|
Packit |
67cb25 |
const double log_x = log(x);
|
|
Packit |
67cb25 |
const double t1 = M_PI*M_PI/3.0;
|
|
Packit |
67cb25 |
const double t2 = ser.val;
|
|
Packit |
67cb25 |
const double t3 = 0.5*log_x*log_x;
|
|
Packit |
67cb25 |
result->val = t1 - t2 - t3;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * fabs(log_x) + ser.err;
|
|
Packit |
67cb25 |
result->err += GSL_DBL_EPSILON * (fabs(t1) + fabs(t2) + fabs(t3));
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return stat_ser;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x > 1.01) {
|
|
Packit |
67cb25 |
gsl_sf_result ser;
|
|
Packit |
67cb25 |
const int stat_ser = dilog_series_2(1.0 - 1.0/x, &ser;;
|
|
Packit |
67cb25 |
const double log_x = log(x);
|
|
Packit |
67cb25 |
const double log_term = log_x * (log(1.0-1.0/x) + 0.5*log_x);
|
|
Packit |
67cb25 |
const double t1 = M_PI*M_PI/6.0;
|
|
Packit |
67cb25 |
const double t2 = ser.val;
|
|
Packit |
67cb25 |
const double t3 = log_term;
|
|
Packit |
67cb25 |
result->val = t1 + t2 - t3;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * fabs(log_x) + ser.err;
|
|
Packit |
67cb25 |
result->err += GSL_DBL_EPSILON * (fabs(t1) + fabs(t2) + fabs(t3));
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return stat_ser;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x > 1.0) {
|
|
Packit |
67cb25 |
/* series around x = 1.0 */
|
|
Packit |
67cb25 |
const double eps = x - 1.0;
|
|
Packit |
67cb25 |
const double lne = log(eps);
|
|
Packit |
67cb25 |
const double c0 = M_PI*M_PI/6.0;
|
|
Packit |
67cb25 |
const double c1 = 1.0 - lne;
|
|
Packit |
67cb25 |
const double c2 = -(1.0 - 2.0*lne)/4.0;
|
|
Packit |
67cb25 |
const double c3 = (1.0 - 3.0*lne)/9.0;
|
|
Packit |
67cb25 |
const double c4 = -(1.0 - 4.0*lne)/16.0;
|
|
Packit |
67cb25 |
const double c5 = (1.0 - 5.0*lne)/25.0;
|
|
Packit |
67cb25 |
const double c6 = -(1.0 - 6.0*lne)/36.0;
|
|
Packit |
67cb25 |
const double c7 = (1.0 - 7.0*lne)/49.0;
|
|
Packit |
67cb25 |
const double c8 = -(1.0 - 8.0*lne)/64.0;
|
|
Packit |
67cb25 |
result->val = c0+eps*(c1+eps*(c2+eps*(c3+eps*(c4+eps*(c5+eps*(c6+eps*(c7+eps*c8)))))));
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x == 1.0) {
|
|
Packit |
67cb25 |
result->val = M_PI*M_PI/6.0;
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * M_PI*M_PI/6.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x > 0.5) {
|
|
Packit |
67cb25 |
gsl_sf_result ser;
|
|
Packit |
67cb25 |
const int stat_ser = dilog_series_2(1.0-x, &ser;;
|
|
Packit |
67cb25 |
const double log_x = log(x);
|
|
Packit |
67cb25 |
const double t1 = M_PI*M_PI/6.0;
|
|
Packit |
67cb25 |
const double t2 = ser.val;
|
|
Packit |
67cb25 |
const double t3 = log_x*log(1.0-x);
|
|
Packit |
67cb25 |
result->val = t1 - t2 - t3;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * fabs(log_x) + ser.err;
|
|
Packit |
67cb25 |
result->err += GSL_DBL_EPSILON * (fabs(t1) + fabs(t2) + fabs(t3));
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return stat_ser;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x > 0.25) {
|
|
Packit |
67cb25 |
return dilog_series_2(x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x > 0.0) {
|
|
Packit |
67cb25 |
return dilog_series_1(x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* x == 0.0 */
|
|
Packit |
67cb25 |
result->val = 0.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Evaluate the series representation for Li2(z):
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Li2(z) = Sum[ |z|^k / k^2 Exp[i k arg(z)], {k,1,Infinity}]
|
|
Packit |
67cb25 |
* |z| = r
|
|
Packit |
67cb25 |
* arg(z) = theta
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Assumes 0 < r < 1.
|
|
Packit |
67cb25 |
* It is used only for small r.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
dilogc_series_1(
|
|
Packit |
67cb25 |
const double r,
|
|
Packit |
67cb25 |
const double x,
|
|
Packit |
67cb25 |
const double y,
|
|
Packit |
67cb25 |
gsl_sf_result * real_result,
|
|
Packit |
67cb25 |
gsl_sf_result * imag_result
|
|
Packit |
67cb25 |
)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double cos_theta = x/r;
|
|
Packit |
67cb25 |
const double sin_theta = y/r;
|
|
Packit |
67cb25 |
const double alpha = 1.0 - cos_theta;
|
|
Packit |
67cb25 |
const double beta = sin_theta;
|
|
Packit |
67cb25 |
double ck = cos_theta;
|
|
Packit |
67cb25 |
double sk = sin_theta;
|
|
Packit |
67cb25 |
double rk = r;
|
|
Packit |
67cb25 |
double real_sum = r*ck;
|
|
Packit |
67cb25 |
double imag_sum = r*sk;
|
|
Packit |
67cb25 |
const int kmax = 50 + (int)(22.0/(-log(r))); /* tuned for double-precision */
|
|
Packit |
67cb25 |
int k;
|
|
Packit |
67cb25 |
for(k=2; k
|
|
Packit |
67cb25 |
double dr, di;
|
|
Packit |
67cb25 |
double ck_tmp = ck;
|
|
Packit |
67cb25 |
ck = ck - (alpha*ck + beta*sk);
|
|
Packit |
67cb25 |
sk = sk - (alpha*sk - beta*ck_tmp);
|
|
Packit |
67cb25 |
rk *= r;
|
|
Packit |
67cb25 |
dr = rk/((double)k*k) * ck;
|
|
Packit |
67cb25 |
di = rk/((double)k*k) * sk;
|
|
Packit |
67cb25 |
real_sum += dr;
|
|
Packit |
67cb25 |
imag_sum += di;
|
|
Packit |
67cb25 |
if(fabs((dr*dr + di*di)/(real_sum*real_sum + imag_sum*imag_sum)) < GSL_DBL_EPSILON*GSL_DBL_EPSILON) break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
real_result->val = real_sum;
|
|
Packit |
67cb25 |
real_result->err = 2.0 * kmax * GSL_DBL_EPSILON * fabs(real_sum);
|
|
Packit |
67cb25 |
imag_result->val = imag_sum;
|
|
Packit |
67cb25 |
imag_result->err = 2.0 * kmax * GSL_DBL_EPSILON * fabs(imag_sum);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* sum_{k=1}{infty} z^k / (k^2 (k+1))
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This is a series which appears in the one-step accelerated
|
|
Packit |
67cb25 |
* method, which splits out one elementary function from the
|
|
Packit |
67cb25 |
* full definition of Li_2.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
series_2_c(
|
|
Packit |
67cb25 |
double r,
|
|
Packit |
67cb25 |
double x,
|
|
Packit |
67cb25 |
double y,
|
|
Packit |
67cb25 |
gsl_sf_result * sum_re,
|
|
Packit |
67cb25 |
gsl_sf_result * sum_im
|
|
Packit |
67cb25 |
)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double cos_theta = x/r;
|
|
Packit |
67cb25 |
const double sin_theta = y/r;
|
|
Packit |
67cb25 |
const double alpha = 1.0 - cos_theta;
|
|
Packit |
67cb25 |
const double beta = sin_theta;
|
|
Packit |
67cb25 |
double ck = cos_theta;
|
|
Packit |
67cb25 |
double sk = sin_theta;
|
|
Packit |
67cb25 |
double rk = r;
|
|
Packit |
67cb25 |
double real_sum = 0.5 * r*ck;
|
|
Packit |
67cb25 |
double imag_sum = 0.5 * r*sk;
|
|
Packit |
67cb25 |
const int kmax = 30 + (int)(18.0/(-log(r))); /* tuned for double-precision */
|
|
Packit |
67cb25 |
int k;
|
|
Packit |
67cb25 |
for(k=2; k
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double dr, di;
|
|
Packit |
67cb25 |
const double ck_tmp = ck;
|
|
Packit |
67cb25 |
ck = ck - (alpha*ck + beta*sk);
|
|
Packit |
67cb25 |
sk = sk - (alpha*sk - beta*ck_tmp);
|
|
Packit |
67cb25 |
rk *= r;
|
|
Packit |
67cb25 |
dr = rk/((double)k*k*(k+1.0)) * ck;
|
|
Packit |
67cb25 |
di = rk/((double)k*k*(k+1.0)) * sk;
|
|
Packit |
67cb25 |
real_sum += dr;
|
|
Packit |
67cb25 |
imag_sum += di;
|
|
Packit |
67cb25 |
if(fabs((dr*dr + di*di)/(real_sum*real_sum + imag_sum*imag_sum)) < GSL_DBL_EPSILON*GSL_DBL_EPSILON) break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
sum_re->val = real_sum;
|
|
Packit |
67cb25 |
sum_re->err = 2.0 * kmax * GSL_DBL_EPSILON * fabs(real_sum);
|
|
Packit |
67cb25 |
sum_im->val = imag_sum;
|
|
Packit |
67cb25 |
sum_im->err = 2.0 * kmax * GSL_DBL_EPSILON * fabs(imag_sum);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute Li_2(z) using the one-step accelerated series.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Li_2(z) = 1 + (1-z)ln(1-z)/z + series_2_c(z)
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* z = r exp(i theta)
|
|
Packit |
67cb25 |
* assumes: r < 1
|
|
Packit |
67cb25 |
* assumes: r > epsilon, so that we take no special care with log(1-z)
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
dilogc_series_2(
|
|
Packit |
67cb25 |
const double r,
|
|
Packit |
67cb25 |
const double x,
|
|
Packit |
67cb25 |
const double y,
|
|
Packit |
67cb25 |
gsl_sf_result * real_dl,
|
|
Packit |
67cb25 |
gsl_sf_result * imag_dl
|
|
Packit |
67cb25 |
)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(r == 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
real_dl->val = 0.0;
|
|
Packit |
67cb25 |
imag_dl->val = 0.0;
|
|
Packit |
67cb25 |
real_dl->err = 0.0;
|
|
Packit |
67cb25 |
imag_dl->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_sf_result sum_re;
|
|
Packit |
67cb25 |
gsl_sf_result sum_im;
|
|
Packit |
67cb25 |
const int stat_s3 = series_2_c(r, x, y, &sum_re, &sum_im);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* t = ln(1-z)/z */
|
|
Packit |
67cb25 |
gsl_sf_result ln_omz_r;
|
|
Packit |
67cb25 |
gsl_sf_result ln_omz_theta;
|
|
Packit |
67cb25 |
const int stat_log = gsl_sf_complex_log_e(1.0-x, -y, &ln_omz_r, &ln_omz_theta);
|
|
Packit |
67cb25 |
const double t_x = ( ln_omz_r.val * x + ln_omz_theta.val * y)/(r*r);
|
|
Packit |
67cb25 |
const double t_y = (-ln_omz_r.val * y + ln_omz_theta.val * x)/(r*r);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* r = (1-z) ln(1-z)/z */
|
|
Packit |
67cb25 |
const double r_x = (1.0 - x) * t_x + y * t_y;
|
|
Packit |
67cb25 |
const double r_y = (1.0 - x) * t_y - y * t_x;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
real_dl->val = sum_re.val + r_x + 1.0;
|
|
Packit |
67cb25 |
imag_dl->val = sum_im.val + r_y;
|
|
Packit |
67cb25 |
real_dl->err = sum_re.err + 2.0*GSL_DBL_EPSILON*(fabs(real_dl->val) + fabs(r_x));
|
|
Packit |
67cb25 |
imag_dl->err = sum_im.err + 2.0*GSL_DBL_EPSILON*(fabs(imag_dl->val) + fabs(r_y));
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_s3, stat_log);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Evaluate a series for Li_2(z) when |z| is near 1.
|
|
Packit |
67cb25 |
* This is uniformly good away from z=1.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Li_2(z) = Sum[ a^n/n! H_n(theta), {n, 0, Infinity}]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* where
|
|
Packit |
67cb25 |
* H_n(theta) = Sum[ e^(i m theta) m^n / m^2, {m, 1, Infinity}]
|
|
Packit |
67cb25 |
* a = ln(r)
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* H_0(t) = Gl_2(t) + i Cl_2(t)
|
|
Packit |
67cb25 |
* H_1(t) = 1/2 ln(2(1-c)) + I atan2(-s, 1-c)
|
|
Packit |
67cb25 |
* H_2(t) = -1/2 + I/2 s/(1-c)
|
|
Packit |
67cb25 |
* H_3(t) = -1/2 /(1-c)
|
|
Packit |
67cb25 |
* H_4(t) = -I/2 s/(1-c)^2
|
|
Packit |
67cb25 |
* H_5(t) = 1/2 (2 + c)/(1-c)^2
|
|
Packit |
67cb25 |
* H_6(t) = I/2 s/(1-c)^5 (8(1-c) - s^2 (3 + c))
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
dilogc_series_3(
|
|
Packit |
67cb25 |
const double r,
|
|
Packit |
67cb25 |
const double x,
|
|
Packit |
67cb25 |
const double y,
|
|
Packit |
67cb25 |
gsl_sf_result * real_result,
|
|
Packit |
67cb25 |
gsl_sf_result * imag_result
|
|
Packit |
67cb25 |
)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double theta = atan2(y, x);
|
|
Packit |
67cb25 |
const double cos_theta = x/r;
|
|
Packit |
67cb25 |
const double sin_theta = y/r;
|
|
Packit |
67cb25 |
const double a = log(r);
|
|
Packit |
67cb25 |
const double omc = 1.0 - cos_theta;
|
|
Packit |
67cb25 |
const double omc2 = omc*omc;
|
|
Packit |
67cb25 |
double H_re[7];
|
|
Packit |
67cb25 |
double H_im[7];
|
|
Packit |
67cb25 |
double an, nfact;
|
|
Packit |
67cb25 |
double sum_re, sum_im;
|
|
Packit |
67cb25 |
gsl_sf_result Him0;
|
|
Packit |
67cb25 |
int n;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
H_re[0] = M_PI*M_PI/6.0 + 0.25*(theta*theta - 2.0*M_PI*fabs(theta));
|
|
Packit |
67cb25 |
gsl_sf_clausen_e(theta, &Him0);
|
|
Packit |
67cb25 |
H_im[0] = Him0.val;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
H_re[1] = -0.5*log(2.0*omc);
|
|
Packit |
67cb25 |
H_im[1] = -atan2(-sin_theta, omc);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
H_re[2] = -0.5;
|
|
Packit |
67cb25 |
H_im[2] = 0.5 * sin_theta/omc;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
H_re[3] = -0.5/omc;
|
|
Packit |
67cb25 |
H_im[3] = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
H_re[4] = 0.0;
|
|
Packit |
67cb25 |
H_im[4] = -0.5*sin_theta/omc2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
H_re[5] = 0.5 * (2.0 + cos_theta)/omc2;
|
|
Packit |
67cb25 |
H_im[5] = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
H_re[6] = 0.0;
|
|
Packit |
67cb25 |
H_im[6] = 0.5 * sin_theta/(omc2*omc2*omc) * (8.0*omc - sin_theta*sin_theta*(3.0 + cos_theta));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
sum_re = H_re[0];
|
|
Packit |
67cb25 |
sum_im = H_im[0];
|
|
Packit |
67cb25 |
an = 1.0;
|
|
Packit |
67cb25 |
nfact = 1.0;
|
|
Packit |
67cb25 |
for(n=1; n<=6; n++) {
|
|
Packit |
67cb25 |
double t;
|
|
Packit |
67cb25 |
an *= a;
|
|
Packit |
67cb25 |
nfact *= n;
|
|
Packit |
67cb25 |
t = an/nfact;
|
|
Packit |
67cb25 |
sum_re += t * H_re[n];
|
|
Packit |
67cb25 |
sum_im += t * H_im[n];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
real_result->val = sum_re;
|
|
Packit |
67cb25 |
real_result->err = 2.0 * 6.0 * GSL_DBL_EPSILON * fabs(sum_re) + fabs(an/nfact);
|
|
Packit |
67cb25 |
imag_result->val = sum_im;
|
|
Packit |
67cb25 |
imag_result->err = 2.0 * 6.0 * GSL_DBL_EPSILON * fabs(sum_im) + Him0.err + fabs(an/nfact);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Calculate complex dilogarithm Li_2(z) in the fundamental region,
|
|
Packit |
67cb25 |
* which we take to be the intersection of the unit disk with the
|
|
Packit |
67cb25 |
* half-space x < MAGIC_SPLIT_VALUE. It turns out that 0.732 is a
|
|
Packit |
67cb25 |
* nice choice for MAGIC_SPLIT_VALUE since then points mapped out
|
|
Packit |
67cb25 |
* of the x > MAGIC_SPLIT_VALUE region and into another part of the
|
|
Packit |
67cb25 |
* unit disk are bounded in radius by MAGIC_SPLIT_VALUE itself.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* If |z| < 0.98 we use a direct series summation. Otherwise z is very
|
|
Packit |
67cb25 |
* near the unit circle, and the series_2 expansion is used; see above.
|
|
Packit |
67cb25 |
* Because the fundamental region is bounded away from z = 1, this
|
|
Packit |
67cb25 |
* works well.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
dilogc_fundamental(double r, double x, double y, gsl_sf_result * real_dl, gsl_sf_result * imag_dl)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(r > 0.98)
|
|
Packit |
67cb25 |
return dilogc_series_3(r, x, y, real_dl, imag_dl);
|
|
Packit |
67cb25 |
else if(r > 0.25)
|
|
Packit |
67cb25 |
return dilogc_series_2(r, x, y, real_dl, imag_dl);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return dilogc_series_1(r, x, y, real_dl, imag_dl);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute Li_2(z) for z in the unit disk, |z| < 1. If z is outside
|
|
Packit |
67cb25 |
* the fundamental region, which means that it is too close to z = 1,
|
|
Packit |
67cb25 |
* then it is reflected into the fundamental region using the identity
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Li2(z) = -Li2(1-z) + zeta(2) - ln(z) ln(1-z).
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
dilogc_unitdisk(double x, double y, gsl_sf_result * real_dl, gsl_sf_result * imag_dl)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
static const double MAGIC_SPLIT_VALUE = 0.732;
|
|
Packit |
67cb25 |
static const double zeta2 = M_PI*M_PI/6.0;
|
|
Packit |
67cb25 |
const double r = hypot(x, y);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x > MAGIC_SPLIT_VALUE)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Reflect away from z = 1 if we are too close. The magic value
|
|
Packit |
67cb25 |
* insures that the reflected value of the radius satisfies the
|
|
Packit |
67cb25 |
* related inequality r_tmp < MAGIC_SPLIT_VALUE.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
const double x_tmp = 1.0 - x;
|
|
Packit |
67cb25 |
const double y_tmp = - y;
|
|
Packit |
67cb25 |
const double r_tmp = hypot(x_tmp, y_tmp);
|
|
Packit |
67cb25 |
/* const double cos_theta_tmp = x_tmp/r_tmp; */
|
|
Packit |
67cb25 |
/* const double sin_theta_tmp = y_tmp/r_tmp; */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_sf_result result_re_tmp;
|
|
Packit |
67cb25 |
gsl_sf_result result_im_tmp;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const int stat_dilog = dilogc_fundamental(r_tmp, x_tmp, y_tmp, &result_re_tmp, &result_im_tmp);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const double lnz = log(r); /* log(|z|) */
|
|
Packit |
67cb25 |
const double lnomz = log(r_tmp); /* log(|1-z|) */
|
|
Packit |
67cb25 |
const double argz = atan2(y, x); /* arg(z) assuming principal branch */
|
|
Packit |
67cb25 |
const double argomz = atan2(y_tmp, x_tmp); /* arg(1-z) */
|
|
Packit |
67cb25 |
real_dl->val = -result_re_tmp.val + zeta2 - lnz*lnomz + argz*argomz;
|
|
Packit |
67cb25 |
real_dl->err = result_re_tmp.err;
|
|
Packit |
67cb25 |
real_dl->err += 2.0 * GSL_DBL_EPSILON * (zeta2 + fabs(lnz*lnomz) + fabs(argz*argomz));
|
|
Packit |
67cb25 |
imag_dl->val = -result_im_tmp.val - argz*lnomz - argomz*lnz;
|
|
Packit |
67cb25 |
imag_dl->err = result_im_tmp.err;
|
|
Packit |
67cb25 |
imag_dl->err += 2.0 * GSL_DBL_EPSILON * (fabs(argz*lnomz) + fabs(argomz*lnz));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return stat_dilog;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return dilogc_fundamental(r, x, y, real_dl, imag_dl);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_dilog_e(const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(x >= 0.0) {
|
|
Packit |
67cb25 |
return dilog_xge0(x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
gsl_sf_result d1, d2;
|
|
Packit |
67cb25 |
int stat_d1 = dilog_xge0( -x, &d1;;
|
|
Packit |
67cb25 |
int stat_d2 = dilog_xge0(x*x, &d2;;
|
|
Packit |
67cb25 |
result->val = -d1.val + 0.5 * d2.val;
|
|
Packit |
67cb25 |
result->err = d1.err + 0.5 * d2.err;
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_d1, stat_d2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_complex_dilog_xy_e(
|
|
Packit |
67cb25 |
const double x,
|
|
Packit |
67cb25 |
const double y,
|
|
Packit |
67cb25 |
gsl_sf_result * real_dl,
|
|
Packit |
67cb25 |
gsl_sf_result * imag_dl
|
|
Packit |
67cb25 |
)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double zeta2 = M_PI*M_PI/6.0;
|
|
Packit |
67cb25 |
const double r2 = x*x + y*y;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(y == 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(x >= 1.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
imag_dl->val = -M_PI * log(x);
|
|
Packit |
67cb25 |
imag_dl->err = 2.0 * GSL_DBL_EPSILON * fabs(imag_dl->val);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
imag_dl->val = 0.0;
|
|
Packit |
67cb25 |
imag_dl->err = 0.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
return gsl_sf_dilog_e(x, real_dl);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(fabs(r2 - 1.0) < GSL_DBL_EPSILON)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Lewin A.2.4.1 and A.2.4.2 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const double theta = atan2(y, x);
|
|
Packit |
67cb25 |
const double term1 = theta*theta/4.0;
|
|
Packit |
67cb25 |
const double term2 = M_PI*fabs(theta)/2.0;
|
|
Packit |
67cb25 |
real_dl->val = zeta2 + term1 - term2;
|
|
Packit |
67cb25 |
real_dl->err = 2.0 * GSL_DBL_EPSILON * (zeta2 + term1 + term2);
|
|
Packit |
67cb25 |
return gsl_sf_clausen_e(theta, imag_dl);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(r2 < 1.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return dilogc_unitdisk(x, y, real_dl, imag_dl);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Reduce argument to unit disk. */
|
|
Packit |
67cb25 |
const double r = sqrt(r2);
|
|
Packit |
67cb25 |
const double x_tmp = x/r2;
|
|
Packit |
67cb25 |
const double y_tmp = -y/r2;
|
|
Packit |
67cb25 |
/* const double r_tmp = 1.0/r; */
|
|
Packit |
67cb25 |
gsl_sf_result result_re_tmp, result_im_tmp;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const int stat_dilog =
|
|
Packit |
67cb25 |
dilogc_unitdisk(x_tmp, y_tmp, &result_re_tmp, &result_im_tmp);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Unwind the inversion.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Li_2(z) + Li_2(1/z) = -zeta(2) - 1/2 ln(-z)^2
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
const double theta = atan2(y, x);
|
|
Packit |
67cb25 |
const double theta_abs = fabs(theta);
|
|
Packit |
67cb25 |
const double theta_sgn = ( theta < 0.0 ? -1.0 : 1.0 );
|
|
Packit |
67cb25 |
const double ln_minusz_re = log(r);
|
|
Packit |
67cb25 |
const double ln_minusz_im = theta_sgn * (theta_abs - M_PI);
|
|
Packit |
67cb25 |
const double lmz2_re = ln_minusz_re*ln_minusz_re - ln_minusz_im*ln_minusz_im;
|
|
Packit |
67cb25 |
const double lmz2_im = 2.0*ln_minusz_re*ln_minusz_im;
|
|
Packit |
67cb25 |
real_dl->val = -result_re_tmp.val - 0.5 * lmz2_re - zeta2;
|
|
Packit |
67cb25 |
real_dl->err = result_re_tmp.err + 2.0*GSL_DBL_EPSILON*(0.5 * fabs(lmz2_re) + zeta2);
|
|
Packit |
67cb25 |
imag_dl->val = -result_im_tmp.val - 0.5 * lmz2_im;
|
|
Packit |
67cb25 |
imag_dl->err = result_im_tmp.err + 2.0*GSL_DBL_EPSILON*fabs(lmz2_im);
|
|
Packit |
67cb25 |
return stat_dilog;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_complex_dilog_e(
|
|
Packit |
67cb25 |
const double r,
|
|
Packit |
67cb25 |
const double theta,
|
|
Packit |
67cb25 |
gsl_sf_result * real_dl,
|
|
Packit |
67cb25 |
gsl_sf_result * imag_dl
|
|
Packit |
67cb25 |
)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double cos_theta = cos(theta);
|
|
Packit |
67cb25 |
const double sin_theta = sin(theta);
|
|
Packit |
67cb25 |
const double x = r * cos_theta;
|
|
Packit |
67cb25 |
const double y = r * sin_theta;
|
|
Packit |
67cb25 |
return gsl_sf_complex_dilog_xy_e(x, y, real_dl, imag_dl);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_complex_spence_xy_e(
|
|
Packit |
67cb25 |
const double x,
|
|
Packit |
67cb25 |
const double y,
|
|
Packit |
67cb25 |
gsl_sf_result * real_sp,
|
|
Packit |
67cb25 |
gsl_sf_result * imag_sp
|
|
Packit |
67cb25 |
)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double oms_x = 1.0 - x;
|
|
Packit |
67cb25 |
const double oms_y = - y;
|
|
Packit |
67cb25 |
return gsl_sf_complex_dilog_xy_e(oms_x, oms_y, real_sp, imag_sp);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "eval.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_dilog(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_dilog_e(x, &result));
|
|
Packit |
67cb25 |
}
|