Merge pull request #77 from niclaurenti/fix-overflow

Guards against numerical overflows that result in `NaN`

Author: Radonirinaunimi

src/Convolution.cc

@@ -126,6 +126,12 @@ double Convolution::regular_integrand(double z, void *p) {
     double x = (params->x);
     int nf = (params->nf);
 
+    // When the integration interval [x, x_max] is very small, a quadrature node can round to
+    // z = x in double precision, giving x/z = 1 exactly. Some splitting functions (e.g. Pgg1reg)
+    // produce NaN at x/z = 1 due to 0*inf in harmonic polylogarithms. And given that the true
+    // integrand is finite there, we can simply return 0.
+    if (x / z >= 1.) return 0.;
+
     return (params->conv)->GetCoeffFunc()->MuIndependentTerms(z, m2Q2, nf)
            * (params->conv)->GetSplitFunc()->Regular(x / z, nf) / z;
 }
@@ -198,6 +204,12 @@ double Convolution::singular_integrand(double z, void *p) {
     double x = (params->x);
     int nf = (params->nf);
 
+    // When the integration interval [x/x_max, 1] is very small, a quadrature node can round to
+    // z = 1 in double precision. At z = 1 the plus-distribution subtraction (f(x/z)/z - f(x))
+    // vanishes while Singular(z) = 1/(1-z) diverges, producing inf * 0 = NaN.  Given that the
+    // true limit is finite, we can simply return 0.
+    if (z >= 1.) return 0.;
+
     return (params->conv)->GetSplitFunc()->Singular(z, nf)
            * ((params->conv)
                       ->GetCoeffFunc()
@@ -468,6 +480,9 @@ double
 
     double z1 = z[0], z2 = z[1];
 
+    // Guard against x/z1 or z1/z2 rounding to 1 (see `regular_integrand`).
+    if (x / z1 >= 1. || z1 / z2 >= 1.) return 0.;
+
     if (z2 > z1) {
         return 1. / (z1 * z2)
                * (params->conv)->GetSplitFunc()->Regular(x / z1, nf)
@@ -494,6 +509,9 @@ double
 
     double z1 = z[0], z2 = z[1];
 
+    // Guard against x/z1 rounding to 1 or z1/z2 rounding to 1.
+    if (x / z1 >= 1. || z1 / z2 >= 1.) return 0.;
+
     if (z2 > z1) {
         return 1. / (z1 * z2)
                * (params->conv)->GetSplitFunc()->Regular(x / z1, nf)
@@ -522,6 +540,9 @@ double DoubleConvolution::regular3_integrand(double z, void *p) {
     double x = (params->x);
     int nf = (params->nf);
 
+    // Guard against x/z rounding to 1.
+    if (x / z >= 1.) return 0.;
+
     double x_max = CoefficientFunction::xMax(m2Q2);
 
     return -1. / z * (params->conv)->GetSplitFunc()->Regular(x / z, nf)
@@ -599,6 +620,11 @@ double DoubleConvolution::singular1_integrand(
 
     double z1 = z[0], z2 = z[1];
 
+    // At z1 = 1 the plus-distribution subtraction vanishes while Singular diverges,
+    // producing inf * 0 = NaN. Given that the true integrand is finite, we can simply
+    // return 0.
+    if (z1 >= 1.) return 0.;
+
     double tmp;
     if (z2 > x / z1) {
         tmp = (params->conv)->GetSplitFunc()->Regular(x / (z1 * z2), nf) / z1;
@@ -630,6 +656,10 @@ double DoubleConvolution::singular2_integrand(
 
     double z1 = z[0], z2 = z[1];
 
+    // At z1 = 1 or z2 = 1 the plus-distribution subtraction vanishes while Singular diverges,
+    // producing inf * 0 = NaN.
+    if (z1 >= 1. || z2 >= 1.) return 0.;
+
     double tmp;
     if (z2 > x / (x_max * z1)) {
         tmp = ((params->conv)
@@ -668,6 +698,10 @@ double DoubleConvolution::singular3_integrand(double z, void *p) {
     double x = (params->x);
     int nf = (params->nf);
 
+    // Guard: at z = 1 the plus-distribution subtraction vanishes while Singular diverges,
+    // producing inf * 0 = NaN.
+    if (z >= 1.) return 0.;
+
     double x_max = CoefficientFunction::xMax(m2Q2);
 
     return -(

src/ExactCoefficientFunction.cc

@@ -436,6 +436,12 @@ double ExactCoefficientFunction::fx(
 double ExactCoefficientFunction::MuIndependentTerms(
     double x, double m2Q2, int nf
 ) const {
+    // The coefficient function vanishes at x = xmax due to a constrained phase-space. We need
+    // to avoid x >= xmax to prevent sqrt of a negative number due to floating-point rounding
+    // in the beta computation: the expression 1 - 4*m2Q2*x/(1-x), which should be 0 at x = xmax,
+    // can become slightly negative (~-1e-9) in floating point.
+    if (x >= xMax(m2Q2)) return 0.;
+
     try {
         return (this->*mu_indep_)(x, m2Q2, nf);
     } catch (NotValidException &e) {
@@ -454,8 +460,7 @@ double ExactCoefficientFunction::MuDependentTerms(
     double x, double m2Q2, double m2mu2, int nf
 ) const {
 
-    if (x <= 0 || x > xMax(m2Q2))
-        return 0.;
+    if (x <= 0 || x >= xMax(m2Q2)) return 0.;
 
     return (this->*mu_dep_)(x, m2Q2, m2mu2, nf);
 }
@@ -480,7 +485,9 @@ Value ExactCoefficientFunction::fxBand(
 double ExactCoefficientFunction::C2_g1(double x, double m2Q2, int /*nf*/)
     const { // m2Q2=m^2/Q^2
 
-    double beta = sqrt(1. - 4. * m2Q2 * x / (1 - x));
+    double beta2 = 1. - 4. * m2Q2 * x / (1. - x);
+    if (beta2 < 0.) beta2 = 0.;
+    double beta = sqrt(beta2);
     double x2 = x * x;
     double m4Q4 = m2Q2 * m2Q2;
     double L = log((1. + beta) / (1. - beta));
@@ -501,7 +508,9 @@ double ExactCoefficientFunction::C2_g1(double x, double m2Q2, int /*nf*/)
 double
     ExactCoefficientFunction::CL_g1(double x, double m2Q2, int /*nf*/) const {
 
-    double beta = sqrt(1. - 4. * m2Q2 * x / (1. - x));
+    double beta2 = 1. - 4. * m2Q2 * x / (1. - x);
+    if (beta2 < 0.) beta2 = 0.;
+    double beta = sqrt(beta2);
     double x2 = x * x;
     double L = log((1. + beta) / (1. - beta));
 

src/ThresholdCoefficientFunction.cc

@@ -331,7 +331,9 @@ double ThresholdCoefficientFunction::BetaIndependentTerms(
 double
     ThresholdCoefficientFunction::C2_g1_threshold(double x, double m2Q2) const {
 
-    double beta = sqrt(1. - 4. * m2Q2 * x / (1. - x));
+    double beta2 = 1. - 4. * m2Q2 * x / (1. - x);
+    if (beta2 < 0.) beta2 = 0.;
+    double beta = sqrt(beta2);
     double xi = 1. / m2Q2;
 
     return xi * TR * beta / (1. + xi / 4.) / x;
@@ -348,7 +350,9 @@ double
 double
     ThresholdCoefficientFunction::CL_g1_threshold(double x, double m2Q2) const {
 
-    double beta = sqrt(1. - 4. * m2Q2 * x / (1. - x));
+    double beta2 = 1. - 4. * m2Q2 * x / (1. - x);
+    if (beta2 < 0.) beta2 = 0.;
+    double beta = sqrt(beta2);
     double beta3 = beta * beta * beta;
 
     double xi = 1. / m2Q2;
@@ -368,7 +372,10 @@ double ThresholdCoefficientFunction::C2_g2_threshold_expansion(
     double x, double m2Q2, double m2mu2, int /*nf*/
 ) const {
 
-    double beta = sqrt(1. - 4. * m2Q2 * x / (1. - x));
+    double beta2 = 1. - 4. * m2Q2 * x / (1. - x);
+    if (beta2 < 0.) beta2 = 0.;
+    double beta = sqrt(beta2);
+    if (beta == 0.) return 0.;
 
     double logb = log(beta);
     double log2b = logb * logb;
@@ -388,7 +395,10 @@ double ThresholdCoefficientFunction::CL_g2_threshold_expansion(
     double x, double m2Q2, double m2mu2, int /*nf*/
 ) const {
 
-    double beta = sqrt(1. - 4. * m2Q2 * x / (1. - x));
+    double beta2 = 1. - 4. * m2Q2 * x / (1. - x);
+    if (beta2 < 0.) beta2 = 0.;
+    double beta = sqrt(beta2);
+    if (beta == 0.) return 0.;
 
     double logb = log(beta);
     double log2b = logb * logb;
@@ -449,7 +459,10 @@ double ThresholdCoefficientFunction::C2_g3_threshold_expansion(
 ) const {
 
     double xi = 1. / m2Q2;
-    double beta = sqrt(1. - 4. * m2Q2 * x / (1. - x));
+    double beta2 = 1. - 4. * m2Q2 * x / (1. - x);
+    if (beta2 < 0.) beta2 = 0.;
+    double beta = sqrt(beta2);
+    if (beta == 0.) return 0.;
 
     double Lm = log(m2mu2);
     double Lm2 = Lm * Lm;
@@ -539,7 +552,10 @@ double ThresholdCoefficientFunction::C2_g3_threshold_expansion_const(
 double ThresholdCoefficientFunction::CL_g3_threshold_expansion(
     double x, double m2Q2, double m2mu2, int nf
 ) const {
-    double beta = sqrt(1. - 4. * m2Q2 * x / (1. - x));
+    double beta2 = 1. - 4. * m2Q2 * x / (1. - x);
+    if (beta2 < 0.) beta2 = 0.;
+    double beta = sqrt(beta2);
+    if (beta == 0.) return 0.;
 
     double rhoq = -4. * m2Q2;
     double betaq = sqrt(1. - rhoq);