spence.c
Go to the documentation of this file.
1 /* spence.c
2  *
3  * Dilogarithm
4  *
5  *
6  *
7  * SYNOPSIS:
8  *
9  * double x, y, spence();
10  *
11  * y = spence( x );
12  *
13  *
14  *
15  * DESCRIPTION:
16  *
17  * Computes the integral
18  *
19  * x
20  * -
21  * | | log t
22  * spence(x) = - | ----- dt
23  * | | t - 1
24  * -
25  * 1
26  *
27  * for x >= 0. A rational approximation gives the integral in
28  * the interval (0.5, 1.5). Transformation formulas for 1/x
29  * and 1-x are employed outside the basic expansion range.
30  *
31  *
32  *
33  * ACCURACY:
34  *
35  * Relative error:
36  * arithmetic domain # trials peak rms
37  * IEEE 0,4 30000 3.9e-15 5.4e-16
38  *
39  *
40  */
41 
42 /* spence.c */
43 
44 
45 /*
46  * Cephes Math Library Release 2.1: January, 1989
47  * Copyright 1985, 1987, 1989 by Stephen L. Moshier
48  * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
49  */
50 
51 #include "mconf.h"
52 
53 static double A[8] = {
54  4.65128586073990045278E-5,
55  7.31589045238094711071E-3,
56  1.33847639578309018650E-1,
57  8.79691311754530315341E-1,
58  2.71149851196553469920E0,
59  4.25697156008121755724E0,
60  3.29771340985225106936E0,
61  1.00000000000000000126E0,
62 };
63 
64 static double B[8] = {
65  6.90990488912553276999E-4,
66  2.54043763932544379113E-2,
67  2.82974860602568089943E-1,
68  1.41172597751831069617E0,
69  3.63800533345137075418E0,
70  5.03278880143316990390E0,
71  3.54771340985225096217E0,
72  9.99999999999999998740E-1,
73 };
74 
75 extern double MACHEP;
76 
77 double spence(double x)
78 {
79  double w, y, z;
80  int flag;
81 
82  if (x < 0.0) {
83  sf_error("spence", SF_ERROR_DOMAIN, NULL);
84  return (NAN);
85  }
86 
87  if (x == 1.0)
88  return (0.0);
89 
90  if (x == 0.0)
91  return (M_PI * M_PI / 6.0);
92 
93  flag = 0;
94 
95  if (x > 2.0) {
96  x = 1.0 / x;
97  flag |= 2;
98  }
99 
100  if (x > 1.5) {
101  w = (1.0 / x) - 1.0;
102  flag |= 2;
103  }
104 
105  else if (x < 0.5) {
106  w = -x;
107  flag |= 1;
108  }
109 
110  else
111  w = x - 1.0;
112 
113 
114  y = -w * polevl(w, A, 7) / polevl(w, B, 7);
115 
116  if (flag & 1)
117  y = (M_PI * M_PI) / 6.0 - log(x) * log(1.0 - x) - y;
118 
119  if (flag & 2) {
120  z = log(x);
121  y = -0.5 * z * z - y;
122  }
123 
124  return (y);
125 }
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
spence
double spence(double x)
Definition: spence.c:77
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
B
Definition: test_numpy_dtypes.cpp:299
log
const EIGEN_DEVICE_FUNC LogReturnType log() const
Definition: ArrayCwiseUnaryOps.h:128
MACHEP
double MACHEP
Definition: const.c:54
A
Definition: test_numpy_dtypes.cpp:298
polevl
static double polevl(double x, const double coef[], int N)
Definition: polevl.h:65
pybind_wrapper_test_script.z
z
Definition: pybind_wrapper_test_script.py:61
y
Scalar * y
Definition: level1_cplx_impl.h:124
mconf.h
sf_error
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Definition: sf_error.c:41
NULL
#define NULL
Definition: ccolamd.c:609
M_PI
#define M_PI
Definition: mconf.h:117
SF_ERROR_DOMAIN
@ SF_ERROR_DOMAIN
Definition: sf_error.h:16


gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:04:59