a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
hpl.h
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2017 HEPfit Collaboration
3  *
4  *
5  * For the licensing terms see doc/COPYING.
6  */
7 
8 /* Code derived from arXiv:0810.4077
9  * authors: Christoph Greub, Volker Pilipp, Christof Schüpbach
10  */
11 
12 #ifndef HPL_H
13 #define HPL_H
14 
15 #include <limits>
16 #include <cmath>
17 #include <complex>
18 
19 typedef std::complex<double> cd;
20 
21 const std::numeric_limits<double> nld = *new std::numeric_limits<double>;
22 
23 inline cd operator*(int i, cd x)
24 {
25  return (double) i*x;
26 }
27 
28 inline cd operator*(cd x, int i)
29 {
30  return (double) i*x;
31 }
32 
33 inline cd hpl_base(int i, cd x)
34 {
35  if (i == 0 && norm(x) < nld.min()) return log(nld.min());
36  if (i == 0) return log(x);
37  if (i == 1 && norm(x) < nld.epsilon()) return x;
38  if (i == 1) return -log(1. - x);
39  else
40  return 0;
41 }
42 
43 inline cd hpl_base(int i1, int i2, cd x)
44 {
45  cd u;
46  if (norm(x) < nld.epsilon()) u = -x;
47  else u = log(1. - x);
48 
49  if (i1 == 0 && i2 == 1) return
50  -1. * u - 0.25 * pow(u, 2) - 0.027777777777777776 * pow(u, 3) +
51  0.0002777777777777778 * pow(u, 5) -
52  4.72411186696901e-6 * pow(u, 7) +
53  9.185773074661964e-8 * pow(u, 9) -
54  1.8978869988971e-9 * pow(u, 11) +
55  4.0647616451442256e-11 * pow(u, 13) -
56  8.921691020456452e-13 * pow(u, 15) +
57  1.9939295860721074e-14 * pow(u, 17) -
58  4.518980029619918e-16 * pow(u, 19) +
59  1.0356517612181247e-17 * pow(u, 21);
60  else
61  return 0;
62 }
63 
64 inline cd hpl_base(int i1, int i2, int i3, cd x)
65 {
66  cd u;
67  if (norm(x) < nld.epsilon()) u = -x;
68  else u = log(1. - x);
69 
70  if (i1 == 0 && i2 == 0 && i3 == 1)
71  return
72  -1. * u - 0.375 * pow(u, 2) - 0.0787037037037037 * pow(u, 3) -
73  0.008680555555555556 * pow(u, 4) -
74  0.00012962962962962963 * pow(u, 5) +
75  0.00008101851851851852 * pow(u, 6) +
76  3.4193571608537595e-6 * pow(u, 7) -
77  1.328656462585034e-6 * pow(u, 8) -
78  8.660871756109851e-8 * pow(u, 9) +
79  2.52608759553204e-8 * pow(u, 10) +
80  2.144694468364065e-9 * pow(u, 11) -
81  5.140110622012979e-10 * pow(u, 12) -
82  5.24958211460083e-11 * pow(u, 13) +
83  1.0887754406636318e-11 * pow(u, 14) +
84  1.2779396094493695e-12 * pow(u, 15) -
85  2.369824177308745e-13 * pow(u, 16) -
86  3.104357887965462e-14 * pow(u, 17) +
87  5.261758629912506e-15 * pow(u, 18) +
88  7.538479549949265e-16 * pow(u, 19) -
89  1.1862322577752286e-16 * pow(u, 20) -
90  1.8316979965491384e-17 * pow(u, 21);
91  if (i1 == 0 && i2 == 1 && i3 == 1)
92  return
93  0.25 * pow(u, 2) + 0.08333333333333333 * pow(u, 3) +
94  0.010416666666666666 * pow(u, 4) -
95  0.00011574074074074075 * pow(u, 6) +
96  2.066798941798942e-6 * pow(u, 8) -
97  4.1335978835978836e-8 * pow(u, 10) +
98  8.698648744945042e-10 * pow(u, 12) -
99  1.887210763816962e-11 * pow(u, 14) +
100  4.182042665838962e-13 * pow(u, 16) -
101  9.415778600896063e-15 * pow(u, 18) +
102  2.146515514069461e-16 * pow(u, 20);
103  else
104  return 0;
105 }
106 
107 inline cd hpl_base(int i1, int i2, int i3, int i4, cd x)
108 {
109  cd u;
110  if (norm(x) < nld.epsilon()) u = -x;
111  else u = log(1. - x);
112 
113  if (i1 == 0 && i2 == 0 && i3 == 0 && i4 == 1)
114  return
115  -1. * u - 0.4375 * pow(u, 2) - 0.11651234567901235 * pow(u, 3) -
116  0.019820601851851853 * pow(u, 4) -
117  0.001927932098765432 * pow(u, 5) -
118  0.000031057098765432096 * pow(u, 6) +
119  0.000015624009114857836 * pow(u, 7) +
120  8.485123546773206e-7 * pow(u, 8) -
121  2.290961660318971e-7 * pow(u, 9) -
122  2.1832614218526917e-8 * pow(u, 10) +
123  3.882824879172015e-9 * pow(u, 11) +
124  5.446292103220332e-10 * pow(u, 12) -
125  6.960805210682725e-11 * pow(u, 13) -
126  1.3375737686445216e-11 * pow(u, 14) +
127  1.2784852685266572e-12 * pow(u, 15) +
128  3.260562858024892e-13 * pow(u, 16) -
129  2.364757116861826e-14 * pow(u, 17) -
130  7.923135122031162e-15 * pow(u, 18) +
131  4.3452915709984186e-16 * pow(u, 19) +
132  1.923627006253592e-16 * pow(u, 20) -
133  7.812414333195955e-18 * pow(u, 21);
134  if (i1 == 0 && i2 == 1 && i3 == 0 && i4 == 1)
135  return
136  0.25 * pow(u, 2) + 0.1111111111111111 * pow(u, 3) +
137  0.022569444444444444 * pow(u, 4) +
138  0.0020833333333333333 * pow(u, 5) -
139  0.000027006172839506174 * pow(u, 6) -
140  0.00001984126984126984 * pow(u, 7) +
141  4.527273872511968e-7 * pow(u, 8) +
142  3.389987682315725e-7 * pow(u, 9) -
143  7.939132443100697e-9 * pow(u, 10) -
144  6.6805622361177916e-9 * pow(u, 11) +
145  1.4490216610627064e-10 * pow(u, 12) +
146  1.39908336457158e-10 * pow(u, 13) -
147  2.7425719106565973e-12 * pow(u, 14) -
148  3.032441227329819e-12 * pow(u, 15) +
149  5.358569182999823e-14 * pow(u, 16) +
150  6.724068599976371e-14 * pow(u, 17) -
151  1.0756816626218996e-15 * pow(u, 18) -
152  1.5158529016922455e-15 * pow(u, 19) +
153  2.208955024323606e-17 * pow(u, 20) +
154  3.460964863954937e-17 * pow(u, 21);
155  if (i1 == 0 && i2 == 1 && i3 == 1 && i4 == 1)
156  return
157  -0.05555555555555555 * pow(u, 3) -
158  0.020833333333333332 * pow(u, 4) -
159  0.002777777777777778 * pow(u, 5) +
160  0.00003306878306878307 * pow(u, 7) -
161  6.123848716441309e-7 * pow(u, 9) +
162  1.252605419272086e-8 * pow(u, 11) -
163  2.6765073061369356e-10 * pow(u, 13) +
164  5.871322376319437e-12 * pow(u, 15) -
165  1.312013385361243e-13 * pow(u, 17) +
166  2.97340376870402e-15 * pow(u, 19) -
167  6.814334965299877e-17 * pow(u, 21);
168  else
169  return 0;
170 }
171 
172 inline cd hpl(int i, cd x)
173 {
174  if (i == 0 && norm(x) < nld.min()) return log(nld.min());
175  if (i == 0) return log(x);
176  if (i == 1 && norm(x) < nld.epsilon()) return x;
177  if (i == 1) return -log(1. - x);
178  else
179  return 0;
180 }
181 
182 inline cd hpl(int i1, int i2, cd x)
183 {
184  const double Pi = M_PI;
185  if (i1 == 0 && i2 == 0) {
186  if (norm(x) > 1) return
187  pow(hpl(0, 1. / x), 2) / 2.;
188  if (real(x) > .5) return
189  pow(hpl(1, 1. - x), 2) / 2.;
190  return
191  pow(hpl_base(0, x), 2) / 2.;
192  }
193  if (i1 == 0 && i2 == 1) {
194  if (norm(x) > 1) return
195  Pi * cd(0, -1) * hpl(0, 1. / x) - hpl(0, 1, 1. / x) + pow(Pi, 2) / 3. -
196  pow(hpl(0, 1. / x), 2) / 2.;
197  if (real(x) > .5) return
198  hpl(0, 1. - x) * hpl(1, 1. - x) - hpl(0, 1, 1. - x) + pow(Pi, 2) / 6.;
199  return
200  hpl_base(0, 1, x);
201  }
202  if (i1 == 1 && i2 == 0) {
203  if (norm(x) > 1) return
204  Pi * cd(0, 1) * hpl(0, 1. / x) -
205  hpl(0, 1. / x)*(Pi * cd(0, 1) + hpl(0, 1. / x) + hpl(1, 1. / x)) +
206  hpl(0, 1, 1. / x) - pow(Pi, 2) / 3. + pow(hpl(0, 1. / x), 2) / 2.;
207  if (real(x) > .5) return
208  hpl(0, 1, 1. - x) - pow(Pi, 2) / 6.;
209  return
210  hpl_base(0, x) * hpl_base(1, x) - hpl_base(0, 1, x);
211  }
212  if (i1 == 1 && i2 == 1) {
213  if (norm(x) > 1) return
214  pow(Pi * cd(0, 1) + hpl(0, 1. / x) + hpl(1, 1. / x), 2) / 2.;
215  if (real(x) > .5) return
216  pow(hpl(0, 1. - x), 2) / 2.;
217  return
218  pow(hpl_base(1, x), 2) / 2.;
219  } else
220  return 0;
221 }
222 
223 inline cd hpl(int i1, int i2, int i3, cd x)
224 {
225  const double Pi = M_PI;
226  if (i1 == 0 && i2 == 0 && i3 == 0) {
227  if (norm(x) > 1) return
228  -pow(hpl(0, 1. / x), 3) / 6.;
229  if (real(x) > .5) return
230  -pow(hpl(1, 1. - x), 3) / 6.;
231  return
232  pow(hpl_base(0, x), 3) / 6.;
233  }
234  if (i1 == 0 && i2 == 0 && i3 == 1) {
235  if (norm(x) > 1) return
236  hpl(0, 0, 1, 1. / x) - (hpl(0, 1. / x) * pow(Pi, 2)) / 3. +
237  Pi * cd(0, 0.5) * pow(hpl(0, 1. / x), 2) + pow(hpl(0, 1. / x), 3) / 6.;
238  if (real(x) > .5) return
239  1.2020569031595942 + hpl(1, 1. - x) * hpl(0, 1, 1. - x) -
240  hpl(0, 1, 1, 1. - x) - (hpl(1, 1. - x) * pow(Pi, 2)) / 6. -
241  (hpl(0, 1. - x) * pow(hpl(1, 1. - x), 2)) / 2.;
242  return
243  hpl_base(0, 0, 1, x);
244  }
245  if (i1 == 0 && i2 == 1 && i3 == 0) {
246  if (norm(x) > 1) return
247  -(hpl(0, 1. / x)*(Pi * cd(0, -1) * hpl(0, 1. / x) - hpl(0, 1, 1. / x) +
248  pow(Pi, 2) / 3. - pow(hpl(0, 1. / x), 2) / 2.)) -
249  2 * (hpl(0, 0, 1, 1. / x) - (hpl(0, 1. / x) * pow(Pi, 2)) / 3. +
250  Pi * cd(0, 0.5) * pow(hpl(0, 1. / x), 2) +
251  pow(hpl(0, 1. / x), 3) / 6.);
252  if (real(x) > .5) return
253  -(hpl(1, 1. - x)*(hpl(0, 1. - x) * hpl(1, 1. - x) -
254  hpl(0, 1, 1. - x) + pow(Pi, 2) / 6.)) -
255  2 * (1.2020569031595942 + hpl(1, 1. - x) * hpl(0, 1, 1. - x) -
256  hpl(0, 1, 1, 1. - x) - (hpl(1, 1. - x) * pow(Pi, 2)) / 6. -
257  (hpl(0, 1. - x) * pow(hpl(1, 1. - x), 2)) / 2.);
258  return
259  hpl_base(0, x) * hpl_base(0, 1, x) - 2 * hpl_base(0, 0, 1, x);
260  }
261  if (i1 == 0 && i2 == 1 && i3 == 1) {
262  if (norm(x) > 1) return
263  1.2020569031595942 + Pi * cd(0, -1) * hpl(0, 1, 1. / x) -
264  hpl(0, 1. / x) * hpl(0, 1, 1. / x) + hpl(0, 0, 1, 1. / x) -
265  hpl(0, 1, 1, 1. / x) + (hpl(0, 1. / x) * pow(Pi, 2)) / 2. +
266  cd(0, 0.16666666666666666) * pow(Pi, 3) +
267  Pi * cd(0, -0.5) * pow(hpl(0, 1. / x), 2) - pow(hpl(0, 1. / x), 3) / 6.;
268  if (real(x) > .5) return
269  1.2020569031595942 + hpl(0, 1. - x) * hpl(0, 1, 1. - x) -
270  hpl(0, 0, 1, 1. - x) - (hpl(1, 1. - x) * pow(hpl(0, 1. - x), 2)) /
271  2.;
272  return
273  hpl_base(0, 1, 1, x);
274  }
275  if (i1 == 1 && i2 == 0 && i3 == 0) {
276  if (norm(x) > 1) return
277  hpl(0, 0, 1, 1. / x) - (hpl(0, 1. / x) * pow(Pi, 2)) / 3. +
278  hpl(0, 1. / x)*(Pi * cd(0, -1) * hpl(0, 1. / x) - hpl(0, 1, 1. / x) +
279  pow(Pi, 2) / 3. - pow(hpl(0, 1. / x), 2) / 2.) +
280  Pi * cd(0, 0.5) * pow(hpl(0, 1. / x), 2) +
281  ((Pi * cd(0, 1) + hpl(0, 1. / x) + hpl(1, 1. / x)) *
282  pow(hpl(0, 1. / x), 2)) / 2. + pow(hpl(0, 1. / x), 3) / 6.;
283  if (real(x) > .5) return
284  1.2020569031595942 + hpl(1, 1. - x) * hpl(0, 1, 1. - x) -
285  hpl(0, 1, 1, 1. - x) + hpl(1, 1. - x)*
286  (hpl(0, 1. - x) * hpl(1, 1. - x) - hpl(0, 1, 1. - x) +
287  pow(Pi, 2) / 6.) - (hpl(1, 1. - x) * pow(Pi, 2)) / 6. -
288  hpl(0, 1. - x) * pow(hpl(1, 1. - x), 2);
289  return
290  -(hpl_base(0, x) * hpl_base(0, 1, x)) + hpl_base(0, 0, 1, x) +
291  (hpl_base(1, x) * pow(hpl_base(0, x), 2)) / 2.;
292  }
293  if (i1 == 1 && i2 == 0 && i3 == 1) {
294  if (norm(x) > 1) return
295  (Pi * cd(0, 1) + hpl(0, 1. / x) + hpl(1, 1. / x))*
296  (Pi * cd(0, -1) * hpl(0, 1. / x) - hpl(0, 1, 1. / x) +
297  pow(Pi, 2) / 3. - pow(hpl(0, 1. / x), 2) / 2.) -
298  2 * (1.2020569031595942 + Pi * cd(0, -1) * hpl(0, 1, 1. / x) -
299  hpl(0, 1. / x) * hpl(0, 1, 1. / x) + hpl(0, 0, 1, 1. / x) -
300  hpl(0, 1, 1, 1. / x) + (hpl(0, 1. / x) * pow(Pi, 2)) / 2. +
301  cd(0, 0.16666666666666666) * pow(Pi, 3) +
302  Pi * cd(0, -0.5) * pow(hpl(0, 1. / x), 2) -
303  pow(hpl(0, 1. / x), 3) / 6.);
304  if (real(x) > .5) return
305  -(hpl(0, 1. - x)*(hpl(0, 1. - x) * hpl(1, 1. - x) -
306  hpl(0, 1, 1. - x) + pow(Pi, 2) / 6.)) -
307  2 * (1.2020569031595942 + hpl(0, 1. - x) * hpl(0, 1, 1. - x) -
308  hpl(0, 0, 1, 1. - x) -
309  (hpl(1, 1. - x) * pow(hpl(0, 1. - x), 2)) / 2.);
310  return
311  hpl_base(1, x) * hpl_base(0, 1, x) - 2 * hpl_base(0, 1, 1, x);
312  }
313  if (i1 == 1 && i2 == 1 && i3 == 0) {
314  if (norm(x) > 1) return
315  1.2020569031595942 + Pi * cd(0, -1) * hpl(0, 1, 1. / x) -
316  hpl(0, 1. / x) * hpl(0, 1, 1. / x) + hpl(0, 0, 1, 1. / x) -
317  hpl(0, 1, 1, 1. / x) + (hpl(0, 1. / x) * pow(Pi, 2)) / 2. +
318  cd(0, 0.16666666666666666) * pow(Pi, 3) -
319  (Pi * cd(0, 1) + hpl(0, 1. / x) + hpl(1, 1. / x))*
320  (Pi * cd(0, -1) * hpl(0, 1. / x) - hpl(0, 1, 1. / x) +
321  pow(Pi, 2) / 3. - pow(hpl(0, 1. / x), 2) / 2.) +
322  Pi * cd(0, -0.5) * pow(hpl(0, 1. / x), 2) -
323  pow(hpl(0, 1. / x), 3) / 6. -
324  (hpl(0, 1. / x) * pow(Pi * cd(0, 1) + hpl(0, 1. / x) + hpl(1, 1. / x),
325  2)) / 2.;
326  if (real(x) > .5) return
327  1.2020569031595942 + hpl(0, 1. - x) * hpl(0, 1, 1. - x) -
328  hpl(0, 0, 1, 1. - x) + hpl(0, 1. - x)*
329  (hpl(0, 1. - x) * hpl(1, 1. - x) - hpl(0, 1, 1. - x) +
330  pow(Pi, 2) / 6.) - hpl(1, 1. - x) * pow(hpl(0, 1. - x), 2);
331  return
332  -(hpl_base(1, x) * hpl_base(0, 1, x)) + hpl_base(0, 1, 1, x) +
333  (hpl_base(0, x) * pow(hpl_base(1, x), 2)) / 2.;
334  }
335  if (i1 == 1 && i2 == 1 && i3 == 1) {
336  if (norm(x) > 1) return
337  pow(Pi * cd(0, 1) + hpl(0, 1. / x) + hpl(1, 1. / x), 3) / 6.;
338  if (real(x) > .5) return
339  -pow(hpl(0, 1. - x), 3) / 6.;
340  return
341  pow(hpl_base(1, x), 3) / 6.;
342  } else
343  return 0;
344 }
345 
346 inline cd hpl(int i1, int i2, int i3, int i4, cd x)
347 {
348  const double Pi = M_PI;
349  if (i1 == 0 && i2 == 0 && i3 == 0 && i4 == 0) {
350  if (norm(x) > 1) {
351  return
352  pow(hpl(0, 1. / x), 4) / 24.;
353  }
354  if (real(x) > .5) {
355  return
356  pow(hpl(1, 1. - x), 4) / 24.;
357  }
358  return
359  pow(hpl_base(0, x), 4) / 24.;
360  }
361  if (i1 == 0 && i2 == 0 && i3 == 0 && i4 == 1) {
362  if (norm(x) > 1) {
363  return
364  -hpl(0, 0, 0, 1, 1. / x) + pow(Pi, 4) / 45. +
365  (pow(Pi, 2) * pow(hpl(0, 1. / x), 2)) / 6. +
366  Pi * cd(0, -0.16666666666666666) * pow(hpl(0, 1. / x), 3) -
367  pow(hpl(0, 1. / x), 4) / 24.;
368  }
369  if (real(x) > .5) {
370  return
371  -1.2020569031595942 * hpl(1, 1. - x) +
372  hpl(1, 1. - x) * hpl(0, 1, 1, 1. - x) - hpl(0, 1, 1, 1, 1. - x) +
373  pow(Pi, 4) / 90. - (hpl(0, 1, 1. - x) * pow(hpl(1, 1. - x), 2)) /
374  2. + (pow(Pi, 2) * pow(hpl(1, 1. - x), 2)) / 12. +
375  (hpl(0, 1. - x) * pow(hpl(1, 1. - x), 3)) / 6.;
376  }
377  return
378  hpl_base(0, 0, 0, 1, x);
379  }
380  if (i1 == 0 && i2 == 0 && i3 == 1 && i4 == 0) {
381  if (norm(x) > 1) {
382  return
383  -(hpl(0, 1. / x)*(hpl(0, 0, 1, 1. / x) -
384  (hpl(0, 1. / x) * pow(Pi, 2)) / 3. +
385  Pi * cd(0, 0.5) * pow(hpl(0, 1. / x), 2) +
386  pow(hpl(0, 1. / x), 3) / 6.)) -
387  3 * (-hpl(0, 0, 0, 1, 1. / x) + pow(Pi, 4) / 45. +
388  (pow(Pi, 2) * pow(hpl(0, 1. / x), 2)) / 6. +
389  Pi * cd(0, -0.16666666666666666) * pow(hpl(0, 1. / x), 3) -
390  pow(hpl(0, 1. / x), 4) / 24.);
391  }
392  if (real(x) > .5) {
393  return
394  -(hpl(1, 1. - x)*(1.2020569031595942 +
395  hpl(1, 1. - x) * hpl(0, 1, 1. - x) - hpl(0, 1, 1, 1. - x) -
396  (hpl(1, 1. - x) * pow(Pi, 2)) / 6. -
397  (hpl(0, 1. - x) * pow(hpl(1, 1. - x), 2)) / 2.)) -
398  3 * (-1.2020569031595942 * hpl(1, 1. - x) +
399  hpl(1, 1. - x) * hpl(0, 1, 1, 1. - x) -
400  hpl(0, 1, 1, 1, 1. - x) + pow(Pi, 4) / 90. -
401  (hpl(0, 1, 1. - x) * pow(hpl(1, 1. - x), 2)) / 2. +
402  (pow(Pi, 2) * pow(hpl(1, 1. - x), 2)) / 12. +
403  (hpl(0, 1. - x) * pow(hpl(1, 1. - x), 3)) / 6.);
404  }
405  return
406  hpl_base(0, x) * hpl_base(0, 0, 1, x) - 3 * hpl_base(0, 0, 0, 1, x);
407  }
408  if (i1 == 0 && i2 == 0 && i3 == 1 && i4 == 1) {
409  if (norm(x) > 1) {
410  return
411  (-2 * (3.7763731361630786 * cd(0, 2) +
412  2.4041138063191885 * hpl(0, 1. / x) +
413  Pi * cd(0, 1) * hpl(0, 1. / x) * hpl(0, 1, 1. / x) +
414  Pi * cd(0, -2) * hpl(0, 0, 1, 1. / x) -
415  2 * hpl(0, 1. / x) * hpl(0, 0, 1, 1. / x) +
416  4 * hpl(0, 0, 0, 1, 1. / x) + hpl(0, 1, 0, 1, 1. / x) -
417  (hpl(0, 1, 1. / x) * pow(Pi, 2)) / 3. + pow(Pi, 4) / 90. +
418  (hpl(0, 1, 1. / x) * pow(hpl(0, 1. / x), 2)) / 2. -
419  (pow(Pi, 2) * pow(hpl(0, 1. / x), 2)) / 6. +
420  Pi * cd(0, 0.16666666666666666) * pow(hpl(0, 1. / x), 3) +
421  pow(hpl(0, 1. / x), 4) / 24.) +
422  pow(Pi * cd(0, -1) * hpl(0, 1. / x) - hpl(0, 1, 1. / x) +
423  pow(Pi, 2) / 3. - pow(hpl(0, 1. / x), 2) / 2., 2)) / 4.;
424  }
425  if (real(x) > .5) {
426  return
427  (-2 * (2.4041138063191885 * hpl(1, 1. - x) +
428  hpl(0, 1. - x) * hpl(1, 1. - x) * hpl(0, 1, 1. - x) -
429  hpl(0, 1, 0, 1, 1. - x) +
430  (hpl(0, 1. - x) * hpl(1, 1. - x) * pow(Pi, 2)) / 6. -
431  (hpl(0, 1, 1. - x) * pow(Pi, 2)) / 6. + pow(Pi, 4) / 120. +
432  pow(hpl(0, 1, 1. - x), 2) -
433  2 * (hpl(1, 1. - x) * hpl(0, 0, 1, 1. - x) +
434  (2 * hpl(0, 1, 0, 1, 1. - x) - pow(hpl(0, 1, 1. - x), 2)) /
435  2. + (-2 * hpl(0, 1, 0, 1, 1. - x) +
436  pow(hpl(0, 1, 1. - x), 2)) / 2.) -
437  2 * (hpl(0, 1. - x) * hpl(0, 1, 1, 1. - x) +
438  (2 * hpl(0, 1, 0, 1, 1. - x) - pow(hpl(0, 1, 1. - x), 2)) /
439  2. + (-2 * hpl(0, 1, 0, 1, 1. - x) +
440  pow(hpl(0, 1, 1. - x), 2)) / 2.)) +
441  pow(hpl(0, 1. - x) * hpl(1, 1. - x) - hpl(0, 1, 1. - x) +
442  pow(Pi, 2) / 6., 2)) / 4.;
443  }
444  return
445  (-2 * hpl_base(0, 1, 0, 1, x) + pow(hpl_base(0, 1, x), 2)) / 4.;
446  }
447  if (i1 == 0 && i2 == 1 && i3 == 0 && i4 == 0) {
448  if (norm(x) > 1) {
449  return
450  ((Pi * cd(0, -1) * hpl(0, 1. / x) - hpl(0, 1, 1. / x) + pow(Pi, 2) / 3. -
451  pow(hpl(0, 1. / x), 2) / 2.) * pow(hpl(0, 1. / x), 2) +
452  4 * hpl(0, 1. / x)*(hpl(0, 0, 1, 1. / x) -
453  (hpl(0, 1. / x) * pow(Pi, 2)) / 3. +
454  Pi * cd(0, 0.5) * pow(hpl(0, 1. / x), 2) +
455  pow(hpl(0, 1. / x), 3) / 6.) +
456  6 * (-hpl(0, 0, 0, 1, 1. / x) + pow(Pi, 4) / 45. +
457  (pow(Pi, 2) * pow(hpl(0, 1. / x), 2)) / 6. +
458  Pi * cd(0, -0.16666666666666666) * pow(hpl(0, 1. / x), 3) -
459  pow(hpl(0, 1. / x), 4) / 24.)) / 2.;
460  }
461  if (real(x) > .5) {
462  return
463  ((hpl(0, 1. - x) * hpl(1, 1. - x) - hpl(0, 1, 1. - x) +
464  pow(Pi, 2) / 6.) * pow(hpl(1, 1. - x), 2) +
465  4 * hpl(1, 1. - x)*(1.2020569031595942 +
466  hpl(1, 1. - x) * hpl(0, 1, 1. - x) - hpl(0, 1, 1, 1. - x) -
467  (hpl(1, 1. - x) * pow(Pi, 2)) / 6. -
468  (hpl(0, 1. - x) * pow(hpl(1, 1. - x), 2)) / 2.) +
469  6 * (-1.2020569031595942 * hpl(1, 1. - x) +
470  hpl(1, 1. - x) * hpl(0, 1, 1, 1. - x) -
471  hpl(0, 1, 1, 1, 1. - x) + pow(Pi, 4) / 90. -
472  (hpl(0, 1, 1. - x) * pow(hpl(1, 1. - x), 2)) / 2. +
473  (pow(Pi, 2) * pow(hpl(1, 1. - x), 2)) / 12. +
474  (hpl(0, 1. - x) * pow(hpl(1, 1. - x), 3)) / 6.)) / 2.;
475  }
476  return
477  (-4 * hpl_base(0, x) * hpl_base(0, 0, 1, x) + 6 * hpl_base(0, 0, 0, 1, x) +
478  hpl_base(0, 1, x) * pow(hpl_base(0, x), 2)) / 2.;
479  }
480  if (i1 == 0 && i2 == 1 && i3 == 0 && i4 == 1) {
481  if (norm(x) > 1) {
482  return
483  3.7763731361630786 * cd(0, 2) +
484  2.4041138063191885 * hpl(0, 1. / x) +
485  Pi * cd(0, 1) * hpl(0, 1. / x) * hpl(0, 1, 1. / x) +
486  Pi * cd(0, -2) * hpl(0, 0, 1, 1. / x) -
487  2 * hpl(0, 1. / x) * hpl(0, 0, 1, 1. / x) + 4 * hpl(0, 0, 0, 1, 1. / x) +
488  hpl(0, 1, 0, 1, 1. / x) - (hpl(0, 1, 1. / x) * pow(Pi, 2)) / 3. +
489  pow(Pi, 4) / 90. + (hpl(0, 1, 1. / x) * pow(hpl(0, 1. / x), 2)) / 2. -
490  (pow(Pi, 2) * pow(hpl(0, 1. / x), 2)) / 6. +
491  Pi * cd(0, 0.16666666666666666) * pow(hpl(0, 1. / x), 3) +
492  pow(hpl(0, 1. / x), 4) / 24.;
493  }
494  if (real(x) > .5) {
495  return
496  2.4041138063191885 * hpl(1, 1. - x) +
497  hpl(0, 1. - x) * hpl(1, 1. - x) * hpl(0, 1, 1. - x) -
498  hpl(0, 1, 0, 1, 1. - x) +
499  (hpl(0, 1. - x) * hpl(1, 1. - x) * pow(Pi, 2)) / 6. -
500  (hpl(0, 1, 1. - x) * pow(Pi, 2)) / 6. + pow(Pi, 4) / 120. +
501  pow(hpl(0, 1, 1. - x), 2) -
502  2 * (hpl(1, 1. - x) * hpl(0, 0, 1, 1. - x) +
503  (2 * hpl(0, 1, 0, 1, 1. - x) - pow(hpl(0, 1, 1. - x), 2)) / 2. +
504  (-2 * hpl(0, 1, 0, 1, 1. - x) + pow(hpl(0, 1, 1. - x), 2)) / 2.)\
505  - 2 * (hpl(0, 1. - x) * hpl(0, 1, 1, 1. - x) +
506  (2 * hpl(0, 1, 0, 1, 1. - x) - pow(hpl(0, 1, 1. - x), 2)) / 2. +
507  (-2 * hpl(0, 1, 0, 1, 1. - x) + pow(hpl(0, 1, 1. - x), 2)) / 2.);
508  }
509  return
510  hpl_base(0, 1, 0, 1, x);
511  }
512  if (i1 == 0 && i2 == 1 && i3 == 1 && i4 == 0) {
513  if (norm(x) > 1) {
514  return
515  3.7763731361630786 * cd(0, -2) -
516  2.4041138063191885 * hpl(0, 1. / x) +
517  Pi * cd(0, -1) * hpl(0, 1. / x) * hpl(0, 1, 1. / x) +
518  Pi * cd(0, 2) * hpl(0, 0, 1, 1. / x) +
519  2 * hpl(0, 1. / x) * hpl(0, 0, 1, 1. / x) - 4 * hpl(0, 0, 0, 1, 1. / x) -
520  hpl(0, 1, 0, 1, 1. / x) + (hpl(0, 1, 1. / x) * pow(Pi, 2)) / 3. -
521  pow(Pi, 4) / 90. - (hpl(0, 1, 1. / x) * pow(hpl(0, 1. / x), 2)) / 2. +
522  (pow(Pi, 2) * pow(hpl(0, 1. / x), 2)) / 6. -
523  hpl(0, 1. / x)*(1.2020569031595942 +
524  Pi * cd(0, -1) * hpl(0, 1, 1. / x) -
525  hpl(0, 1. / x) * hpl(0, 1, 1. / x) + hpl(0, 0, 1, 1. / x) -
526  hpl(0, 1, 1, 1. / x) + (hpl(0, 1. / x) * pow(Pi, 2)) / 2. +
527  cd(0, 0.16666666666666666) * pow(Pi, 3) +
528  Pi * cd(0, -0.5) * pow(hpl(0, 1. / x), 2) -
529  pow(hpl(0, 1. / x), 3) / 6.) +
530  Pi * cd(0, -0.16666666666666666) * pow(hpl(0, 1. / x), 3) -
531  pow(hpl(0, 1. / x), 4) / 24. +
532  (2 * (3.7763731361630786 * cd(0, 2) +
533  2.4041138063191885 * hpl(0, 1. / x) +
534  Pi * cd(0, 1) * hpl(0, 1. / x) * hpl(0, 1, 1. / x) +
535  Pi * cd(0, -2) * hpl(0, 0, 1, 1. / x) -
536  2 * hpl(0, 1. / x) * hpl(0, 0, 1, 1. / x) +
537  4 * hpl(0, 0, 0, 1, 1. / x) + hpl(0, 1, 0, 1, 1. / x) -
538  (hpl(0, 1, 1. / x) * pow(Pi, 2)) / 3. + pow(Pi, 4) / 90. +
539  (hpl(0, 1, 1. / x) * pow(hpl(0, 1. / x), 2)) / 2. -
540  (pow(Pi, 2) * pow(hpl(0, 1. / x), 2)) / 6. +
541  Pi * cd(0, 0.16666666666666666) * pow(hpl(0, 1. / x), 3) +
542  pow(hpl(0, 1. / x), 4) / 24.) -
543  pow(Pi * cd(0, -1) * hpl(0, 1. / x) - hpl(0, 1, 1. / x) +
544  pow(Pi, 2) / 3. - pow(hpl(0, 1. / x), 2) / 2., 2)) / 2.;
545  }
546  if (real(x) > .5) {
547  return
548  -2.4041138063191885 * hpl(1, 1. - x) -
549  hpl(0, 1. - x) * hpl(1, 1. - x) * hpl(0, 1, 1. - x) +
550  hpl(0, 1, 0, 1, 1. - x) -
551  (hpl(0, 1. - x) * hpl(1, 1. - x) * pow(Pi, 2)) / 6. +
552  (hpl(0, 1, 1. - x) * pow(Pi, 2)) / 6. - pow(Pi, 4) / 120. -
553  hpl(1, 1. - x)*(1.2020569031595942 +
554  hpl(0, 1. - x) * hpl(0, 1, 1. - x) - hpl(0, 0, 1, 1. - x) -
555  (hpl(1, 1. - x) * pow(hpl(0, 1. - x), 2)) / 2.) -
556  pow(hpl(0, 1, 1. - x), 2) +
557  2 * (hpl(1, 1. - x) * hpl(0, 0, 1, 1. - x) +
558  (2 * hpl(0, 1, 0, 1, 1. - x) - pow(hpl(0, 1, 1. - x), 2)) / 2. +
559  (-2 * hpl(0, 1, 0, 1, 1. - x) + pow(hpl(0, 1, 1. - x), 2)) / 2.)\
560  + 2 * (hpl(0, 1. - x) * hpl(0, 1, 1, 1. - x) +
561  (2 * hpl(0, 1, 0, 1, 1. - x) - pow(hpl(0, 1, 1. - x), 2)) / 2. +
562  (-2 * hpl(0, 1, 0, 1, 1. - x) + pow(hpl(0, 1, 1. - x), 2)) / 2.)\
563  + (2 * (2.4041138063191885 * hpl(1, 1. - x) +
564  hpl(0, 1. - x) * hpl(1, 1. - x) * hpl(0, 1, 1. - x) -
565  hpl(0, 1, 0, 1, 1. - x) +
566  (hpl(0, 1. - x) * hpl(1, 1. - x) * pow(Pi, 2)) / 6. -
567  (hpl(0, 1, 1. - x) * pow(Pi, 2)) / 6. + pow(Pi, 4) / 120. +
568  pow(hpl(0, 1, 1. - x), 2) -
569  2 * (hpl(1, 1. - x) * hpl(0, 0, 1, 1. - x) +
570  (2 * hpl(0, 1, 0, 1, 1. - x) -
571  pow(hpl(0, 1, 1. - x), 2)) / 2. +
572  (-2 * hpl(0, 1, 0, 1, 1. - x) +
573  pow(hpl(0, 1, 1. - x), 2)) / 2.) -
574  2 * (hpl(0, 1. - x) * hpl(0, 1, 1, 1. - x) +
575  (2 * hpl(0, 1, 0, 1, 1. - x) -
576  pow(hpl(0, 1, 1. - x), 2)) / 2. +
577  (-2 * hpl(0, 1, 0, 1, 1. - x) +
578  pow(hpl(0, 1, 1. - x), 2)) / 2.)) -
579  pow(hpl(0, 1. - x) * hpl(1, 1. - x) - hpl(0, 1, 1. - x) +
580  pow(Pi, 2) / 6., 2)) / 2.;
581  }
582  return
583  hpl_base(0, x) * hpl_base(0, 1, 1, x) - hpl_base(0, 1, 0, 1, x) +
584  (2 * hpl_base(0, 1, 0, 1, x) - pow(hpl_base(0, 1, x), 2)) / 2.;
585  }
586  if (i1 == 0 && i2 == 1 && i3 == 1 && i4 == 1) {
587  if (norm(x) > 1) {
588  return
589  Pi * cd(0, -1) * hpl(0, 1. / x) * hpl(0, 1, 1. / x) +
590  Pi * cd(0, 1) * hpl(0, 0, 1, 1. / x) +
591  hpl(0, 1. / x) * hpl(0, 0, 1, 1. / x) +
592  Pi * cd(0, -1) * hpl(0, 1, 1, 1. / x) -
593  hpl(0, 1. / x) * hpl(0, 1, 1, 1. / x) - hpl(0, 0, 0, 1, 1. / x) -
594  hpl(0, 1, 0, 1, 1. / x) / 2. - hpl(0, 1, 1, 1, 1. / x) +
595  (hpl(0, 1, 1. / x) * pow(Pi, 2)) / 2. +
596  cd(0, 0.16666666666666666) * hpl(0, 1. / x) * pow(Pi, 3) -
597  (19 * pow(Pi, 4)) / 360. -
598  (hpl(0, 1, 1. / x) * pow(hpl(0, 1. / x), 2)) / 2. +
599  (pow(Pi, 2) * pow(hpl(0, 1. / x), 2)) / 4. +
600  Pi * cd(0, -0.16666666666666666) * pow(hpl(0, 1. / x), 3) -
601  pow(hpl(0, 1. / x), 4) / 24. +
602  (2 * hpl(0, 1, 0, 1, 1. / x) - pow(hpl(0, 1, 1. / x), 2)) / 2. +
603  pow(hpl(0, 1, 1. / x), 2) / 4. +
604  (-2 * hpl(0, 1, 0, 1, 1. / x) + pow(hpl(0, 1, 1. / x), 2)) / 2.;
605  }
606  if (real(x) > .5) {
607  return
608  hpl(0, 1. - x) * hpl(0, 0, 1, 1. - x) - hpl(0, 0, 0, 1, 1. - x) +
609  pow(Pi, 4) / 90. - (hpl(0, 1, 1. - x) * pow(hpl(0, 1. - x), 2)) /
610  2. + (hpl(1, 1. - x) * pow(hpl(0, 1. - x), 3)) / 6.;
611  }
612  return
613  hpl_base(0, 1, 1, 1, x);
614  }
615  if (i1 == 1 && i2 == 0 && i3 == 0 && i4 == 0) {
616  if (norm(x) > 1) {
617  return
618  (-3 * (Pi * cd(0, -1) * hpl(0, 1. / x) - hpl(0, 1, 1. / x) +
619  pow(Pi, 2) / 3. - pow(hpl(0, 1. / x), 2) / 2.) *
620  pow(hpl(0, 1. / x), 2) -
621  6 * hpl(0, 1. / x)*(hpl(0, 0, 1, 1. / x) -
622  (hpl(0, 1. / x) * pow(Pi, 2)) / 3. +
623  Pi * cd(0, 0.5) * pow(hpl(0, 1. / x), 2) +
624  pow(hpl(0, 1. / x), 3) / 6.) -
625  (Pi * cd(0, 1) + hpl(0, 1. / x) + hpl(1, 1. / x)) *
626  pow(hpl(0, 1. / x), 3) -
627  6 * (-hpl(0, 0, 0, 1, 1. / x) + pow(Pi, 4) / 45. +
628  (pow(Pi, 2) * pow(hpl(0, 1. / x), 2)) / 6. +
629  Pi * cd(0, -0.16666666666666666) * pow(hpl(0, 1. / x), 3) -
630  pow(hpl(0, 1. / x), 4) / 24.)) / 6.;
631  }
632  if (real(x) > .5) {
633  return
634  (-3 * (hpl(0, 1. - x) * hpl(1, 1. - x) - hpl(0, 1, 1. - x) +
635  pow(Pi, 2) / 6.) * pow(hpl(1, 1. - x), 2) -
636  6 * hpl(1, 1. - x)*(1.2020569031595942 +
637  hpl(1, 1. - x) * hpl(0, 1, 1. - x) - hpl(0, 1, 1, 1. - x) -
638  (hpl(1, 1. - x) * pow(Pi, 2)) / 6. -
639  (hpl(0, 1. - x) * pow(hpl(1, 1. - x), 2)) / 2.) +
640  hpl(0, 1. - x) * pow(hpl(1, 1. - x), 3) -
641  6 * (-1.2020569031595942 * hpl(1, 1. - x) +
642  hpl(1, 1. - x) * hpl(0, 1, 1, 1. - x) -
643  hpl(0, 1, 1, 1, 1. - x) + pow(Pi, 4) / 90. -
644  (hpl(0, 1, 1. - x) * pow(hpl(1, 1. - x), 2)) / 2. +
645  (pow(Pi, 2) * pow(hpl(1, 1. - x), 2)) / 12. +
646  (hpl(0, 1. - x) * pow(hpl(1, 1. - x), 3)) / 6.)) / 6.;
647  }
648  return
649  (6 * hpl_base(0, x) * hpl_base(0, 0, 1, x) - 6 * hpl_base(0, 0, 0, 1, x) -
650  3 * hpl_base(0, 1, x) * pow(hpl_base(0, x), 2) + hpl_base(1, x) * pow(hpl_base(0, x), 3)
651  ) / 6.;
652  }
653  if (i1 == 1 && i2 == 0 && i3 == 0 && i4 == 1) {
654  if (norm(x) > 1) {
655  return
656  3.7763731361630786 * cd(0, -2) -
657  2.4041138063191885 * hpl(0, 1. / x) +
658  Pi * cd(0, -1) * hpl(0, 1. / x) * hpl(0, 1, 1. / x) +
659  Pi * cd(0, 2) * hpl(0, 0, 1, 1. / x) +
660  2 * hpl(0, 1. / x) * hpl(0, 0, 1, 1. / x) - 4 * hpl(0, 0, 0, 1, 1. / x) -
661  hpl(0, 1, 0, 1, 1. / x) + (hpl(0, 1, 1. / x) * pow(Pi, 2)) / 3. -
662  pow(Pi, 4) / 90. - (hpl(0, 1, 1. / x) * pow(hpl(0, 1. / x), 2)) / 2. +
663  (pow(Pi, 2) * pow(hpl(0, 1. / x), 2)) / 6. +
664  (Pi * cd(0, 1) + hpl(0, 1. / x) + hpl(1, 1. / x))*
665  (hpl(0, 0, 1, 1. / x) - (hpl(0, 1. / x) * pow(Pi, 2)) / 3. +
666  Pi * cd(0, 0.5) * pow(hpl(0, 1. / x), 2) +
667  pow(hpl(0, 1. / x), 3) / 6.) +
668  Pi * cd(0, -0.16666666666666666) * pow(hpl(0, 1. / x), 3) -
669  pow(hpl(0, 1. / x), 4) / 24. +
670  (2 * (3.7763731361630786 * cd(0, 2) +
671  2.4041138063191885 * hpl(0, 1. / x) +
672  Pi * cd(0, 1) * hpl(0, 1. / x) * hpl(0, 1, 1. / x) +
673  Pi * cd(0, -2) * hpl(0, 0, 1, 1. / x) -
674  2 * hpl(0, 1. / x) * hpl(0, 0, 1, 1. / x) +
675  4 * hpl(0, 0, 0, 1, 1. / x) + hpl(0, 1, 0, 1, 1. / x) -
676  (hpl(0, 1, 1. / x) * pow(Pi, 2)) / 3. + pow(Pi, 4) / 90. +
677  (hpl(0, 1, 1. / x) * pow(hpl(0, 1. / x), 2)) / 2. -
678  (pow(Pi, 2) * pow(hpl(0, 1. / x), 2)) / 6. +
679  Pi * cd(0, 0.16666666666666666) * pow(hpl(0, 1. / x), 3) +
680  pow(hpl(0, 1. / x), 4) / 24.) -
681  pow(Pi * cd(0, -1) * hpl(0, 1. / x) - hpl(0, 1, 1. / x) +
682  pow(Pi, 2) / 3. - pow(hpl(0, 1. / x), 2) / 2., 2)) / 2.;
683  }
684  if (real(x) > .5) {
685  return
686  -2.4041138063191885 * hpl(1, 1. - x) -
687  hpl(0, 1. - x) * hpl(1, 1. - x) * hpl(0, 1, 1. - x) +
688  hpl(0, 1, 0, 1, 1. - x) -
689  (hpl(0, 1. - x) * hpl(1, 1. - x) * pow(Pi, 2)) / 6. +
690  (hpl(0, 1, 1. - x) * pow(Pi, 2)) / 6. - pow(Pi, 4) / 120. -
691  hpl(0, 1. - x)*(1.2020569031595942 +
692  hpl(1, 1. - x) * hpl(0, 1, 1. - x) - hpl(0, 1, 1, 1. - x) -
693  (hpl(1, 1. - x) * pow(Pi, 2)) / 6. -
694  (hpl(0, 1. - x) * pow(hpl(1, 1. - x), 2)) / 2.) -
695  pow(hpl(0, 1, 1. - x), 2) +
696  2 * (hpl(1, 1. - x) * hpl(0, 0, 1, 1. - x) +
697  (2 * hpl(0, 1, 0, 1, 1. - x) - pow(hpl(0, 1, 1. - x), 2)) / 2. +
698  (-2 * hpl(0, 1, 0, 1, 1. - x) + pow(hpl(0, 1, 1. - x), 2)) / 2.)\
699  + 2 * (hpl(0, 1. - x) * hpl(0, 1, 1, 1. - x) +
700  (2 * hpl(0, 1, 0, 1, 1. - x) - pow(hpl(0, 1, 1. - x), 2)) / 2. +
701  (-2 * hpl(0, 1, 0, 1, 1. - x) + pow(hpl(0, 1, 1. - x), 2)) / 2.)\
702  + (2 * (2.4041138063191885 * hpl(1, 1. - x) +
703  hpl(0, 1. - x) * hpl(1, 1. - x) * hpl(0, 1, 1. - x) -
704  hpl(0, 1, 0, 1, 1. - x) +
705  (hpl(0, 1. - x) * hpl(1, 1. - x) * pow(Pi, 2)) / 6. -
706  (hpl(0, 1, 1. - x) * pow(Pi, 2)) / 6. + pow(Pi, 4) / 120. +
707  pow(hpl(0, 1, 1. - x), 2) -
708  2 * (hpl(1, 1. - x) * hpl(0, 0, 1, 1. - x) +
709  (2 * hpl(0, 1, 0, 1, 1. - x) -
710  pow(hpl(0, 1, 1. - x), 2)) / 2. +
711  (-2 * hpl(0, 1, 0, 1, 1. - x) +
712  pow(hpl(0, 1, 1. - x), 2)) / 2.) -
713  2 * (hpl(0, 1. - x) * hpl(0, 1, 1, 1. - x) +
714  (2 * hpl(0, 1, 0, 1, 1. - x) -
715  pow(hpl(0, 1, 1. - x), 2)) / 2. +
716  (-2 * hpl(0, 1, 0, 1, 1. - x) +
717  pow(hpl(0, 1, 1. - x), 2)) / 2.)) -
718  pow(hpl(0, 1. - x) * hpl(1, 1. - x) - hpl(0, 1, 1. - x) +
719  pow(Pi, 2) / 6., 2)) / 2.;
720  }
721  return
722  hpl_base(1, x) * hpl_base(0, 0, 1, x) - hpl_base(0, 1, 0, 1, x) +
723  (2 * hpl_base(0, 1, 0, 1, x) - pow(hpl_base(0, 1, x), 2)) / 2.;
724  }
725  if (i1 == 1 && i2 == 0 && i3 == 1 && i4 == 0) {
726  if (norm(x) > 1) {
727  return
728  3.7763731361630786 * cd(0, -2) -
729  2.4041138063191885 * hpl(0, 1. / x) +
730  Pi * cd(0, -1) * hpl(0, 1. / x) * hpl(0, 1, 1. / x) +
731  Pi * cd(0, 2) * hpl(0, 0, 1, 1. / x) +
732  2 * hpl(0, 1. / x) * hpl(0, 0, 1, 1. / x) - 4 * hpl(0, 0, 0, 1, 1. / x) -
733  hpl(0, 1, 0, 1, 1. / x) + (hpl(0, 1, 1. / x) * pow(Pi, 2)) / 3. -
734  pow(Pi, 4) / 90. - hpl(0, 1. / x)*
735  (Pi * cd(0, 1) + hpl(0, 1. / x) + hpl(1, 1. / x))*
736  (Pi * cd(0, -1) * hpl(0, 1. / x) - hpl(0, 1, 1. / x) +
737  pow(Pi, 2) / 3. - pow(hpl(0, 1. / x), 2) / 2.) -
738  (hpl(0, 1, 1. / x) * pow(hpl(0, 1. / x), 2)) / 2. +
739  (pow(Pi, 2) * pow(hpl(0, 1. / x), 2)) / 6. +
740  2 * hpl(0, 1. / x)*(1.2020569031595942 +
741  Pi * cd(0, -1) * hpl(0, 1, 1. / x) -
742  hpl(0, 1. / x) * hpl(0, 1, 1. / x) + hpl(0, 0, 1, 1. / x) -
743  hpl(0, 1, 1, 1. / x) + (hpl(0, 1. / x) * pow(Pi, 2)) / 2. +
744  cd(0, 0.16666666666666666) * pow(Pi, 3) +
745  Pi * cd(0, -0.5) * pow(hpl(0, 1. / x), 2) -
746  pow(hpl(0, 1. / x), 3) / 6.) -
747  2 * (Pi * cd(0, 1) + hpl(0, 1. / x) + hpl(1, 1. / x))*
748  (hpl(0, 0, 1, 1. / x) - (hpl(0, 1. / x) * pow(Pi, 2)) / 3. +
749  Pi * cd(0, 0.5) * pow(hpl(0, 1. / x), 2) +
750  pow(hpl(0, 1. / x), 3) / 6.) +
751  Pi * cd(0, -0.16666666666666666) * pow(hpl(0, 1. / x), 3) -
752  pow(hpl(0, 1. / x), 4) / 24. +
753  pow(Pi * cd(0, -1) * hpl(0, 1. / x) - hpl(0, 1, 1. / x) +
754  pow(Pi, 2) / 3. - pow(hpl(0, 1. / x), 2) / 2., 2);
755  }
756  if (real(x) > .5) {
757  return
758  -2.4041138063191885 * hpl(1, 1. - x) -
759  hpl(0, 1. - x) * hpl(1, 1. - x) * hpl(0, 1, 1. - x) +
760  hpl(0, 1, 0, 1, 1. - x) +
761  hpl(0, 1. - x) * hpl(1, 1. - x)*
762  (hpl(0, 1. - x) * hpl(1, 1. - x) - hpl(0, 1, 1. - x) +
763  pow(Pi, 2) / 6.) - (hpl(0, 1. - x) * hpl(1, 1. - x) *
764  pow(Pi, 2)) / 6. + (hpl(0, 1, 1. - x) * pow(Pi, 2)) / 6. -
765  pow(Pi, 4) / 120. + 2 * hpl(1, 1. - x)*
766  (1.2020569031595942 + hpl(0, 1. - x) * hpl(0, 1, 1. - x) -
767  hpl(0, 0, 1, 1. - x) -
768  (hpl(1, 1. - x) * pow(hpl(0, 1. - x), 2)) / 2.) +
769  2 * hpl(0, 1. - x)*(1.2020569031595942 +
770  hpl(1, 1. - x) * hpl(0, 1, 1. - x) - hpl(0, 1, 1, 1. - x) -
771  (hpl(1, 1. - x) * pow(Pi, 2)) / 6. -
772  (hpl(0, 1. - x) * pow(hpl(1, 1. - x), 2)) / 2.) -
773  pow(hpl(0, 1, 1. - x), 2) +
774  2 * (hpl(1, 1. - x) * hpl(0, 0, 1, 1. - x) +
775  (2 * hpl(0, 1, 0, 1, 1. - x) - pow(hpl(0, 1, 1. - x), 2)) / 2. +
776  (-2 * hpl(0, 1, 0, 1, 1. - x) + pow(hpl(0, 1, 1. - x), 2)) / 2.)\
777  + 2 * (hpl(0, 1. - x) * hpl(0, 1, 1, 1. - x) +
778  (2 * hpl(0, 1, 0, 1, 1. - x) - pow(hpl(0, 1, 1. - x), 2)) / 2. +
779  (-2 * hpl(0, 1, 0, 1, 1. - x) + pow(hpl(0, 1, 1. - x), 2)) / 2.)\
780  + pow(hpl(0, 1. - x) * hpl(1, 1. - x) - hpl(0, 1, 1. - x) +
781  pow(Pi, 2) / 6., 2);
782  }
783  return
784  hpl_base(0, x) * hpl_base(1, x) * hpl_base(0, 1, x) - 2 * hpl_base(1, x) * hpl_base(0, 0, 1, x) -
785  2 * hpl_base(0, x) * hpl_base(0, 1, 1, x) - hpl_base(0, 1, 0, 1, x) +
786  pow(hpl_base(0, 1, x), 2);
787  }
788  if (i1 == 1 && i2 == 0 && i3 == 1 && i4 == 1) {
789  if (norm(x) > 1) {
790  return
791  (Pi * cd(0, 1) + hpl(0, 1. / x) + hpl(1, 1. / x))*
792  (1.2020569031595942 + Pi * cd(0, -1) * hpl(0, 1, 1. / x) -
793  hpl(0, 1. / x) * hpl(0, 1, 1. / x) + hpl(0, 0, 1, 1. / x) -
794  hpl(0, 1, 1, 1. / x) + (hpl(0, 1. / x) * pow(Pi, 2)) / 2. +
795  cd(0, 0.16666666666666666) * pow(Pi, 3) +
796  Pi * cd(0, -0.5) * pow(hpl(0, 1. / x), 2) -
797  pow(hpl(0, 1. / x), 3) / 6.) -
798  3 * (Pi * cd(0, -1) * hpl(0, 1. / x) * hpl(0, 1, 1. / x) +
799  Pi * cd(0, 1) * hpl(0, 0, 1, 1. / x) +
800  hpl(0, 1. / x) * hpl(0, 0, 1, 1. / x) +
801  Pi * cd(0, -1) * hpl(0, 1, 1, 1. / x) -
802  hpl(0, 1. / x) * hpl(0, 1, 1, 1. / x) - hpl(0, 0, 0, 1, 1. / x) -
803  hpl(0, 1, 0, 1, 1. / x) / 2. - hpl(0, 1, 1, 1, 1. / x) +
804  (hpl(0, 1, 1. / x) * pow(Pi, 2)) / 2. +
805  cd(0, 0.16666666666666666) * hpl(0, 1. / x) * pow(Pi, 3) -
806  (19 * pow(Pi, 4)) / 360. -
807  (hpl(0, 1, 1. / x) * pow(hpl(0, 1. / x), 2)) / 2. +
808  (pow(Pi, 2) * pow(hpl(0, 1. / x), 2)) / 4. +
809  Pi * cd(0, -0.16666666666666666) * pow(hpl(0, 1. / x), 3) -
810  pow(hpl(0, 1. / x), 4) / 24. +
811  (2 * hpl(0, 1, 0, 1, 1. / x) - pow(hpl(0, 1, 1. / x), 2)) / 2. +
812  pow(hpl(0, 1, 1. / x), 2) / 4. +
813  (-2 * hpl(0, 1, 0, 1, 1. / x) + pow(hpl(0, 1, 1. / x), 2)) / 2.);
814  }
815  if (real(x) > .5) {
816  return
817  -(hpl(0, 1. - x)*(1.2020569031595942 +
818  hpl(0, 1. - x) * hpl(0, 1, 1. - x) - hpl(0, 0, 1, 1. - x) -
819  (hpl(1, 1. - x) * pow(hpl(0, 1. - x), 2)) / 2.)) -
820  3 * (hpl(0, 1. - x) * hpl(0, 0, 1, 1. - x) -
821  hpl(0, 0, 0, 1, 1. - x) + pow(Pi, 4) / 90. -
822  (hpl(0, 1, 1. - x) * pow(hpl(0, 1. - x), 2)) / 2. +
823  (hpl(1, 1. - x) * pow(hpl(0, 1. - x), 3)) / 6.);
824  }
825  return
826  hpl_base(1, x) * hpl_base(0, 1, 1, x) - 3 * hpl_base(0, 1, 1, 1, x);
827  }
828  if (i1 == 1 && i2 == 1 && i3 == 0 && i4 == 0) {
829  if (norm(x) > 1) {
830  return
831  hpl(0, 1. / x)*(Pi * cd(0, 1) + hpl(0, 1. / x) + hpl(1, 1. / x))*
832  (Pi * cd(0, -1) * hpl(0, 1. / x) - hpl(0, 1, 1. / x) +
833  pow(Pi, 2) / 3. - pow(hpl(0, 1. / x), 2) / 2.) -
834  hpl(0, 1. / x)*(1.2020569031595942 +
835  Pi * cd(0, -1) * hpl(0, 1, 1. / x) -
836  hpl(0, 1. / x) * hpl(0, 1, 1. / x) + hpl(0, 0, 1, 1. / x) -
837  hpl(0, 1, 1, 1. / x) + (hpl(0, 1. / x) * pow(Pi, 2)) / 2. +
838  cd(0, 0.16666666666666666) * pow(Pi, 3) +
839  Pi * cd(0, -0.5) * pow(hpl(0, 1. / x), 2) -
840  pow(hpl(0, 1. / x), 3) / 6.) +
841  (Pi * cd(0, 1) + hpl(0, 1. / x) + hpl(1, 1. / x))*
842  (hpl(0, 0, 1, 1. / x) - (hpl(0, 1. / x) * pow(Pi, 2)) / 3. +
843  Pi * cd(0, 0.5) * pow(hpl(0, 1. / x), 2) +
844  pow(hpl(0, 1. / x), 3) / 6.) +
845  (pow(hpl(0, 1. / x), 2) * pow(Pi * cd(0, 1) + hpl(0, 1. / x) +
846  hpl(1, 1. / x), 2)) / 4. +
847  (2 * (3.7763731361630786 * cd(0, 2) +
848  2.4041138063191885 * hpl(0, 1. / x) +
849  Pi * cd(0, 1) * hpl(0, 1. / x) * hpl(0, 1, 1. / x) +
850  Pi * cd(0, -2) * hpl(0, 0, 1, 1. / x) -
851  2 * hpl(0, 1. / x) * hpl(0, 0, 1, 1. / x) +
852  4 * hpl(0, 0, 0, 1, 1. / x) + hpl(0, 1, 0, 1, 1. / x) -
853  (hpl(0, 1, 1. / x) * pow(Pi, 2)) / 3. + pow(Pi, 4) / 90. +
854  (hpl(0, 1, 1. / x) * pow(hpl(0, 1. / x), 2)) / 2. -
855  (pow(Pi, 2) * pow(hpl(0, 1. / x), 2)) / 6. +
856  Pi * cd(0, 0.16666666666666666) * pow(hpl(0, 1. / x), 3) +
857  pow(hpl(0, 1. / x), 4) / 24.) -
858  pow(Pi * cd(0, -1) * hpl(0, 1. / x) - hpl(0, 1, 1. / x) +
859  pow(Pi, 2) / 3. - pow(hpl(0, 1. / x), 2) / 2., 2)) / 4.;
860  }
861  if (real(x) > .5) {
862  return
863  -(hpl(0, 1. - x) * hpl(1, 1. - x)*
864  (hpl(0, 1. - x) * hpl(1, 1. - x) - hpl(0, 1, 1. - x) +
865  pow(Pi, 2) / 6.)) -
866  hpl(1, 1. - x)*(1.2020569031595942 +
867  hpl(0, 1. - x) * hpl(0, 1, 1. - x) - hpl(0, 0, 1, 1. - x) -
868  (hpl(1, 1. - x) * pow(hpl(0, 1. - x), 2)) / 2.) +
869  (pow(hpl(0, 1. - x), 2) * pow(hpl(1, 1. - x), 2)) / 4. -
870  hpl(0, 1. - x)*(1.2020569031595942 +
871  hpl(1, 1. - x) * hpl(0, 1, 1. - x) - hpl(0, 1, 1, 1. - x) -
872  (hpl(1, 1. - x) * pow(Pi, 2)) / 6. -
873  (hpl(0, 1. - x) * pow(hpl(1, 1. - x), 2)) / 2.) +
874  (2 * (2.4041138063191885 * hpl(1, 1. - x) +
875  hpl(0, 1. - x) * hpl(1, 1. - x) * hpl(0, 1, 1. - x) -
876  hpl(0, 1, 0, 1, 1. - x) +
877  (hpl(0, 1. - x) * hpl(1, 1. - x) * pow(Pi, 2)) / 6. -
878  (hpl(0, 1, 1. - x) * pow(Pi, 2)) / 6. + pow(Pi, 4) / 120. +
879  pow(hpl(0, 1, 1. - x), 2) -
880  2 * (hpl(1, 1. - x) * hpl(0, 0, 1, 1. - x) +
881  (2 * hpl(0, 1, 0, 1, 1. - x) -
882  pow(hpl(0, 1, 1. - x), 2)) / 2. +
883  (-2 * hpl(0, 1, 0, 1, 1. - x) +
884  pow(hpl(0, 1, 1. - x), 2)) / 2.) -
885  2 * (hpl(0, 1. - x) * hpl(0, 1, 1, 1. - x) +
886  (2 * hpl(0, 1, 0, 1, 1. - x) -
887  pow(hpl(0, 1, 1. - x), 2)) / 2. +
888  (-2 * hpl(0, 1, 0, 1, 1. - x) +
889  pow(hpl(0, 1, 1. - x), 2)) / 2.)) -
890  pow(hpl(0, 1. - x) * hpl(1, 1. - x) - hpl(0, 1, 1. - x) +
891  pow(Pi, 2) / 6., 2)) / 4.;
892  }
893  return
894  -(hpl_base(0, x) * hpl_base(1, x) * hpl_base(0, 1, x)) + hpl_base(1, x) * hpl_base(0, 0, 1, x) +
895  hpl_base(0, x) * hpl_base(0, 1, 1, x) +
896  (pow(hpl_base(0, x), 2) * pow(hpl_base(1, x), 2)) / 4. +
897  (2 * hpl_base(0, 1, 0, 1, x) - pow(hpl_base(0, 1, x), 2)) / 4.;
898  }
899  if (i1 == 1 && i2 == 1 && i3 == 0 && i4 == 1) {
900  if (norm(x) > 1) {
901  return
902  (-4 * (Pi * cd(0, 1) + hpl(0, 1. / x) + hpl(1, 1. / x))*
903  (1.2020569031595942 + Pi * cd(0, -1) * hpl(0, 1, 1. / x) -
904  hpl(0, 1. / x) * hpl(0, 1, 1. / x) + hpl(0, 0, 1, 1. / x) -
905  hpl(0, 1, 1, 1. / x) + (hpl(0, 1. / x) * pow(Pi, 2)) / 2. +
906  cd(0, 0.16666666666666666) * pow(Pi, 3) +
907  Pi * cd(0, -0.5) * pow(hpl(0, 1. / x), 2) -
908  pow(hpl(0, 1. / x), 3) / 6.) +
909  (Pi * cd(0, -1) * hpl(0, 1. / x) - hpl(0, 1, 1. / x) +
910  pow(Pi, 2) / 3. - pow(hpl(0, 1. / x), 2) / 2.) *
911  pow(Pi * cd(0, 1) + hpl(0, 1. / x) + hpl(1, 1. / x), 2) +
912  6 * (Pi * cd(0, -1) * hpl(0, 1. / x) * hpl(0, 1, 1. / x) +
913  Pi * cd(0, 1) * hpl(0, 0, 1, 1. / x) +
914  hpl(0, 1. / x) * hpl(0, 0, 1, 1. / x) +
915  Pi * cd(0, -1) * hpl(0, 1, 1, 1. / x) -
916  hpl(0, 1. / x) * hpl(0, 1, 1, 1. / x) - hpl(0, 0, 0, 1, 1. / x) -
917  hpl(0, 1, 0, 1, 1. / x) / 2. - hpl(0, 1, 1, 1, 1. / x) +
918  (hpl(0, 1, 1. / x) * pow(Pi, 2)) / 2. +
919  cd(0, 0.16666666666666666) * hpl(0, 1. / x) * pow(Pi, 3) -
920  (19 * pow(Pi, 4)) / 360. -
921  (hpl(0, 1, 1. / x) * pow(hpl(0, 1. / x), 2)) / 2. +
922  (pow(Pi, 2) * pow(hpl(0, 1. / x), 2)) / 4. +
923  Pi * cd(0, -0.16666666666666666) * pow(hpl(0, 1. / x), 3) -
924  pow(hpl(0, 1. / x), 4) / 24. +
925  (2 * hpl(0, 1, 0, 1, 1. / x) - pow(hpl(0, 1, 1. / x), 2)) / 2. +
926  pow(hpl(0, 1, 1. / x), 2) / 4. +
927  (-2 * hpl(0, 1, 0, 1, 1. / x) + pow(hpl(0, 1, 1. / x), 2)) / 2.)) /
928  2.;
929  }
930  if (real(x) > .5) {
931  return
932  ((hpl(0, 1. - x) * hpl(1, 1. - x) - hpl(0, 1, 1. - x) +
933  pow(Pi, 2) / 6.) * pow(hpl(0, 1. - x), 2) +
934  4 * hpl(0, 1. - x)*(1.2020569031595942 +
935  hpl(0, 1. - x) * hpl(0, 1, 1. - x) - hpl(0, 0, 1, 1. - x) -
936  (hpl(1, 1. - x) * pow(hpl(0, 1. - x), 2)) / 2.) +
937  6 * (hpl(0, 1. - x) * hpl(0, 0, 1, 1. - x) -
938  hpl(0, 0, 0, 1, 1. - x) + pow(Pi, 4) / 90. -
939  (hpl(0, 1, 1. - x) * pow(hpl(0, 1. - x), 2)) / 2. +
940  (hpl(1, 1. - x) * pow(hpl(0, 1. - x), 3)) / 6.)) / 2.;
941  }
942  return
943  (-4 * hpl_base(1, x) * hpl_base(0, 1, 1, x) + 6 * hpl_base(0, 1, 1, 1, x) +
944  hpl_base(0, 1, x) * pow(hpl_base(1, x), 2)) / 2.;
945  }
946  if (i1 == 1 && i2 == 1 && i3 == 1 && i4 == 0) {
947  if (norm(x) > 1) {
948  return
949  (6 * (Pi * cd(0, 1) + hpl(0, 1. / x) + hpl(1, 1. / x))*
950  (1.2020569031595942 + Pi * cd(0, -1) * hpl(0, 1, 1. / x) -
951  hpl(0, 1. / x) * hpl(0, 1, 1. / x) + hpl(0, 0, 1, 1. / x) -
952  hpl(0, 1, 1, 1. / x) + (hpl(0, 1. / x) * pow(Pi, 2)) / 2. +
953  cd(0, 0.16666666666666666) * pow(Pi, 3) +
954  Pi * cd(0, -0.5) * pow(hpl(0, 1. / x), 2) -
955  pow(hpl(0, 1. / x), 3) / 6.) -
956  3 * (Pi * cd(0, -1) * hpl(0, 1. / x) - hpl(0, 1, 1. / x) +
957  pow(Pi, 2) / 3. - pow(hpl(0, 1. / x), 2) / 2.) *
958  pow(Pi * cd(0, 1) + hpl(0, 1. / x) + hpl(1, 1. / x), 2) -
959  hpl(0, 1. / x) * pow(Pi * cd(0, 1) + hpl(0, 1. / x) + hpl(1, 1. / x),
960  3) - 6 * (Pi * cd(0, -1) * hpl(0, 1. / x) * hpl(0, 1, 1. / x) +
961  Pi * cd(0, 1) * hpl(0, 0, 1, 1. / x) +
962  hpl(0, 1. / x) * hpl(0, 0, 1, 1. / x) +
963  Pi * cd(0, -1) * hpl(0, 1, 1, 1. / x) -
964  hpl(0, 1. / x) * hpl(0, 1, 1, 1. / x) - hpl(0, 0, 0, 1, 1. / x) -
965  hpl(0, 1, 0, 1, 1. / x) / 2. - hpl(0, 1, 1, 1, 1. / x) +
966  (hpl(0, 1, 1. / x) * pow(Pi, 2)) / 2. +
967  cd(0, 0.16666666666666666) * hpl(0, 1. / x) * pow(Pi, 3) -
968  (19 * pow(Pi, 4)) / 360. -
969  (hpl(0, 1, 1. / x) * pow(hpl(0, 1. / x), 2)) / 2. +
970  (pow(Pi, 2) * pow(hpl(0, 1. / x), 2)) / 4. +
971  Pi * cd(0, -0.16666666666666666) * pow(hpl(0, 1. / x), 3) -
972  pow(hpl(0, 1. / x), 4) / 24. +
973  (2 * hpl(0, 1, 0, 1, 1. / x) - pow(hpl(0, 1, 1. / x), 2)) / 2. +
974  pow(hpl(0, 1, 1. / x), 2) / 4. +
975  (-2 * hpl(0, 1, 0, 1, 1. / x) + pow(hpl(0, 1, 1. / x), 2)) / 2.)) /
976  6.;
977  }
978  if (real(x) > .5) {
979  return
980  (-3 * (hpl(0, 1. - x) * hpl(1, 1. - x) - hpl(0, 1, 1. - x) +
981  pow(Pi, 2) / 6.) * pow(hpl(0, 1. - x), 2) -
982  6 * hpl(0, 1. - x)*(1.2020569031595942 +
983  hpl(0, 1. - x) * hpl(0, 1, 1. - x) - hpl(0, 0, 1, 1. - x) -
984  (hpl(1, 1. - x) * pow(hpl(0, 1. - x), 2)) / 2.) +
985  hpl(1, 1. - x) * pow(hpl(0, 1. - x), 3) -
986  6 * (hpl(0, 1. - x) * hpl(0, 0, 1, 1. - x) -
987  hpl(0, 0, 0, 1, 1. - x) + pow(Pi, 4) / 90. -
988  (hpl(0, 1, 1. - x) * pow(hpl(0, 1. - x), 2)) / 2. +
989  (hpl(1, 1. - x) * pow(hpl(0, 1. - x), 3)) / 6.)) / 6.;
990  }
991  return
992  (6 * hpl_base(1, x) * hpl_base(0, 1, 1, x) - 6 * hpl_base(0, 1, 1, 1, x) -
993  3 * hpl_base(0, 1, x) * pow(hpl_base(1, x), 2) + hpl_base(0, x) * pow(hpl_base(1, x), 3)
994  ) / 6.;
995  }
996  if (i1 == 1 && i2 == 1 && i3 == 1 && i4 == 1) {
997  if (norm(x) > 1) {
998  return
999  pow(Pi * cd(0, 1) + hpl(0, 1. / x) + hpl(1, 1. / x), 4) / 24.;
1000  }
1001  if (real(x) > .5) {
1002  return
1003  pow(hpl(0, 1. - x), 4) / 24.;
1004  }
1005  return
1006  pow(hpl_base(1, x), 4) / 24.;
1007  } else
1008  return 0;
1009 }
1010 
1011 inline cd Li2(cd x)
1012 {
1013  return hpl(0, 1, x);
1014 }
1015 
1016 inline cd Li3(cd x)
1017 {
1018  return hpl(0, 0, 1, x);
1019 }
1020 
1021 inline cd Li4(cd x)
1022 {
1023  return hpl(0, 0, 0, 1, x);
1024 }
1025 
1026 inline double Cl2(double x)
1027 {
1028  return imag(Li2(exp(cd(0, x))));
1029 }
1030 
1031 inline double Cl3(double x)
1032 {
1033  return real(Li3(exp(cd(0, x))));
1034 }
1035 
1036 #endif /* HPL_H */
Li3
cd Li3(cd x)
Definition: hpl.h:1016
Li2
cd Li2(cd x)
Definition: hpl.h:1011
Li4
cd Li4(cd x)
Definition: hpl.h:1021
gslpp::log
complex log(const complex &z)
Definition: gslpp_complex.cpp:342
operator*
cd operator*(int i, cd x)
Definition: hpl.h:23
gslpp::pow
complex pow(const complex &z1, const complex &z2)
Definition: gslpp_complex.cpp:395
hpl
cd hpl(int i, cd x)
Definition: hpl.h:172
cd
std::complex< double > cd
Definition: hpl.h:19
nld
const std::numeric_limits< double > nld
Definition: hpl.h:21
Cl2
double Cl2(double x)
Definition: hpl.h:1026
Cl3
double Cl3(double x)
Definition: hpl.h:1031
gslpp::exp
complex exp(const complex &z)
Definition: gslpp_complex.cpp:333
hpl_base
cd hpl_base(int i, cd x)
Definition: hpl.h:33