Hyperiso 1.0.3
Modular flavour-physics calculations, Wilson coefficients and statistical inference
Loading...
Searching...
No Matches
polylog.cpp
Go to the documentation of this file.
1#include <cmath>
2#include <complex>
3#include <array>
4#include <vector>
5#include "special_generic.h"
6#include "gsl/gsl_sf_dilog.h"
7#include "gsl/gsl_sf_clausen.h"
8
9scalar_t cd(double x, double y) {
10 return {x, y};
11}
12
13scalar_t hpl_base1(int i, scalar_t x) {
14 if (i == 0) return std::log(x);
15 if (i == 1) return -std::log(1.0 - x);
16 return {0.0, 0.0};
17}
18
19scalar_t hpl_base2(int i1, int i2, scalar_t x) {
20 scalar_t u = std::log(1.0 - x);
21
22 if (i1 == 0 && i2 == 1) {
23
24 return -u - 0.25 * std::pow(u, 2) - 0.027777777777777776 * std::pow(u, 3) +
25 0.0002777777777777778 * std::pow(u, 5) -
26 4.72411186696901e-6 * std::pow(u, 7) +
27 9.185773074661964e-8 * std::pow(u, 9) -
28 1.8978869988971e-9 * std::pow(u, 11) +
29 4.0647616451442256e-11 * std::pow(u, 13) -
30 8.921691020456452e-13 * std::pow(u, 15) +
31 1.9939295860721074e-14 * std::pow(u, 17) -
32 4.518980029619918e-16 * std::pow(u, 19) +
33 1.0356517612181247e-17 * std::pow(u, 21);
34 }
35
36 return {0.0, 0.0};
37}
38
39
40scalar_t hpl_base3(int i1, int i2, int i3, scalar_t x) {
41 scalar_t u = std::log(1.0 - x);
42
43 if (i1 == 0 && i2 == 0 && i3 == 1) {
44
45 return -u - 0.375 * std::pow(u, 2) - 0.0787037037037037 * std::pow(u, 3) -
46 0.008680555555555556 * std::pow(u, 4) -
47 0.00012962962962962963 * std::pow(u, 5) +
48 0.00008101851851851852 * std::pow(u, 6) +
49 3.4193571608537595e-6 * std::pow(u, 7) -
50 1.328656462585034e-6 * std::pow(u, 8) -
51 8.660871756109851e-8 * std::pow(u, 9) +
52 2.52608759553204e-8 * std::pow(u, 10) +
53 2.144694468364065e-9 * std::pow(u, 11) -
54 5.140110622012979e-10 * std::pow(u, 12) -
55 5.24958211460083e-11 * std::pow(u, 13) +
56 1.0887754406636318e-11 * std::pow(u, 14) +
57 1.2779396094493695e-12 * std::pow(u, 15) -
58 2.369824177308745e-13 * std::pow(u, 16) -
59 3.104357887965462e-14 * std::pow(u, 17) +
60 5.261758629912506e-15 * std::pow(u, 18) +
61 7.538479549949265e-16 * std::pow(u, 19) -
62 1.1862322577752286e-16 * std::pow(u, 20) -
63 1.8316979965491384e-17 * std::pow(u, 21);
64 }
65
66 if (i1 == 0 && i2 == 1 && i3 == 1) {
67 return 0.25 * std::pow(u, 2) + 0.08333333333333333 * std::pow(u, 3) +
68 0.010416666666666666 * std::pow(u, 4) -
69 0.00011574074074074075 * std::pow(u, 6) +
70 2.066798941798942e-6 * std::pow(u, 8) -
71 4.1335978835978836e-8 * std::pow(u, 10) +
72 8.698648744945042e-10 * std::pow(u, 12) -
73 1.887210763816962e-11 * std::pow(u, 14) +
74 4.182042665838962e-13 * std::pow(u, 16) -
75 9.415778600896063e-15 * std::pow(u, 18) +
76 2.146515514069461e-16 * std::pow(u, 20);
77 }
78
79 return {0.0, 0.0};
80}
81
82scalar_t hpl_base4(int i1, int i2, int i3, int i4, scalar_t x) {
83 scalar_t u = std::log(1.0 - x);
84
85 if (i1 == 0 && i2 == 0 && i3 == 0 && i4 == 1) {
86 return -1.0 * u - 0.4375 * std::pow(u, 2) - 0.11651234567901235 * std::pow(u, 3) -
87 0.019820601851851853 * std::pow(u, 4) - 0.001927932098765432 * std::pow(u, 5) -
88 0.000031057098765432096 * std::pow(u, 6) + 0.000015624009114857836 * std::pow(u, 7) +
89 8.485123546773206e-7 * std::pow(u, 8) - 2.290961660318971e-7 * std::pow(u, 9) -
90 2.1832614218526917e-8 * std::pow(u, 10) + 3.882824879172015e-9 * std::pow(u, 11) +
91 5.446292103220332e-10 * std::pow(u, 12) - 6.960805210682725e-11 * std::pow(u, 13) -
92 1.3375737686445216e-11 * std::pow(u, 14) + 1.2784852685266572e-12 * std::pow(u, 15) +
93 3.260562858024892e-13 * std::pow(u, 16) - 2.364757116861826e-14 * std::pow(u, 17) -
94 7.923135122031162e-15 * std::pow(u, 18) + 4.3452915709984186e-16 * std::pow(u, 19) +
95 1.923627006253592e-16 * std::pow(u, 20) - 7.812414333195955e-18 * std::pow(u, 21);
96 }
97
98 if (i1 == 0 && i2 == 1 && i3 == 0 && i4 == 1) {
99 return 0.25 * std::pow(u, 2) + 0.1111111111111111 * std::pow(u, 3) +
100 0.022569444444444444 * std::pow(u, 4) + 0.0020833333333333333 * std::pow(u, 5) -
101 0.000027006172839506174 * std::pow(u, 6) - 0.00001984126984126984 * std::pow(u, 7) +
102 4.527273872511968e-7 * std::pow(u, 8) + 3.389987682315725e-7 * std::pow(u, 9) -
103 7.939132443100697e-9 * std::pow(u, 10) - 6.6805622361177916e-9 * std::pow(u, 11) +
104 1.4490216610627064e-10 * std::pow(u, 12) + 1.39908336457158e-10 * std::pow(u, 13) -
105 2.7425719106565973e-12 * std::pow(u, 14) - 3.032441227329819e-12 * std::pow(u, 15) +
106 5.358569182999823e-14 * std::pow(u, 16) + 6.724068599976371e-14 * std::pow(u, 17) -
107 1.0756816626218996e-15 * std::pow(u, 18) - 1.5158529016922455e-15 * std::pow(u, 19) +
108 2.208955024323606e-17 * std::pow(u, 20) + 3.460964863954937e-17 * std::pow(u, 21);
109 }
110
111 if (i1 == 0 && i2 == 1 && i3 == 1 && i4 == 1) {
112 return -0.05555555555555555 * std::pow(u, 3) - 0.020833333333333332 * std::pow(u, 4) -
113 0.002777777777777778 * std::pow(u, 5) + 0.00003306878306878307 * std::pow(u, 7) -
114 6.123848716441309e-7 * std::pow(u, 9) + 1.252605419272086e-8 * std::pow(u, 11) -
115 2.6765073061369356e-10 * std::pow(u, 13) + 5.871322376319437e-12 * std::pow(u, 15) -
116 1.312013385361243e-13 * std::pow(u, 17) + 2.97340376870402e-15 * std::pow(u, 19) -
117 6.814334965299877e-17 * std::pow(u, 21);
118 }
119
120 return {0.0, 0.0};
121}
122
123scalar_t hpl1(int i, scalar_t x) {
124 if (i == 0) return std::log(x);
125 if (i == 1) return -std::log(1.0 - x);
126 return {0.0, 0.0};
127}
128
129scalar_t hpl2(int i1, int i2, scalar_t x) {
130 const scalar_t pi = std::acos(-1);
131 const scalar_t i(0.0, 1.0);
132
133 if (i1 == 0 && i2 == 0) {
134 if (std::abs(x) > 1) return std::pow(hpl1(0, 1.0 / x), 2) / 2.0;
135 if (std::real(x) > 0.5) return std::pow(hpl1(1, 1.0 - x), 2) / 2.0;
136 return std::pow(hpl_base1(0, x), 2) / 2.0;
137 }
138
139 if (i1 == 0 && i2 == 1) {
140 if (std::abs(x) > 1) return
141 pi * i * hpl1(0, 1.0 / x) - hpl2(0, 1, 1.0 / x) + std::pow(pi, 2) / 3.0 -
142 std::pow(hpl1(0, 1.0 / x), 2) / 2.0;
143
144 if (std::real(x) > 0.5) return
145 hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) - hpl2(0, 1, 1.0 - x) + std::pow(pi, 2) / 6.0;
146
147 return hpl_base2(0, 1, x);
148 }
149
150 if (i1 == 1 && i2 == 0) {
151 if (std::abs(x) > 1) return
152 pi * -i * hpl1(0, 1.0 / x) -
153 hpl1(0, 1.0 / x) * (pi * -i + hpl1(0, 1.0 / x) + hpl1(1, 1.0 / x)) +
154 hpl2(0, 1, 1.0 / x) - std::pow(pi, 2) / 3.0 + std::pow(hpl1(0, 1.0 / x), 2) / 2.0;
155
156 if (std::real(x) > 0.5) return hpl2(0, 1, 1.0 - x) - std::pow(pi, 2) / 6.0;
157
158 return hpl_base1(0, x) * hpl_base1(1, x) - hpl_base2(0, 1, x);
159 }
160
161 if (i1 == 1 && i2 == 1) {
162 if (std::abs(x) > 1) return std::pow(pi * i + hpl1(0, 1.0 / x) + hpl1(1, 1.0 / x), 2) / 2.0;
163 if (std::real(x) > 0.5) return std::pow(hpl1(0, 1.0 - x), 2) / 2.0;
164 return std::pow(hpl_base1(1, x), 2) / 2.0;
165 }
166
167 return {0.0, 0.0};
168}
169
170scalar_t hpl3(int i1, int i2, int i3, scalar_t x) {
171 const scalar_t pi = std::acos(-1);
172 const scalar_t i(0.0, 1.0);
173 const double zeta3 = 1.2020569031595942; // Apery's constant
174
175 if (i1 == 0 && i2 == 0 && i3 == 0) {
176 if (std::abs(x) > 1) return -std::pow(hpl1(0, 1.0 / x), 3) / 6.0;
177 if (std::real(x) > 0.5) return -std::pow(hpl1(1, 1.0 - x), 3) / 6.0;
178 return std::pow(hpl_base1(0, x), 3) / 6.0;
179 }
180
181 if (i1 == 0 && i2 == 0 && i3 == 1) {
182 if (std::abs(x) > 1) return hpl3(0, 0, 1, 1.0 / x) - (hpl1(0, 1.0 / x) * std::pow(pi, 2)) / 3.0 + pi * i * 0.5 * std::pow(hpl1(0, 1.0 / x), 2) + std::pow(hpl1(0, 1.0 / x), 3) / 6.0;
183 if (std::real(x) > 0.5) return zeta3 + hpl1(1, 1.0 - x) * hpl2(0, 1, 1.0 - x) - hpl3(0, 1, 1, 1.0 - x) - (hpl1(1, 1.0 - x) * std::pow(pi, 2)) / 6.0 - (hpl1(0, 1.0 - x) * std::pow(hpl1(1, 1.0 - x), 2)) / 2.0;
184 return hpl_base3(0, 0, 1, x);
185 }
186
187 if (i1 == 0 && i2 == 1 && i3 == 0) {
188 if (abs(x) > 1) return -(hpl1(0, 1.0 / x) * (-pi * i * hpl1(0, 1.0 / x) - hpl2(0, 1, 1.0 / x) + pow(pi, 2) / 3.0 - pow(hpl1(0, 1.0 / x), 2) / 2.0)) - 2.0 * (hpl3(0, 0, 1, 1.0 / x) - (hpl1(0, 1.0 / x) * pow(pi, 2)) / 3.0 + pi * i * 0.5 * pow(hpl1(0, 1.0 / x), 2) + pow(hpl1(0, 1.0 / x), 3) / 6.0);
189 if (real(x) > 0.5) return -(hpl1(1, 1.0 - x) * (hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) - hpl2(0, 1, 1.0 - x) + pow(pi, 2) / 6.0)) - 2.0 * (1.2020569031595942 + hpl1(1, 1.0 - x) * hpl2(0, 1, 1.0 - x) - hpl3(0, 1, 1, 1.0 - x) - (hpl1(1, 1.0 - x) * pow(pi, 2)) / 6.0 - (hpl1(0, 1.0 - x) * pow(hpl1(1, 1.0 - x), 2)) / 2.0);
190 return hpl_base1(0, x) * hpl_base2(0, 1, x) - 2.0 * hpl_base3(0, 0, 1, x);
191}
192
193 if (i1 == 0 && i2 == 1 && i3 == 1) {
194 if (abs(x) > 1) return 1.2020569031595942 - pi * i * hpl2(0, 1, 1.0 / x) - hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) + hpl3(0, 0, 1, 1.0 / x) - hpl3(0, 1, 1, 1.0 / x) + (hpl1(0, 1.0 / x) * pow(pi, 2)) / 2.0 + i * 0.16666666666666666 * pow(pi, 3) + pi * i * -0.5 * pow(hpl1(0, 1.0 / x), 2) - pow(hpl1(0, 1.0 / x), 3) / 6.0;
195 if (real(x) > 0.5) return 1.2020569031595942 + hpl1(0, 1.0 - x) * hpl2(0, 1, 1.0 - x) - hpl3(0, 0, 1, 1.0 - x) - (hpl1(1, 1.0 - x) * pow(hpl1(0, 1.0 - x), 2)) / 2.0;
196 return hpl_base3(0, 1, 1, x);
197}
198 if (i1 == 1 && i2 == 0 && i3 == 0) {
199 if (abs(x) > 1) return hpl3(0, 0, 1, 1.0 / x) - (hpl1(0, 1.0 / x) * pow(pi, 2)) / 3.0 + hpl1(0, 1.0 / x) * (-pi * i * hpl1(0, 1.0 / x) - hpl2(0, 1, 1.0 / x) + pow(pi, 2) / 3.0 - pow(hpl1(0, 1.0 / x), 2) / 2.0) + pi * i * 0.5 * pow(hpl1(0, 1.0 / x), 2) + ((pi * i + hpl1(0, 1.0 / x) + hpl1(1, 1.0 / x)) * pow(hpl1(0, 1.0 / x), 2)) / 2.0 + pow(hpl1(0, 1.0 / x), 3) / 6.0;
200 if (real(x) > 0.5) return 1.2020569031595942 + hpl1(1, 1.0 - x) * hpl2(0, 1, 1.0 - x) - hpl3(0, 1, 1, 1.0 - x) + hpl1(1, 1.0 - x) * (hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) - hpl2(0, 1, 1.0 - x) + pow(pi, 2) / 6.0) - (hpl1(1, 1.0 - x) * pow(pi, 2)) / 6.0 - hpl1(0, 1.0 - x) * pow(hpl1(1, 1.0 - x), 2);
201 return -(hpl_base1(0, x) * hpl_base2(0, 1, x)) + hpl_base3(0, 0, 1, x) + (hpl_base1(1, x) * pow(hpl_base1(0, x), 2)) / 2.0;
202}
203
204 if (i1 == 1 && i2 == 0 && i3 == 1) {
205 if (abs(x) > 1) return (pi * i + hpl1(0, 1.0 / x) + hpl1(1, 1.0 / x)) * (pi * i * -1. * hpl1(0, 1.0 / x) - hpl2(0, 1, 1.0 / x) + pow(pi, 2) / 3.0 - pow(hpl1(0, 1.0 / x), 2) / 2.0) - 2.0 * (1.2020569031595942 + pi * i * -1. * hpl2(0, 1, 1.0 / x) - hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) + hpl3(0, 0, 1, 1.0 / x) - hpl3(0, 1, 1, 1.0 / x) + (hpl1(0, 1.0 / x) * pow(pi, 2)) / 2.0 + i * 0.16666666666666666 * pow(pi, 3) + pi * i * -0.5 * pow(hpl1(0, 1.0 / x), 2) - pow(hpl1(0, 1.0 / x), 3) / 6.0);
206 if (real(x) > 0.5) return -(hpl1(0, 1.0 - x) * (hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) - hpl2(0, 1, 1.0 - x) + pow(pi, 2) / 6.0)) - 2.0 * (1.2020569031595942 + hpl1(0, 1.0 - x) * hpl2(0, 1, 1.0 - x) - hpl3(0, 0, 1, 1.0 - x) - (hpl1(1, 1.0 - x) * pow(hpl1(0, 1.0 - x), 2)) / 2.0);
207 return hpl_base3(0, 0, 1, x) - (hpl_base1(0, x) * pow(pi, 2)) / 3.0 + pi * i * 0.5 * pow(hpl_base1(0, x), 2) + pow(hpl_base1(0, x), 3) / 6.0;
208}
209 if (i1 == 1 && i2 == 1 && i3 == 0) {
210 if (abs(x) > 1) return -hpl3(0, 1, 1, 1.0 / x) + hpl2(0, 1, 1.0 / x) * (pi * i + hpl1(0, 1.0 / x) + hpl1(1, 1.0 / x)) - hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) + (hpl1(0, 1.0 / x) * pow(pi, 2)) / 2.0 - pi * i * 0.5 * pow(hpl1(0, 1.0 / x), 2) - pow(hpl1(0, 1.0 / x), 3) / 6.0;
211 if (real(x) > 0.5) return -hpl3(0, 1, 1, 1.0 - x) + hpl2(0, 1, 1.0 - x) * hpl1(1, 1.0 - x) - hpl1(0, 1.0 - x) * hpl2(0, 1, 1.0 - x) + (hpl1(0, 1.0 - x) * pow(pi, 2)) / 6.0 - (hpl1(1, 1.0 - x) * pow(hpl1(0, 1.0 - x), 2)) / 2.0;
212 return hpl_base3(0, 1, 0, x);
213}
214
215 if (i1 == 1 && i2 == 1 && i3 == 1) {
216 if (abs(x) > 1) return -hpl3(0, 1, 1, 1.0 / x) + (pi * i + hpl1(0, 1.0 / x) + hpl1(1, 1.0 / x)) * hpl2(0, 1, 1.0 / x) - hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) + (hpl1(0, 1.0 / x) * pow(pi, 2)) / 2.0 + i * 0.16666666666666666 * pow(pi, 3) + pi * i * -0.5 * pow(hpl1(0, 1.0 / x), 2) - pow(hpl1(0, 1.0 / x), 3) / 6.0;
217 if (real(x) > 0.5) return -hpl3(0, 1, 1, 1.0 - x) + hpl1(0, 1.0 - x) * hpl2(0, 1, 1.0 - x) - (hpl1(1, 1.0 - x) * pow(hpl1(0, 1.0 - x), 2)) / 2.0 + 1.2020569031595942;
218 return hpl_base3(0, 1, 1, x);
219}
220
221 return {0.0, 0.0}; // Default return for unhandled cases
222}
223
224scalar_t hpl4(int i1, int i2, int i3, int i4, scalar_t x) {
225 const scalar_t i(0, 1);
226 const double pi = acos(-1);
227
228 if (i1 == 0 && i2 == 0 && i3 == 0 && i4 == 0) {
229 if (abs(x) > 1) return pow(hpl1(0, 1.0 / x), 4) / 24.0;
230 if (real(x) > 0.5) return pow(hpl1(1, 1.0 - x), 4) / 24.0;
231 return pow(hpl_base1(0, x), 4) / 24.0;
232 }
233
234 if (i1 == 0 && i2 == 0 && i3 == 0 && i4 == 1) {
235 if (abs(x) > 1) return -hpl4(0, 0, 0, 1, 1.0 / x) + pow(pi, 4) / 45.0 + (pow(pi, 2) * pow(hpl1(0, 1.0 / x), 2)) / 6.0 + pi * i * -0.16666666666666666 * pow(hpl1(0, 1.0 / x), 3) - pow(hpl1(0, 1.0 / x), 4) / 24.0;
236 if (real(x) > 0.5) return -1.2020569031595942 * hpl1(1, 1.0 - x) + hpl1(1, 1.0 - x) * hpl3(0, 1, 1, 1.0 - x) - hpl4(0, 1, 1, 1, 1.0 - x) + pow(pi, 4) / 90.0 - (hpl2(0, 1, 1.0 - x) * pow(hpl1(1, 1.0 - x), 2)) / 2.0 + (pow(pi, 2) * pow(hpl1(1, 1.0 - x), 2)) / 12.0 + (hpl1(0, 1.0 - x) * pow(hpl1(1, 1.0 - x), 3)) / 6.0;
237 return hpl_base4(0, 0, 0, 1, x);
238 }
239
240 if (i1 == 0 && i2 == 0 && i3 == 1 && i4 == 0) {
241 if (abs(x) > 1) return -(hpl1(0, 1.0 / x) * (hpl3(0, 0, 1, 1.0 / x) - (hpl1(0, 1.0 / x) * pow(pi, 2)) / 3.0 + pi * i * 0.5 * pow(hpl1(0, 1.0 / x), 2) + pow(hpl1(0, 1.0 / x), 3) / 6.0)) - 3. * (-hpl4(0, 0, 0, 1, 1.0 / x) + pow(pi, 4) / 45.0 + (pow(pi, 2) * pow(hpl1(0, 1.0 / x), 2)) / 6.0 + pi * i * -0.16666666666666666 * pow(hpl1(0, 1.0 / x), 3) - pow(hpl1(0, 1.0 / x), 4) / 24.0);
242 if (real(x) > 0.5) return -(hpl1(1, 1.0 - x) * (1.2020569031595942 + hpl1(1, 1.0 - x) * hpl2(0, 1, 1.0 - x) - hpl3(0, 1, 1, 1.0 - x) - (hpl1(1, 1.0 - x) * pow(pi, 2)) / 6.0 - (hpl1(0, 1.0 - x) * pow(hpl1(1, 1.0 - x), 2)) / 2.0)) - 3. * (-1.2020569031595942 * hpl1(1, 1.0 - x) + hpl1(1, 1.0 - x) * hpl3(0, 1, 1, 1.0 - x) - hpl4(0, 1, 1, 1, 1.0 - x) + pow(pi, 4) / 90.0 - (hpl2(0, 1, 1.0 - x) * pow(hpl1(1, 1.0 - x), 2)) / 2.0 + (pow(pi, 2) * pow(hpl1(1, 1.0 - x), 2)) / 12.0 + (hpl1(0, 1.0 - x) * pow(hpl1(1, 1.0 - x), 3)) / 6.0);
243 return hpl_base1(0, x) * hpl_base3(0, 0, 1, x) - 3. * hpl_base4(0, 0, 0, 1, x);
244 }
245
246 if (i1 == 0 && i2 == 0 && i3 == 1 && i4 == 1) {
247 if (abs(x) > 1) {
248 return (-2. * (3.7763731361630786 * scalar_t(0, 2) +
249 2.4041138063191885 * hpl1(0, 1.0 / x) +
250 pi * scalar_t(0, 1) * hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) +
251 pi * scalar_t(0, -2) * hpl3(0, 0, 1, 1.0 / x) -
252 2.0 * hpl1(0, 1.0 / x) * hpl3(0, 0, 1, 1.0 / x) +
253 4. * hpl4(0, 0, 0, 1, 1.0 / x) + hpl4(0, 1, 0, 1, 1.0 / x) -
254 (hpl2(0, 1, 1.0 / x) * pow(pi, 2)) / 3.0 + pow(pi, 4) / 90.0 +
255 (hpl2(0, 1, 1.0 / x) * pow(hpl1(0, 1.0 / x), 2)) / 2.0 -
256 (pow(pi, 2) * pow(hpl1(0, 1.0 / x), 2)) / 6.0 +
257 pi * scalar_t(0, 0.16666666666666666) * pow(hpl1(0, 1.0 / x), 3) +
258 pow(hpl1(0, 1.0 / x), 4) / 24.0) +
259 pow(pi * scalar_t(0, -1) * hpl1(0, 1.0 / x) - hpl2(0, 1, 1.0 / x) +
260 pow(pi, 2) / 3.0 - pow(hpl1(0, 1.0 / x), 2) / 2.0, 2)) / 4.0;
261 } else if (real(x) > 0.5) {
262 return (-2. * (2.4041138063191885 * hpl1(1, 1.0 - x) +
263 hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) * hpl2(0, 1, 1.0 - x) -
264 hpl4(0, 1, 0, 1, 1.0 - x) +
265 (hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) * pow(pi, 2)) / 6.0 -
266 (hpl2(0, 1, 1.0 - x) * pow(pi, 2)) / 6.0 +
267 pow(pi, 4) / 120.0 +
268 pow(hpl2(0, 1, 1.0 - x), 2) -
269 2. * (hpl1(1, 1.0 - x) * hpl3(0, 0, 1, 1.0 - x) +
270 (2. * hpl4(0, 1, 0, 1, 1.0 - x) - pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0 +
271 (-2. * hpl4(0, 1, 0, 1, 1.0 - x) + pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0) -
272 2. * (hpl1(0, 1.0 - x) * hpl3(0, 1, 1, 1.0 - x) +
273 (2. * hpl4(0, 1, 0, 1, 1.0 - x) - pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0 +
274 (-2. * hpl4(0, 1, 0, 1, 1.0 - x) + pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0)) +
275 pow(hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) - hpl2(0, 1, 1.0 - x) +
276 pow(pi, 2) / 6.0, 2)) / 4.0;
277 } else {
278 return (-2. * hpl_base4(0, 1, 0, 1, x) + pow(hpl_base2(0, 1, x), 2)) / 4.0;
279 }
280 }
281
282 if (i1 == 0 && i2 == 1 && i3 == 0 && i4 == 0) {
283 if (abs(x) > 1) {
284 return (((pi * scalar_t(0, -1) * hpl1(0, 1.0 / x) - hpl2(0, 1, 1.0 / x) + pow(pi, 2) / 3.0 -
285 pow(hpl1(0, 1.0 / x), 2) / 2.0) * pow(hpl1(0, 1.0 / x), 2) +
286 4.0 * hpl1(0, 1.0 / x) * (hpl3(0, 0, 1, 1.0 / x) -
287 (hpl1(0, 1.0 / x) * pow(pi, 2)) / 3.0 +
288 pi * scalar_t(0, 0.5) * pow(hpl1(0, 1.0 / x), 2) +
289 pow(hpl1(0, 1.0 / x), 3) / 6.0) +
290 6.0 * (-hpl4(0, 0, 0, 1, 1.0 / x) + pow(pi, 4) / 45.0 +
291 (pow(pi, 2) * pow(hpl1(0, 1.0 / x), 2)) / 6.0 +
292 pi * scalar_t(0, -0.16666666666666666) * pow(hpl1(0, 1.0 / x), 3) -
293 pow(hpl1(0, 1.0 / x), 4) / 24.0)) / 2.0);
294 } else if (real(x) > 0.5) {
295 return ((hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) - hpl2(0, 1, 1.0 - x) + pow(pi, 2) / 6.0) * pow(hpl1(1, 1.0 - x), 2) +
296 4.0 * hpl1(1, 1.0 - x) * (1.2020569031595942 + hpl1(1, 1.0 - x) * hpl2(0, 1, 1.0 - x) - hpl3(0, 1, 1, 1.0 - x) -
297 (hpl1(1, 1.0 - x) * pow(pi, 2)) / 6.0 - (hpl1(0, 1.0 - x) * pow(hpl1(1, 1.0 - x), 2)) / 2.0) +
298 6.0 * (-1.2020569031595942 * hpl1(1, 1.0 - x) + hpl1(1, 1.0 - x) * hpl3(0, 1, 1, 1.0 - x) -
299 hpl4(0, 1, 1, 1, 1.0 - x) + pow(pi, 4) / 90.0 - (hpl2(0, 1, 1.0 - x) * pow(hpl1(1, 1.0 - x), 2)) / 2.0 +
300 (pow(pi, 2) * pow(hpl1(1, 1.0 - x), 2)) / 12.0 + (hpl1(0, 1.0 - x) * pow(hpl1(1, 1.0 - x), 3)) / 6.0)) / 2.0;
301 } else {
302
303 return
304 (-4.*hpl_base1(0,x)*hpl_base3(0,0,1,x) + 6.*hpl_base4(0,0,0,1,x) +
305 hpl_base2(0,1,x)*pow(hpl_base1(0,x),2))/2.;
306 }
307}
308 if (i1 == 0 && i2 == 1 && i3 == 0 && i4 == 1) {
309 if (abs(x) > 1) {
310 return 3.7763731361630786 * scalar_t(0, 2) +
311 2.4041138063191885 * hpl1(0, 1.0 / x) +
312 pi * scalar_t(0, 1) * hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) +
313 pi * scalar_t(0, -2) * hpl3(0, 0, 1, 1.0 / x) -
314 2.0 * hpl1(0, 1.0 / x) * hpl3(0, 0, 1, 1.0 / x) +
315 4.0 * hpl4(0, 0, 0, 1, 1.0 / x) +
316 hpl4(0, 1, 0, 1, 1.0 / x) -
317 (hpl2(0, 1, 1.0 / x) * pow(pi, 2)) / 3.0 +
318 pow(pi, 4) / 90.0 +
319 (hpl2(0, 1, 1.0 / x) * pow(hpl1(0, 1.0 / x), 2)) / 2.0 -
320 (pow(pi, 2) * pow(hpl1(0, 1.0 / x), 2)) / 6.0 +
321 pi * scalar_t(0, 0.16666666666666666) * pow(hpl1(0, 1.0 / x), 3) +
322 pow(hpl1(0, 1.0 / x), 4) / 24.0;
323 } else if (real(x) > 0.5) {
324 return 2.4041138063191885 * hpl1(1, 1.0 - x) +
325 hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) * hpl2(0, 1, 1.0 - x) -
326 hpl4(0, 1, 0, 1, 1.0 - x) +
327 (hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) * pow(pi, 2)) / 6.0 -
328 (hpl2(0, 1, 1.0 - x) * pow(pi, 2)) / 6.0 +
329 pow(pi, 4) / 120.0 +
330 pow(hpl2(0, 1, 1.0 - x), 2) -
331 2.0 * (hpl1(1, 1.0 - x) * hpl3(0, 0, 1, 1.0 - x) +
332 (2.0 * hpl4(0, 1, 0, 1, 1.0 - x) - pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0 +
333 (-2.0 * hpl4(0, 1, 0, 1, 1.0 - x) + pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0) -
334 2.0 * (hpl1(0, 1.0 - x) * hpl3(0, 1, 1, 1.0 - x) +
335 (2.0 * hpl4(0, 1, 0, 1, 1.0 - x) - pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0 +
336 (-2.0 * hpl4(0, 1, 0, 1, 1.0 - x) + pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0);
337 } else {
338 return hpl_base4(0,1,0,1,x);
339}
340 }
341
342 if (i1 == 0 && i2 == 1 && i3 == 1 && i4 == 0) {
343 if (abs(x) > 1) {
344 return 3.7763731361630786 * scalar_t(0, -2) -
345 2.4041138063191885 * hpl1(0, 1.0 / x) +
346 pi * scalar_t(0, -1) * hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) +
347 pi * scalar_t(0, 2) * hpl3(0, 0, 1, 1.0 / x) +
348 2.0 * hpl1(0, 1.0 / x) * hpl3(0, 0, 1, 1.0 / x) -
349 4. * hpl4(0, 0, 0, 1, 1.0 / x) -
350 hpl4(0, 1, 0, 1, 1.0 / x) +
351 (hpl2(0, 1, 1.0 / x) * pow(pi, 2)) / 3.0 -
352 pow(pi, 4) / 90.0 -
353 (hpl2(0, 1, 1.0 / x) * pow(hpl1(0, 1.0 / x), 2)) / 2.0 +
354 (pow(pi, 2) * pow(hpl1(0, 1.0 / x), 2)) / 6.0 -
355 hpl1(0, 1.0 / x) * (1.2020569031595942 +
356 pi * scalar_t(0, -1) * hpl2(0, 1, 1.0 / x) -
357 hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) +
358 hpl3(0, 0, 1, 1.0 / x) -
359 hpl3(0, 1, 1, 1.0 / x) +
360 (hpl1(0, 1.0 / x) * pow(pi, 2)) / 2.0 +
361 scalar_t(0, 0.16666666666666666) * pow(pi, 3) +
362 pi * scalar_t(0, -0.5) * pow(hpl1(0, 1.0 / x), 2) -
363 pow(hpl1(0, 1.0 / x), 3) / 6.0) +
364 pi * scalar_t(0, -0.16666666666666666) * pow(hpl1(0, 1.0 / x), 3) -
365 pow(hpl1(0, 1.0 / x), 4) / 24.0 +
366 (2.0 * (3.7763731361630786 * scalar_t(0, 2) +
367 2.4041138063191885 * hpl1(0, 1.0 / x) +
368 pi * scalar_t(0, 1) * hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) +
369 pi * scalar_t(0, -2) * hpl3(0, 0, 1, 1.0 / x) -
370 2.0 * hpl1(0, 1.0 / x) * hpl3(0, 0, 1, 1.0 / x) +
371 4. * hpl4(0, 0, 0, 1, 1.0 / x) +
372 hpl4(0, 1, 0, 1, 1.0 / x) -
373 (hpl2(0, 1, 1.0 / x) * pow(pi, 2)) / 3.0 +
374 pow(pi, 4) / 90.0 +
375 (hpl2(0, 1, 1.0 / x) * pow(hpl1(0, 1.0 / x), 2)) / 2.0 -
376 (pow(pi, 2) * pow(hpl1(0, 1.0 / x), 2)) / 6.0 +
377 pi * scalar_t(0, 0.16666666666666666) * pow(hpl1(0, 1.0 / x), 3) +
378 pow(hpl1(0, 1.0 / x), 4) / 24.0) -
379 pow(pi * scalar_t(0, -1) * hpl1(0, 1.0 / x) - hpl2(0, 1, 1.0 / x) +
380 pow(pi, 2) / 3.0 - pow(hpl1(0, 1.0 / x), 2) / 2.0, 2)) / 2.0;
381 } else if (real(x) > 0.5) {
382 return -2.4041138063191885 * hpl1(1, 1.0 - x) -
383 hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) * hpl2(0, 1, 1.0 - x) +
384 hpl4(0, 1, 0, 1, 1.0 - x) -
385 (hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) * pow(pi, 2)) / 6.0 +
386 (hpl2(0, 1, 1.0 - x) * pow(pi, 2)) / 6.0 -
387 pow(pi, 4) / 120.0 -
388 hpl1(1, 1.0 - x) * (1.2020569031595942 +
389 hpl1(0, 1.0 - x) * hpl2(0, 1, 1.0 - x) -
390 hpl3(0, 0, 1, 1.0 - x) -
391 (hpl1(1, 1.0 - x) * pow(hpl1(0, 1.0 - x), 2)) / 2.0) -
392 pow(hpl2(0, 1, 1.0 - x), 2) +
393 2.0 * (hpl1(1, 1.0 - x) * hpl3(0, 0, 1, 1.0 - x) +
394 (2.0 * hpl4(0, 1, 0, 1, 1.0 - x) - pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0 +
395 (-2.0 * hpl4(0, 1, 0, 1, 1.0 - x) + pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0) +
396 2.0 * (hpl1(0, 1.0 - x) * hpl3(0, 1, 1, 1.0 - x) +
397 (2.0 * hpl4(0, 1, 0, 1, 1.0 - x) - pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0 +
398 (-2.0 * hpl4(0, 1, 0, 1, 1.0 - x) + pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0) +
399 (2.0 * (2.4041138063191885 * hpl1(1, 1.0 - x) +
400 hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) * hpl2(0, 1, 1.0 - x) -
401 hpl4(0, 1, 0, 1, 1.0 - x) +
402 (hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) * pow(pi, 2)) / 6.0 -
403 (hpl2(0, 1, 1.0 - x) * pow(pi, 2)) / 6.0 +
404 pow(pi, 4) / 120.0 +
405 pow(hpl2(0, 1, 1.0 - x), 2)) -
406 2.0 * (hpl1(1, 1.0 - x) * hpl3(0, 0, 1, 1.0 - x) +
407 (2.0 * hpl4(0, 1, 0, 1, 1.0 - x) - pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0 +
408 (-2.0 * hpl4(0, 1, 0, 1, 1.0 - x) + pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0) -
409 2.0 * (hpl1(0, 1.0 - x) * hpl3(0, 1, 1, 1.0 - x) +
410 (2.0 * hpl4(0, 1, 0, 1, 1.0 - x) - pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0 +
411 (-2.0 * hpl4(0, 1, 0, 1, 1.0 - x) + pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0) -
412 pow(hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) - hpl2(0, 1, 1.0 - x) +
413 pow(pi, 2) / 6.0, 2)) / 2.0;
414 }
415 else {
416 hpl_base1(0,x)*hpl_base3(0,1,1,x) - hpl_base4(0,1,0,1,x) +
417 (2.*hpl_base4(0,1,0,1,x) - pow(hpl_base2(0,1,x),2))/2.;
418 }
419 }
420
421 if (i1 == 0 && i2 == 1 && i3 == 1 && i4 == 1) {
422 if (abs(x) > 1) {
423 return pi * scalar_t(0, -1) * hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) +
424 pi * scalar_t(0, 1) * hpl3(0, 0, 1, 1.0 / x) +
425 hpl1(0, 1.0 / x) * hpl3(0, 0, 1, 1.0 / x) +
426 pi * scalar_t(0, -1) * hpl3(0, 1, 1, 1.0 / x) -
427 hpl1(0, 1.0 / x) * hpl3(0, 1, 1, 1.0 / x) -
428 hpl4(0, 0, 0, 1, 1.0 / x) -
429 hpl4(0, 1, 0, 1, 1.0 / x) / 2.0 -
430 hpl4(0, 1, 1, 1, 1.0 / x) +
431 (hpl2(0, 1, 1.0 / x) * pow(pi, 2)) / 2.0 +
432 scalar_t(0, 0.16666666666666666) * hpl1(0, 1.0 / x) * pow(pi, 3) -
433 (19 * pow(pi, 4)) / 360.0 -
434 (hpl2(0, 1, 1.0 / x) * pow(hpl1(0, 1.0 / x), 2)) / 2.0 +
435 (pow(pi, 2) * pow(hpl1(0, 1.0 / x), 2)) / 4.0 +
436 pi * scalar_t(0, -0.16666666666666666) * pow(hpl1(0, 1.0 / x), 3) -
437 pow(hpl1(0, 1.0 / x), 4) / 24.0 +
438 (2.0 * hpl4(0, 1, 0, 1, 1.0 / x) - pow(hpl2(0, 1, 1.0 / x), 2)) / 2.0 +
439 pow(hpl2(0, 1, 1.0 / x), 2) / 4.0 +
440 (-2.0 * hpl4(0, 1, 0, 1, 1.0 / x) + pow(hpl2(0, 1, 1.0 / x), 2)) / 2.0;
441 } else if (real(x) > 0.5) {
442 return hpl1(0, 1.0 - x) * hpl3(0, 0, 1, 1.0 - x) -
443 hpl4(0, 0, 0, 1, 1.0 - x) +
444 pow(pi, 4) / 90.0 -
445 (hpl2(0, 1, 1.0 - x) * pow(hpl1(0, 1.0 - x), 2)) / 2.0 +
446 (hpl1(1, 1.0 - x) * pow(hpl1(0, 1.0 - x), 3)) / 6.0;
447 } else {
448 return hpl_base4(0, 1, 1, 1, x);
449 }
450 }
451
452 if (i1 == 1 && i2 == 0 && i3 == 0 && i4 == 0) {
453 if (abs(x) > 1) {
454 return (-3. * (pi * scalar_t(0, -1) * hpl1(0, 1.0 / x) - hpl2(0, 1, 1.0 / x) +
455 pow(pi, 2) / 3.0 - pow(hpl1(0, 1.0 / x), 2) / 2.0) * pow(hpl1(0, 1.0 / x), 2) -
456 6. * hpl1(0, 1.0 / x) * (hpl3(0, 0, 1, 1.0 / x) -
457 (hpl1(0, 1.0 / x) * pow(pi, 2)) / 3.0 +
458 pi * scalar_t(0, 0.5) * pow(hpl1(0, 1.0 / x), 2) +
459 pow(hpl1(0, 1.0 / x), 3) / 6.0) -
460 (pi * scalar_t(0, 1) + hpl1(0, 1.0 / x) + hpl1(1, 1.0 / x)) * pow(hpl1(0, 1.0 / x), 3) -
461 6. * (-hpl4(0, 0, 0, 1, 1.0 / x) + pow(pi, 4) / 45.0 +
462 (pow(pi, 2) * pow(hpl1(0, 1.0 / x), 2)) / 6.0 +
463 pi * scalar_t(0, -0.16666666666666666) * pow(hpl1(0, 1.0 / x), 3) -
464 pow(hpl1(0, 1.0 / x), 4) / 24.0)) / 6.0;
465 } else if (real(x) > 0.5) {
466 return (-3. * (hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) - hpl2(0, 1, 1.0 - x) +
467 pow(pi, 2) / 6.0) * pow(hpl1(1, 1.0 - x), 2) -
468 6. * hpl1(1, 1.0 - x) * (1.2020569031595942 +
469 hpl1(1, 1.0 - x) * hpl2(0, 1, 1.0 - x) - hpl3(0, 1, 1, 1.0 - x) -
470 (hpl1(1, 1.0 - x) * pow(pi, 2)) / 6.0 -
471 (hpl1(0, 1.0 - x) * pow(hpl1(1, 1.0 - x), 2)) / 2.0) +
472 hpl1(0, 1.0 - x) * pow(hpl1(1, 1.0 - x), 3) -
473 6. * (-1.2020569031595942 * hpl1(1, 1.0 - x) +
474 hpl1(1, 1.0 - x) * hpl3(0, 1, 1, 1.0 - x) -
475 hpl4(0, 1, 1, 1, 1.0 - x) + pow(pi, 4) / 90.0 -
476 (hpl2(0, 1, 1.0 - x) * pow(hpl1(1, 1.0 - x), 2)) / 2.0 +
477 (pow(pi, 2) * pow(hpl1(1, 1.0 - x), 2)) / 12.0 +
478 (hpl1(0, 1.0 - x) * pow(hpl1(1, 1.0 - x), 3)) / 6.0)) / 6.0;
479 } else {
480 return (6. * hpl_base1(0, x) * hpl_base3(0, 0, 1, x) - 6. * hpl_base4(0, 0, 0, 1, x) -
481 3. * hpl_base2(0, 1, x) * pow(hpl_base1(0, x), 2) + hpl_base1(1, x) * pow(hpl_base1(0, x), 3)) / 6.0;
482 }
483 }
484
485 if (i1 == 1 && i2 == 0 && i3 == 0 && i4 == 1) {
486 if (abs(x) > 1) {
487 return 3.7763731361630786 * scalar_t(0, -2) -
488 2.4041138063191885 * hpl1(0, 1.0 / x) +
489 pi * scalar_t(0, -1) * hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) +
490 pi * scalar_t(0, 2) * hpl3(0, 0, 1, 1.0 / x) +
491 2.0 * hpl1(0, 1.0 / x) * hpl3(0, 0, 1, 1.0 / x) -
492 4. * hpl4(0, 0, 0, 1, 1.0 / x) -
493 hpl4(0, 1, 0, 1, 1.0 / x) +
494 (hpl2(0, 1, 1.0 / x) * pow(pi, 2)) / 3.0 -
495 pow(pi, 4) / 90.0 -
496 (hpl2(0, 1, 1.0 / x) * pow(hpl1(0, 1.0 / x), 2)) / 2.0 +
497 (pow(pi, 2) * pow(hpl1(0, 1.0 / x), 2)) / 6.0 +
498 (pi * scalar_t(0, 1) + hpl1(0, 1.0 / x) + hpl1(1, 1.0 / x)) *
499 (hpl3(0, 0, 1, 1.0 / x) - (hpl1(0, 1.0 / x) * pow(pi, 2)) / 3.0 +
500 pi * scalar_t(0, 0.5) * pow(hpl1(0, 1.0 / x), 2) +
501 pow(hpl1(0, 1.0 / x), 3) / 6.0) +
502 pi * scalar_t(0, -0.16666666666666666) * pow(hpl1(0, 1.0 / x), 3) -
503 pow(hpl1(0, 1.0 / x), 4) / 24.0 +
504 (2.0 * (3.7763731361630786 * scalar_t(0, 2) +
505 2.4041138063191885 * hpl1(0, 1.0 / x) +
506 pi * scalar_t(0, 1) * hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) +
507 pi * scalar_t(0, -2) * hpl3(0, 0, 1, 1.0 / x) -
508 2.0 * hpl1(0, 1.0 / x) * hpl3(0, 0, 1, 1.0 / x) +
509 4. * hpl4(0, 0, 0, 1, 1.0 / x) + hpl4(0, 1, 0, 1, 1.0 / x) -
510 (hpl2(0, 1, 1.0 / x) * pow(pi, 2)) / 3.0 + pow(pi, 4) / 90.0 +
511 (hpl2(0, 1, 1.0 / x) * pow(hpl1(0, 1.0 / x), 2)) / 2.0 -
512 (pow(pi, 2) * pow(hpl1(0, 1.0 / x), 2)) / 6.0 +
513 pi * scalar_t(0, 0.16666666666666666) * pow(hpl1(0, 1.0 / x), 3) +
514 pow(hpl1(0, 1.0 / x), 4) / 24.0) -
515 pow(pi * scalar_t(0, -1) * hpl1(0, 1.0 / x) - hpl2(0, 1, 1.0 / x) +
516 pow(pi, 2) / 3.0 - pow(hpl1(0, 1.0 / x), 2) / 2.0, 2)) / 2.0;
517 } else if (real(x) > 0.5) {
518 return -2.4041138063191885 * hpl1(1, 1.0 - x) -
519 hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) * hpl2(0, 1, 1.0 - x) +
520 hpl4(0, 1, 0, 1, 1.0 - x) -
521 (hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) * pow(pi, 2)) / 6.0 +
522 (hpl2(0, 1, 1.0 - x) * pow(pi, 2)) / 6.0 - pow(pi, 4) / 120.0 -
523 hpl1(0, 1.0 - x) * (1.2020569031595942 +
524 hpl1(1, 1.0 - x) * hpl2(0, 1, 1.0 - x) - hpl3(0, 1, 1, 1.0 - x) -
525 (hpl1(1, 1.0 - x) * pow(pi, 2)) / 6.0 -
526 (hpl1(0, 1.0 - x) * pow(hpl1(1, 1.0 - x), 2)) / 2.0) -
527 pow(hpl2(0, 1, 1.0 - x), 2) +
528 2.0 * (hpl1(1, 1.0 - x) * hpl3(0, 0, 1, 1.0 - x) +
529 (2.0 * hpl4(0, 1, 0, 1, 1.0 - x) - pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0 +
530 (-2.0 * hpl4(0, 1, 0, 1, 1.0 - x) + pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0) +
531 2.0 * (hpl1(0, 1.0 - x) * hpl3(0, 1, 1, 1.0 - x) +
532 (2.0 * hpl4(0, 1, 0, 1, 1.0 - x) - pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0 +
533 (-2.0 * hpl4(0, 1, 0, 1, 1.0 - x) + pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0) +
534 (2.0 * (2.4041138063191885 * hpl1(1, 1.0 - x) +
535 hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) * hpl2(0, 1, 1.0 - x) -
536 hpl4(0, 1, 0, 1, 1.0 - x) +
537 (hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) * pow(pi, 2)) / 6.0 -
538 (hpl2(0, 1, 1.0 - x) * pow(pi, 2)) / 6.0 + pow(pi, 4) / 120.0 +
539 pow(hpl2(0, 1, 1.0 - x), 2)) -
540 2.0 * (hpl1(1, 1.0 - x) * hpl3(0, 0, 1, 1.0 - x) +
541 (2.0 * hpl4(0, 1, 0, 1, 1.0 - x) - pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0 +
542 (-2.0 * hpl4(0, 1, 0, 1, 1.0 - x) + pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0) -
543 pow(hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) - hpl2(0, 1, 1.0 - x) +
544 pow(pi, 2) / 6.0, 2)) / 2.0;
545 } else {
546 return hpl_base1(1, x) * hpl_base3(0, 0, 1, x) - hpl_base4(0, 1, 0, 1, x) +
547 (2.0 * hpl_base4(0, 1, 0, 1, x) - pow(hpl_base2(0, 1, x), 2)) / 2.0;
548 }
549 }
550
551 if (i1 == 1 && i2 == 0 && i3 == 1 && i4 == 0) {
552 if (abs(x) > 1) {
553 return 3.7763731361630786 * scalar_t(0, -2) -
554 2.4041138063191885 * hpl1(0, 1.0 / x) +
555 pi * scalar_t(0, -1) * hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) +
556 pi * scalar_t(0, 2) * hpl3(0, 0, 1, 1.0 / x) +
557 2.0 * hpl1(0, 1.0 / x) * hpl3(0, 0, 1, 1.0 / x) -
558 4. * hpl4(0, 0, 0, 1, 1.0 / x) -
559 hpl4(0, 1, 0, 1, 1.0 / x) +
560 (hpl2(0, 1, 1.0 / x) * pow(pi, 2)) / 3.0 -
561 pow(pi, 4) / 90.0 - hpl1(0, 1.0 / x) *
562 (pi * scalar_t(0, 1) + hpl1(0, 1.0 / x) + hpl1(1, 1.0 / x)) *
563 (pi * scalar_t(0, -1) * hpl1(0, 1.0 / x) - hpl2(0, 1, 1.0 / x) +
564 pow(pi, 2) / 3.0 - pow(hpl1(0, 1.0 / x), 2) / 2.0) -
565 (hpl2(0, 1, 1.0 / x) * pow(hpl1(0, 1.0 / x), 2)) / 2.0 +
566 (pow(pi, 2) * pow(hpl1(0, 1.0 / x), 2)) / 6.0 +
567 2.0 * hpl1(0, 1.0 / x) * (1.2020569031595942 +
568 pi * scalar_t(0, -1) * hpl2(0, 1, 1.0 / x) -
569 hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) + hpl3(0, 0, 1, 1.0 / x) -
570 hpl3(0, 1, 1, 1.0 / x) + (hpl1(0, 1.0 / x) * pow(pi, 2)) / 2.0 +
571 scalar_t(0, 0.16666666666666666) * pow(pi, 3) +
572 pi * scalar_t(0, -0.5) * pow(hpl1(0, 1.0 / x), 2) -
573 pow(hpl1(0, 1.0 / x), 3) / 6.0) -
574 2.0 * (pi * scalar_t(0, 1) + hpl1(0, 1.0 / x) + hpl1(1, 1.0 / x)) *
575 (hpl3(0, 0, 1, 1.0 / x) - (hpl1(0, 1.0 / x) * pow(pi, 2)) / 3.0 +
576 pi * scalar_t(0, 0.5) * pow(hpl1(0, 1.0 / x), 2) +
577 pow(hpl1(0, 1.0 / x), 3) / 6.0) +
578 pi * scalar_t(0, -0.16666666666666666) * pow(hpl1(0, 1.0 / x), 3) -
579 pow(hpl1(0, 1.0 / x), 4) / 24.0 +
580 pow(pi * scalar_t(0, -1) * hpl1(0, 1.0 / x) - hpl2(0, 1, 1.0 / x) +
581 pow(pi, 2) / 3.0 - pow(hpl1(0, 1.0 / x), 2) / 2.0, 2);
582 } else if (real(x) > 0.5) {
583 return -2.4041138063191885 * hpl1(1, 1.0 - x) -
584 hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) * hpl2(0, 1, 1.0 - x) +
585 hpl4(0, 1, 0, 1, 1.0 - x) +
586 hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) *
587 (hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) - hpl2(0, 1, 1.0 - x) +
588 pow(pi, 2) / 6.0) - (hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) *
589 pow(pi, 2)) / 6.0 + (hpl2(0, 1, 1.0 - x) * pow(pi, 2)) / 6.0 -
590 pow(pi, 4) / 120.0 + 2. * hpl1(1, 1.0 - x) *
591 (1.2020569031595942 + hpl1(0, 1.0 - x) * hpl2(0, 1, 1.0 - x) -
592 hpl3(0, 0, 1, 1.0 - x) -
593 (hpl1(1, 1.0 - x) * pow(hpl1(0, 1.0 - x), 2)) / 2.0) +
594 2.0 * hpl1(0, 1.0 - x) * (1.2020569031595942 +
595 hpl1(1, 1.0 - x) * hpl2(0, 1, 1.0 - x) - hpl3(0, 1, 1, 1.0 - x) -
596 (hpl1(1, 1.0 - x) * pow(pi, 2)) / 6.0 -
597 (hpl1(0, 1.0 - x) * pow(hpl1(1, 1.0 - x), 2)) / 2.0) -
598 pow(hpl2(0, 1, 1.0 - x), 2) +
599 2.0 * (hpl1(1, 1.0 - x) * hpl3(0, 0, 1, 1.0 - x) +
600 (2.0 * hpl4(0, 1, 0, 1, 1.0 - x) - pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0 +
601 (-2.0 * hpl4(0, 1, 0, 1, 1.0 - x) + pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0) +
602 2.0 * (hpl1(0, 1.0 - x) * hpl3(0, 1, 1, 1.0 - x) +
603 (2.0 * hpl4(0, 1, 0, 1, 1.0 - x) - pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0 +
604 (-2.0 * hpl4(0, 1, 0, 1, 1.0 - x) + pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0) +
605 pow(hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) - hpl2(0, 1, 1.0 - x) +
606 pow(pi, 2) / 6.0, 2);
607 } else {
608 return hpl_base1(0, x) * hpl_base1(1, x) * hpl_base2(0, 1, x) -
609 2. * hpl_base1(1, x) * hpl_base3(0, 0, 1, x) -
610 2.0 * hpl_base1(0, x) * hpl_base3(0, 1, 1, x) -
611 hpl_base4(0, 1, 0, 1, x) +
612 pow(hpl_base2(0, 1, x), 2);
613 }
614 }
615
616 if (i1 == 1 && i2 == 0 && i3 == 1 && i4 == 1) {
617 if (abs(x) > 1) {
618 return (pi * scalar_t(0, 1) + hpl1(0, 1.0 / x) + hpl1(1, 1.0 / x)) *
619 (1.2020569031595942 + pi * scalar_t(0, -1) * hpl2(0, 1, 1.0 / x) -
620 hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) + hpl3(0, 0, 1, 1.0 / x) -
621 hpl3(0, 1, 1, 1.0 / x) + (hpl1(0, 1.0 / x) * pow(pi, 2)) / 2.0 +
622 scalar_t(0, 0.16666666666666666) * pow(pi, 3) +
623 pi * scalar_t(0, -0.5) * pow(hpl1(0, 1.0 / x), 2) -
624 pow(hpl1(0, 1.0 / x), 3) / 6.0) -
625 3. * (pi * scalar_t(0, -1) * hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) +
626 pi * scalar_t(0, 1) * hpl3(0, 0, 1, 1.0 / x) +
627 hpl1(0, 1.0 / x) * hpl3(0, 0, 1, 1.0 / x) +
628 pi * scalar_t(0, -1) * hpl3(0, 1, 1, 1.0 / x) -
629 hpl1(0, 1.0 / x) * hpl3(0, 1, 1, 1.0 / x) - hpl4(0, 0, 0, 1, 1.0 / x) -
630 hpl4(0, 1, 0, 1, 1.0 / x) / 2.0 - hpl4(0, 1, 1, 1, 1.0 / x) +
631 (hpl2(0, 1, 1.0 / x) * pow(pi, 2)) / 2.0 +
632 scalar_t(0, 0.16666666666666666) * hpl1(0, 1.0 / x) * pow(pi, 3) -
633 (19 * pow(pi, 4)) / 360.0 -
634 (hpl2(0, 1, 1.0 / x) * pow(hpl1(0, 1.0 / x), 2)) / 2.0 +
635 (pow(pi, 2) * pow(hpl1(0, 1.0 / x), 2)) / 4.0 +
636 pi * scalar_t(0, -0.16666666666666666) * pow(hpl1(0, 1.0 / x), 3) -
637 pow(hpl1(0, 1.0 / x), 4) / 24.0 +
638 (2.0 * hpl4(0, 1, 0, 1, 1.0 / x) - pow(hpl2(0, 1, 1.0 / x), 2)) / 2.0 +
639 pow(hpl2(0, 1, 1.0 / x), 2) / 4.0 +
640 (-2.0 * hpl4(0, 1, 0, 1, 1.0 / x) + pow(hpl2(0, 1, 1.0 / x), 2)) / 2.0);
641 } else if (real(x) > 0.5) {
642 return -(hpl1(0, 1.0 - x) * (1.2020569031595942 +
643 hpl1(0, 1.0 - x) * hpl2(0, 1, 1.0 - x) - hpl3(0, 0, 1, 1.0 - x) -
644 (hpl1(1, 1.0 - x) * pow(hpl1(0, 1.0 - x), 2)) / 2.0)) -
645 3. * (hpl1(0, 1.0 - x) * hpl3(0, 0, 1, 1.0 - x) -
646 hpl4(0, 0, 0, 1, 1.0 - x) + pow(pi, 4) / 90.0 -
647 (hpl2(0, 1, 1.0 - x) * pow(hpl1(0, 1.0 - x), 2)) / 2.0 +
648 (hpl1(1, 1.0 - x) * pow(hpl1(0, 1.0 - x), 3)) / 6.0);
649 } else {
650 return hpl_base1(1, x) * hpl_base3(0, 1, 1, x) - 3. * hpl_base4(0, 1, 1, 1, x);
651 }
652 }
653
654 if (i1 == 1 && i2 == 1 && i3 == 0 && i4 == 0) {
655 if (abs(x) > 1) {
656 return hpl1(0, 1.0 / x) * (pi * scalar_t(0, 1) + hpl1(0, 1.0 / x) + hpl1(1, 1.0 / x)) *
657 (pi * scalar_t(0, -1) * hpl1(0, 1.0 / x) - hpl2(0, 1, 1.0 / x) +
658 pow(pi, 2) / 3.0 - pow(hpl1(0, 1.0 / x), 2) / 2.0) -
659 hpl1(0, 1.0 / x) * (1.2020569031595942 +
660 pi * scalar_t(0, -1) * hpl2(0, 1, 1.0 / x) -
661 hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) + hpl3(0, 0, 1, 1.0 / x) -
662 hpl3(0, 1, 1, 1.0 / x) + (hpl1(0, 1.0 / x) * pow(pi, 2)) / 2.0 +
663 scalar_t(0, 0.16666666666666666) * pow(pi, 3) +
664 pi * scalar_t(0, -0.5) * pow(hpl1(0, 1.0 / x), 2) -
665 pow(hpl1(0, 1.0 / x), 3) / 6.0) +
666 (pi * scalar_t(0, 1) + hpl1(0, 1.0 / x) + hpl1(1, 1.0 / x)) *
667 (hpl3(0, 0, 1, 1.0 / x) - (hpl1(0, 1.0 / x) * pow(pi, 2)) / 3.0 +
668 pi * scalar_t(0, 0.5) * pow(hpl1(0, 1.0 / x), 2) +
669 pow(hpl1(0, 1.0 / x), 3) / 6.0) +
670 (pow(hpl1(0, 1.0 / x), 2) * pow(pi * scalar_t(0, 1) + hpl1(0, 1.0 / x) +
671 hpl1(1, 1.0 / x), 2)) / 4.0 +
672 (2.0 * (3.7763731361630786 * scalar_t(0, 2) +
673 2.4041138063191885 * hpl1(0, 1.0 / x) +
674 pi * scalar_t(0, 1) * hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) +
675 pi * scalar_t(0, -2) * hpl3(0, 0, 1, 1.0 / x) -
676 2.0 * hpl1(0, 1.0 / x) * hpl3(0, 0, 1, 1.0 / x) +
677 4. * hpl4(0, 0, 0, 1, 1.0 / x) + hpl4(0, 1, 0, 1, 1.0 / x) -
678 (hpl2(0, 1, 1.0 / x) * pow(pi, 2)) / 3.0 + pow(pi, 4) / 90.0 +
679 (hpl2(0, 1, 1.0 / x) * pow(hpl1(0, 1.0 / x), 2)) / 2.0 -
680 (pow(pi, 2) * pow(hpl1(0, 1.0 / x), 2)) / 6.0 +
681 pi * scalar_t(0, 0.16666666666666666) * pow(hpl1(0, 1.0 / x), 3) +
682 pow(hpl1(0, 1.0 / x), 4) / 24.0) -
683 pow(pi * scalar_t(0, -1) * hpl1(0, 1.0 / x) - hpl2(0, 1, 1.0 / x) +
684 pow(pi, 2) / 3.0 - pow(hpl1(0, 1.0 / x), 2) / 2.0, 2)) / 4.0;
685 } else if (real(x) > 0.5) {
686 return -(hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) *
687 (hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) - hpl2(0, 1, 1.0 - x) +
688 pow(pi, 2) / 6.0)) -
689 hpl1(1, 1.0 - x) * (1.2020569031595942 +
690 hpl1(0, 1.0 - x) * hpl2(0, 1, 1.0 - x) - hpl3(0, 0, 1, 1.0 - x) -
691 (hpl1(1, 1.0 - x) * pow(hpl1(0, 1.0 - x), 2)) / 2.0) +
692 (pow(hpl1(0, 1.0 - x), 2) * pow(hpl1(1, 1.0 - x), 2)) / 4.0 -
693 hpl1(0, 1.0 - x) * (1.2020569031595942 +
694 hpl1(1, 1.0 - x) * hpl2(0, 1, 1.0 - x) - hpl3(0, 1, 1, 1.0 - x) -
695 (hpl1(1, 1.0 - x) * pow(pi, 2)) / 6.0 -
696 (hpl1(0, 1.0 - x) * pow(hpl1(1, 1.0 - x), 2)) / 2.0) +
697 (2.0 * (2.4041138063191885 * hpl1(1, 1.0 - x) +
698 hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) * hpl2(0, 1, 1.0 - x) -
699 hpl4(0, 1, 0, 1, 1.0 - x) +
700 (hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) * pow(pi, 2)) / 6.0 -
701 (hpl2(0, 1, 1.0 - x) * pow(pi, 2)) / 6.0 + pow(pi, 4) / 120.0 +
702 pow(hpl2(0, 1, 1.0 - x), 2) -
703 2.0 * (hpl1(1, 1.0 - x) * hpl3(0, 0, 1, 1.0 - x) +
704 (2.0 * hpl4(0, 1, 0, 1, 1.0 - x) -
705 pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0 +
706 (-2.0 * hpl4(0, 1, 0, 1, 1.0 - x) +
707 pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0) -
708 2.0 * (hpl1(0, 1.0 - x) * hpl3(0, 1, 1, 1.0 - x) +
709 (2.0 * hpl4(0, 1, 0, 1, 1.0 - x) -
710 pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0 +
711 (-2.0 * hpl4(0, 1, 0, 1, 1.0 - x) +
712 pow(hpl2(0, 1, 1.0 - x), 2)) / 2.0)) -
713 pow(hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) - hpl2(0, 1, 1.0 - x) +
714 pow(pi, 2) / 6.0, 2)) / 4.0;
715 } else {
716 return -(hpl_base1(0, x) * hpl_base1(1, x) * hpl_base2(0, 1, x)) + hpl_base1(1, x) * hpl_base3(0, 0, 1, x) +
717 hpl_base1(0, x) * hpl_base3(0, 1, 1, x) +
718 (pow(hpl_base1(0, x), 2) * pow(hpl_base1(1, x), 2)) / 4.0 +
719 (2.0 * hpl_base4(0, 1, 0, 1, x) - pow(hpl_base2(0, 1, x), 2)) / 4.0;
720 }
721 }
722
723 if (i1 == 1 && i2 == 1 && i3 == 0 && i4 == 1) {
724 if (abs(x) > 1) {
725 return (-4. * (pi * scalar_t(0, 1) + hpl1(0, 1.0 / x) + hpl1(1, 1.0 / x)) *
726 (1.2020569031595942 + pi * scalar_t(0, -1) * hpl2(0, 1, 1.0 / x) -
727 hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) + hpl3(0, 0, 1, 1.0 / x) -
728 hpl3(0, 1, 1, 1.0 / x) + (hpl1(0, 1.0 / x) * pow(pi, 2)) / 2.0 +
729 scalar_t(0, 0.16666666666666666) * pow(pi, 3) +
730 pi * scalar_t(0, -0.5) * pow(hpl1(0, 1.0 / x), 2) -
731 pow(hpl1(0, 1.0 / x), 3) / 6.0) +
732 (pi * scalar_t(0, -1) * hpl1(0, 1.0 / x) - hpl2(0, 1, 1.0 / x) +
733 pow(pi, 2) / 3.0 - pow(hpl1(0, 1.0 / x), 2) / 2.0) *
734 pow(pi * scalar_t(0, 1) + hpl1(0, 1.0 / x) + hpl1(1, 1.0 / x), 2) +
735 6. * (pi * scalar_t(0, -1) * hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) +
736 pi * scalar_t(0, 1) * hpl3(0, 0, 1, 1.0 / x) +
737 hpl1(0, 1.0 / x) * hpl3(0, 0, 1, 1.0 / x) +
738 pi * scalar_t(0, -1) * hpl3(0, 1, 1, 1.0 / x) -
739 hpl1(0, 1.0 / x) * hpl3(0, 1, 1, 1.0 / x) - hpl4(0, 0, 0, 1, 1.0 / x) -
740 hpl4(0, 1, 0, 1, 1.0 / x) / 2.0 - hpl4(0, 1, 1, 1, 1.0 / x) +
741 (hpl2(0, 1, 1.0 / x) * pow(pi, 2)) / 2.0 +
742 scalar_t(0, 0.16666666666666666) * hpl1(0, 1.0 / x) * pow(pi, 3) -
743 (19 * pow(pi, 4)) / 360.0 -
744 (hpl2(0, 1, 1.0 / x) * pow(hpl1(0, 1.0 / x), 2)) / 2.0 +
745 (pow(pi, 2) * pow(hpl1(0, 1.0 / x), 2)) / 4.0 +
746 pi * scalar_t(0, -0.16666666666666666) * pow(hpl1(0, 1.0 / x), 3) -
747 pow(hpl1(0, 1.0 / x), 4) / 24.0 +
748 (2.0 * hpl4(0, 1, 0, 1, 1.0 / x) - pow(hpl2(0, 1, 1.0 / x), 2)) / 2.0 +
749 pow(hpl2(0, 1, 1.0 / x), 2) / 4.0 +
750 (-2.0 * hpl4(0, 1, 0, 1, 1.0 / x) + pow(hpl2(0, 1, 1.0 / x), 2)) / 2.0)) / 2.0;
751 } else if (real(x) > 0.5) {
752 return ((hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) - hpl2(0, 1, 1.0 - x) +
753 pow(pi, 2) / 6.0) * pow(hpl1(0, 1.0 - x), 2) +
754 4. * hpl1(0, 1.0 - x) * (1.2020569031595942 +
755 hpl1(0, 1.0 - x) * hpl2(0, 1, 1.0 - x) - hpl3(0, 0, 1, 1.0 - x) -
756 (hpl1(1, 1.0 - x) * pow(hpl1(0, 1.0 - x), 2)) / 2.0) +
757 6. * (hpl1(0, 1.0 - x) * hpl3(0, 0, 1, 1.0 - x) -
758 hpl4(0, 0, 0, 1, 1.0 - x) + pow(pi, 4) / 90.0 -
759 (hpl2(0, 1, 1.0 - x) * pow(hpl1(0, 1.0 - x), 2)) / 2.0 +
760 (hpl1(1, 1.0 - x) * pow(hpl1(0, 1.0 - x), 3)) / 6.0)) / 2.0;
761 }
762 return (-4. * hpl_base1(1, x) * hpl_base3(0, 1, 1, x) + 6. * hpl_base4(0, 1, 1, 1, x) +
763 hpl_base2(0, 1, x) * pow(hpl_base1(1, x), 2)) / 2.0;
764 }
765
766
767 if (i1 == 1 && i2 == 1 && i3 == 1 && i4 == 0) {
768 if (abs(x) > 1) {
769 return (6. * (pi * scalar_t(0, 1) + hpl1(0, 1.0 / x) + hpl1(1, 1.0 / x)) *
770 (1.2020569031595942 + pi * scalar_t(0, -1) * hpl2(0, 1, 1.0 / x) -
771 hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) + hpl3(0, 0, 1, 1.0 / x) -
772 hpl3(0, 1, 1, 1.0 / x) + (hpl1(0, 1.0 / x) * pow(pi, 2)) / 2.0 +
773 scalar_t(0, 0.16666666666666666) * pow(pi, 3) +
774 pi * scalar_t(0, -0.5) * pow(hpl1(0, 1.0 / x), 2) -
775 pow(hpl1(0, 1.0 / x), 3) / 6.0) -
776 3. * (pi * scalar_t(0, -1) * hpl1(0, 1.0 / x) - hpl2(0, 1, 1.0 / x) +
777 pow(pi, 2) / 3.0 - pow(hpl1(0, 1.0 / x), 2) / 2.0) *
778 pow(pi * scalar_t(0, 1) + hpl1(0, 1.0 / x) + hpl1(1, 1.0 / x), 2) -
779 hpl1(0, 1.0 / x) * pow(pi * scalar_t(0, 1) + hpl1(0, 1.0 / x) + hpl1(1, 1.0 / x), 3) -
780 6. * (pi * scalar_t(0, -1) * hpl1(0, 1.0 / x) * hpl2(0, 1, 1.0 / x) +
781 pi * scalar_t(0, 1) * hpl3(0, 0, 1, 1.0 / x) +
782 hpl1(0, 1.0 / x) * hpl3(0, 0, 1, 1.0 / x) +
783 pi * scalar_t(0, -1) * hpl3(0, 1, 1, 1.0 / x) -
784 hpl1(0, 1.0 / x) * hpl3(0, 1, 1, 1.0 / x) - hpl4(0, 0, 0, 1, 1.0 / x) -
785 hpl4(0, 1, 0, 1, 1.0 / x) / 2.0 - hpl4(0, 1, 1, 1, 1.0 / x) +
786 (hpl2(0, 1, 1.0 / x) * pow(pi, 2)) / 2.0 +
787 scalar_t(0, 0.16666666666666666) * hpl1(0, 1.0 / x) * pow(pi, 3) -
788 (19 * pow(pi, 4)) / 360.0 -
789 (hpl2(0, 1, 1.0 / x) * pow(hpl1(0, 1.0 / x), 2)) / 2.0 +
790 (pow(pi, 2) * pow(hpl1(0, 1.0 / x), 2)) / 4.0 +
791 pi * scalar_t(0, -0.16666666666666666) * pow(hpl1(0, 1.0 / x), 3) -
792 pow(hpl1(0, 1.0 / x), 4) / 24.0 +
793 (2.0 * hpl4(0, 1, 0, 1, 1.0 / x) - pow(hpl2(0, 1, 1.0 / x), 2)) / 2.0 +
794 pow(hpl2(0, 1, 1.0 / x), 2) / 4.0 +
795 (-2.0 * hpl4(0, 1, 0, 1, 1.0 / x) + pow(hpl2(0, 1, 1.0 / x), 2)) / 2.0)) / 6.0;
796 } else if (real(x) > 0.5) {
797 return (-3. * (hpl1(0, 1.0 - x) * hpl1(1, 1.0 - x) - hpl2(0, 1, 1.0 - x) +
798 pow(pi, 2) / 6.0) * pow(hpl1(0, 1.0 - x), 2) -
799 6. * hpl1(0, 1.0 - x) * (1.2020569031595942 +
800 hpl1(0, 1.0 - x) * hpl2(0, 1, 1.0 - x) - hpl3(0, 0, 1, 1.0 - x) -
801 (hpl1(1, 1.0 - x) * pow(hpl1(0, 1.0 - x), 2)) / 2.0) +
802 hpl1(1, 1.0 - x) * pow(hpl1(0, 1.0 - x), 3) -
803 6. * (hpl1(0, 1.0 - x) * hpl3(0, 0, 1, 1.0 - x) -
804 hpl4(0, 0, 0, 1, 1.0 - x) + pow(pi, 4) / 90.0 -
805 (hpl2(0, 1, 1.0 - x) * pow(hpl1(0, 1.0 - x), 2)) / 2.0 +
806 (hpl1(1, 1.0 - x) * pow(hpl1(0, 1.0 - x), 3)) / 6.0)) / 6.0;
807 }
808 return (6. * hpl_base1(1, x) * hpl_base3(0, 1, 1, x) - 6. * hpl_base4(0, 1, 1, 1, x) -
809 3. * hpl_base2(0, 1, x) * pow(hpl_base1(1, x), 2) + hpl_base1(0, x) * pow(hpl_base1(1, x), 3)) / 6.0;
810 }
811
812
813 if(i1==1&&i2==1&&i3==1&&i4==1)
814 {
815 if(abs(x)>1) return pow(pi*cd(0.,1.) + hpl1(0,1./x) + hpl1(1,1./x),4)/24.;
816 if(real(x)>.5) return pow(hpl1(0,1. - x),4)/24.;
817 return pow(hpl_base1(1,x),4)/24.;
818 }
819
820
821 return {0.,0.};
822
823}
824
825
826template <size_t MAX_N, size_t MAX_M, size_t MAX_LE, size_t MAX_NC>
827void initialize_constants(std::array<std::array<double, MAX_M>, MAX_N> &s1, std::array<std::array<double, MAX_M>, MAX_N> &c, std::array<std::array<double, MAX_NC>, MAX_LE> &a) {
828 // Initialisation des constantes s1
829 s1[1][1] = 1.6449340668482;
830 s1[1][2] = 1.2020569031596;
831 s1[1][3] = 1.0823232337111;
832 s1[1][4] = 1.0369277551434;
833 s1[2][1] = 1.2020569031596;
834 s1[2][2] = 2.7058080842778e-1;
835 s1[2][3] = 9.6551159989444e-2;
836 s1[3][1] = 1.0823232337111;
837 s1[3][2] = 9.6551159989444e-2;
838 s1[4][1] = 1.0369277551434;
839
840 // Initialisation des constantes c
841 c[1][1] = 1.6449340668482;
842 c[1][2] = 1.2020569031596;
843 c[1][3] = 1.0823232337111;
844 c[1][4] = 1.0369277551434;
845 c[2][1] = 0.0;
846 c[2][2] = -1.8940656589945;
847 c[2][3] = -3.0142321054407;
848 c[3][1] = 1.8940656589945;
849 c[3][2] = 3.0142321054407;
850 c[4][1] = 0.0;
851
852
853
854 a[0][1]=.96753215043498;
855 a[1][1]=.16607303292785;
856 a[2][1]=.02487932292423;
857 a[3][1]=.00468636195945;
858 a[4][1]=.00100162749616;
859 a[5][1]=.00023200219609;
860 a[6][1]=.00005681782272;
861 a[7][1]=.00001449630056;
862 a[8][1]=.00000381632946;
863 a[9][1]=.00000102990426;
864 a[10][1]=.00000028357538;
865 a[11][1]=.00000007938705;
866 a[12][1]=.00000002253670;
867 a[13][1]=.00000000647434;
868 a[14][1]=.00000000187912;
869 a[15][1]=.00000000055029;
870 a[16][1]=.00000000016242;
871 a[17][1]=.00000000004827;
872 a[18][1]=.00000000001444;
873 a[19][1]=.00000000000434;
874 a[20][1]=.00000000000131;
875 a[21][1]=.00000000000040;
876 a[22][1]=.00000000000012;
877 a[23][1]=.00000000000004;
878 a[24][1]=.00000000000001;
879 a[0][2]=.95180889127832;
880 a[1][2]=.43131131846532;
881 a[2][2]=.10002250714905;
882 a[3][2]=.02442415595220;
883 a[4][2]=.00622512463724;
884 a[5][2]=.00164078831235;
885 a[6][2]=.00044407920265;
886 a[7][2]=.00012277494168;
887 a[8][2]=.00003453981284;
888 a[9][2]=.00000985869565;
889 a[10][2]=.00000284856995;
890 a[11][2]=.00000083170847;
891 a[12][2]=.00000024503950;
892 a[13][2]=.00000007276496;
893 a[14][2]=.00000002175802;
894 a[15][2]=.00000000654616;
895 a[16][2]=.00000000198033;
896 a[17][2]=.00000000060204;
897 a[18][2]=.00000000018385;
898 a[19][2]=.00000000005637;
899 a[20][2]=.00000000001735;
900 a[21][2]=.00000000000536;
901 a[22][2]=.00000000000166;
902 a[23][2]=.00000000000052;
903 a[24][2]=.00000000000016;
904 a[25][2]=.00000000000005;
905 a[26][2]=.00000000000002;
906 a[0][3]=.98161027991365;
907 a[1][3]=.72926806320726;
908 a[2][3]=.22774714909321;
909 a[3][3]=.06809083296197;
910 a[4][3]=.02013701183064;
911 a[5][3]=.00595478480197;
912 a[6][3]=.00176769013959;
913 a[7][3]=.00052748218502;
914 a[8][3]=.00015827461460;
915 a[9][3]=.00004774922076;
916 a[10][3]=.00001447920408;
917 a[11][3]=.00000441154886;
918 a[12][3]=.00000135003870;
919 a[13][3]=.00000041481779;
920 a[14][3]=.00000012793307;
921 a[15][3]=.00000003959070;
922 a[16][3]=.00000001229055;
923 a[17][3]=.00000000382658;
924 a[18][3]=.00000000119459;
925 a[19][3]=.00000000037386;
926 a[20][3]=.00000000011727;
927 a[21][3]=.00000000003687;
928 a[22][3]=.00000000001161;
929 a[23][3]=.00000000000366;
930 a[24][3]=.00000000000116;
931 a[25][3]=.00000000000037;
932 a[26][3]=.00000000000012;
933 a[27][3]=.00000000000004;
934 a[28][3]=.00000000000001;
935 a[0][4]=1.0640521184614;
936 a[1][4]=1.0691720744981;
937 a[2][4]=.41527193251768;
938 a[3][4]=.14610332936222;
939 a[4][4]=.04904732648784;
940 a[5][4]=.01606340860396;
941 a[6][4]=.00518889350790;
942 a[7][4]=.00166298717324;
943 a[8][4]=.00053058279969;
944 a[9][4]=.00016887029251;
945 a[10][4]=.00005368328059;
946 a[11][4]=.00001705923313;
947 a[12][4]=.00000542174374;
948 a[13][4]=.00000172394082;
949 a[14][4]=.00000054853275;
950 a[15][4]=.00000017467795;
951 a[16][4]=.00000005567550;
952 a[17][4]=.00000001776234;
953 a[18][4]=.00000000567224;
954 a[19][4]=.00000000181313;
955 a[20][4]=.00000000058012;
956 a[21][4]=.00000000018579;
957 a[22][4]=.00000000005955;
958 a[23][4]=.00000000001911;
959 a[24][4]=.00000000000614;
960 a[25][4]=.00000000000197;
961 a[26][4]=.00000000000063;
962 a[27][4]=.00000000000020;
963 a[28][4]=.00000000000007;
964 a[29][4]=.00000000000002;
965 a[30][4]=.00000000000001;
966 a[0][5]=.97920860669175;
967 a[1][5]=.08518813148683;
968 a[2][5]=.00855985222013;
969 a[3][5]=.00121177214413;
970 a[4][5]=.00020722768531;
971 a[5][5]=.00003996958691;
972 a[6][5]=.00000838064065;
973 a[7][5]=.00000186848945;
974 a[8][5]=.00000043666087;
975 a[9][5]=.00000010591733;
976 a[10][5]=.00000002647892;
977 a[11][5]=.00000000678700;
978 a[12][5]=.00000000177654;
979 a[13][5]=.00000000047342;
980 a[14][5]=.00000000012812;
981 a[15][5]=.00000000003514;
982 a[16][5]=.00000000000975;
983 a[17][5]=.00000000000274;
984 a[18][5]=.00000000000077;
985 a[19][5]=.00000000000022;
986 a[20][5]=.00000000000006;
987 a[21][5]=.00000000000002;
988 a[22][5]=.00000000000001;
989 a[0][6]=.95021851963952;
990 a[1][6]=.29052529161433;
991 a[2][6]=.05081774061716;
992 a[3][6]=.00995543767280;
993 a[4][6]=.00211733895031;
994 a[5][6]=.00047859470550;
995 a[6][6]=.00011334321308;
996 a[7][6]=.00002784733104;
997 a[8][6]=.00000704788108;
998 a[9][6]=.00000182788740;
999 a[10][6]=.00000048387492;
1000 a[11][6]=.00000013033842;
1001 a[12][6]=.00000003563769;
1002 a[13][6]=.00000000987174;
1003 a[14][6]=.00000000276586;
1004 a[15][6]=.00000000078279;
1005 a[16][6]=.00000000022354;
1006 a[17][6]=.00000000006435;
1007 a[18][6]=.00000000001866;
1008 a[19][6]=.00000000000545;
1009 a[20][6]=.00000000000160;
1010 a[21][6]=.00000000000047;
1011 a[22][6]=.00000000000014;
1012 a[23][6]=.00000000000004;
1013 a[24][6]=.00000000000001;
1014 a[0][7]=.95064032186777;
1015 a[1][7]=.54138285465171;
1016 a[2][7]=.13649979590321;
1017 a[3][7]=.03417942328207;
1018 a[4][7]=.00869027883583;
1019 a[5][7]=.00225284084155;
1020 a[6][7]=.00059516089806;
1021 a[7][7]=.00015995617766;
1022 a[8][7]=.00004365213096;
1023 a[9][7]=.00001207474688;
1024 a[10][7]=.00000338018176;
1025 a[11][7]=.00000095632476;
1026 a[12][7]=.00000027313129;
1027 a[13][7]=.00000007866968;
1028 a[14][7]=.00000002283195;
1029 a[15][7]=.00000000667205;
1030 a[16][7]=.00000000196191;
1031 a[17][7]=.00000000058018;
1032 a[18][7]=.00000000017246;
1033 a[19][7]=.00000000005151;
1034 a[20][7]=.00000000001545;
1035 a[21][7]=.00000000000465;
1036 a[22][7]=.00000000000141;
1037 a[23][7]=.00000000000043;
1038 a[24][7]=.00000000000013;
1039 a[25][7]=.00000000000004;
1040 a[26][7]=.00000000000001;
1041 a[0][8]=.98800011672229;
1042 a[1][8]=.04364067609601;
1043 a[2][8]=.00295091178278;
1044 a[3][8]=.00031477809720;
1045 a[4][8]=.00004314846029;
1046 a[5][8]=.00000693818230;
1047 a[6][8]=.00000124640350;
1048 a[7][8]=.00000024293628;
1049 a[8][8]=.00000005040827;
1050 a[9][8]=.00000001099075;
1051 a[10][8]=.00000000249467;
1052 a[11][8]=.00000000058540;
1053 a[12][8]=.00000000014127;
1054 a[13][8]=.00000000003492;
1055 a[14][8]=.00000000000881;
1056 a[15][8]=.00000000000226;
1057 a[16][8]=.00000000000059;
1058 a[17][8]=.00000000000016;
1059 a[18][8]=.00000000000004;
1060 a[19][8]=.00000000000001;
1061 a[0][9]=.95768506546350;
1062 a[1][9]=.19725249679534;
1063 a[2][9]=.02603370313918;
1064 a[3][9]=.00409382168261;
1065 a[4][9]=.00072681707110;
1066 a[5][9]=.00014091879261;
1067 a[6][9]=.00002920458914;
1068 a[7][9]=.00000637631144;
1069 a[8][9]=.00000145167850;
1070 a[9][9]=.00000034205281;
1071 a[10][9]=.00000008294302;
1072 a[11][9]=.00000002060784;
1073 a[12][9]=.00000000522823;
1074 a[13][9]=.00000000135066;
1075 a[14][9]=.00000000035451;
1076 a[15][9]=.00000000009436;
1077 a[16][9]=.00000000002543;
1078 a[17][9]=.00000000000693;
1079 a[18][9]=.00000000000191;
1080 a[19][9]=.00000000000053;
1081 a[20][9]=.00000000000015;
1082 a[21][9]=.00000000000004;
1083 a[22][9]=.00000000000001;
1084 a[0][10]=.99343651671347;
1085 a[1][10]=.02225770126826;
1086 a[2][10]=.00101475574703;
1087 a[3][10]=.00008175156250;
1088 a[4][10]=.00000899973547;
1089 a[5][10]=.00000120823987;
1090 a[6][10]=.00000018616913;
1091 a[7][10]=.00000003174723;
1092 a[8][10]=.00000000585215;
1093 a[9][10]=.00000000114739;
1094 a[10][10]=.00000000023652;
1095 a[11][10]=.00000000005082;
1096 a[12][10]=.00000000001131;
1097 a[13][10]=.00000000000259;
1098 a[14][10]=.00000000000061;
1099 a[15][10]=.00000000000015;
1100 a[16][10]=.00000000000004;
1101 a[17][10]=.00000000000001;
1102}
1103
1104
1105
1106
1107
1108
1109scalar_t polylog(size_t n, size_t m, double x) {
1110
1111 const int MAX_N = 5;
1112 const int MAX_M = 5;
1113 const int MAX_LE = 31;
1114 const int MAX_NC = 11;
1115
1116 std::array<std::array<double, MAX_M>, MAX_N> s1, c;
1117 std::array<std::array<double, MAX_NC>, MAX_LE> a;
1118
1119 initialize_constants<MAX_N, MAX_M, MAX_LE, MAX_NC>(s1, c, a);
1120
1121
1122 // double u[5],s1[5][5],c[5][5],a[31][11];
1123 double u[5];
1124 const int fct[5] = {1, 1, 2, 6, 24};
1125 const int sgn[5] = {1, -1, 1, -1, 1};
1126 const size_t index[32] = {0, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 5, 6, 7, 0, 0, 0, 0, 0, 0, 0, 8, 9, 0, 0, 0, 0, 0, 0, 0, 0, 10};
1127 const size_t nc[11] = {0, 24, 26, 28, 30, 22, 24, 26, 19, 22, 17};
1128
1129 scalar_t z = 0.0;
1130 scalar_t v[6] = {0.0};
1131 scalar_t sk;
1132 scalar_t sj;
1133
1134 double z1=1.;
1135 double hf=z1/2.;
1136 double b0{0.};
1137
1138 if ((n < 1) || (n > 4) || (m < 1) || (m > 4) || ((n + m) > 5)) {
1139 return 0.0;
1140 }
1141
1142 if (x == 1.0) {
1143 z = s1[n][m];
1144 } else if ((x > 2.0) || (x < -1.0)) {
1145 // Branch for x > 2 or x < -1
1146 double x1 = 1.0 / x;
1147 double h = 4.0 / 3.0 * x1 + 1.0 / 3.0;
1148 double alfa = h + h;
1149 v[0] = 1.0;
1150 v[1] = log(-x);
1151
1152 for (size_t le = 2; le <= n + m; le++) {
1153 v[le] = v[1] * v[le - 1] / (le*1.);
1154 }
1155
1156
1157 for (size_t ke = 0; ke <= m - 1; ke++) {
1158 size_t m1 = m - ke;
1159 double r = pow(x1, m1) / ((fct[m1] * fct[n - 1])*1.);
1160 sj = 0.0;
1161
1162 for (size_t je = 0; je <= ke; je++) {
1163 size_t n1 = n + ke - je;
1164 size_t le = index[10 * n1 + m1 - 10];
1165 double b1 = 0.0;
1166 double b2 = 0.0;
1167
1168 auto nc_le = nc[le]; // size_t
1169 for (std::ptrdiff_t it = static_cast<std::ptrdiff_t>(nc_le); it >= 0; --it) {
1170 b0 = a[le][static_cast<std::size_t>(it)] + alfa * b1 - b2;
1171 b2 = b1;
1172 b1 = b0;
1173 }
1174
1175 double q = (fct[n1 - 1] / fct[ke - je]) * (b0 - h * b2) * r / pow(m1, n1);
1176 sj += v[je] * q;
1177 }
1178
1179 sk += sgn[ke]*1. * sj;
1180 }
1181
1182
1183
1184 for (size_t je = 0; je <= n - 1; je++) {
1185 sj += v[je] * c[n - je][m];
1186 }
1187
1188 z = sgn[n] *1.* sk + sgn[m]*1. * (sj + v[n + m]);
1189 } else if (x>hf) {
1190 // Branch for x > 0.5
1191 double x1 = 1.0 - x;
1192 double h = 4.0 / 3.0 * x1 + 1.0 / 3.0;
1193 double alfa = h + h;
1194 v[0] = 1.0;
1195 u[0] = 1.0;
1196 v[1] = log(x1);
1197 u[1] = log(x);
1198
1199 for (size_t le = 2; le <= m; le++) {
1200 v[le] = v[1] * v[le - 1] / (1.*le);
1201 }
1202
1203 for (size_t le = 2; le <= n; le++) {
1204 u[le] = u[1] * u[le - 1] / le;
1205 }
1206
1207 sk = 0.0;
1208
1209 for (size_t ke = 0; ke <= n - 1; ke++) {
1210 size_t m1 = n - ke;
1211 double r = pow(x1, m1) / fct[m1];
1212
1213 for (size_t je = 0; je <= m - 1; je++) {
1214 size_t n1 = m - je;
1215 size_t le = index[10 * n1 + m1 - 10];
1216 double b1 = 0.0;
1217 double b2 = 0.0;
1218
1219 auto nc_le = nc[le]; // size_t
1220 for (std::ptrdiff_t it = static_cast<std::ptrdiff_t>(nc_le); it >= 0; --it) {
1221 b0 = a[le][static_cast<std::size_t>(it)] + alfa * b1 - b2;
1222 b2 = b1;
1223 b1 = b0;
1224 }
1225
1226 double q = sgn[je] * (b0 - h * b2) * r / pow(m1, n1);
1227 sj += v[je] * q;
1228 }
1229
1230 sk = sk + u[ke] * (s1[m1][m] - sj);
1231 }
1232
1233 z = sk + sgn[m] * u[n] * v[m];
1234 } else {
1235 // Branch for the remaining cases
1236 size_t le = index[10 * n + m - 10];
1237 double h = 4.0 / 3.0 * x + 1.0 / 3.0;
1238 double alfa = h + h;
1239 double b1 = 0.0;
1240 double b2 = 0.0;
1241
1242 auto nc_le = nc[le]; // size_t
1243 for (std::ptrdiff_t it = static_cast<std::ptrdiff_t>(nc_le); it >= 0; --it) {
1244 b0 = a[le][static_cast<std::size_t>(it)] + alfa * b1 - b2;
1245 b2 = b1;
1246 b1 = b0;
1247 }
1248
1249 z = (b0 - h * b2) * pow(x, m) / (fct[m] * pow(m, n));
1250 }
1251
1252 return z;
1253}
1254
1255// LI2 ! -------------------------------------------------------------------------------------------------------------------
1256
1257double Li2(double x) {
1258 return gsl_sf_dilog(x);
1259}
1260
1261
1262// LI3 ! -------------------------------------------------------------------------------------------------------------------
1263
1264
1265const double zeta3 = 1.202056903159594; // Apery's constant
1266
1267double Li3(double x) {
1268 const double pisq6 = M_PI * M_PI / 6.0;
1269 const double x_0 = -1.0;
1270 const double x_1 = -0.85;
1271 const double x_2 = 0.25;
1272 const double x_3 = 0.63;
1273 const double x_4 = 1.0;
1274
1275 if (fpeq(x, 1.0)) return zeta3;
1276 if (fpeq(x, -1.0)) return -0.75 * zeta3;
1277
1278 if (x <= x_0) {
1279 double lnx = log(-x);
1280 return Li3(1.0 / x) - pisq6 * lnx - lnx * lnx * lnx / 6.0;
1281 } else if (x < x_1) {
1282 return Li3(x * x) / 4.0 - Li3(-x);
1283 } else if (x < x_2) {
1284 double z = -log(1.0 - x);
1285 double temp = z * (1.0 - 3.0 * z / 8.0 * (1.0 - 17.0 * z / 81.0 * (1.0 - 15.0 * z / 136.0 *
1286 (1.0 - 28.0 * z / 1875.0 * (1.0 + 5.0 * z / 8.0 * (1.0 - 304.0 * z / 7203.0 *
1287 (1.0 + 945.0 * z / 2432.0 * (1.0 - 44.0 * z / 675.0 * (1.0 + 7.0 * z / 24.0 *
1288 (1.0 - 26104.0 * z / 307461.0 * (1.0 + 1925.0 * z / 8023.0 *
1289 (1.0 - 53598548.0 * z / 524808375.0 *
1290 (1.0 + 22232925.0 * z / 107197096.0)))))))))))));
1291 return temp;
1292} else if (x < x_3) {
1293 return Li3(x * x) / 4.0 - Li3(-x);
1294 } else if (x < x_4) {
1295 double ln1x = log(1.0 - x);
1296 return -Li3(1.0 - x) - Li3(-x / (1.0 - x)) + zeta3 + pisq6 * ln1x - log(x) * ln1x * ln1x / 2.0 + ln1x * ln1x * ln1x / 6.0;
1297 } else {
1298 double lnx = log(x);
1299 return Li3(1.0 / x) + 2.0 * pisq6 * lnx - lnx * lnx * lnx / 6.0;
1300 }
1301}
1302
1303scalar_t Li4(double x)
1304/* calculates the quadrilogarithm function of x */
1305{
1306 return polylog(3,1,x);
1307}
1308
1309/* calculates the dilogarithm function of x, extended to complex numbers */
1310/*-------------------------------------------------------------*/
1311
1312scalar_t CLi2(scalar_t x) {
1313 gsl_sf_result res_r;
1314 gsl_sf_result res_i;
1315 gsl_sf_complex_dilog_xy_e(x.real(), x.imag(), &res_r, &res_i);
1316 return scalar_t(res_r.val, res_i.val);
1317}
1318
1319
1320scalar_t CLi3(scalar_t x)
1321/* calculates the trilogarithm function of x, extended to complex numbers */
1322{
1323 return hpl3(0,0,1,x);
1324}
1325
1326scalar_t CLi4(scalar_t x)
1327/* calculates the quadrilogarithm function of x, extended to complex numbers */
1328{
1329 return hpl4(0,0,0,1,x);
1330}
1331
1332
1333double Cl2(double x) {
1334 return gsl_sf_clausen(x);
1335}
1336
1337double Cl3(double x) {
1338 scalar_t z = std::cos(x) + scalar_t(0, 1.) * std::sin(x);
1339 return std::real(CLi3(z));
1340}
1341
1342
1343
1344
scalar_t
Definition scalar.h:11
CLi2
scalar_t CLi2(scalar_t x)
Computes the complex dilogarithm function.
Definition polylog.cpp:1312
Cl2
double Cl2(double x)
Computes the Clausen function Cl2(x).
Definition polylog.cpp:1333
hpl1
scalar_t hpl1(int i, scalar_t x)
Definition polylog.cpp:123
cd
scalar_t cd(double x, double y)
Definition polylog.cpp:9
polylog
scalar_t polylog(size_t n, size_t m, double x)
Definition polylog.cpp:1109
hpl3
scalar_t hpl3(int i1, int i2, int i3, scalar_t x)
Definition polylog.cpp:170
hpl_base2
scalar_t hpl_base2(int i1, int i2, scalar_t x)
Definition polylog.cpp:19
hpl4
scalar_t hpl4(int i1, int i2, int i3, int i4, scalar_t x)
Definition polylog.cpp:224
hpl2
scalar_t hpl2(int i1, int i2, scalar_t x)
Definition polylog.cpp:129
zeta3
const double zeta3
Definition polylog.cpp:1265
Cl3
double Cl3(double x)
Computes the Clausen function Cl3(x).
Definition polylog.cpp:1337
Li3
double Li3(double x)
Computes the trilogarithm function Li3(x).
Definition polylog.cpp:1267
Li4
scalar_t Li4(double x)
Definition polylog.cpp:1303
CLi3
scalar_t CLi3(scalar_t x)
Computes the complex trilogarithm function.
Definition polylog.cpp:1320
hpl_base4
scalar_t hpl_base4(int i1, int i2, int i3, int i4, scalar_t x)
Definition polylog.cpp:82
CLi4
scalar_t CLi4(scalar_t x)
Computes the complex quadrilogarithm function.
Definition polylog.cpp:1326
hpl_base3
scalar_t hpl_base3(int i1, int i2, int i3, scalar_t x)
Definition polylog.cpp:40
initialize_constants
void initialize_constants(std::array< std::array< double, MAX_M >, MAX_N > &s1, std::array< std::array< double, MAX_M >, MAX_N > &c, std::array< std::array< double, MAX_NC >, MAX_LE > &a)
Definition polylog.cpp:827
Li2
double Li2(double x)
Computes the dilogarithm function Li2(x).
Definition polylog.cpp:1257
hpl_base1
scalar_t hpl_base1(int i, scalar_t x)
Definition polylog.cpp:13
pow
scalar_t pow(const scalar_t &base, const scalar_t &exp)
Definition scalar.cpp:75
real
double real(const scalar_t &z)
Definition scalar.cpp:83
fpeq
std::enable_if_t< not std::numeric_limits< T >::is_integer, bool > fpeq(T, T, std::size_t n=10)
Compares two floating point numbers with a given precision.
Definition special_generic.h:107
special_generic.h
sgn
int sgn(T val)
Returns the sign of a value.
Definition special_generic.h:93