Merge branch 'development' into feature-SplitMSSM

addons/gm2calc/ffunctions.cpp

@@ -17,7 +17,7 @@
 // ====================================================================
 
 #include "ffunctions.hpp"
-#include "dilog.h"
+#include "dilog.hpp"
 #include "numerics2.hpp"
 
 #include <iostream>
@@ -30,91 +30,10 @@ namespace gm2calc {
 
 using namespace flexiblesusy;
 
-double dilog(double x) {
-   // The DiLogarithm function
-   // Code translated by R.Brun from CERNLIB DILOG function C332
-
-   const double PI = M_PI;
-   const double HF  = 0.5;
-   const double PI2 = PI*PI;
-   const double PI3 = PI2/3;
-   const double PI6 = PI2/6;
-   const double PI12 = PI2/12;
-   const double C[20] = {0.42996693560813697, 0.40975987533077105,
-     -0.01858843665014592, 0.00145751084062268,-0.00014304184442340,
-      0.00001588415541880,-0.00000190784959387, 0.00000024195180854,
-     -0.00000003193341274, 0.00000000434545063,-0.00000000060578480,
-      0.00000000008612098,-0.00000000001244332, 0.00000000000182256,
-     -0.00000000000027007, 0.00000000000004042,-0.00000000000000610,
-      0.00000000000000093,-0.00000000000000014, 0.00000000000000002};
-
-   double T,H,Y,S,A,ALFA,B1,B2,B0;
-
-   if (x == 1) {
-       H = PI6;
-   } else if (x == -1) {
-       H = -PI12;
-   } else {
-       T = -x;
-       if (T <= -2) {
-           Y = -1/(1+T);
-           S = 1;
-           B1= log(-T);
-           B2= log(1+1/T);
-           A = -PI3+HF*(B1*B1-B2*B2);
-       } else if (T < -1) {
-           Y = -1-T;
-           S = -1;
-           A = log(-T);
-           A = -PI6+A*(A+log(1+1/T));
-       } else if (T <= -0.5) {
-           Y = -(1+T)/T;
-           S = 1;
-           A = log(-T);
-           A = -PI6+A*(-HF*A+log(1+T));
-       } else if (T < 0) {
-           Y = -T/(1+T);
-           S = -1;
-           B1= log(1+T);
-           A = HF*B1*B1;
-       } else if (T <= 1) {
-           Y = T;
-           S = 1;
-           A = 0;
-       } else {
-           Y = 1/T;
-           S = -1;
-           B1= log(T);
-           A = PI6+HF*B1*B1;
-       }
-       H    = Y+Y-1;
-       ALFA = H+H;
-       B1   = 0;
-       B2   = 0;
-       for (int i=19;i>=0;i--){
-          B0 = C[i] + ALFA*B1-B2;
-          B2 = B1;
-          B1 = B0;
-       }
-       H = -(S*(B0-H*B2)+A);
-    }
-    return H;
-}
-
-std::complex<double> dilog(const std::complex<double>& x) {
-   double a = x.real(), b = x.imag();
-   double ansreal = 0., ansimag = 0.;
-
-   dilogc_(&a, &b, &ansreal, &ansimag);
-
-   return std::complex<double>(ansreal, ansimag);
-}
-
-double sqr(double x) { return x*x; }
-double cube(double x) { return x*x*x; }
-double quad(double x) { return sqr(x)*sqr(x); }
 int sign(double x) { return x < 0 ? -1 : 1; }
+
 double signed_sqr(double x) { return sign(x) * x * x; }
+
 double signed_abs_sqrt(double x) {
    return sign(x) * std::sqrt(std::abs(x));
 }

addons/gm2calc/ffunctions.hpp

@@ -25,13 +25,17 @@ namespace gm2calc {
 double F1C(double);
 /// \f$F_2^C(x)\f$, Eq (55) arXiv:hep-ph/0609168
 double F2C(double);
+/// \f$F_3^C(x)\f$, Eq (37) arXiv:1003.5820
 double F3C(double);
+/// \f$F_4^C(x)\f$, Eq (38) arXiv:1003.5820
 double F4C(double);
 /// \f$F_1^N(x)\f$, Eq (52) arXiv:hep-ph/0609168
 double F1N(double);
 /// \f$F_2^N(x)\f$, Eq (53) arXiv:hep-ph/0609168
 double F2N(double);
+/// \f$F_3^N(x)\f$, Eq (39) arXiv:1003.5820
 double F3N(double);
+/// \f$F_4^N(x)\f$, Eq (40) arXiv:1003.5820
 double F4N(double);
 /// \f$F_a(x)\f$, Eq (6.3a) arXiv:1311.1775
 double Fa(double, double);
@@ -43,6 +47,7 @@ double G3(double);
 double G4(double);
 /// \f$H_2(x,y)\f$, Eq (34) of arxiv:0901.2065
 double H2(double, double);
+/// \f$I_{abc}(a,b,c)\f$ (arguments are interpreted as unsquared)
 double Iabc(double, double, double);
 /// \f$f_{PS}(z)\f$, Eq (70) arXiv:hep-ph/0609168
 double f_PS(double);
@@ -51,12 +56,18 @@ double f_S(double);
 /// \f$f_{\tilde{f}}(z)\f$, Eq (72) arXiv:hep-ph/0609168
 double f_sferm(double);
 
-double cube(double);
-double quad(double);
+/// returns number squared
+template <typename T> T sqr(T x) { return x*x; }
+/// returns number cubed
+template <typename T> T cube(T x) { return x*x*x; }
+/// returns number to the power 4
+template <typename T> T quad(T x) { return x*x*x*x; }
+/// returns sign of real number
 int sign(double);
+/// returns square root of absolute of number, times sign
 double signed_abs_sqrt(double);
+/// returns square of number, times sign
 double signed_sqr(double);
-double sqr(double);
 
 } // namespace gm2calc
 

src/wrappers.hpp

@@ -296,7 +296,7 @@ inline int Sign(int x)
 template <typename T>
 T PolyLog(int n, T z) {
    if (n == 2)
-      return dilog(z);
+      return gm2calc::dilog(z);
    assert(false && "PolyLog(n!=2) not implemented");
 }