#include "dilog.h"
#include "Wrappers.H"
#include "ATOOLS/Math/MathTools.H"
using namespace ATOOLS;
// routines only valid for real values of the argument;
// imaginary parts not presently given.
// this version uses 't Hooft and Veltman's change of variable
// good to ~ 10^(-16)
double li2(double x){
double x_0 = -0.30;
double x_1 = 0.25;
double x_2 = 0.51;
if (x == 1.) return PISQ6;
if (x <= x_0){
double temp = log(Abs(1.0-x));
return -li2(-x/(1.0-x)) - temp*temp/2 ; }
else if (x < x_1){
double z = - log(1.0-x);
double temp = z*(1.0-z/4.0*(1.0-z/9.0*(1.0-z*z/100.0
*(1.0-5.0*z*z/294.0*(1.0-7.0*z*z/360.0
*(1.0-5.0*z*z/242.0*(1.0-7601.0*z*z/354900.0
*(1.0-91.0*z*z/4146.0*(1.0-3617.0*z*z/161840.0)
))))))));
return temp; }
else if (x < x_2) return - li2(-x) + li2(x*x)/2.0 ;
else { return PISQ6 - li2(1.0-x)
- log(Abs(x))*log(Abs(1.0-x)) ; }
}
// COMPLEX version using 't Hooft and Veltman's change of variable
Complex CLi2(Complex x){
double x_0 = -0.30;
double x_1 = 0.25;
double x_2 = 0.51;
if (x == Complex(1.0,0.0)) return PISQ6;
if (real(x) >= x_2) { return PISQ6 - CLi2(1.-x) - log(x)*log(1.-x) ; }
if ((Abs(imag(x)) > 1.) || (real(x)*real(x) + imag(x)*imag(x) > 1.2))
return - CLi2(1./x) - 0.5 * log(-x) * log(-x) - PISQ6 ;
if (real(x) <= x_0){
Complex zz = log(1.-x);
return -CLi2(-x/(1.-x)) - zz*zz/2. ; }
else if (real(x) < x_1){
Complex z = - log(1.-x);
Complex temp = z*(1.-z/4.*(1.-z/9.*(1.-z*z/100.
*(1.-5.*z*z/294.*(1.-7.*z*z/360.
*(1.-5.*z*z/242.*(1.-7601.*z*z/354900.
*(1.-91.*z*z/4146.*(1.-3617.*z*z/161840.)
))))))));
return temp; }
else return - CLi2(-x) + CLi2(x*x)/2. ;
}
double ReLi2(Complex x){ return real(CLi2(x)) ; }
double ImLi2(Complex x){ return imag(CLi2(x)) ; }
// maple's definition
double dilog(double x){
return li2(1.-x);
}
// the trilog, li3.
// good to ~ 10^(-14)
double li3(double x){
double x_0 = -1.0;
double x_1 = -0.85;
double x_2 = 0.25;
double x_3 = 0.63;
double x_4 = 1.0;
if (x == 1.) return ZETA3;
if (x == -1.) return - 0.75 * ZETA3;
if (x <= x_0){
double lnx = log(-x);
return li3(1.0/x) - PISQ6*lnx - lnx*lnx*lnx/6.0; }
else if (x < x_1){
return li3(x*x)/4.0 - li3(-x); }
else if (x < x_2){
double z = - log(1.0-x);
double temp = z*(1.0-3.0*z/8.0*(1.0-17.0*z/81.0*(1.0-15*z/136.0
*(1.0-28.0*z/1875.0*(1.0+5.0*z/8.0*(1.0-304.0*z/7203.0
*(1.0+945.0*z/2432.0*(1.0-44.0*z/675.0*(1.0+7.0*z/24.0
*(1.0-26104.0*z/307461.0*(1.0+1925.0*z/8023.0
*(1.0-53598548.0*z/524808375.0
*(1.0+22232925.0*z/107197096.0
)))))))))))));
return temp; }
else if (x < x_3){
return li3(x*x)/4.0 - li3(-x); }
else if (x < x_4){
double ln1x = log(1.0-x);
return -li3(1.0-x) - li3(-x/(1.0-x)) + ZETA3 + PISQ6*ln1x
- log(x)*ln1x*ln1x/2.0 + ln1x*ln1x*ln1x/6.0; }
else {
double lnx = log(x);
return li3(1./x) + 2.0*PISQ6*lnx - lnx*lnx*lnx/6.0; }
}
// S12(x) -- only gets the real part correct for x > 1:
double S12(double x){
double c;
if (x > 1){
return - li3(1.-x) + ZETA3 + log(x-1.) * li2(1.-x)
+ 0.5 * log(x) * ( log(x-1.)*log(x-1.) - PISQ ); }
else if (x==1){ return ZETA3; }
else if ((0 < x) && (x < 1)){
return - li3(1.-x) + ZETA3 + log(1.-x) * li2(1.-x)
+ 0.5 * log(x) * log(1.-x)*log(1.-x); }
else if (x==0){ return 0.; }
else if (x<0){
c = 1./(1.-x);
return - li3(c) + ZETA3 + log(c) * li2(c)
+ 0.5 * log(c)*log(c) * log(1.-c)
- 1./6. * log(c)*log(c)*log(c); }
else { return 0.; }
}
// the tetralog, li4.
// good to ~ 1.7 * 10^(-12) (worst areas: x = 0.9501, also x = - 0.9701)
double li4(double x){
double x_0 = -1.0;
double x_1 = -0.97;
double x_2 = 0.25;
double x_3 = 0.95;
double x_4 = 1.0;
if (x == -1) return -0.35 * PISQ6*PISQ6;
if (x == 1) return 0.4 * PISQ6*PISQ6;
if (x <= x_0){
double lnx = log(-x);
return - li4(1./x) - 0.5 * PISQ6*lnx*lnx
- 1./24. * lnx*lnx*lnx*lnx - 0.7 * PISQ6*PISQ6; }
else if (x < x_1){
return li4(x*x)/8. - li4(-x); }
else if (x < x_2){
double z = - log(1.-x);
double temp = z*(1.-7.*z/16.*(1.-151.*z/567.*(1.0-411.*z/2416.
*(1.-24986.*z/256875.*(1.-805.*z/49972.
*(1.+583406.*z/1159683.*(1.-7455.*z/137272.
*(1.+659921.*z/2444175.*(1.-251559.*z/2639684.
*(1.+24259894.*z/136410197.*(1.-30625595.*z/218339046.
*(1.+2134239258113.*z/16698772722450.
*(1.-1640443805715.*z/8536957032452.
)))))))))))));
return temp; }
else if (x < x_3){
return li4(x*x)/8. - li4(-x); }
else if (x < x_4){
double y = 1.-x;
double lny = log(y);
return 0.4*PISQ6*PISQ6 + ZETA3 * log(1.-y) + 0.5*PISQ6 * y*y
+ (-11./36.+0.5*PISQ6+1./6.*lny) * y*y*y
+ (11./24.*PISQ6+1./4.*lny-19./48.) * y*y*y*y
+ (7./24.*lny+5./12.*PISQ6-599./1440.) * y*y*y*y*y
+ (137./360.*PISQ6+5./16.*lny-79./192.) * y*y*y*y*y*y
+ (7./20.*PISQ6-3343./8400.+29./90.*lny) * y*y*y*y*y*y*y
+ (363./1120*PISQ6-21977./57600.+469./1440.*lny) * y*y*y*y*y*y*y*y;
}
else {
double lnx = log(x);
return - li4(1/x) + PISQ6*lnx*lnx
- 1./24. * lnx*lnx*lnx*lnx + 0.8 * PISQ6*PISQ6; }
}
// the li2 below does not use 't Hooft and Veltman's change of variable
// -> not quite as good
double myli2(double x){
int n = 25;
double x_0 = -0.30;
double x_1 = 0.25;
double x_2 = 0.51;
if (x == 1) return PISQ6;
if (x <= x_0){
double temp = log(Abs(1.0-x));
return -myli2(-x/(1.0-x)) - temp*temp/2 ; }
else if (x < x_1){
double temp = 0.0;
double xk = 1.0;
for (int k=1; k<=n; k++){
xk = xk*x;
temp = temp + xk/k/k; }
return temp; }
else if (x < x_2) return - myli2(-x) + myli2(x*x)/2.0 ;
else { return PISQ6 - myli2(1.0-x)
- log(Abs(x))*log(Abs(1.0-x)) ; }
}
// the 3-mass triangle, i3_3m(x,y,z), presently only for the pure
// Minkowski region, x,y,z > 0:
Complex i3_3m(double s1, double s2, double s3){
double x, y, lambda, rho, mylamsq, mylam, eta1, eta2, eta3, rho2;
x = s1/s3;
y = s2/s3;
lambda = sqrt((1.0-x-y)*(1.0-x-y) - 4.0*x*y);
rho = 2.0/(1.0-x-y+lambda);
mylamsq = s1*s1+s2*s2+s3*s3-2.*(s1*s2+s2*s3+s3*s1);
mylam = sqrt(mylamsq);
// Euclidean region:
if (s1<0 && s2<0 && s3<0) {
if (mylamsq > 0) {
double temp = 2.0 * ( li2(-rho*x) + li2(-rho*y) )
+ log(Abs(rho*x))*log(Abs(rho*y))
+ log(y/x)*log(Abs((1.0+rho*y)/(1.0+rho*x)))
+ PISQ/3. ;
if (rho<0) temp = temp - PISQ;
return -temp/s3/lambda;
}
else {
eta1 = (s1-s2-s3)/2.0;
eta2 = (s2-s3-s1)/2.0;
eta3 = (s3-s1-s2)/2.0;
rho2 = sqrt(eta1*eta2 + eta2*eta3 + eta3*eta1);
return 1.0/rho2 * ( fastCl(2.0*atan(rho2/eta1))
+ fastCl(2.0*atan(rho2/eta2))
+ fastCl(2.0*atan(rho2/eta3)) ); }
}
else if (s1>0 && s2>0 && s3>0){ return -i3_3m(-s1,-s2,-s3); }
else if (s1<0 && s2>0 && s3<0){
double xplus = (s1+s2-s3+mylam)/2./s2;
double xminus = (s1+s2-s3-mylam)/2./s2;
double temp = 2.0 * ( ReLi2(-rho*x) + ReLi2(-rho*y) )
+ log(Abs(rho*x))*log(Abs(rho*y))
+ log(Abs(y/x))*log(Abs((1.0+rho*y)/(1.0+rho*x)))
+ PISQ/3. ;
if (rho<0){ temp = temp - PISQ; }
return - temp/s3/lambda
- PI_DEF/mylam * I * log(Abs((1.-xplus)*xminus/xplus/(1.-xminus))) ;
}
else if (s1<0 && s2<0 && s3>0){ return i3_3m(s2,s3,s1); }
else if (s1>0 && s2<0 && s3<0){ return i3_3m(s3,s1,s2); }
else return - conj(i3_3m(-s1,-s2,-s3));
}
// Auxiliary function for the 3-mass triangle:
double fastCl(double x){
double xb, x_lower, x_upper, x_switch;
x_switch = 0.75*PI_DEF;
x_lower = -x_switch;
x_upper = 2.0*PI_DEF+x_lower;
// put x in the range [x_lower,x_upper]:
if (x > x_upper) return fastCl(x-2.0*PI_DEF);
if (x < x_lower) return fastCl(x+2.0*PI_DEF);
// dividing line for the two expansions is x_switch:
if (x < x_switch){ // expansion around 0:
return -x*log(Abs(x)) + x * ( 1.0 + x*x/72.0 * ( 1.0 + x*x/200.0
* ( 1.0 + x*x*5.0/441.0 * ( 1.0 + x*x*7.0/480.0
* ( 1.0 + x*x*2.0/121.0 ) ) ) ) ); }
else {
xb = x-PI_DEF;
return -xb*log(2.0) + xb*xb*xb/24.0 * ( 1.0 + xb*xb/40.0
* ( 1.0 + xb*xb/21.0 * ( 1.0 + xb*xb/288.0
* ( 1.0 + xb*xb*62.0/55.0 * ( 1.0 + xb*xb*691.0/9672.0
* ( 1.0 + xb*xb*5461.0/72555.0 ) ) ) ) ) );
}
}
//==========================================================
// Basic functions of s_{ij}:
// The L_i and Ls_i functions.
// get ln((-s1)/(-s2))
Complex Clog(double s1, double s2){ // Complex log of a ratio, s1/s2
Complex temp = log(Abs(s1/s2));
if (s1 > 0 && s2 < 0) { temp = temp - I * PI; }
if (s2 > 0 && s1 < 0) { temp = temp + I * PI; }
return temp;
}
// get ln(-s/mu_sq)
Complex Clog1(double s){ // Complex log of a single s
Complex temp = log(Abs(s/mu_sq));
if (s > 0) { temp = temp - I * PI; }
return temp;
}
Complex L0(double s1, double s2){
double r = s1/s2;
return Clog(s1,s2)/(1.-r);
}
Complex L1(double s1, double s2){
double r = s1/s2;
return (Clog(s1,s2)+1.-r)/(1.-r)/(1.-r);
}
Complex L2(double s1, double s2){
double r = s1/s2;
return (Clog(s1,s2)-(r-1./r)/2.)/(1.-r)/(1.-r)/(1.-r);
}
// Li_2(1-r), for r = s1/s2:
Complex CLi2r(double s1, double s2){
double r = s1/s2;
return li2(1.-r) - log(Abs(1.-r)) * I * imag(Clog(s1,s2));
}
Complex Lsm1(double s1, double s2, double s3, double s4){
return CLi2r(s1,s2) + CLi2r(s3,s4) + Clog(s1,s2)*Clog(s3,s4) - PISQ6;
}
Complex Ls0(double s1, double s2, double s3, double s4){
double r1 = s1/s2;
double r2 = s3/s4;
return Lsm1(s1,s2,s3,s4)/(1.-r1-r2);
}
Complex Ls1(double s1, double s2, double s3, double s4){
double r1 = s1/s2;
double r2 = s3/s4;
return ( Ls0(s1,s2,s3,s4) + L0(s1,s2) + L0(s3,s4) )/(1.-r1-r2);
}
Complex Ls2(double s1, double s2, double s3, double s4){
double r1 = s1/s2;
double r2 = s3/s4;
return ( Ls1(s1,s2,s3,s4) + 0.5 * ( L1(s1,s2) + L1(s3,s4) ) )/(1.-r1-r2);
}
Complex Ls3(double s1, double s2, double s3, double s4){
double r1 = s1/s2;
double r2 = s3/s4;
return ( Ls2(s1,s2,s3,s4) + 1./3. * ( L2(s1,s2) + L2(s3,s4) )
- 1./6. * (1./r1+1./r2) )/(1.-r1-r2);
}