Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 2 additions & 4 deletions src/bubble.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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]) +
Expand Down
76 changes: 42 additions & 34 deletions src/tools.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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():
Expand Down Expand Up @@ -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;
}
Expand Down Expand Up @@ -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;
}
Expand Down Expand Up @@ -819,7 +819,7 @@ namespace ql {
{
const TOutput arg4 = (y0-_cone)/y0;
res += extra*cLn(arg4, Sign(Imag(arg4)));
}
}

return res;
}
Expand Down Expand Up @@ -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
{
Expand Down