Skip to content

Commit f0578cf

Browse files
ExpanderExpander
authored
Improved implementation for Clausen functions (#19)
1 parent bad77ed commit f0578cf

18 files changed

+3518
-149
lines changed

‎src/cpp/Cl2.cpp

+158-24
Original file line numberDiff line numberDiff line change
@@ -6,62 +6,196 @@
66

77
#include "Cl2.hpp"
88
#include "Li2.hpp"
9-
#include <limits>
109

1110
namespace polylogarithm {
1211

1312
/**
1413
* @brief Clausen function \f$\mathrm{Cl}_2(\theta) = \mathrm{Im}(\mathrm{Li}_2(e^{i\theta}))\f$
1514
* @param x real angle
1615
* @return \f$\mathrm{Cl}_2(\theta)\f$
16+
* @author K.S. Kölbig
17+
* @note Implementation translated by Alexander Voigt from CERNLIB DCLAUS function C326
18+
*
19+
* Journal of Computational and Applied Mathematics 64 (1995) 295-297.
1720
*/
1821
double Cl2(double x) noexcept
1922
{
20-
const double PI = 3.141592653589793;
21-
const std::complex<double> i(0.0, 1.0);
23+
const double PI = 3.14159265358979324;
24+
const double PI2 = 2*PI, PIH = PI/2, RPIH = 2/PI;
25+
const double A[9] = {
26+
0.0279528319735756613,
27+
0.0001763088743898116,
28+
0.0000012662741461157,
29+
0.0000000117171818134,
30+
0.0000000001230064129,
31+
0.0000000000013952729,
32+
0.0000000000000166908,
33+
0.0000000000000002076,
34+
0.0000000000000000027
35+
};
36+
const double B[14] = {
37+
0.639097088857265341,
38+
-0.054980569301851716,
39+
-0.000961261945950606,
40+
-0.000032054686822550,
41+
-0.000001329461695426,
42+
-0.000000062093601824,
43+
-0.000000003129600656,
44+
-0.000000000166351954,
45+
-0.000000000009196527,
46+
-0.000000000000524004,
47+
-0.000000000000030580,
48+
-0.000000000000001820,
49+
-0.000000000000000110,
50+
-0.000000000000000007,
51+
};
2252

23-
while (x >= 2*PI) {
24-
x -= 2*PI;
25-
}
53+
double h = 0;
54+
double v = std::fmod(std::abs(x), PI2);
55+
double sgn = x >= 0 ? 1 : -1;
2656

27-
while (x < 0.0) {
28-
x += 2*PI;
57+
if (v > PI) {
58+
const double p0 = 6.28125;
59+
const double p1 = 0.0019353071795864769253;
60+
v = (p0 - v) + p1;
61+
sgn = -sgn;
2962
}
3063

31-
if (std::abs(x) < std::numeric_limits<double>::epsilon() ||
32-
std::abs(x - PI) < std::numeric_limits<double>::epsilon() ||
33-
std::abs(x - 2*PI) < std::numeric_limits<double>::epsilon()) {
34-
return 0.0;
64+
if (v == 0 || v == PI) {
65+
h = 0;
66+
} else if (v < PIH) {
67+
const double u = RPIH*v;
68+
h = 2*u*u - 1;
69+
const double alfa = h + h;
70+
double b0 = 0, b1 = 0, b2 = 0;
71+
for (int i = 8; i >= 0; i--) {
72+
b0 = A[i] + alfa*b1 - b2;
73+
b2 = b1;
74+
b1 = b0;
75+
}
76+
h = v*(1 - std::log(v) + v*v*(b0 - h*b2)/2);
77+
} else {
78+
const double u = RPIH*v - 2;
79+
h = 2*u*u - 1;
80+
const double alfa = h + h;
81+
double b0 = 0, b1 = 0, b2 = 0;
82+
for (int i = 13; i >= 0; i--) {
83+
b0 = B[i] + alfa*b1 - b2;
84+
b2 = b1;
85+
b1 = b0;
86+
}
87+
h = (PI - v)*(b0 - h*b2);
3588
}
3689

37-
return std::imag(Li2(std::exp(i*x)));
90+
return sgn*h;
3891
}
3992

4093
/**
4194
* @brief Clausen function \f$\mathrm{Cl}_2(\theta) = \mathrm{Im}(\mathrm{Li}_2(e^{i\theta}))\f$ with long double precision
4295
* @param x real angle
4396
* @return \f$\mathrm{Cl}_2(\theta)\f$
97+
* @author K.S. Kölbig
98+
* @note Implementation translated from CERNLIB DCLAUS function C326
99+
* and extended to long double precision by Alexander Voigt.
100+
*
101+
* Journal of Computational and Applied Mathematics 64 (1995) 295-297.
44102
*/
45103
long double Cl2(long double x) noexcept
46104
{
47105
const long double PI = 3.14159265358979323846264338327950288L;
48-
const std::complex<long double> i(0.0L, 1.0L);
106+
const long double PI2 = 2*PI, PIH = PI/2, RPIH = 2/PI;
107+
const long double A[19] = {
108+
0.0279528319735756613494585924765551791L,
109+
0.0001763088743898115653057636473920103L,
110+
0.0000012662741461156530021975187159184L,
111+
0.0000000117171818134392379295428166866L,
112+
0.0000000001230064128833746922855709386L,
113+
0.0000000000013952728970012911958374309L,
114+
0.0000000000000166907761628567345146740L,
115+
0.0000000000000002076091315145432983502L,
116+
0.0000000000000000026609198306058056092L,
117+
0.0000000000000000000349249563561378275L,
118+
0.0000000000000000000004673313082962865L,
119+
0.0000000000000000000000063542322337428L,
120+
0.0000000000000000000000000875698871820L,
121+
0.0000000000000000000000000012208003299L,
122+
0.0000000000000000000000000000171890569L,
123+
0.0000000000000000000000000000002441331L,
124+
0.0000000000000000000000000000000034940L,
125+
0.0000000000000000000000000000000000503L,
126+
0.0000000000000000000000000000000000007L
127+
};
128+
const double B[30] = {
129+
0.6390970888572653413071869135953864197L,
130+
-0.0549805693018517156397035696498958507L,
131+
-0.0009612619459506064293859076874070709L,
132+
-0.0000320546868225504765586825318112711L,
133+
-0.0000013294616954255450141343828695514L,
134+
-0.0000000620936018243975194590942773212L,
135+
-0.0000000031296006563911126723262365339L,
136+
-0.0000000001663519538192669775933926077L,
137+
-0.0000000000091965272507194254496027281L,
138+
-0.0000000000005240037738758450093649037L,
139+
-0.0000000000000305803841873659454134183L,
140+
-0.0000000000000018196918249487950988000L,
141+
-0.0000000000000001100398263196261522324L,
142+
-0.0000000000000000067451775715424687730L,
143+
-0.0000000000000000004182784651572477035L,
144+
-0.0000000000000000000261987180876106127L,
145+
-0.0000000000000000000016553211620337322L,
146+
-0.0000000000000000000001053943580037580L,
147+
-0.0000000000000000000000067562822456482L,
148+
-0.0000000000000000000000004357492971838L,
149+
-0.0000000000000000000000000282576222954L,
150+
-0.0000000000000000000000000018415141361L,
151+
-0.0000000000000000000000000001205472536L,
152+
-0.0000000000000000000000000000079233873L,
153+
-0.0000000000000000000000000000005227403L,
154+
-0.0000000000000000000000000000000346060L,
155+
-0.0000000000000000000000000000000022982L,
156+
-0.0000000000000000000000000000000001531L,
157+
-0.0000000000000000000000000000000000102L,
158+
-0.0000000000000000000000000000000000007L
159+
};
49160

50-
while (x >= 2*PI) {
51-
x -= 2*PI;
52-
}
161+
long double h = 0;
162+
long double v = std::fmod(std::abs(x), PI2);
163+
long double sgn = x >= 0 ? 1 : -1;
53164

54-
while (x < 0.0L) {
55-
x += 2*PI;
165+
if (v > PI) {
166+
const long double p0 = 6.28125L;
167+
const long double p1 = 0.0019353071795864769252867665590057683943L;
168+
v = (p0 - v) + p1;
169+
sgn = -sgn;
56170
}
57171

58-
if (std::abs(x) < std::numeric_limits<long double>::epsilon() ||
59-
std::abs(x - PI) < std::numeric_limits<long double>::epsilon() ||
60-
std::abs(x - 2*PI) < std::numeric_limits<long double>::epsilon()) {
61-
return 0.0L;
172+
if (v == 0 || v == PI) {
173+
h = 0;
174+
} else if (v < PIH) {
175+
const long double u = RPIH*v;
176+
h = 2*u*u - 1;
177+
const long double alfa = h + h;
178+
long double b0 = 0, b1 = 0, b2 = 0;
179+
for (int i = 18; i >= 0; i--) {
180+
b0 = A[i] + alfa*b1 - b2;
181+
b2 = b1;
182+
b1 = b0;
183+
}
184+
h = v*(1 - std::log(v) + v*v*(b0 - h*b2)/2);
185+
} else {
186+
const long double u = RPIH*v - 2;
187+
h = 2*u*u - 1;
188+
const long double alfa = h + h;
189+
long double b0 = 0, b1 = 0, b2 = 0;
190+
for (int i = 29; i >= 0; i--) {
191+
b0 = B[i] + alfa*b1 - b2;
192+
b2 = b1;
193+
b1 = b0;
194+
}
195+
h = (PI - v)*(b0 - h*b2);
62196
}
63197

64-
return std::imag(Li2(std::exp(i*x)));
198+
return sgn*h;
65199
}
66200

67201
} // namespace polylogarithm

‎src/cpp/Cl3.cpp

+4-25
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,6 @@
77
#include "Cl3.hpp"
88
#include "Li3.hpp"
99
#include <cmath>
10-
#include <limits>
1110

1211
namespace polylogarithm {
1312

@@ -18,18 +17,8 @@ namespace polylogarithm {
1817
*/
1918
double Cl3(double x) noexcept
2019
{
21-
const double PI = 3.141592653589793;
22-
const std::complex<double> i(0.0, 1.0);
23-
24-
while (x >= 2*PI) {
25-
x -= 2*PI;
26-
}
27-
28-
while (x < 0.0) {
29-
x += 2*PI;
30-
}
31-
32-
return std::real(Li3(std::exp(i*x)));
20+
const std::complex<double> ix(0.0, x);
21+
return std::real(Li3(std::exp(ix)));
3322
}
3423

3524
/**
@@ -39,18 +28,8 @@ double Cl3(double x) noexcept
3928
*/
4029
long double Cl3(long double x) noexcept
4130
{
42-
const long double PI = 3.14159265358979323846264338327950288L;
43-
const std::complex<long double> i(0.0L, 1.0L);
44-
45-
while (x >= 2*PI) {
46-
x -= 2*PI;
47-
}
48-
49-
while (x < 0.0L) {
50-
x += 2*PI;
51-
}
52-
53-
return std::real(Li3(std::exp(i*x)));
31+
const std::complex<long double> ix(0.0L, x);
32+
return std::real(Li3(std::exp(ix)));
5433
}
5534

5635
} // namespace polylogarithm

‎src/cpp/Cl4.cpp

+4-37
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,6 @@
77
#include "Cl4.hpp"
88
#include "Li4.hpp"
99
#include <cmath>
10-
#include <limits>
1110

1211
namespace polylogarithm {
1312

@@ -18,24 +17,8 @@ namespace polylogarithm {
1817
*/
1918
double Cl4(double x) noexcept
2019
{
21-
const double PI = 3.141592653589793;
22-
const std::complex<double> i(0.0, 1.0);
23-
24-
while (x >= 2*PI) {
25-
x -= 2*PI;
26-
}
27-
28-
while (x < 0.0) {
29-
x += 2*PI;
30-
}
31-
32-
if (std::abs(x) < std::numeric_limits<double>::epsilon() ||
33-
std::abs(x - PI) < std::numeric_limits<double>::epsilon() ||
34-
std::abs(x - 2*PI) < std::numeric_limits<double>::epsilon()) {
35-
return 0.0;
36-
}
37-
38-
return std::imag(Li4(std::exp(i*x)));
20+
const std::complex<double> ix(0.0, x);
21+
return std::imag(Li4(std::exp(ix)));
3922
}
4023

4124
/**
@@ -45,24 +28,8 @@ double Cl4(double x) noexcept
4528
*/
4629
long double Cl4(long double x) noexcept
4730
{
48-
const long double PI = 3.14159265358979323846264338327950288L;
49-
const std::complex<long double> i(0.0L, 1.0L);
50-
51-
while (x >= 2*PI) {
52-
x -= 2*PI;
53-
}
54-
55-
while (x < 0.0L) {
56-
x += 2*PI;
57-
}
58-
59-
if (std::abs(x) < std::numeric_limits<long double>::epsilon() ||
60-
std::abs(x - PI) < std::numeric_limits<long double>::epsilon() ||
61-
std::abs(x - 2*PI) < std::numeric_limits<long double>::epsilon()) {
62-
return 0.0L;
63-
}
64-
65-
return std::imag(Li4(std::exp(i*x)));
31+
const std::complex<long double> ix(0.0L, x);
32+
return std::imag(Li4(std::exp(ix)));
6633
}
6734

6835
} // namespace polylogarithm

‎src/cpp/Cl5.cpp

+4-25
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,6 @@
77
#include "Cl5.hpp"
88
#include "Li5.hpp"
99
#include <cmath>
10-
#include <limits>
1110

1211
namespace polylogarithm {
1312

@@ -18,18 +17,8 @@ namespace polylogarithm {
1817
*/
1918
double Cl5(double x) noexcept
2019
{
21-
const double PI = 3.141592653589793;
22-
const std::complex<double> i(0.0, 1.0);
23-
24-
while (x >= 2*PI) {
25-
x -= 2*PI;
26-
}
27-
28-
while (x < 0.0) {
29-
x += 2*PI;
30-
}
31-
32-
return std::real(Li5(std::exp(i*x)));
20+
const std::complex<double> ix(0.0, x);
21+
return std::real(Li5(std::exp(ix)));
3322
}
3423

3524
/**
@@ -39,18 +28,8 @@ double Cl5(double x) noexcept
3928
*/
4029
long double Cl5(long double x) noexcept
4130
{
42-
const long double PI = 3.14159265358979323846264338327950288L;
43-
const std::complex<long double> i(0.0L, 1.0L);
44-
45-
while (x >= 2*PI) {
46-
x -= 2*PI;
47-
}
48-
49-
while (x < 0.0L) {
50-
x += 2*PI;
51-
}
52-
53-
return std::real(Li5(std::exp(i*x)));
31+
const std::complex<long double> ix(0.0L, x);
32+
return std::real(Li5(std::exp(ix)));
5433
}
5534

5635
} // namespace polylogarithm

0 commit comments

Comments
 (0)