Hyperiso 1.0.3
Modular flavour-physics calculations, Wilson coefficients and statistical inference
Loading...
Searching...
No Matches
special_SM.cpp
Go to the documentation of this file.
1#include "special_SM.h"
2#include "special_SUSY.h"
3#include "special_THDM.h"
4
5// MAJ : Complete separation between SM and THDM/SUSY according to headers
6
7/*--------------------------------------------------------------------*/
8
9double Bplus(double x, double y)
10{
11 return y/(x-y)*(log(y)/(y-1.)-log(x)/(x-1.));
12}
13
14/*---------------------------------------------------------------------*/
15
16
17double D3(double x)
18{
19 if(fabs(x)<1.e-5) return 0.;
20 if(fabs(x-1.)<1.e-5) return -1.;
21
22 return x*log(x)/(1.-x);
23}
24
25
26/*--------------------------------------------------------------------*/
27
28double D2(double x, double y)
29{
30 if(fabs(x-y)<1.e-5)
31 {
32 if(fabs(x-1.)<1.e-5) return -0.5;
33 return (1.-x+log(x))/(1.-x)/(1.-x);
34 }
35
36 return (D3(x)-D3(y))/(x-y);
37}
38
39/*---------------------------------------------------------------------*/
40
41
42double h10(double x)
43{
44 if(fabs(1.-x)<1.e-5) return h10(0.9999);
45
46 return (3.*x*x-2.*x)/3./pow(x-1.,4.)*log(x)
47 -(8.*x*x+5.*x-7.)/18./pow(x-1.,3.);
48}
49
50/*--------------------------------------------------------------------*/
51
52double h20(double x)
53{
54 if(fabs(1.-x)<1.e-5) return h20(0.9999);
55
56 return (-6.*x*x+4.*x)/3./pow(x-1.,3.)*log(x)
57 +(7.*x-5.)/3./pow(x-1.,2.);
58}
59
60/*--------------------------------------------------------------------*/
61
62double h30(double x)
63{
64 if(fabs(1.-x)<1.e-5) return h30(0.9999);
65
66 return (-6.*x*x*x+9.*x*x-2)/9./pow(x-1.,4.)*log(x)
67 +(52.*x*x-101.*x+43.)/54./pow(x-1.,3.);
68}
69
70/*----------------------------------------------------------------------*/
71
72double h40(double x)
73{
74 if(fabs(1.-x)<1.e-5) return h40(0.9999);
75
76 return -log(x)/3./pow(x-1.,4.)+(2.*x*x-7.*x+11.)/18./pow(x-1.,3.);
77}
78
79/*----------------------------------------------------------------------*/
80
81double h50(double x)
82{
83 if(fabs(1.-x)<1.e-5) return h50(0.9999);
84
85 return -x/pow(x-1.,4.)*log(x)+(-x*x+5.*x+2.)/6./pow(x-1.,3.);
86}
87
88/*----------------------------------------------------------------------*/
89
90double h60(double x)
91{
92 if(fabs(1.-x)<1.e-5) return h60(0.9999);
93
94 return 2.*x/pow(x-1.,3.)*log(x)-(x+1.)/pow(x-1.,2.);
95}
96
97/*----------------------------------------------------------------------*/
98
99double f20(double x)
100{
101 if(fabs(1.-x)<1.e-5) return f20(0.9999);
102
103 return -x/(x-1.)*(1.-1./(x-1.)*log(x));
104}
105
106/*----------------------------------------------------------------------*/
107
108double f30(double x, double y)
109{
110 if((fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)) return f30(0.9998,1.0002);
111
112 if(fabs(1.-x)<1.e-5) return f30(0.9999,y);
113 if(fabs(1.-y)<1.e-5) return f30(x,0.9999);
114
115 if(fabs(1.-x/y)<1.e-5) return f30(y*0.9998,y);
116
117 return x*log(x)/(x-1.)/(x-y)+y*log(y)/(y-1.)/(y-x);
118}
119
120/*----------------------------------------------------------------------*/
121
122double f40(double x, double y)
123{
124 if((fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)) return f40(0.9998,1.0002);
125
126 if(fabs(1.-x)<1.e-5) return f40(0.9999,y);
127 if(fabs(1.-y)<1.e-5) return f40(x,0.9999);
128
129 if(fabs(1.-x/y)<1.e-5) return f40(y*0.9998,y);
130
131 return x*x*log(x)/(x-1.)/(x-y)+y*y*log(y)/(y-1.)/(y-x);
132}
133
134/*----------------------------------------------------------------------*/
135
136double f50(double x, double y, double z)
137{
138 if((fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)&&(fabs(1.-z)<1.e-5)) return f50(0.9996,0.9998,1.0002);
139
140 if((fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)) return f50(0.9998,1.0002,z);
141 if((fabs(1.-x)<1.e-5)&&(fabs(1.-z)<1.e-5)) return f50(0.9998,y,1.0002);
142 if((fabs(1.-y)<1.e-5)&&(fabs(1.-z)<1.e-5)) return f50(x,0.9998,1.0002);
143
144 if(fabs(1.-x)<1.e-5) return f50(0.9999,y,z);
145 if(fabs(1.-y)<1.e-5) return f50(x,0.9999,z);
146 if(fabs(1.-z)<1.e-5) return f50(x,y,0.9999);
147
148 if(fabs(1.-x/y)<1.e-5) return f50(y*0.9998,y,z);
149 if(fabs(1.-x/y)<1.e-5) return f50(y*0.9998,y,z);
150 if(fabs(1.-y/z)<1.e-5) return f50(x,z*0.9998,z);
151 if(fabs(1.-x/z)<1.e-5) return f50(x,y,x*0.9998);
152
153 return x*x*log(x)/(x-1.)/(x-y)/(x-z)+y*y*log(y)/(y-1.)/(y-x)/(y-z)
154 +z*z*log(z)/(z-1.)/(z-x)/(z-y);
155}
156
157/*----------------------------------------------------------------------*/
158
159double f60(double x, double y, double z)
160{
161 if((fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)&&(fabs(1.-z)<1.e-5)) return f60(0.9996,0.9998,1.0002);
162
163 if((fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)) return f60(0.9998,1.0002,z);
164 if((fabs(1.-x)<1.e-5)&&(fabs(1.-z)<1.e-5)) return f60(0.9998,y,1.0002);
165 if((fabs(1.-y)<1.e-5)&&(fabs(1.-z)<1.e-5)) return f60(x,0.9998,1.0002);
166
167 if(fabs(1.-x)<1.e-5) return f60(0.9999,y,z);
168 if(fabs(1.-y)<1.e-5) return f60(x,0.9999,z);
169 if(fabs(1.-z)<1.e-5) return f60(x,y,0.9999);
170
171 if(fabs(1.-x/y)<1.e-5) return f60(y*0.9998,y,z);
172 if(fabs(1.-y/z)<1.e-5) return f60(x,z*0.9998,z);
173 if(fabs(1.-x/z)<1.e-5) return f60(x,y,x*0.9998);
174
175 return x*log(x)/(x-1.)/(x-y)/(x-z)+y*log(y)/(y-1.)/(y-x)/(y-z)
176 +z*log(z)/(z-1.)/(z-x)/(z-y);
177}
178
179/*----------------------------------------------------------------------*/
180
181double f70(double x, double y)
182{
183 if((fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)) return f70(0.9998,1.0002);
184
185 if(fabs(1.-x)<1.e-5) return f70(0.9999,y);
186 if(fabs(1.-y)<1.e-5) return f70(x,0.9999);
187
188 if(fabs(1.-x/y)<1.e-5) return f70(y*0.9998,y);
189
190 return x*log(x)/(x-1.)/(x-y)+x*log(y)/(y-1.)/(y-x);
191}
192
193/*----------------------------------------------------------------------*/
194
195double f80(double x)
196{
197 if(fabs(x-1.)<1.e-5) return 1.;
198
199 return x*log(x)/(x-1.);
200}
201
202/*----------------------------------------------------------------------*/
203
204double f90(double w, double x, double y, double z)
205{
206 if((fabs(1.-w)<1.e-5)&&(fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)&&(fabs(1.-z)<1.e-5)) return f90(0.9996,0.9998,1.0002,1.0004);
207
208 if((fabs(1.-w)<1.e-5)&&(fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)) return f90(0.9996,0.9998,1.0002,z);
209 if((fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)&&(fabs(1.-z)<1.e-5)) return f90(w,0.9996,0.9998,1.0002);
210 if((fabs(1.-w)<1.e-5)&&(fabs(1.-y)<1.e-5)&&(fabs(1.-z)<1.e-5)) return f90(0.9996,x,0.9998,1.0002);
211 if((fabs(1.-w)<1.e-5)&&(fabs(1.-x)<1.e-5)&&(fabs(1.-z)<1.e-5)) return f90(0.9996,0.9998,y,1.0002);
212
213 if((fabs(1.-w)<1.e-5)&&(fabs(1.-x)<1.e-5)) return f90(0.9998,1.0002,y,z);
214 if((fabs(1.-w)<1.e-5)&&(fabs(1.-y)<1.e-5)) return f90(0.9998,x,1.0002,z);
215 if((fabs(1.-w)<1.e-5)&&(fabs(1.-z)<1.e-5)) return f90(0.9998,x,y,1.0002);
216 if((fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)) return f90(w,0.9998,1.0002,z);
217 if((fabs(1.-x)<1.e-5)&&(fabs(1.-z)<1.e-5)) return f90(w,0.9998,y,1.0002);
218 if((fabs(1.-y)<1.e-5)&&(fabs(1.-z)<1.e-5)) return f90(w,x,0.9998,1.0002);
219
220 if(fabs(1.-w)<1.e-5) return f90(0.9999,x,y,z);
221 if(fabs(1.-x)<1.e-5) return f90(w,0.9999,y,z);
222 if(fabs(1.-y)<1.e-5) return f90(w,x,0.9999,z);
223 if(fabs(1.-z)<1.e-5) return f90(w,x,y,0.9999);
224
225 if(fabs(1.-w/x)<1.e-5) return f90(x*0.9998,x,y,z);
226 if(fabs(1.-w/y)<1.e-5) return f90(y*0.9998,x,y,z);
227 if(fabs(1.-w/z)<1.e-5) return f90(z*0.9998,x,y,z);
228 if(fabs(1.-x/y)<1.e-5) return f90(w,y*0.9998,y,z);
229 if(fabs(1.-y/z)<1.e-5) return f90(w,x,z*0.9998,z);
230 if(fabs(1.-x/z)<1.e-5) return f90(w,x,y,x*0.9998);
231
232
233 return w*w*log(w)/(w-1.)/(w-x)/(w-y)/(w-z) +x*x*log(x)/(x-1.)/(x-w)/(x-y)/(x-z)
234 +y*y*log(y)/(y-1.)/(y-x)/(y-w)/(y-z) +z*z*log(z)/(z-1.)/(z-x)/(z-y)/(z-w);
235}
236
237/*----------------------------------------------------------------------*/
238
239double f100(double w, double x, double y, double z)
240{
241 if((fabs(1.-w)<1.e-5)&&(fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)&&(fabs(1.-z)<1.e-5)) return f100(0.9996,0.9998,1.0002,1.0004);
242
243 if((fabs(1.-w)<1.e-5)&&(fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)) return f100(0.9996,0.9998,1.0002,z);
244 if((fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)&&(fabs(1.-z)<1.e-5)) return f100(w,0.9996,0.9998,1.0002);
245 if((fabs(1.-w)<1.e-5)&&(fabs(1.-y)<1.e-5)&&(fabs(1.-z)<1.e-5)) return f100(0.9996,x,0.9998,1.0002);
246 if((fabs(1.-w)<1.e-5)&&(fabs(1.-x)<1.e-5)&&(fabs(1.-z)<1.e-5)) return f100(0.9996,0.9998,y,1.0002);
247
248 if((fabs(1.-w)<1.e-5)&&(fabs(1.-x)<1.e-5)) return f100(0.9998,1.0002,y,z);
249 if((fabs(1.-w)<1.e-5)&&(fabs(1.-y)<1.e-5)) return f100(0.9998,x,1.0002,z);
250 if((fabs(1.-w)<1.e-5)&&(fabs(1.-z)<1.e-5)) return f100(0.9998,x,y,1.0002);
251 if((fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)) return f100(w,0.9998,1.0002,z);
252 if((fabs(1.-x)<1.e-5)&&(fabs(1.-z)<1.e-5)) return f100(w,0.9998,y,1.0002);
253 if((fabs(1.-y)<1.e-5)&&(fabs(1.-z)<1.e-5)) return f100(w,x,0.9998,1.0002);
254
255 if(fabs(1.-w)<1.e-5) return f100(0.9999,x,y,z);
256 if(fabs(1.-x)<1.e-5) return f100(w,0.9999,y,z);
257 if(fabs(1.-y)<1.e-5) return f100(w,x,0.9999,z);
258 if(fabs(1.-z)<1.e-5) return f100(w,x,y,0.9999);
259
260 if(fabs(1.-w/x)<1.e-5) return f100(x*0.9998,x,y,z);
261 if(fabs(1.-w/y)<1.e-5) return f100(y*0.9998,x,y,z);
262 if(fabs(1.-w/z)<1.e-5) return f100(z*0.9998,x,y,z);
263 if(fabs(1.-x/y)<1.e-5) return f100(w,y*0.9998,y,z);
264 if(fabs(1.-y/z)<1.e-5) return f100(w,x,z*0.9998,z);
265 if(fabs(1.-x/z)<1.e-5) return f100(w,x,y,x*0.9998);
266
267 return w*log(w)/(w-1.)/(w-x)/(w-y)/(w-z) +x*log(x)/(x-1.)/(x-w)/(x-y)/(x-z)
268 +y*log(y)/(y-1.)/(y-x)/(y-w)/(y-z) +z*log(z)/(z-1.)/(z-x)/(z-y)/(z-w);
269}
270
271/*----------------------------------------------------------------------*/
272
273double f110(double x, double y)
274{
275 if((fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)) return f110(0.9998,1.0002);
276
277 if(fabs(1.-x)<1.e-5) return f110(0.9999,y);
278 if(fabs(1.-y)<1.e-5) return f110(x,0.9999);
279
280 if(fabs(1.-x/y)<1.e-5) return f110(y*0.9998,y);
281
282 return x*log(x)/(x-y)+x*log(y)/(y-x);
283}
284
285/*----------------------------------------------------------------------*/
286
287double h11(double x, double y)
288{
289 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return h11(0.98,1.02);
290
291 if(fabs(1.-x)<1.e-3) return h11(0.99,y);
292 if(fabs(1.-y)<1.e-3) return h11(x,0.99);
293
294 if(fabs(1.-x/y)<5.e-3) return h11(y*0.98,y);
295
296 return ((-48.*x*x*x-104.*x*x+64.*x)*Li2(1.-1./x)
297 +(-378.*x*x*x-1566.*x*x+850.*x+86.)/9./(x-1.)*log(x)
298 +(2060.*x*x*x+3798.*x*x-2664.*x-170.)/27.
299 +((12.*x*x*x-124.*x*x+64.*x)/(x-1.)*log(x)+(-56.*x*x*x+258.*x*x+24.*x-82.)/3.)*y)/9./pow(x-1.,4.);
300}
301
302/*----------------------------------------------------------------------*/
303
304double h21(double x, double y)
305{
306 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return h21(0.98,1.02);
307
308 if(fabs(1.-x)<1.e-3) return h21(0.99,y);
309 if(fabs(1.-y)<1.e-3) return h21(x,0.99);
310
311 if(fabs(1.-x/y)<5.e-3) return h21(y*0.98,y);
312
313 return ((224.*x*x-96.*x)*Li2(1.-1./x)
314 +(-24.*x*x*x+352.*x*x-128.*x-32.)*log(x)/(x-1.)
315 +(-340.*x*x+132.*x+40.)
316 +((-24.*x*x*x+176.*x*x-80.*x)*log(x)/(x-1.)+(-28.*x*x-108.*x+64.))*y)/9./pow(x-1.,3.);
317}
318
319/*----------------------------------------------------------------------*/
320
321double h31(double x, double y)
322{
323 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return h31(0.98,1.02);
324
325 if(fabs(1.-x)<1.e-3) return h31(0.99,y);
326 if(fabs(1.-y)<1.e-3) return h31(x,0.99);
327
328 if(fabs(1.-x/y)<5.e-3) return h31(y*0.98,y);
329
330 return (32.*x*x*x+120.*x*x-384.*x+128.)*Li2(1.-1./x)/81./pow(x-1.,4.)
331 +(-108.*x*x*x*x+1058.*x*x*x-898.*x*x-1098.*x+710.)/81./pow(x-1.,5.)*log(x)
332 +(-304.*x*x*x-13686.*x*x+29076.*x-12062.)/729./pow(x-1.,4.)
333 +((540.*x*x*x-972.*x*x+232.*x+56.)/81./pow(x-1.,5.)*log(x)
334 +(-664.*x*x*x+54.*x*x+1944.*x-902.)/243./pow(x-1.,4.))*y;
335}
336
337/*----------------------------------------------------------------------*/
338
339double h41(double x, double y)
340{
341 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return h41(0.98,1.02);
342
343 if(fabs(1.-x)<1.e-3) return h41(0.99,y);
344 if(fabs(1.-y)<1.e-3) return h41(x,0.99);
345
346 if(fabs(1.-x/y)<5.e-3) return h41(y*0.98,y);
347
348 return (-562.*x*x*x+1101.*x*x-420.*x+101.)/54./pow(x-1.,4.)*Li2(1.-1./x)
349 +(-562.*x*x*x+1604.*x*x-799.*x+429.)/54./pow(x-1.,5.)*log(x)
350 +(17470.*x*x*x-47217.*x*x+31098.*x-13447.)/972./pow(x-1.,4.)
351 +((89.*x+55.)*log(x)/27./pow(x-1.,5.)+(38.*x*x*x-135.*x*x+54.*x-821.)/162./pow(x-1.,4.))*y;
352}
353
354/*----------------------------------------------------------------------*/
355
356double h51(double x, double y)
357{
358 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return h51(0.98,1.02);
359
360 if(fabs(1.-x)<1.e-3) return h51(0.99,y);
361 if(fabs(1.-y)<1.e-3) return h51(x,0.99);
362
363 if(fabs(1.-x/y)<5.e-3) return h51(y*0.98,y);
364
365 return ((9.*x*x*x+46.*x*x+49.*x)*Li2(1.-1./x)/2.
366 +(81.*x*x*x+594.*x*x+1270.*x+71.)*log(x)/18./(x-1.)
367 +(-923.*x*x*x-3042.*x*x-6921.*x-1210.)/108.
368 +((10.*x*x+38.*x)/(x-1.)*log(x)+(-7.*x*x*x+30.*x*x-141.*x-26.)/3.)*y)/3./pow(x-1.,4.);
369}
370
371/*----------------------------------------------------------------------*/
372
373double h61(double x, double y)
374{
375 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return h61(0.98,1.02);
376
377 if(fabs(1.-x)<1.e-3) return h61(0.99,y);
378 if(fabs(1.-y)<1.e-3) return h61(x,0.99);
379
380 if(fabs(1.-x/y)<5.e-3) return h61(y*0.98,y);
381
382 return ((-32.*x*x-24.*x)*Li2(1.-1./x)
383 +(-52.*x*x-109.*x-7.)*log(x)/(x-1.)
384 +(95.*x*x+180.*x+61.)/2.
385 +((-20.*x*x-52.*x)/(x-1.)*log(x)+(-2.*x*x+60.*x+14.))*y)/3./pow(x-1.,3.);
386}
387
388/*----------------------------------------------------------------------*/
389
390double h71(double x, double y)
391{
392 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return h71(0.98,1.02);
393
394 if(fabs(1.-x)<1.e-3) return h71(0.99,y);
395 if(fabs(1.-y)<1.e-3) return h71(x,0.99);
396
397 if(fabs(1.-x/y)<5.e-3) return h71(y*0.98,y);
398
399 return (-20.*x*x*x+60.*x*x-60.*x-20.)*Li2(1.-1./x)/27./pow(x-1.,4.)
400 +(-60.*x*x+240.*x+4.)/81./pow(x-1.,4.)*log(x)
401 +(132.*x*x-382.*x+186.)/81./pow(x-1.,3.)
402 +(20.*log(x)/27./pow(x-1.,4.)+(-20.*x*x+70.*x-110.)/81./pow(x-1.,3.))*y;
403}
404
405/*----------------------------------------------------------------------*/
406
407double f31(double x, double y)
408{
409 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return f31(0.98,1.02);
410
411 if(fabs(1.-x)<1.e-3) return f31(0.99,y);
412 if(fabs(1.-y)<1.e-3) return f31(x,0.99);
413
414 if(fabs(1.-x/y)<5.e-3) return f31(y*0.98,y);
415
416 return -28.*y/3./(y-1.)/(x-y)+2.*x*(11.*x+3.*y)/3./(x-1.)/(x-y)/(x-y)*log(x)
417 +2.*y*(25.*x-11.*x*y-11.*y-3.*y*y)/3./(y-1.)/(y-1.)/(x-y)/(x-y)*log(y)
418 +4.*(1.+y)/(x-1.)/(y-1)*Li2(1.-1./y)+4.*(x+y)/(x-1.)/(x-y)*Li2(1.-x/y);
419}
420
421/*----------------------------------------------------------------------*/
422
423double f41(double x, double y)
424{
425 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return f41(0.98,1.02);
426
427 if(fabs(1.-x)<1.e-3) return f41(0.99,y);
428 if(fabs(1.-y)<1.e-3) return f41(x,0.99);
429
430 if(fabs(1.-x/y)<5.e-3) return f41(y*0.98,y);
431
432 return (59.*x*(1.-y)-y*(59.-3.*y))/6./(y-1.)/(x-y)
433 +4.*x*(7.*x*x-3.*x*y+3.*y*y)/3./(x-1.)/(x-y)/(x-y)*log(x)+2.*log(y)*log(y)
434 +4.*y*y*(18.*x-11.*x*y-11.*y+4.*y*y)/3./(y-1.)/(y-1.)/(x-y)/(x-y)*log(y)
435 +4.*(1.+y*y)/(x-1.)/(y-1)*Li2(1.-1./y)+4.*(x*x+y*y)/(x-1.)/(x-y)*Li2(1.-x/y);
436}
437
438/*----------------------------------------------------------------------*/
439
440double f51(double x, double y)
441{
442 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return f51(0.98,1.02);
443
444 if(fabs(1.-x)<1.e-3) return f51(0.99,y);
445 if(fabs(1.-y)<1.e-3) return f51(x,0.99);
446
447 if(fabs(1.-x/y)<5.e-3) return f51(y*0.98,y);
448
449 return (-83.-27.*x*(y-1.)+27.*y)/6./(x-1.)/(y-1.)
450
451 -4.*x*(1.+x*(12.+y)-y-6.*x*x)/3./(x-1.)/(x-1.)/(x-y)*log(x)
452 +2.*(1.+6.*x*x*(y-1.)-3.*x*x*x*(y-1.)+x*(3.*y-4.))/3./(x-1.)/(x-1.)/(y-1.)/(x-y)*log(x)*log(x)
453 -4.*y*(3.*x*x*(y-1.)+x*y*(3.-2.*y)+y*y*(y-2.))/3./(x-1.)/(y-1)/(x-y)/(x-y)*Li2(1.-x/y)
454 -4.*(1.-3.*x-x*x*(3.-6.*y)-x*x*x)/3./(x-1.)/(y-1)/(x-y)*Li2(1.-1./x)
455
456 -4.*y*(1.+y*(12.+x)-x-6.*y*y)/3./(y-1.)/(y-1.)/(y-x)*log(y)
457 +2.*(1.+6.*y*y*(x-1.)-3.*y*y*y*(x-1.)+y*(3.*x-4.))/3./(y-1.)/(y-1.)/(x-1.)/(y-x)*log(y)*log(y)
458 -4.*x*(3.*y*y*(x-1.)+y*x*(3.-2.*x)+x*x*(x-2.))/3./(y-1.)/(x-1)/(y-x)/(y-x)*Li2(1.-y/x)
459 -4.*(1.-3.*y-y*y*(3.-6.*x)-y*y*y)/3./(y-1.)/(x-1)/(y-x)*Li2(1.-1./y)
460 +4.*log(x)*(f40(x,y)+(f40(1.0001*x,y)-f40(0.9999*x,y))/0.0002+(f40(x,1.0001*y)-f40(x,0.9999*y))/0.0002);
461}
462
463/*----------------------------------------------------------------------*/
464
465double f81(double x, double y, double z)
466{
467 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)&&(fabs(1.-z)<1.e-3)) return f81(0.96,0.98,1.02);
468
469 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return f81(0.98,1.02,z);
470 if((fabs(1.-x)<1.e-3)&&(fabs(1.-z)<1.e-3)) return f81(0.98,y,1.02);
471 if((fabs(1.-y)<1.e-3)&&(fabs(1.-z)<1.e-3)) return f81(x,0.98,1.02);
472
473 if(fabs(1.-x)<1.e-3) return f81(0.99,y,z);
474 if(fabs(1.-y)<1.e-3) return f81(x,0.99,z);
475 if(fabs(1.-z)<1.e-3) return f81(x,y,0.99);
476
477 if(fabs(1.-x/y)<5.e-3) return f81(y*0.98,y,z);
478 if(fabs(1.-y/z)<5.e-3) return f81(x,z*0.98,z);
479 if(fabs(1.-x/z)<5.e-3) return f81(x,y,x*0.98);
480
481 return -28.*y*y/3./(y-1.)/(x-y)/(y-z)
482 +4.*x*(7.*x*x-3.*x*y+3.*y*y)/3./(x-1.)/(x-y)/(x-y)/(x-z)*log(x)
483 +4.*z*(7.*z*z-3.*z*y+3.*y*y)/3./(z-1.)/(z-y)/(z-y)/(z-x)*log(z)
484 -4.*y*y*(x*(4.*y*y+18.*z-11.*y*(1.+z))+y*(3.*y*y-11.*z+4.*y*(1.+z)))/3./(y-1.)/(y-1.)/(x-y)/(x-y)/(y-z)/(y-z)*log(y)
485 -4.*(1.+y*y)/(x-1.)/(y-1.)/(z-1.)*Li2(1.-1./y)
486 +4.*(x*x+y*y)/(x-1.)/(x-y)/(x-z)*Li2(1.-x/y)
487 +4.*(z*z+y*y)/(z-1.)/(z-y)/(z-x)*Li2(1.-z/y);
488}
489
490/*----------------------------------------------------------------------*/
491
492double f91(double x, double y, double z)
493{
494 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)&&(fabs(1.-z)<1.e-3)) return f91(0.96,0.98,1.02);
495
496 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return f91(0.98,1.02,z);
497 if((fabs(1.-x)<1.e-3)&&(fabs(1.-z)<1.e-3)) return f91(0.98,y,1.02);
498 if((fabs(1.-y)<1.e-3)&&(fabs(1.-z)<1.e-3)) return f91(x,0.98,1.02);
499
500 if(fabs(1.-x)<1.e-3) return f91(0.99,y,z);
501 if(fabs(1.-y)<1.e-3) return f91(x,0.99,z);
502 if(fabs(1.-z)<1.e-3) return f91(x,y,0.99);
503
504 if(fabs(1.-x/y)<5.e-3) return f91(y*0.98,y,z);
505 if(fabs(1.-y/z)<5.e-3) return f91(x,z*0.98,z);
506 if(fabs(1.-x/z)<5.e-3) return f91(x,y,x*0.98);
507
508 return -28.*y/3./(y-1.)/(x-y)/(y-z)
509 +2.*x*(11.*x+3.*y)/3./(x-1.)/(x-y)/(x-y)/(x-z)*log(x)
510 +2.*z*(11.*z+3.*y)/3./(z-1.)/(z-y)/(z-y)/(z-x)*log(z)
511 +2.*y*(x*(3.*y*y-25.*z+11.*y*(1.+x))+y*(-17.*y*y+11.*z+3.*y*(1.+z)))/3./(y-1.)/(y-1.)/(x-y)/(x-y)/(y-z)/(y-z)*log(y)
512 -4.*(1.+y)/(x-1.)/(y-1.)/(z-1.)*Li2(1.-1./y)
513 +4.*(x+y)/(x-1.)/(x-y)/(x-z)*Li2(1.-x/y)
514 +4.*(z+y)/(z-1.)/(z-y)/(z-x)*Li2(1.-z/y);
515}
516
517/*----------------------------------------------------------------------*/
518
519double f111(double x, double y)
520{
521 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return f111(0.98,1.02);
522
523 if(fabs(1.-x)<1.e-3) return f111(0.99,y);
524 if(fabs(1.-y)<1.e-3) return f111(x,0.99);
525
526 if(fabs(1.-x/y)<5.e-3) return f111(y*0.98,y);
527
528 return 4.*x*(8.*y+(x-1.)*(x-y)*PI*PI)/3./y/(x-1.)/(x-y)
529 -8.*x*(x*x-7.*y+3.*x*(1.+y))/3./(x-y)/(x-y)/(x-1.)/(x-1.)*log(x)
530 -8.*x*(3.*x-7.*y)/3./(x-y)/(x-y)/(y-1.)*log(y)
531 -8.*x/(y-1.)*Li2(1.-1./x)
532 +8.*x/y/(y-1.)*Li2(1.-y/x);
533}
534
535/*----------------------------------------------------------------------*/
536
537double f121(double x, double y, double z)
538{
539 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)&&(fabs(1.-z)<1.e-3)) return f121(0.96,0.98,1.02);
540
541 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return f121(0.98,1.02,z);
542 if((fabs(1.-x)<1.e-3)&&(fabs(1.-z)<1.e-3)) return f121(0.98,y,1.02);
543 if((fabs(1.-y)<1.e-3)&&(fabs(1.-z)<1.e-3)) return f121(x,0.98,1.02);
544
545 if(fabs(1.-x)<1.e-3) return f121(0.99,y,z);
546 if(fabs(1.-y)<1.e-3) return f121(x,0.99,z);
547 if(fabs(1.-z)<1.e-3) return f121(x,y,0.99);
548
549 if(fabs(1.-x/y)<5.e-3) return f121(y*0.98,y,z);
550 if(fabs(1.-y/z)<5.e-3) return f121(x,z*0.98,z);
551 if(fabs(1.-x/z)<5.e-3) return f121(x,y,x*0.98);
552
553 return -28.*y*y/3./(x-y)/(y-1.)/(y-z)
554 +4.*x*x*(6.*x+y)/3./(x-1.)/(x-y)/(x-y)/(x-z)*log(x)
555 +4.*z*z*(6.*z+y)/3./(z-1.)/(z-y)/(z-y)/(z-x)*log(z)
556 -4.*y*y*(x*(6.*y*y+20.*z-13.*y*(1.+z))+y*(y*y-13.*z+6.*y*(1.+z)))/3./(x-y)/(x-y)/(y-1.)/(y-1.)/(y-z)/(y-z)*log(y);
557
558}
559
560/*----------------------------------------------------------------------*/
561
562double f131(double x, double y, double z)
563{
564 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)&&(fabs(1.-z)<1.e-3)) return f131(0.96,0.98,1.02);
565
566 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return f131(0.98,1.02,z);
567 if((fabs(1.-x)<1.e-3)&&(fabs(1.-z)<1.e-3)) return f131(0.98,y,1.02);
568 if((fabs(1.-y)<1.e-3)&&(fabs(1.-z)<1.e-3)) return f131(x,0.98,1.02);
569
570 if(fabs(1.-x)<1.e-3) return f131(0.99,y,z);
571 if(fabs(1.-y)<1.e-3) return f131(x,0.99,z);
572 if(fabs(1.-z)<1.e-3) return f131(x,y,0.99);
573
574 if(fabs(1.-x/y)<5.e-3) return f131(y*0.98,y,z);
575 if(fabs(1.-y/z)<5.e-3) return f131(x,z*0.98,z);
576 if(fabs(1.-x/z)<5.e-3) return f131(x,y,x*0.98);
577
578 return -28.*y/3./(x-y)/(y-1.)/(y-z)
579 +4.*x*(6.*x+y)/3./(x-1.)/(x-y)/(x-y)/(x-z)*log(x)
580 +4.*z*(6.*z+y)/3./(z-1.)/(z-y)/(z-y)/(z-x)*log(z)
581 +4.*y*(x*(y*y-13.*z+6.*y*(1.+z))+y*(y-8.*y*y+6.*z+y*z))/3./(x-y)/(x-y)/(y-1.)/(y-1.)/(y-z)/(y-z)*log(y);
582}
583
584/*----------------------------------------------------------------------*/
585
586double f141(double x, double y)
587{
588 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return f141(0.98,1.02);
589
590 if(fabs(1.-x)<1.e-3) return f141(0.99,y);
591 if(fabs(1.-y)<1.e-3) return f141(x,0.99);
592
593 if(fabs(1.-x/y)<5.e-3) return f141(y*0.98,y);
594
595 return 32.*x*x/3./(x-1.)/(x-y)
596 -8.*x*x*(7.*x*(1.+y)-11.*y-3.*x*x)/3./(x-1.)/(x-1.)/(x-y)/(x-y)*log(x)
597 -8.*x*y*(3.*x-7.*y)/3./(x-y)/(x-y)/(y-1.)*log(y)
598 -8.*x/(y-1.)*Li2(1.-1./x)
599 +8.*x/(y-1.)*Li2(1.-y/x);
600}
601
602/*----------------------------------------------------------------------*/
603
604double f151(double x)
605{
606 if(fabs(1.-x)<1.e-3) return f151(0.99);
607
608 return (1.-3.*x)/(x-1.)+2.*x/(x-1.)/(x-1.)*log(x)+2.*x/(x-1.)*Li2(1.-1./x);
609}
610
611
612/*----------------------------------------------------------------------*/
613
614double f161(double x)
615{
616 if(fabs(1.-x)<1.e-3) return f161(0.99);
617
618 return 28./3./(x-1.)-4.*x*(13.-6.*x)/3./(x-1.)/(x-1.)*log(x);
619}
620
621/*----------------------------------------------------------------------*/
622
623double f171(double x, double y)
624{
625 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return f171(0.98,1.02);
626
627 if(fabs(1.-x)<1.e-3) return f171(0.99,y);
628 if(fabs(1.-y)<1.e-3) return f171(x,0.99);
629
630 if(fabs(1.-x/y)<5.e-3) return f171(y*0.98,y);
631
632 return -28./3./(x-1.)/(y-1.)
633 +4.*y*(10.-3.*y)/3./(x-y)/(y-1.)/(y-1.)*log(y)
634 -4.*y/(x-y)/(y-1.)/(y-1.)*log(y)*log(y)
635 +(4.*(13.*x-6.*x*x-3.*y-7.*x*y+3.*x*x*y)/3./(x-1.)/(x-1.)/(x-y)/(y-1.)+4.*y*log(y)/(x-y)/(y-1.)/(y-1.))*log(x);
636}
637
638/*----------------------------------------------------------------------*/
639
640double f181(double x, double y)
641{
642 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return f181(0.98,1.02);
643
644 if(fabs(1.-x)<1.e-3) return f181(0.99,y);
645 if(fabs(1.-y)<1.e-3) return f181(x,0.99);
646
647 if(fabs(1.-x/y)<5.e-3) return f181(y*0.98,y);
648
649 return -28.*y/3./(x-y)/(y-1.)
650 +4.*x*(6.*x+y)/3./(x-1.)/(x-y)/(x-y)*log(x)
651 -4.*y*(y*(6.+y)-x*(13.-6.*y))/3./(x-y)/(x-y)/(y-1.)/(y-1.)*log(y);
652}
653
654/*----------------------------------------------------------------------*/
655
656double f191(double x, double y)
657{
658 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return f191(0.98,1.02);
659
660 if(fabs(1.-x)<1.e-3) return f191(0.99,y);
661 if(fabs(1.-y)<1.e-3) return f191(x,0.99);
662
663 if(fabs(1.-x/y)<5.e-3) return f191(y*0.98,y);
664
665 return -28.*(x*(y-1.)+y)/3./(x-y)/(y-1.)
666 +4.*x*x*(6.*x+y)/3./(x-1.)/(x-y)/(x-y)*log(x)
667 +4.*y*y*(x*(20.-13.*y)-y*(13.-6.*y))/3./(x-y)/(x-y)/(y-1.)/(y-1.)*log(y);
668}
669
670/*----------------------------------------------------------------------*/
671
672double q11(double x, double y)
673{
674 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return q11(0.98,1.02);
675
676 if(fabs(1.-x)<1.e-3) return q11(0.99,y);
677 if(fabs(1.-y)<1.e-3) return q11(x,0.99);
678
679 if(fabs(1.-x/y)<5.e-3) return q11(y*0.98,y);
680
681 return 4./3./(x-y)*(x*x*log(x)/pow(x-1.,4.)-y*y*log(y)/pow(y-1.,4.)) +(4.*x*x*y*y+10.*x*y*y-2.*y*y+10.*x*x*y-44.*x*y+10.*y-2.*x*x+10.*x+4.)/9./pow(x-1.,3.)/pow(y-1.,3.);
682}
683
684/*----------------------------------------------------------------------*/
685
686double q21(double x, double y)
687{
688 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return q21(0.98,1.02);
689
690 if(fabs(1.-x)<1.e-3) return q21(0.99,y);
691 if(fabs(1.-y)<1.e-3) return q21(x,0.99);
692
693 if(fabs(1.-x/y)<5.e-3) return q21(y*0.98,y);
694
695 return 4./3./(x-y)*(x*log(x)/pow(x-1.,4.)-y*log(y)/pow(y-1.,4.)) +(-2.*x*x*y*y+10.*x*y*y+4.*y*y+10.*x*x*y-20.*x*y-14.*y+4.*x*x-14.*x+22.)/9./pow(x-1.,3.)/pow(y-1.,3.);
696}
697
698/*----------------------------------------------------------------------*/
699
700double q31(double x, double y)
701{
702 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return q31(0.98,1.02);
703
704 if(fabs(1.-x)<1.e-3) return q31(0.99,y);
705 if(fabs(1.-y)<1.e-3) return q31(x,0.99);
706
707 if(fabs(1.-x/y)<5.e-3) return q31(y*0.98,y);
708
709 return 8./3./(x-y)*(-x*x*log(x)/pow(x-1.,3.)+y*y*log(y)/pow(y-1.,3.)) +(-12.*x*y+4.*y+4.*x+4.)/3./pow(x-1.,2.)/pow(y-1.,2.);
710}
711
712/*----------------------------------------------------------------------*/
713
714double q41(double x, double y)
715{
716 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return q41(0.98,1.02);
717
718 if(fabs(1.-x)<1.e-3) return q41(0.99,y);
719 if(fabs(1.-y)<1.e-3) return q41(x,0.99);
720
721 if(fabs(1.-x/y)<5.e-3) return q41(y*0.98,y);
722
723 return 8./3./(x-y)*(-x*log(x)/pow(x-1.,3.)+y*log(y)/pow(y-1.,3.)) +(-4.*x*y-4.*y-4.*x+12.)/3./pow(x-1.,2.)/pow(y-1.,2.);
724}
725
726/*----------------------------------------------------------------------*/
727
728double q51(double x, double y)
729{
730 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return q51(0.98,1.02);
731
732 if(fabs(1.-x)<1.e-3) return q51(0.99,y);
733 if(fabs(1.-y)<1.e-3) return q51(x,0.99);
734
735 if(fabs(1.-x/y)<5.e-3) return q51(y*0.98,y);
736
737 return 4./27./(x-y)*((6.*x*x*x-9.*x*x+2.)*log(x)/pow(x-1.,4.)-(6.*y*y*y-9.*y*y+2.)*log(y)/pow(y-1.,4.))
738 +(104.*x*x*y*y-202.*x*y*y+86.*y*y-202.*x*x*y+380.*x*y-154.*y+86.*x*x-154.*x+56.)/81./pow(x-1.,3.)/pow(y-1.,3.);
739}
740
741/*----------------------------------------------------------------------*/
742
743double q61(double x, double y)
744{
745 if((fabs(1.-x)<1.e-3)&&(fabs(1.-y)<1.e-3)) return q61(0.98,1.02);
746
747 if(fabs(1.-x)<1.e-3) return q61(0.99,y);
748 if(fabs(1.-y)<1.e-3) return q61(x,0.99);
749
750 if(fabs(1.-x/y)<5.e-3) return q61(y*0.98,y);
751
752 return 4./9./(x-y)*(log(x)/pow(x-1.,4.)-log(y)/pow(y-1.,4.))
753 +(4.*x*x*y*y-14.*x*y*y+22.*y*y-14.*x*x*y+52.*x*y-62.*y+22.*x*x-62.*x+52.)/27./pow(x-1.,3.)/pow(y-1.,3.);
754}
755
756/*----------------------------------------------------------------------*/
757
758double A0t(double x)
759{
760 return (-46.+159.*x-153.*x*x+22.*x*x*x)/36./pow(1.-x,3.)+x*x*(-3.*x+2.)/2./pow(1.-x,4.)*log(x);
761}
762
763/*----------------------------------------------------------------------*/
764
765double A1t(double x, double l)
766{
767 return (32.*pow(x,4.)+244.*pow(x,3.)-160.*x*x+16.*x)/9./pow(1.-x,4.)*Li2(1.-1./x)
768 +(-774.*x*x*x*x-2826.*x*x*x+1994.*x*x-130.*x+8.)/81./pow(1.-x,5.)*log(x)
769 +(-94.*x*x*x*x-18665.*x*x*x+20682.*x*x-9113.*x+2006.)/243./pow(1.-x,4.)
770 +((-12.*x*x*x*x-92.*x*x*x+56.*x*x)/3./pow(1.-x,5.)*log(x)+(-68.*x*x*x*x-202.*x*x*x-804.*x*x+794.*x-152.)/27./pow(1.-x,4.))*l;
771}
772
773/*--------------------------------------------------------------------*/
774
775double B0t(double x)
776{
777 return x/4./pow(1.-x,2.)*log(x)+1./4./(1.-x);
778}
779
780/*----------------------------------------------------------------------*/
781
782double C0t(double x)
783{
784 return (3.*x+2.)*x/8./pow(1.-x,2.)*log(x)+(-x+6.)*x/8./(1.-x);
785}
786
787/*----------------------------------------------------------------------*/
788
789double D0t(double x)
790{
791 return (-3.*pow(x,4.)+30.*pow(x,3.)-54.*x*x+32.*x-8.)/18./pow(1.-x,4.)*log(x)+(-47.*pow(x,3.)+237.*x*x-312.*x+104.)/108./pow(1.-x,3.);
792}
793
794/*----------------------------------------------------------------------*/
795
796double B1t(double x, double l)
797{
798 if(fabs(1.-x)<1.e-5) return B1t(0.9999,l);
799
800 return -2.*x/(1.-x)/(1.-x)*Li2(1.-1./x) +(-x+17.)*x/3./pow(1.-x,3.)*log(x) +(13.*x+3)/3./(1.-x)/(1.-x) +((2.*x+2)*x/pow(1.-x,3.)*log(x)+4.*x/(1.-x)/(1.-x))*l;
801}
802
803/*----------------------------------------------------------------------*/
804
805double C1t(double x, double l)
806{
807 return (-x*x-4.)*x/(1.-x)/(1.-x)*Li2(1.-1./x) +(3.*x*x+14.*x+23.)*x/3./pow(1.-x,3.)*log(x) +(4.*x*x+7.*x+29.)*x/3./(1.-x)/(1.-x) +((8.*x+2.)*x/pow(1.-x,3.)*log(x)+(x*x+x+8.)*x/(1.-x)/(1.-x))*l;
808}
809
810/*----------------------------------------------------------------------*/
811
812double D1t(double x, double l)
813{
814 return (380.*pow(x,4.)-1352.*pow(x,3.)+1656.*x*x-784.*x+256.)/81./pow(1.-x,4.)*Li2(1.-1./x) +(304.*pow(x,4.)+1716.*pow(x,3.)-4644.*x*x+2768.*x-720.)/81./pow(1.-x,5.)*log(x) +(-6175.*pow(x,4.)+41608.*pow(x,3.)-66723.*x*x+33106.*x-7000.)/729./pow(1.-x,4.) +((648.*pow(x,4.)-720.*pow(x,3.)-232.*x*x-160.*x+32.)/81./pow(1.-x,5.)*log(x)
815 +(-352.*pow(x,4.)+4912.*pow(x,3.)-8280.*x*x+3304.*x-880.)/243./pow(1.-x,4.))*l;
816}
817
818/*----------------------------------------------------------------------*/
819
820double F0t(double x)
821{
822 if (fabs(1.-x)<1.e-5) return F0t(0.9999);
823 return (5.*x*x*x-9.*x*x+30.*x-8.)/12./pow(1.-x,3.)+3.*x*x/2./pow(1.-x,4.)*log(x);
824}
825
826/*----------------------------------------------------------------------*/
827
828double F1t(double x,double l)
829{
830 return (4.*pow(x,4.)-40.*pow(x,3.)-41.*x*x-x)/3./pow(1.-x,4.)*Li2(1.-1./x)
831 +(-144.*x*x*x*x+3177.*x*x*x+3661.*x*x+250.*x-32.)/108./pow(1.-x,5.)*log(x)
832 +(-247.*x*x*x*x+11890.*x*x*x+31779.*x*x-2966.*x+1016.)/648./pow(1.-x,4.)
833 +((17.*x*x*x+31.*x*x)/pow(1.-x,5.)*log(x)+(-35.*x*x*x*x+170.*x*x*x+447.*x*x+338.*x-56.)/18./pow(1.-x,4.))*l;
834}
835
836/*----------------------------------------------------------------------*/
837
838double E0t(double x)
839{
840 return (-9.*x*x+16.*x-4.)/6./pow(1.-x,4.)*log(x)+(-7.*x*x*x-21.*x*x+42.*x+4.)/36./pow(1.-x,3.);
841}
842
843/*----------------------------------------------------------------------*/
844
845double G1t(double x, double l)
846{
847 return (10.*pow(x,4.)-100.*pow(x,3.)+30.*x*x+160.*x-40.)/27./pow(1.-x,4.)*Li2(1.-1./x)
848 +(30.*pow(x,3.)-42.*x*x-332.*x+68.)/81./pow(1.-x,4.)*log(x)
849 +(-6.*pow(x,3.)-293.*x*x+161.*x+42.)/81./pow(1.-x,3.)
850 +((90.*x*x-160.*x+40.)/27./pow(1.-x,4.)*log(x)+(35.*pow(x,3.)+105.*x*x-210.*x-20.)/81./pow(1.-x,3.))*l;
851}
852
853/*----------------------------------------------------------------------*/
854
855double E1t(double x, double l)
856{
857 return (515.*pow(x,4.)-614.*pow(x,3.)-81.*x*x-190.*x+40.)/54./pow(1.-x,4.)*Li2(1.-1./x)
858 +(-1030.*pow(x,4.)+435.*pow(x,3.)+1373.*x*x+1950.*x-424.)/108./pow(1.-x,5.)*log(x)
859 +(-29467.*pow(x,4.)+45604.*pow(x,3.)-30237.*x*x+66532.*x-10960.)/1944./pow(1.-x,4.)
860 +((-1125.*pow(x,3.)+1685.*x*x+380.*x-76.)/54./pow(1.-x,5.)*log(x)
861 +(133.*pow(x,4.)-2758.*pow(x,3.)-2061.*x*x+11522.*x-1652.)/324./pow(1.-x,4.))*l;
862}
863
864/*----------------------------------------------------------------------*/
865#include <iostream>
866double T(double x)
867{
868 return -(16.*x+8.)*sqrt(4.*x-1.)*Cl2(2.*asin(0.5/sqrt(x)))+(16.*x+20./3.)*log(x)+32.*x+112./9.;
869}
870
871/*----------------------------------------------------------------------*/
872
873double F7_1(double x)
874{
875 return x*(7.-5.*x-8.*x*x)/24./pow(x-1.,3.)+x*x*(3.*x-2.)/4./pow(x-1.,4.)*log(x);
876}
877
878/*----------------------------------------------------------------------*/
879
880double F7_2(double x)
881{
882 return x*(3.-5.*x)/12./pow(x-1.,2.)+x*(3.*x-2.)/6./pow(x-1.,3.)*log(x);
883}
884
885/*----------------------------------------------------------------------*/
886
887double F8_1(double x)
888{
889 return x*(2.+5.*x-x*x)/8./pow(x-1.,3.)-3.*x*x/4./pow(x-1.,4.)*log(x);
890}
891
892/*----------------------------------------------------------------------*/
893
894double F8_2(double x)
895{
896 return x*(3.-x)/4./pow(x-1.,2.)-x/2./pow(x-1.,3.)*log(x);
897}
898
899/*----------------------------------------------------------------------*/
900
901double H2(double x, double y)
902{
903 return D2(x,y);
904}
905
906/*----------------------------------------------------------------------*/
907
908double B(double m1, double m2, double Q)
909{
910 double x=pow(m2/m1,2.);
911
912 if(fabs(x-1.)<1.e-5) return -0.5*log(m2*m2/Q/Q);
913
914 return 0.5*(0.5+1./(1.-x)+log(x)/pow(1.-x,2.)-log(m2*m2/Q/Q));
915}
916
917/*----------------------------------------------------------------------*/
918
919double G7H(double x, double lu, double ld)
920{
921 return lu*(ld*4.*x*(4.*(-3.+7.*x-2.*x*x)*Li2(1.-1./x)+(8.-14.*x-3.*x*x)*log(x)*log(x)/(x-1.)
922 +2.*(-3.-x+12.*x*x-2.*x*x*x)*log(x)/(x-1.)+3.*(7.-13.*x+2.*x*x))/9./pow((x-1.),3.)
923 +lu*2.*x*(x*(18.-37.*x+8.*x*x)*Li2(1.-1./x)+x*(-14.+23.*x+3.*x*x)*log(x)*log(x)/(x-1.)+
924 (-50.+251.*x-174.*x*x-192.*x*x*x+21.*x*x*x*x)*log(x)/9./(x-1.)+(797.-5436.*x+7569.*x*x-1202.*x*x*x)/108.)/pow((x-1.),4.)/9.);
925}
926
927/*----------------------------------------------------------------------*/
928
929double Delta7H(double x, double lu, double ld)
930{
931 return 2.*x/9./pow((x-1.),4.)*lu*
932 (ld*((x-1.)*(21.-47.*x+8.*x*x)+2.*(-8.+14.*x+3.*x*x)*log(x))
933 +lu*((-31.-18.*x+135.*x*x-14.*x*x*x)/6.+x*(14.-23.*x-3.*x*x)*log(x)/(x-1.)));
934}
935
936/*----------------------------------------------------------------------*/
937
938double G8H(double x, double lu, double ld)
939{
940 return lu*(ld*x*(0.5*(-36.+25.*x-17.*x*x)*Li2(1.-1./x)+(19.+17.*x)*log(x)*log(x)/(x-1.)
941 +0.25*(-3.-187.*x+12.*x*x-14.*x*x*x)*log(x)/(x-1.)+3.*(143.-44.*x+29.*x*x)/8.)/3./pow((x-1.),3.)
942 + lu*x*(x*(30.-17.*x+13.*x*x)*Li2(1.-1./x)-x*(31.+17.*x)*log(x)*log(x)/(x-1.)+
943 (-226.+817.*x+1353.*x*x+318.*x*x*x+42.*x*x*x*x)*log(x)/36./(x-1.)+(1130.-18153.*x+7650.*x*x-4451.*x*x*x)/216.)/pow((x-1.),4.)/6.);
944}
945
946/*----------------------------------------------------------------------*/
947
948double Delta8H(double x, double lu, double ld)
949{
950 return (x/6./pow((x-1.),4.))*lu*(ld*((x-1.)*(81.-16.*x+7.*x*x)-2.*(19.+17.*x)*log(x))
951 +lu*((-38.-261.*x+18.*x*x-7.*x*x*x)/6.+x*(31.+17.*x)*log(x)/(x-1.)));
952}
953
954/*----------------------------------------------------------------------*/
955
956double EH(double x, double lu)
957{
958 return lu*lu*x*((x-1.)*(16.-29.*x+7.*x*x)+6.*(3.*x-2.)*log(x))/36./pow((x-1.),4.);
959}
960
961/*----------------------------------------------------------------------*/
962
963double G4H(double x, double lu)
964{
965 return lu*lu*((515.*x*x*x-906.*x*x+99.*x+182.)*x*Li2(1.-1./x)/54./pow(x-1.,4.)
966 +(1030.*x*x*x-2763.*x*x-15.*x+980.)*x*log(x)/108./pow(x-1.,5.)
967 +(-29467.*x*x*x+68142.*x*x-6717.*x-18134.)*x/1944./pow(x-1.,4.));
968}
969
970/*----------------------------------------------------------------------*/
971
972double Delta4H(double x, double lu)
973{
974 return -lu*lu*((-375.*x*x-95.*x+182.)*x*log(x)/54./pow(x-1.,5.)
975 +(133.*x*x*x-108.*x*x+4023.*x-2320.)*x/324./pow(x-1.,4.));
976}
977
978/*----------------------------------------------------------------------*/
979
980double G3H(double x, double lu)
981{
982 return lu*lu*((10.*x*x*x+30.*x-20.)*x*Li2(1.-1./x)/27./pow(x-1.,4.)
983 +(30.*x*x-66.*x-56.)*x*log(x)/81./pow(x-1.,4.)
984 +(6.*x*x-187.*x+213.)*x/81./pow(x-1.,3.));
985}
986
987/*----------------------------------------------------------------------*/
988
989double Delta3H(double x, double lu)
990{
991 return -lu*lu*((-30.*x+20.)*x*log(x)/27./pow(x-1.,4.)
992 +(-35.*x*x+145.*x-80.)*x/81./pow(x-1.,3.));
993}
994
995/*----------------------------------------------------------------------*/
996
997double C9llH0(double x, double y, double lu)
998{
999 if(fabs(1.-y)<1.e-5) {
1000 return C9llH0(x,0.9999,lu);
1001 }
1002 return x/y/8.*lu*lu*(-log(y)/(y-1.)+1.)*y*y/(y-1.);
1003}
1004
1005/*----------------------------------------------------------------------*/
1006
1007double D9H0(double x, double lu)
1008{
1009 if(fabs(1.-x)<1.e-5) {
1010 return D9H0(0.9999,lu);
1011 }
1012 return lu*lu*((-3.*x*x*x+6.*x-4.)*x/18./pow(x-1.,4.)*log(x)+(47.*x*x-79.*x+38.)*x/108./pow(x-1.,3.));
1013}
1014
1015/*----------------------------------------------------------------------*/
1016
1017double C9llH1(double x, double y, double lu, double L)
1018{
1019 if(fabs(y-1.)<1.e-5) return C9llH1(x,0.9999,lu,L);
1020
1021 return x/y/8.*lu*lu*((-8.*y*y*y+16.*y*y)/pow(y-1.,2.)*Li2(1.-1./y)
1022 +(-24.*y*y*y+88.*y*y)/3./pow(y-1.,3.)*log(y)
1023 +(32.*y*y*y-96.*y*y)/3./(y-1.)/(y-1.)
1024 +(16.*y*y*log(y)/pow(y-1.,3.)+(8.*y*y*y-24.*y*y)/(y-1.)/(y-1.))*L);
1025}
1026
1027/*----------------------------------------------------------------------*/
1028
1029double D9H1(double x, double lu, double L)
1030{
1031 if(fabs(x-1.)<1.e-5) return D9H1(0.9999,lu,L);
1032
1033 return lu*lu*((380.*x*x*x-528.*x*x+72.*x+128.)*x/81.*pow(1.-x,4.)*Li2(1.-1./x)
1034 +(596.*x*x*x-672.*x*x+64.*x+204.)*x*log(x)/81./pow(x-1.,5.)
1035 +(-6175.*x*x*x+9138.*x*x-3927.*x-764.)*x/729./pow(x-1.,4.)
1036 +((432.*x*x*x-456.*x*x+40.*x+128.)*x*log(x)/81./pow(x-1.,5.)
1037 +(-352.*x*x*x-972.*x*x+1944.*x-1052.)*x/243./pow(x-1.,4.))*L);
1038}
1039
1040/*----------------------------------------------------------------------*/
1041
1042double C7t2mt(double x)
1043{
1044 double z=1./x;
1045 double w=1.-z;
1046 double y=sqrt(z);
1047
1048 if(y<0.4) return 12.06+12.93*z+3.013*z*log(z)+96.71*z*z+52.73*z*z*log(z)+147.9*pow(z,3.)
1049 +187.7*pow(z,3.)*log(z)-144.9*pow(z,4.)+236.1*pow(z,4.)*log(z);
1050
1051 else return 11.74+0.3642*w+0.1155*w*w-0.003145*pow(w,3.)-0.03263*pow(w,4.)-0.03528*pow(w,5.)
1052 -0.03076*pow(w,6.)-0.02504*pow(w,7.)-0.01985*pow(w,8.);
1053}
1054
1055/*----------------------------------------------------------------------*/
1056
1057double C7c2MW(double x)
1058{
1059 double z=1./x;
1060 double w=1.-z;
1061 double y=sqrt(z);
1062
1063 if(y<0.4) return 1.525-0.1165*z+0.01975*z*log(z)+0.06283*z*z+0.005349*z*z*log(z)
1064 +0.01005*pow(z*log(z),2.)-0.04202*pow(z,3.)+0.01535*pow(z,3.)*log(z)-0.00329*z*pow(z*log(z),2.)
1065 +0.002372*pow(z,4.)-0.0007910*pow(z,4.)*log(z);
1066
1067 else return 1.432+0.06709*w+0.01257*w*w+0.004710*pow(w,3.)+0.002373*pow(w,4.)
1068 +0.001406*pow(w,5.)+0.0009216*pow(w,6.)+0.00064730*pow(w,7.)+0.0004779*pow(w,8.);
1069}
1070
1071/*----------------------------------------------------------------------*/
1072
1073double C8t2mt(double x)
1074{
1075 double z=1./x;
1076 double w=1.-z;
1077 double y=sqrt(z);
1078
1079 if(y<0.35) return -0.8954-7.043*z-98.34*z*z-46.21*z*z*log(z)-127.1*pow(z,3.)
1080 -181.6*pow(z,3.)*log(z)+535.8*pow(z,4.)-76.76*pow(z,4.)*log(z);
1081
1082 else return -0.614093-0.897491*w-0.0349242*w*w+0.0679089*pow(w,3.)+0.0796552*pow(w,4.)
1083 +0.0722587 * pow(w,5.0) + 0.0613235*pow(w,6.)+0.0509632*pow(w,7.)+0.0421555*pow(w,8.) + 0.0349448 * pow(w,9.)
1084 +0.0291171 * pow(w,10.) + 0.0244182 * pow(w,11.) + 0.0206190 * pow(w,12.) + 0.0175317 * pow(w, 13)
1085 +0.0150070 * pow(w,14.) + 0.0129285 * pow(w,15.) + 0.0112057 * pow(w,16.);
1086}
1087
1088/*----------------------------------------------------------------------*/
1089
1090double C8c2MW(double x)
1091{
1092 double z=1./x;
1093 double w=1.-z;
1094 double y=sqrt(z);
1095
1096 if(y<0.35) return -1.870+0.1010*z-0.1218*z*log(z)+0.1045*z*z-0.03748*z*z*log(z)
1097 +0.01151*pow(z*log(z),2.)-0.01023*pow(z,3.)+0.004342*pow(z,3.)*log(z)+0.0003031*z*pow(z*log(z),2.)
1098 -0.001537*pow(z,4.)+0.0007532*pow(z,4.)*log(z);
1099
1100 else return -1.676-0.1179*w-0.02926*w*w-0.01297*pow(w,3.)-0.007296*pow(w,4.)
1101 -0.004672*pow(w,5.)-0.003248*pow(w,6.)-0.002389*pow(w,7.)-0.001831*pow(w,8.);
1102}
1103
1104double C10Wt2mt(double x)
1105{
1106 double z=1./x;
1107 double w=1.-z;
1108 double y=sqrt(z);
1109
1110 if(y<0.35) return 2.7096*z-8.15633*z*z-0.539366*pow(z,3.)+35.3187*pow(z,4.)+103.909*pow(z,5.)+207.722*pow(z,6.)+log(y)*(6.00974*z-1.13112*z*z-13.9722*pow(z,3.)+15.6381*pow(z,4.)+149.164*pow(z,5.)+ 454.776*pow(z,6.));
1111
1112 else return -0.449493-0.584517*w+0.132958*w*w+0.156329*pow(w,3.)+0.123277*pow(w,4.)+0.0933327*pow(w,5.)+0.0713398*pow(w,6.)+0.0556133*pow(w,7.)+0.044254*pow(w,8.)+0.0358898*pow(w,9.)+0.0296021*pow(w,10.)+0.0247807*pow(w,11.)+0.0210161*pow(w,12.)+0.018028*pow(w,13.)+0.0156213*pow(w,14.)+0.0136572*pow(w,15.)+0.0120353*pow(w,16.);
1113}
1114
1115
1116double C10Wc2MW(double x)
1117{
1118 double z=1./x;
1119 double w=1.-z;
1120 double y=sqrt(z);
1121
1122 if(y<0.35) return -5.22222-0.221454*z+0.04146*z*z-0.00109247*pow(z,3.)-0.0000428586*pow(z,4.)-3.10949e-6*pow(z,5.)-3.00894e-7*pow(z,6.)+log(y)*(0.124444*z-0.0295465*z*z+0.000634921*pow(z,3.)+0.0000320667*pow(z,4.)+2.64286e-6*pow(z,5.)+2.775e-7*pow(z,6.)+log(y)*(-0.0888889*z+0.00952381*z*z));
1123
1124 else return -5.40335+0.0942163*w+0.0278556*w*w+0.0135525*pow(w,3.)+0.00812898*pow(w,4.)+0.00546931*pow(w,5.)+0.00395656*pow(w,6.)+0.0030088*pow(w,7.)+0.00237313*pow(w,8.)+0.00192462*pow(w,9.)+0.00159554*pow(w,10.)+0.00134643*pow(w,11.)+0.00115299*pow(w,12.)+0.000999582*pow(w,13.)+0.000875719*pow(w,14.)+0.000774173*pow(w,15.)+0.000689814*pow(w,16.);
1125}
1126
1127double C10Zt2mt(double x)
1128{
1129 double z=1./x;
1130 double w=1.-z;
1131 double y=sqrt(z);
1132
1133 if(y<0.35) return 2.13944+0.189707/z+28.593*z+28.0108*z*z-31.409*pow(z,3.)-166.975*pow(z,4.)-387.429*pow(z,5.)-697.854*pow(z,6.)+log(y)*(33.8492*z+97.9845*z*z+106.223*pow(z,3.)-78.5852*pow(z,4.)-618.344*pow(z,5.)-1688.32*pow(z,6.));
1134
1135 else return -1.93409+0.896643*w+0.739905*w*w+0.605819*pow(w,3.)+0.511284*pow(w,4.)+0.443902*pow(w,5.)+0.394799*pow(w,6.)+0.358192*pow(w,7.)+0.330311*pow(w,8.)+0.30866*pow(w,9.)+0.291553*pow(w,10.)+0.277824*pow(w,11.)+0.266654*pow(w,12.)+0.257453*pow(w,13.)+0.249788*pow(w,14.)+0.243341*pow(w,15.)+0.237868*pow(w,16.);
1136}
1137
1138double C10Z2tri(double x)
1139{
1140 double z=1./x;
1141 double w=1.-z;
1142 double y=sqrt(z);
1143
1144 if(y<0.35) return -1.13764-0.987062/z-1.09363*z-1.94448*z*z-2.0342*pow(z,3.)-2.20973*pow(z,4.)-2.35256*pow(z,5.)-2.47336*pow(z,6.)+log(y)*(-1.5-3.79364*z-6.87688*z*z-10.8345*pow(z,3.)-15.0857*pow(z,4.)-19.6523*pow(z,5.)-24.4813*pow(z,6.)+log(y)*(0.222222*z-0.0277778*z*z));
1145
1146 else return -0.746008-0.800433*w-0.836781*w*w-0.861648*pow(w,3.)-0.879561*pow(w,4.)-0.893037*pow(w,5.)- 0.903527*pow(w,6.)-0.911919*pow(w,7.)-0.918783*pow(w,8.)-0.924499*pow(w,9.)-0.929333*pow(w,10.)-0.933474*pow(w,11.)-0.937061*pow(w,12.)-0.940198*pow(w,13.)-0.942964*pow(w,14.)-0.945422*pow(w,15.)-0.947621*pow(w,16.);
1147}
1148
1149
1150
1151double F0SP(double xt)
1152{
1153 if(fabs(1.-xt)<1.e-5) return F0SP(0.9999);
1154
1155 return 1./8./(xt-1.)/(xt-1.)*((xt-3.)/2.-xt*(xt-2.)/(xt-1.)*log(xt));
1156}
1157
1158
1159double S0(double x) {
1160 if(abs(1.-x)<1.e-5) return S0(0.9999);
1161
1162 return (4.*x-11.*x*x+x*x*x)/4./pow(1.-x,2.) - 3.*x*x*x*log(x)/2./pow(1.-x,3.);
1163}
1164
1165double D0(double w, double x, double y, double z)
1166{
1167 if((fabs(1.-w)<1.e-5)&&(fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)&&(fabs(1.-z)<1.e-5)) return D0(0.9996,0.9998,1.0002,1.0004);
1168 if((fabs(1.-w)<1.e-5)&&(fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)) return D0(0.9996,0.9998,1.0002,z);
1169 if((fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)&&(fabs(1.-z)<1.e-5)) return D0(w,0.9996,0.9998,1.0002);
1170 if((fabs(1.-w)<1.e-5)&&(fabs(1.-y)<1.e-5)&&(fabs(1.-z)<1.e-5)) return D0(0.9996,x,0.9998,1.0002);
1171 if((fabs(1.-w)<1.e-5)&&(fabs(1.-x)<1.e-5)&&(fabs(1.-z)<1.e-5)) return D0(0.9996,0.9998,y,1.0002);
1172 if((fabs(1.-w)<1.e-5)&&(fabs(1.-x)<1.e-5)) return D0(0.9998,1.0002,y,z);
1173 if((fabs(1.-w)<1.e-5)&&(fabs(1.-y)<1.e-5)) return D0(0.9998,x,1.0002,z);
1174 if((fabs(1.-w)<1.e-5)&&(fabs(1.-z)<1.e-5)) return D0(0.9998,x,y,1.0002);
1175 if((fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)) return D0(w,0.9998,1.0002,z);
1176 if((fabs(1.-x)<1.e-5)&&(fabs(1.-z)<1.e-5)) return D0(w,0.9998,y,1.0002);
1177 if((fabs(1.-y)<1.e-5)&&(fabs(1.-z)<1.e-5)) return D0(w,x,0.9998,1.0002);
1178 if(fabs(1.-w)<1.e-5) return D0(0.9999,x,y,z);
1179 if(fabs(1.-x)<1.e-5) return D0(w,0.9999,y,z);
1180 if(fabs(1.-y)<1.e-5) return D0(w,x,0.9999,z);
1181 if(fabs(1.-z)<1.e-5) return D0(w,x,y,0.9999);
1182 if(fabs(1.-w/x)<1.e-5) return D0(x*0.9998,x,y,z);
1183 if(fabs(1.-w/y)<1.e-5) return D0(y*0.9998,x,y,z);
1184 if(fabs(1.-w/z)<1.e-5) return D0(z*0.9998,x,y,z);
1185 if(fabs(1.-x/y)<1.e-5) return D0(w,y*0.9998,y,z);
1186 if(fabs(1.-y/z)<1.e-5) return D0(w,x,z*0.9998,z);
1187 if(fabs(1.-x/z)<1.e-5) return D0(w,x,y,x*0.9998);
1188
1189 return (w*log(w)/((z-w)*(y-w)*(x-w)))+(x*log(x)/((z-x)*(y-x)*(w-x)))+(y*log(y)/((z-y)*(w-y)*(x-y)))+(z*log(z)/((w-z)*(y-z)*(x-z)));
1190}
1191
1192double D2p(double w, double x, double y, double z)
1193{
1194 if((fabs(1.-w)<1.e-5)&&(fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)&&(fabs(1.-z)<1.e-5)) return D2p(0.9996,0.9998,1.0002,1.0004);
1195 if((fabs(1.-w)<1.e-5)&&(fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)) return D2p(0.9996,0.9998,1.0002,z);
1196 if((fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)&&(fabs(1.-z)<1.e-5)) return D2p(w,0.9996,0.9998,1.0002);
1197 if((fabs(1.-w)<1.e-5)&&(fabs(1.-y)<1.e-5)&&(fabs(1.-z)<1.e-5)) return D2p(0.9996,x,0.9998,1.0002);
1198 if((fabs(1.-w)<1.e-5)&&(fabs(1.-x)<1.e-5)&&(fabs(1.-z)<1.e-5)) return D2p(0.9996,0.9998,y,1.0002);
1199 if((fabs(1.-w)<1.e-5)&&(fabs(1.-x)<1.e-5)) return D2p(0.9998,1.0002,y,z);
1200 if((fabs(1.-w)<1.e-5)&&(fabs(1.-y)<1.e-5)) return D2p(0.9998,x,1.0002,z);
1201 if((fabs(1.-w)<1.e-5)&&(fabs(1.-z)<1.e-5)) return D2p(0.9998,x,y,1.0002);
1202 if((fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)) return D2p(w,0.9998,1.0002,z);
1203 if((fabs(1.-x)<1.e-5)&&(fabs(1.-z)<1.e-5)) return D2p(w,0.9998,y,1.0002);
1204 if((fabs(1.-y)<1.e-5)&&(fabs(1.-z)<1.e-5)) return D2p(w,x,0.9998,1.0002);
1205 if(fabs(1.-w)<1.e-5) return D2p(0.9999,x,y,z);
1206 if(fabs(1.-x)<1.e-5) return D2p(w,0.9999,y,z);
1207 if(fabs(1.-y)<1.e-5) return D2p(w,x,0.9999,z);
1208 if(fabs(1.-z)<1.e-5) return D2p(w,x,y,0.9999);
1209 if(fabs(1.-w/x)<1.e-5) return D2p(x*0.9998,x,y,z);
1210 if(fabs(1.-w/y)<1.e-5) return D2p(y*0.9998,x,y,z);
1211 if(fabs(1.-w/z)<1.e-5) return D2p(z*0.9998,x,y,z);
1212 if(fabs(1.-x/y)<1.e-5) return D2p(w,y*0.9998,y,z);
1213 if(fabs(1.-y/z)<1.e-5) return D2p(w,x,z*0.9998,z);
1214 if(fabs(1.-x/z)<1.e-5) return D2p(w,x,y,x*0.9998);
1215
1216 return 0.25*((w*w*log(w)/((z-w)*(y-w)*(x-w)))+(x*x*log(x)/((z-x)*(y-x)*(w-x)))+(y*y*log(y)/((z-y)*(w-y)*(x-y)))+(z*z*log(z)/((w-z)*(y-z)*(x-z))));
1217}
1218
1219void getDelta(scalar_t delta[6][6],scalar_t Z[6][6],double M[6],double m_av,scalar_t delta_LL[3][3],scalar_t delta_LR[3][3],scalar_t delta_RL[3][3],scalar_t delta_RR[3][3]) {
1220 scalar_t mass_matrix[6][6];
1221
1222
1223 for(int i=0; i<6; ++i) for(int j=0; j<6; ++j)
1224 {
1225 if (i==j) mass_matrix[i][j]=pow(M[i],2.);
1226 else mass_matrix[i][j]=0.;
1227 }
1228
1229 scalar_t ZMZdag[6][6]; /* NM: ZMZdag = Z*M*Z^dagger */
1230 for(int i=0; i<6; ++i) for(int j=0; j<6; ++j)
1231 {
1232 ZMZdag[i][j]=0.;
1233 for(int k=0; k<6; ++k) for(int l=0; l<6; ++l) ZMZdag[i][j]+=Z[i][k]*mass_matrix[k][l]*conj(Z[j][l]); /* NM: matrix product routine removed and replaced locally */
1234 }
1235
1236 for(int i=0; i<6; ++i) for(int j=0; j<6; ++j) delta[i][j] = (ZMZdag[i][j]-pow(m_av,2)*(i==j))/(pow(m_av,2));
1237
1238 for(int i=0; i<6; ++i) for(int j=0; j<6; ++j)
1239 {
1240 delta_LL[i][j]=delta[i][j];
1241 delta_LR[i][j]=delta[i][j+3];
1242 delta_RL[i][j]=delta[i+3][j];
1243 delta_RR[i][j]=delta[i+3][j+3];
1244 }
1245
1246 return;
1247}
1248
1249double h3(double x)
1250{
1251 if(fabs(x-1.)<1.e-5) return -1./4;
1252 return -0.5/(1.-x)-0.5*x*log(x)*pow(1.-x,-2.);
1253}
1254
1255double h1(double x)
1256{
1257 if(fabs(x-1.)<1.e-5) return 4./9;
1258 return 4.*(1.+x)/(3.*pow(1.-x,2.)) + 8.*x*log(x)/(3.*pow((1.-x),3.));
1259}
1260
1261double h4(double x, double y)
1262{
1263 if((fabs(1.-x)<1.e-5)&&(fabs(1.-y)<1.e-5)) return 1./6;
1264 if(fabs(1.-x)<1.e-5) return h4(0.9999,y);
1265 if(fabs(1.-y)<1.e-5) return h4(x,0.9999);
1266 if(fabs(1.-x/y)<1.e-5) return h4(y*0.9998,y);
1267
1268 return -1./((1.-x)*(1.-y))+x*log(x)/((y-x)*pow((1.-x),2.))+y*log(y)/((x-y)*pow((1.-y),2.));
1269}
1270
1271double f(double x)
1272{
1273 if(fabs(x-1.)<1.e-5) return 0.5;
1274 return 1./(1.-x)+x*log(x)/pow(1.-x,2.);
1275}
1276
1277
1278/* hep-ph/9901288 */
1279double Y0(double xt)
1280{
1281 return xt/8. * ((4.-xt)/(1.-xt) + (3.*xt)/pow(1.-xt,2.) * log(xt));
1282}
1283
1284/* hep-ph/9901288 */
1285double Y1(double xt, double mu, double mass_W)
1286{
1287 return (10.*xt + 10*xt*xt + 4.*pow(xt,3.))/(3.*pow(1.-xt,2.)) - (2.*xt - 8.*xt*xt - pow(xt,3.) - pow(xt,4.))/(pow(1.-xt,3.)) * log(xt)
1288 + (2.*xt - 14.*xt*xt + xt*xt*xt - pow(xt,4.))/(2*pow(1.-xt,3)) * log(xt)*log(xt) + (2.*xt + pow(xt,3.))/(pow(1.-xt,2.)) * Li2(1.-xt)
1289 + 8.*xt*log(mu*mu/(mass_W*mass_W))*( (-4. + 3.*xt + xt*xt*xt -6.*xt*log(xt))/(8.*pow(xt-1.,3.)) );
1290}
1291
1292double X0(double xt)
1293{
1294 return xt/8. * ((xt + 2.)/(xt - 1.) + (3.*xt - 6.)/pow(xt - 1.,2.) * log(xt));
1295}
1296
1297/*----------------------------------------------------------------------------*/
1298
1299double X1(double xt, double mu, double mass_W)
1300{
1301 return -(29.*xt - xt*xt - 4.*pow(xt,3.))/(3.*pow(1.-xt,2.)) - (xt + 9.*xt*xt - pow(xt,3.) - pow(xt,4.))/(pow(1.-xt,3.)) * log(xt)
1302 + (8.*xt + 4.*xt*xt + xt*xt*xt - pow(xt,4.))/(2*pow(1.-xt,3)) * log(xt)*log(xt) - (4.*xt - pow(xt,3.))/(pow(1.-xt,2.)) * Li2(1.-xt)
1303 + 8.*xt*( (8. - 9.*xt + pow(xt,3.) +6.*log(xt))/(8.*pow(-1+xt,3.)) )*log(mu*mu/(mass_W*mass_W));
1304}
1305
1306double psi(int n) {
1307 double sum = 0;
1308 for (int k = 1; k < n; k++) {
1309 sum += 1. / k;
1310 }
1311 return sum - GAMMA;
1312}
scalar_t
Definition scalar.h:11
PI
constexpr double PI
Definition constants.h:7
GAMMA
constexpr double GAMMA
Definition constants.h:17
Cl2
double Cl2(double x)
Computes the Clausen function Cl2(x).
Definition polylog.cpp:1333
Li2
double Li2(double x)
Computes the dilogarithm function Li2(x).
Definition polylog.cpp:1257
pow
scalar_t pow(const scalar_t &base, const scalar_t &exp)
Definition scalar.cpp:75
F7_2
double F7_2(double x)
Wilson coefficient F7_2.
Definition special_SM.cpp:880
B1t
double B1t(double x, double l)
Wilson coefficient B1t depending on x and scale l.
Definition special_SM.cpp:796
Delta3H
double Delta3H(double x, double lu)
Wilson coefficient Delta3H depending on x and lu.
Definition special_SM.cpp:989
h51
double h51(double x, double y)
Definition special_SM.cpp:356
q61
double q61(double x, double y)
Definition special_SM.cpp:743
C8c2MW
double C8c2MW(double x)
Wilson coefficient C8c at M_W scale.
Definition special_SM.cpp:1090
T
double T(double x)
Wilson coefficient T(x).
Definition special_SM.cpp:866
q41
double q41(double x, double y)
Definition special_SM.cpp:714
D2
double D2(double x, double y)
Definition special_SM.cpp:28
f41
double f41(double x, double y)
Definition special_SM.cpp:423
D9H1
double D9H1(double x, double lu, double L)
Wilson coefficient D9H1 depending on x, lu, and L.
Definition special_SM.cpp:1029
q11
double q11(double x, double y)
Definition special_SM.cpp:672
C10Wt2mt
double C10Wt2mt(double x)
Wilson coefficient C10Wt at 2m_t scale.
Definition special_SM.cpp:1104
f141
double f141(double x, double y)
Definition special_SM.cpp:586
S0
double S0(double x)
Wilson special function S0 depending on xt.
Definition special_SM.cpp:1159
G1t
double G1t(double x, double l)
Wilson coefficient G1t depending on x and scale l.
Definition special_SM.cpp:845
f90
double f90(double w, double x, double y, double z)
Definition special_SM.cpp:204
F7_1
double F7_1(double x)
Wilson coefficient F7_1.
Definition special_SM.cpp:873
f110
double f110(double x, double y)
Definition special_SM.cpp:273
H2
double H2(double x, double y)
Computes the two-variable function H2(x, y).
Definition special_SM.cpp:901
f81
double f81(double x, double y, double z)
Definition special_SM.cpp:465
f40
double f40(double x, double y)
Definition special_SM.cpp:122
h40
double h40(double x)
Definition special_SM.cpp:72
A1t
double A1t(double x, double l)
Wilson coefficient A1t depending on x and scale l.
Definition special_SM.cpp:765
C0t
double C0t(double x)
Scalar one-loop Wilson coefficient C0t.
Definition special_SM.cpp:782
X0
double X0(double xt)
Definition special_SM.cpp:1292
h11
double h11(double x, double y)
Definition special_SM.cpp:287
h41
double h41(double x, double y)
Definition special_SM.cpp:339
f
double f(double x)
Wilson special function f depending on x.
Definition special_SM.cpp:1271
getDelta
void getDelta(scalar_t delta[6][6], scalar_t Z[6][6], double M[6], double m_av, scalar_t delta_LL[3][3], scalar_t delta_LR[3][3], scalar_t delta_RL[3][3], scalar_t delta_RR[3][3])
Wilson special function getDelta.
Definition special_SM.cpp:1219
G7H
double G7H(double x, double lu, double ld)
Wilson coefficient G7H depending on x, lu, and ld.
Definition special_SM.cpp:919
f111
double f111(double x, double y)
Definition special_SM.cpp:519
q21
double q21(double x, double y)
Definition special_SM.cpp:686
f30
double f30(double x, double y)
Definition special_SM.cpp:108
D2p
double D2p(double w, double x, double y, double z)
Wilson special function D2p depending on 4 parameters.
Definition special_SM.cpp:1192
q51
double q51(double x, double y)
Definition special_SM.cpp:728
h4
double h4(double x, double y)
Wilson special function h4 depending on x and y.
Definition special_SM.cpp:1261
D1t
double D1t(double x, double l)
Wilson coefficient D1t depending on x and scale l.
Definition special_SM.cpp:812
f70
double f70(double x, double y)
Definition special_SM.cpp:181
h1
double h1(double x)
Wilson special function h1 depending on x.
Definition special_SM.cpp:1255
Delta7H
double Delta7H(double x, double lu, double ld)
Wilson coefficient Delta7H depending on x, lu, and ld.
Definition special_SM.cpp:929
C1t
double C1t(double x, double l)
Wilson coefficient C1t depending on x and scale l.
Definition special_SM.cpp:805
f151
double f151(double x)
Definition special_SM.cpp:604
EH
double EH(double x, double lu)
Wilson coefficient EH depending on x and lu.
Definition special_SM.cpp:956
h3
double h3(double x)
Wilson special function h3 depending on x.
Definition special_SM.cpp:1249
C7c2MW
double C7c2MW(double x)
Wilson coefficient C7c at M_W scale.
Definition special_SM.cpp:1057
f60
double f60(double x, double y, double z)
Definition special_SM.cpp:159
E0t
double E0t(double x)
Scalar one-loop Wilson coefficient E0t.
Definition special_SM.cpp:838
C8t2mt
double C8t2mt(double x)
Wilson coefficient C8t at 2m_t scale.
Definition special_SM.cpp:1073
A0t
double A0t(double x)
Scalar one-loop Wilson coefficient A0t.
Definition special_SM.cpp:758
C10Z2tri
double C10Z2tri(double x)
Wilson coefficient C10Z from triangle diagrams.
Definition special_SM.cpp:1138
Delta8H
double Delta8H(double x, double lu, double ld)
Wilson coefficient Delta8H depending on x, lu, and ld.
Definition special_SM.cpp:948
h60
double h60(double x)
Definition special_SM.cpp:90
h30
double h30(double x)
Definition special_SM.cpp:62
Y1
double Y1(double xt, double mu, double mass_W)
Definition special_SM.cpp:1285
G3H
double G3H(double x, double lu)
Wilson coefficient G3H depending on x and lu.
Definition special_SM.cpp:980
F8_1
double F8_1(double x)
Wilson coefficient F8_1.
Definition special_SM.cpp:887
f181
double f181(double x, double y)
Definition special_SM.cpp:640
Bplus
double Bplus(double x, double y)
Definition special_SM.cpp:9
f121
double f121(double x, double y, double z)
Definition special_SM.cpp:537
h50
double h50(double x)
Definition special_SM.cpp:81
D9H0
double D9H0(double x, double lu)
Wilson coefficient D9H0 depending on x and lu.
Definition special_SM.cpp:1007
h71
double h71(double x, double y)
Definition special_SM.cpp:390
f100
double f100(double w, double x, double y, double z)
Definition special_SM.cpp:239
h21
double h21(double x, double y)
Definition special_SM.cpp:304
D0
double D0(double w, double x, double y, double z)
Wilson special function D0 depending on 4 parameters.
Definition special_SM.cpp:1165
q31
double q31(double x, double y)
Definition special_SM.cpp:700
B0t
double B0t(double x)
Scalar one-loop Wilson coefficient B0t.
Definition special_SM.cpp:775
h10
double h10(double x)
Definition special_SM.cpp:42
C9llH0
double C9llH0(double x, double y, double lu)
Wilson coefficient C9llH0 depending on x, y, and lu.
Definition special_SM.cpp:997
f51
double f51(double x, double y)
Definition special_SM.cpp:440
Delta4H
double Delta4H(double x, double lu)
Wilson coefficient Delta4H depending on x and lu.
Definition special_SM.cpp:972
f80
double f80(double x)
Definition special_SM.cpp:195
X1
double X1(double xt, double mu, double mass_W)
Definition special_SM.cpp:1299
f161
double f161(double x)
Definition special_SM.cpp:614
F0SP
double F0SP(double xt)
Wilson coefficient F0SP depending on xt.
Definition special_SM.cpp:1151
F0t
double F0t(double x)
Scalar one-loop Wilson coefficient F0t.
Definition special_SM.cpp:820
D0t
double D0t(double x)
Scalar one-loop Wilson coefficient D0t.
Definition special_SM.cpp:789
f50
double f50(double x, double y, double z)
Definition special_SM.cpp:136
D3
double D3(double x)
Definition special_SM.cpp:17
psi
double psi(int n)
Definition special_SM.cpp:1306
f31
double f31(double x, double y)
Definition special_SM.cpp:407
Y0
double Y0(double xt)
Definition special_SM.cpp:1279
h20
double h20(double x)
Definition special_SM.cpp:52
F1t
double F1t(double x, double l)
Wilson coefficient F1t depending on x and scale l.
Definition special_SM.cpp:828
f20
double f20(double x)
Definition special_SM.cpp:99
C9llH1
double C9llH1(double x, double y, double lu, double L)
Wilson coefficient C9llH1 depending on x, y, lu, and L.
Definition special_SM.cpp:1017
f131
double f131(double x, double y, double z)
Definition special_SM.cpp:562
G4H
double G4H(double x, double lu)
Wilson coefficient G4H depending on x and lu.
Definition special_SM.cpp:963
f191
double f191(double x, double y)
Definition special_SM.cpp:656
G8H
double G8H(double x, double lu, double ld)
Wilson coefficient G8H depending on x, lu, and ld.
Definition special_SM.cpp:938
C10Zt2mt
double C10Zt2mt(double x)
Wilson coefficient C10Zt at 2m_t scale.
Definition special_SM.cpp:1127
f91
double f91(double x, double y, double z)
Definition special_SM.cpp:492
h31
double h31(double x, double y)
Definition special_SM.cpp:321
f171
double f171(double x, double y)
Definition special_SM.cpp:623
F8_2
double F8_2(double x)
Wilson coefficient F8_2.
Definition special_SM.cpp:894
h61
double h61(double x, double y)
Definition special_SM.cpp:373
C10Wc2MW
double C10Wc2MW(double x)
Wilson coefficient C10Wc at M_W scale.
Definition special_SM.cpp:1116
C7t2mt
double C7t2mt(double x)
Wilson coefficient C7t at 2m_t scale.
Definition special_SM.cpp:1042
E1t
double E1t(double x, double l)
Wilson coefficient E1t depending on x and scale l.
Definition special_SM.cpp:855
special_SM.h
special_SUSY.h
special_THDM.h