Merge pull request #22 from scarrazza/fixgz

importing changes from GZ
scarrazza · May 20, 2022 · 6226972 · 6226972
2 parents ba82bf6 + d22e433
commit 6226972
Show file tree

Hide file tree

Showing 2 changed files with 44 additions and 38 deletions.
diff --git a/src/bubble.cc b/src/bubble.cc
@@ -111,8 +111,7 @@ namespace ql
     const TMass sqm1 = Sqrt(m1);
     const TOutput bb = TOutput(m0+m1-s);
     const TOutput rtt= Sqrt(bb*bb - this->_cfour*TOutput(m1*m0));
-    const TOutput sgn = Sign(Real(Conjg(bb)*rtt));
-    const TOutput x1 = this->_chalf*(bb+sgn*rtt)/(sqm0*sqm1);
+    const TOutput x1 = this->_chalf*(bb+rtt)/(sqm0*sqm1);
     const TOutput x2 = this->_cone/x1;
     res[0] = this->_ctwo - Log(sqm0*sqm1/mu2) + (m0-m1)/s*Log(sqm1/sqm0) - sqm0*sqm1/s*(x2-x1)*this->cLn(x1, Sign(Real(x1-x2)));
     res[1] = this->_cone;
@@ -279,8 +278,7 @@ namespace ql
         const TMass c = a;
         const TOutput b = m[0] + m[1] - TMass(p[0]) - this->_ieps;
         const TOutput root = Sqrt(Pow(b, 2) - this->_cfour*a*c);
-        const TOutput sgn = TOutput(Sign(Real( Conjg(b) * root)));
-        const TOutput q = this->_chalf * (b + sgn*root);
+        const TOutput q = this->_chalf * (b + root);
         const TOutput rm = q / a;
         const TOutput r = rm;
         res[0] = - this->_chalf * (m[0] - m[1]) / Pow(p[0], 2) * Log(m[1]/m[0]) +

diff --git a/src/tools.cc b/src/tools.cc
@@ -15,7 +15,7 @@ using std::cout;
 using std::endl;
 using std::is_same;
 
-namespace ql {    
+namespace ql {
 
   template<typename TOutput, typename TMass, typename TScale>
   Tools<TOutput,TMass,TScale>::Tools():
@@ -269,21 +269,21 @@ namespace ql {
   template<typename TOutput, typename TMass, typename TScale>
   TOutput Tools<TOutput,TMass,TScale>::fndd(int const& n, TOutput const& x, TScale const& iep) const
   {
-    const int infty = 16;    
-    TOutput res = _czero;    
+    const int infty = 16;
+    TOutput res = _czero;
     if (Abs(x) < _ten)
       {
         if (!iszero(Abs(x-_cone)))
           res = (_cone-Pow(x,n+1))*(cLn(x-_cone, iep) - cLn(x, iep));
 
         for (int j = 0; j <= n; j++)
           res -= Pow(x, n-j)/(j+_one);
-      }    
+      }
     else
-      {        
+      {
         res = cLn(_cone-_cone/x, iep);
         for (int j = n+1; j <= n+infty; j++)
-          res += Pow(x, n-j)/(j+_one);          
+          res += Pow(x, n-j)/(j+_one);
       }
     return res;
   }
@@ -481,18 +481,18 @@ namespace ql {
 
     if (Real(z12) > _half)
       {
-        cspence = ltspence(1, z12, _zero);       
+        cspence = ltspence(1, z12, _zero);
         const int etas = eta(z1, im1, z2, im2, im12);
         if (etas != 0) cspence += TOutput(etas)*cLn(_cone-z12, -im12)*_2ipi;
       }
     else if (Abs(z12) < _eps4)
       {
-        cspence = TOutput(_pi2o6);        
+        cspence = TOutput(_pi2o6);
         if (Abs(z12) > _eps14)
           cspence += -ltspence(0, z12, _zero) + (cLn(z1,im1) + cLn(z2,im2))*z12*(_cone + z12*(_chalf + z12*(_cone/_cthree + z12/_cfour)));
       }
     else
-      cspence = TOutput(_pi2o6) - ltspence(0, z12, _zero) - (cLn(z1, im1) + cLn(z2, im2))*cLn(_cone-z12,_zero);      
+      cspence = TOutput(_pi2o6) - ltspence(0, z12, _zero) - (cLn(z1, im1) + cLn(z2, im2))*cLn(_cone-z12,_zero);
 
     return cspence;
   }
@@ -819,7 +819,7 @@ namespace ql {
       {
         const TOutput arg4 = (y0-_cone)/y0;
         res += extra*cLn(arg4, Sign(Imag(arg4)));
-      }      
+      }
 
     return res;
   }
@@ -1221,30 +1221,38 @@ namespace ql {
     if (iszero(Imag(discr)))
       {
         const TMass sgnb = Sign(Real(b));
-        if (Real(discr) > 0)
-            {
-              const TMass q = -_half*(b+sgnb*Sqrt(discr));
-              if (Real(b) > 0)
-                {
-                  z[0] = TOutput(c/q);
-                  z[1] = TOutput(q/a);
-                }
-              else
-                {
-                  z[0] = TOutput(q/a);
-                  z[1] = TOutput(c/q);
-                }
-            }
-          else
-            {
-              z[1] = -(TOutput(b)+sgnb*Sqrt(TOutput(discr)))/(_ctwo*a);
-              z[0] = Conjg(z[1]);
-              if (Real(b) < 0)
-                {
-                  z[0] = z[1];
-                  z[1] = Conjg(z[0]);
-                }
-            }
+        if (iszero(Real(b)))
+        {
+          z[0] = -(TOutput(b)-Sqrt(TOutput(discr)))/(_ctwo*a);
+          z[1] = -(TOutput(b)+Sqrt(TOutput(discr)))/(_ctwo*a);
+        }
+        else
+        {
+          if (Real(discr) > 0)
+              {
+                const TMass q = -_half*(b+sgnb*Sqrt(discr));
+                if (Real(b) > 0)
+                  {
+                    z[0] = TOutput(c/q);
+                    z[1] = TOutput(q/a);
+                  }
+                else
+                  {
+                    z[0] = TOutput(q/a);
+                    z[1] = TOutput(c/q);
+                  }
+              }
+            else
+              {
+                z[1] = -(TOutput(b)+sgnb*Sqrt(TOutput(discr)))/(_ctwo*a);
+                z[0] = Conjg(z[1]);
+                if (Real(b) < 0)
+                  {
+                    z[0] = z[1];
+                    z[1] = Conjg(z[0]);
+                  }
+              }
+        }
       }
     else
       {