diff --git a/examples/test.cpp b/examples/test.cpp
index 03f58c8a..ca2b4c1f 100644
--- a/examples/test.cpp
+++ b/examples/test.cpp
@@ -19,16 +19,16 @@ int main() {
     int nf = 4;
 
     ApproximateCoefficientFunction F2g(
-        3, '2', 'g', true, "abmp", 1e-3, 1e-3, 1000, false, 25000
+        3, '2', 'g', true, "abmp", 1e-3, 1e-3, 1000, "analytical", 25000
     );
     ApproximateCoefficientFunction FLg(
-        3, 'L', 'g', true, "abmp", 1e-3, 1e-3, 1000, false, 25000
+        3, 'L', 'g', true, "abmp", 1e-3, 1e-3, 1000, "analytical", 25000
     );
     ApproximateCoefficientFunction F2q(
-        3, '2', 'q', true, "exact", 1e-3, 1e-3, 1000, false, 25000
+        3, '2', 'q', true, "exact", 1e-3, 1e-3, 1000
     );
     ApproximateCoefficientFunction FLq(
-        3, 'L', 'q', true, "exact", 1e-3, 1e-3, 1000, false, 25000
+        3, 'L', 'q', true, "exact", 1e-3, 1e-3, 1000
     );
 
     for (int i = 0; i < N; i++) {
diff --git a/examples/test.py b/examples/test.py
index 5c7787bd..fe963827 100644
--- a/examples/test.py
+++ b/examples/test.py
@@ -4,10 +4,10 @@
 nf = 4
 m = 4.92
 
-F2g = ad.ApproximateCoefficientFunction(3, '2', 'g', True, "abmp", 1e-3, 1e-3, 1000, False, 25000)
-FLg = ad.ApproximateCoefficientFunction(3, 'L', 'g', True, "abmp", 1e-3, 1e-3, 1000, False, 25000)
-F2q = ad.ApproximateCoefficientFunction(3, '2', 'q', True, "exact", 1e-3, 1e-3, 1000, False, 25000)
-FLq = ad.ApproximateCoefficientFunction(3, 'L', 'q', True, "exact", 1e-3, 1e-3, 1000, False, 25000)
+F2g = ad.ApproximateCoefficientFunction(3, '2', 'g', True, "abmp", 1e-3, 1e-3, 1000, "analytical", 25000)
+FLg = ad.ApproximateCoefficientFunction(3, 'L', 'g', True, "abmp", 1e-3, 1e-3, 1000, "analytical", 25000)
+F2q = ad.ApproximateCoefficientFunction(3, '2', 'q', True, "exact", 1e-3, 1e-3, 1000)
+FLq = ad.ApproximateCoefficientFunction(3, 'L', 'q', True, "exact", 1e-3, 1e-3, 1000)
 
 for x in np.geomspace(1e-4, 1, 5):
     for Q in np.geomspace(10, 100, 5):
diff --git a/inc/adani/ApproximateCoefficientFunction.h b/inc/adani/ApproximateCoefficientFunction.h
index ed71e59a..6eacd192 100644
--- a/inc/adani/ApproximateCoefficientFunction.h
+++ b/inc/adani/ApproximateCoefficientFunction.h
@@ -64,7 +64,8 @@ class AbstractApproximate : public CoefficientFunction {
         AbstractApproximate(
             const int &order, const char &kind, const char &channel,
             const double &abserr = 1e-3, const double &relerr = 1e-3,
-            const int &dim = 1000, const bool &analytical_PggxPgg = true, const bool &MCintegral = false,
+            const int &dim = 1000,
+            const string &double_int_method = "analytical",
             const int &MCcalls = 25000
         );
         ~AbstractApproximate();
@@ -91,7 +92,8 @@ class ApproximateCoefficientFunction : public AbstractApproximate {
             const int &order, const char &kind, const char &channel,
             const bool &NLL = true, const string &highscale_version = "klmv",
             const double &abserr = 1e-3, const double &relerr = 1e-3,
-            const int &dim = 1000, const bool &analytical_PggxPgg = true, const bool &MCintegral = false,
+            const int &dim = 1000,
+            const string &double_int_method = "analytical",
             const int &MCcalls = 25000
         );
         ~ApproximateCoefficientFunction() override;
@@ -124,7 +126,8 @@ class ApproximateCoefficientFunctionKLMV : public AbstractApproximate {
             const int &order, const char &kind, const char &channel,
             const string &highscale_version = "klmv", const bool &lowxi = false,
             const double &abserr = 1e-3, const double &relerr = 1e-3,
-            const int &dim = 1000, const bool &analytical_PggxPgg = true, const bool &MCintegral = false,
+            const int &dim = 1000,
+            const string &double_int_method = "analytical",
             const int &MCcalls = 25000
         );
         ~ApproximateCoefficientFunctionKLMV() override;
diff --git a/inc/adani/ExactCoefficientFunction.h b/inc/adani/ExactCoefficientFunction.h
index 74404bd7..ef345e67 100644
--- a/inc/adani/ExactCoefficientFunction.h
+++ b/inc/adani/ExactCoefficientFunction.h
@@ -22,8 +22,11 @@
 #include "adani/Convolution.h"
 #include "adani/SplittingFunction.h"
 
+#include <string>
 #include <vector>
 
+using std::string;
+
 //==========================================================================================//
 //  forward declaration of class ThresholdCoefficientFunction to avoid circular
 //  dependencies
@@ -41,8 +44,8 @@ class ExactCoefficientFunction : public CoefficientFunction {
         ExactCoefficientFunction(
             const int &order, const char &kind, const char &channel,
             const double &abserr = 1e-3, const double &relerr = 1e-3,
-            const int &dim = 1000, const bool &analytical_PggxPgg = true,
-            const bool &MCintegral = false,
+            const int &dim = 1000,
+            const string &double_int_method = "analytical",
             const int &MCcalls = 25000
         );
         ~ExactCoefficientFunction() override;
@@ -85,7 +88,7 @@ class ExactCoefficientFunction : public CoefficientFunction {
         SplittingFunction *Pgg1_;
         SplittingFunction *Pqg0_;
         ConvolutedSplittingFunctions *Pgq0Pqg0_;
-        ConvolutedSplittingFunctions * Pgg0Pgg0_;
+        ConvolutedSplittingFunctions *Pgg0Pgg0_;
 
         Delta *delta_;
 
diff --git a/inc/adani/SplittingFunction.h b/inc/adani/SplittingFunction.h
index 96624d73..19cc20c9 100644
--- a/inc/adani/SplittingFunction.h
+++ b/inc/adani/SplittingFunction.h
@@ -151,10 +151,10 @@ class ConvolutedSplittingFunctions : public AbstractSplittingFunction {
 
         // Components of the Convoluted Splitting Function
         double Regular(double x, int nf) const override;
-        double Singular(double x, int nf) const override ;
-        double Local(int nf) const override ;
+        double Singular(double x, int nf) const override;
+        double Local(int nf) const override;
         // Integral from 0 to x of the Singular part
-        double SingularIntegrated(double x, int nf) const override ;
+        double SingularIntegrated(double x, int nf) const override;
 
         void SetFunctions();
 
@@ -179,7 +179,7 @@ class ConvolutedSplittingFunctions : public AbstractSplittingFunction {
         double (ConvolutedSplittingFunctions::*sing_int_)(double, int) const;
         double (ConvolutedSplittingFunctions::*loc_)(int) const;
 
-        SplittingFunction* Pgg0_ ;
+        SplittingFunction *Pgg0_;
 
         //==========================================================================================//
         //  Convolution between the splitting functions Pgg0/Pqq0
diff --git a/output/output_grid.py b/output/output_grid.py
index 2618ee60..817cd19c 100644
--- a/output/output_grid.py
+++ b/output/output_grid.py
@@ -31,7 +31,7 @@
     1e-3,
     1e-3,
     1000,
-    False,
+    "analytical",
     25000,
 )
 
diff --git a/pywrap/pywrap.cc b/pywrap/pywrap.cc
index 07de9ca7..03553037 100644
--- a/pywrap/pywrap.cc
+++ b/pywrap/pywrap.cc
@@ -56,11 +56,11 @@ PYBIND11_MODULE(_core, m) {
             py::init<
                 const int &, const char &, const char &, const bool &,
                 const string &, const double &, const double &, const int &,
-                const bool &, const bool &, const int &>(),
+                const string &, const int &>(),
             py::arg("order"), py::arg("kind"), py::arg("channel"),
             py::arg("NLL") = true, py::arg("highscale_version") = "klmv",
             py::arg("abserr") = 1e-3, py::arg("relerr") = 1e-3,
-            py::arg("dim") = 1000, py::arg("analytical_PggxPgg") = true, py::arg("MCintegral") = false,
+            py::arg("dim") = 1000, py::arg("double_int_method") = "analytical",
             py::arg("MCcalls") = 25000
         )
         .def(
@@ -100,11 +100,11 @@ PYBIND11_MODULE(_core, m) {
             py::init<
                 const int &, const char &, const char &, const string &,
                 const bool &, const double &, const double &, const int &,
-                const bool &, const bool &, const int &>(),
+                const string &, const int &>(),
             py::arg("order"), py::arg("kind"), py::arg("channel"),
             py::arg("highscale_version") = "klmv", py::arg("lowxi") = false,
             py::arg("abserr") = 1e-3, py::arg("relerr") = 1e-3,
-            py::arg("dim") = 1000, py::arg("analytical_PggxPgg") = true, py::arg("MCintegral") = false,
+            py::arg("dim") = 1000, py::arg("double_int_method") = "analytical",
             py::arg("MCcalls") = 25000
         )
         .def(
@@ -183,10 +183,10 @@ PYBIND11_MODULE(_core, m) {
         .def(
             py::init<
                 const int &, const char &, const char &, const double &,
-                const double &, const int &, const bool &, const bool &, const int &>(),
+                const double &, const int &, const string &, const int &>(),
             py::arg("order"), py::arg("kind"), py::arg("channel"),
             py::arg("abserr") = 1e-3, py::arg("relerr") = 1e-3,
-            py::arg("dim") = 1000, py::arg("analytical_PggxPgg") = true, py::arg("MCintegral") = false,
+            py::arg("dim") = 1000, py::arg("double_int_method") = "analytical",
             py::arg("MCcalls") = 25000
         )
         .def(
diff --git a/src/ApproximateCoefficientFunction.cc b/src/ApproximateCoefficientFunction.cc
index 15717513..b5812532 100644
--- a/src/ApproximateCoefficientFunction.cc
+++ b/src/ApproximateCoefficientFunction.cc
@@ -19,13 +19,13 @@ using std::vector;
 
 AbstractApproximate::AbstractApproximate(
     const int &order, const char &kind, const char &channel,
-    const double &abserr, const double &relerr, const int &dim, const bool &analytical_PggxPgg,
-    const bool &MCintegral, const int &MCcalls
+    const double &abserr, const double &relerr, const int &dim,
+    const string &double_int_method, const int &MCcalls
 )
     : CoefficientFunction(order, kind, channel) {
 
     muterms_ = new ExactCoefficientFunction(
-        order, kind, channel, abserr, relerr, dim, analytical_PggxPgg, MCintegral, MCcalls
+        order, kind, channel, abserr, relerr, dim, double_int_method, MCcalls
     );
 }
 
@@ -100,10 +100,10 @@ struct variation_parameters CL_var = { 0.2, 2. };
 ApproximateCoefficientFunction::ApproximateCoefficientFunction(
     const int &order, const char &kind, const char &channel, const bool &NLL,
     const string &highscale_version, const double &abserr, const double &relerr,
-    const int &dim, const bool &analytical_PggxPgg, const bool &MCintegral, const int &MCcalls
+    const int &dim, const string &double_int_method, const int &MCcalls
 )
     : AbstractApproximate(
-          order, kind, channel, abserr, relerr, dim, analytical_PggxPgg, MCintegral, MCcalls
+          order, kind, channel, abserr, relerr, dim, double_int_method, MCcalls
       ) {
 
     threshold_ = new ThresholdCoefficientFunction(order, kind, channel);
@@ -266,11 +266,11 @@ struct klmv_params klmv_C2g3B_lowxi = { 0.8, 10.7, 0.055125, 2, 0.3825 };
 ApproximateCoefficientFunctionKLMV::ApproximateCoefficientFunctionKLMV(
     const int &order, const char &kind, const char &channel,
     const string &highscale_version, const bool &lowxi, const double &abserr,
-    const double &relerr, const int &dim, const bool &analytical_PggxPgg, const bool &MCintegral,
+    const double &relerr, const int &dim, const string &double_int_method,
     const int &MCcalls
 )
     : AbstractApproximate(
-          order, kind, channel, abserr, relerr, dim, analytical_PggxPgg, MCintegral, MCcalls
+          order, kind, channel, abserr, relerr, dim, double_int_method, MCcalls
       ) {
     if (GetOrder() == 1) {
         cout << "Error: KLMV approximation is not implemented at O(as)!"
diff --git a/src/ExactCoefficientFunction.cc b/src/ExactCoefficientFunction.cc
index d9529f78..aa162e4d 100644
--- a/src/ExactCoefficientFunction.cc
+++ b/src/ExactCoefficientFunction.cc
@@ -15,8 +15,8 @@ using std::endl;
 
 ExactCoefficientFunction::ExactCoefficientFunction(
     const int &order, const char &kind, const char &channel,
-    const double &abserr, const double &relerr, const int &dim, const bool &analytical_PggxPgg,
-    const bool &MCintegral, const int &MCcalls
+    const double &abserr, const double &relerr, const int &dim,
+    const string &double_int_method, const int &MCcalls
 )
     : CoefficientFunction(order, kind, channel) {
 
@@ -37,6 +37,16 @@ ExactCoefficientFunction::ExactCoefficientFunction(
 
     delta_ = nullptr;
 
+    // check double_int_method
+    if (double_int_method != "analytical"
+        && double_int_method != "double_numerical"
+        && double_int_method != "monte_carlo") {
+        cout << "Error in ExactCoefficientFunction: double_int_method must be "
+                "'analytical', 'double_numerical' or 'monte_carlo'! Got "
+             << double_int_method << endl;
+        exit(-1);
+    }
+
     if (GetOrder() == 2) {
         asy_ = new AsymptoticCoefficientFunction(order, kind, channel);
         thr_ = new ThresholdCoefficientFunction(order, kind, channel);
@@ -63,6 +73,7 @@ ExactCoefficientFunction::ExactCoefficientFunction(
         Pgg1_ = new SplittingFunction(1, 'g', 'g');
         Pqg0_ = new SplittingFunction(0, 'q', 'g');
         Pgq0Pqg0_ = new ConvolutedSplittingFunctions(0, 'g', 'q', 0, 'q', 'g');
+        Pgg0Pgg0_ = new ConvolutedSplittingFunctions(0, 'g', 'g', 0, 'g', 'g');
     }
 
     if (GetOrder() == 2) {
@@ -111,12 +122,18 @@ ExactCoefficientFunction::ExactCoefficientFunction(
             );
             convolutions_lmu1_.push_back(new Convolution(gluon_as2_, delta_));
 
-            if (!analytical_PggxPgg) {
+            if (double_int_method == "monte_carlo") {
+                convolutions_lmu2_.push_back(new DoubleConvolution(
+                    gluon_as1_, Pgg0_, abserr, relerr, dim, true, MCcalls
+                ));
+            } else if (double_int_method == "double_numerical") {
                 convolutions_lmu2_.push_back(new DoubleConvolution(
-                    gluon_as1_, Pgg0_, abserr, relerr, dim, MCintegral, MCcalls
+                    gluon_as1_, Pgg0_, abserr, relerr, dim, false, MCcalls
                 ));
             } else {
-                convolutions_lmu2_.push_back(new Convolution(gluon_as1_, Pgg0Pgg0_, abserr, relerr, dim));
+                convolutions_lmu2_.push_back(
+                    new Convolution(gluon_as1_, Pgg0Pgg0_, abserr, relerr, dim)
+                );
             }
             convolutions_lmu2_.push_back(
                 new Convolution(gluon_as1_, Pgq0Pqg0_, abserr, relerr, dim)
@@ -158,6 +175,7 @@ ExactCoefficientFunction::~ExactCoefficientFunction() {
     delete Pgg1_;
     delete Pqg0_;
     delete Pgq0Pqg0_;
+    delete Pgg0Pgg0_;
 
     delete delta_;
 
diff --git a/src/SplittingFunction.cc b/src/SplittingFunction.cc
index 394b3e5b..71ffd45b 100644
--- a/src/SplittingFunction.cc
+++ b/src/SplittingFunction.cc
@@ -231,7 +231,8 @@ ConvolutedSplittingFunctions::ConvolutedSplittingFunctions(
         exit(-1);
     }
 
-    if (order1_ == 0 && entry1_ == 'g' && entry2_ == 'g' && order2_ == 0 && entry3_ == 'g' && entry4_ == 'g') {
+    if (order1_ == 0 && entry1_ == 'g' && entry2_ == 'g' && order2_ == 0
+        && entry3_ == 'g' && entry4_ == 'g') {
         Pgg0_ = new SplittingFunction(0, 'g', 'g');
     } else {
         Pgg0_ = nullptr;
@@ -240,10 +241,7 @@ ConvolutedSplittingFunctions::ConvolutedSplittingFunctions(
     SetFunctions();
 }
 
-ConvolutedSplittingFunctions::~ConvolutedSplittingFunctions() {
-
-    delete Pgg0_;
-}
+ConvolutedSplittingFunctions::~ConvolutedSplittingFunctions() { delete Pgg0_; }
 
 //==========================================================================================//
 //  ConvolutedSplittingFunctions: regular part
@@ -265,7 +263,8 @@ double ConvolutedSplittingFunctions::Singular(double x, int nf) const {
 //  ConvolutedSplittingFunctions: integral from 0 to x of the singular part
 //------------------------------------------------------------------------------------------//
 
-double ConvolutedSplittingFunctions::SingularIntegrated(double x, int nf) const {
+double
+    ConvolutedSplittingFunctions::SingularIntegrated(double x, int nf) const {
     return GetMultFact() * (this->*sing_int_)(x, nf);
 }
 
@@ -332,22 +331,23 @@ void ConvolutedSplittingFunctions::SetFunctions() {
             sing_int_ = &ConvolutedSplittingFunctions::ZeroFunction_x_nf;
             loc_ = &ConvolutedSplittingFunctions::ZeroFunction_nf;
         } else if (entry1_ == 'g' && entry2_ == 'g' && entry3_ == 'g'
-                 && entry4_ == 'q') {
+                   && entry4_ == 'q') {
             reg_ = &ConvolutedSplittingFunctions::Pgg0_x_Pgq0;
             sing_ = &ConvolutedSplittingFunctions::ZeroFunction_x_nf;
             sing_int_ = &ConvolutedSplittingFunctions::ZeroFunction_x_nf;
             loc_ = &ConvolutedSplittingFunctions::ZeroFunction_nf;
         } else if (entry1_ == 'q' && entry2_ == 'q' && entry3_ == 'g'
-                 && entry4_ == 'q') {
+                   && entry4_ == 'q') {
             reg_ = &ConvolutedSplittingFunctions::Pqq0_x_Pgq0;
             sing_ = &ConvolutedSplittingFunctions::ZeroFunction_x_nf;
             sing_int_ = &ConvolutedSplittingFunctions::ZeroFunction_x_nf;
             loc_ = &ConvolutedSplittingFunctions::ZeroFunction_nf;
         } else if (entry1_ == 'g' && entry2_ == 'g' && entry3_ == 'g'
-                 && entry4_ == 'g') {
+                   && entry4_ == 'g') {
             reg_ = &ConvolutedSplittingFunctions::Pgg0_x_Pgg0_reg;
             sing_ = &ConvolutedSplittingFunctions::Pgg0_x_Pgg0_sing;
-            sing_int_ = &ConvolutedSplittingFunctions::Pgg0_x_Pgg0_sing_integrated;
+            sing_int_ =
+                &ConvolutedSplittingFunctions::Pgg0_x_Pgg0_sing_integrated;
             loc_ = &ConvolutedSplittingFunctions::Pgg0_x_Pgg0_loc;
         } else {
             cout << "Error: P" << entry1_ << entry2_ << order1_ << " x P"
@@ -684,7 +684,10 @@ double ConvolutedSplittingFunctions::Pgq0_x_Pqg0(double x, int nf) const {
 
 double ConvolutedSplittingFunctions::Pgg0reg_x_Pgg0reg(double x) const {
 
-    return 16 * CA * CA *((-1. + x) * (11. + x * (5. + 2. * x)) - 3. * (1. + x * (4. + x + x * x)) * log(x)) / (3. * x);
+    return 16 * CA * CA
+           * ((-1. + x) * (11. + x * (5. + 2. * x))
+              - 3. * (1. + x * (4. + x + x * x)) * log(x))
+           / (3. * x);
 }
 
 //==========================================================================================//
@@ -693,7 +696,10 @@ double ConvolutedSplittingFunctions::Pgg0reg_x_Pgg0reg(double x) const {
 
 double ConvolutedSplittingFunctions::Pgg0reg_x_Pgg0sing(double x) const {
 
-    return 16. * CA * CA * (0.5 * (1. - 4. * x + 3. * x * x) + (-2. + 1. / x + x - x * x) * log(1. - x) + (2. - x + x * x) * log(x));
+    return 16. * CA * CA
+           * (0.5 * (1. - 4. * x + 3. * x * x)
+              + (-2. + 1. / x + x - x * x) * log(1. - x)
+              + (2. - x + x * x) * log(x));
 }
 
 //==========================================================================================//
@@ -702,7 +708,7 @@ double ConvolutedSplittingFunctions::Pgg0reg_x_Pgg0sing(double x) const {
 
 double ConvolutedSplittingFunctions::Pgg0sing_x_Pgg0sing_reg(double x) const {
 
-    return - 16 * CA * CA * log(x) / (1. - x);
+    return -16 * CA * CA * log(x) / (1. - x);
 }
 
 //==========================================================================================//
@@ -718,10 +724,12 @@ double ConvolutedSplittingFunctions::Pgg0sing_x_Pgg0sing_sing(double x) const {
 //  Analytical convolution between the splitting functions Pgg0reg and Pgg0sing
 //------------------------------------------------------------------------------------------//
 
-double ConvolutedSplittingFunctions::Pgg0sing_x_Pgg0sing_sing_integrated(double x) const {
+double ConvolutedSplittingFunctions::Pgg0sing_x_Pgg0sing_sing_integrated(
+    double x
+) const {
 
     double L1 = log(1. - x);
-    return - 16. * CA * CA * L1 * L1 ;
+    return -16. * CA * CA * L1 * L1;
 }
 
 //==========================================================================================//
@@ -730,46 +738,52 @@ double ConvolutedSplittingFunctions::Pgg0sing_x_Pgg0sing_sing_integrated(double
 
 double ConvolutedSplittingFunctions::Pgg0sing_x_Pgg0sing_loc() const {
 
-    return - 16 * CA * CA * zeta2;
+    return -16 * CA * CA * zeta2;
 }
 
 //==========================================================================================//
-//  Regular part of the analytical convolution between the splitting functions Pgg0 and Pgg0
+//  Regular part of the analytical convolution between the splitting functions
+//  Pgg0 and Pgg0
 //------------------------------------------------------------------------------------------//
 
 double ConvolutedSplittingFunctions::Pgg0_x_Pgg0_reg(double x, int nf) const {
 
-    return Pgg0reg_x_Pgg0reg(x)
-           + 2 * Pgg0reg_x_Pgg0sing(x)
-           + 2 * Pgg0_ -> Regular(x, nf) * Pgg0_ -> Local(nf)
+    return Pgg0reg_x_Pgg0reg(x) + 2 * Pgg0reg_x_Pgg0sing(x)
+           + 2 * Pgg0_->Regular(x, nf) * Pgg0_->Local(nf)
            + Pgg0sing_x_Pgg0sing_reg(x);
 }
 
 //==========================================================================================//
-//  Regular part of the analytical convolution between the splitting functions Pgg0 and Pgg0
+//  Regular part of the analytical convolution between the splitting functions
+//  Pgg0 and Pgg0
 //------------------------------------------------------------------------------------------//
 
 double ConvolutedSplittingFunctions::Pgg0_x_Pgg0_sing(double x, int nf) const {
 
     return Pgg0sing_x_Pgg0sing_sing(x)
-           + 2 * Pgg0_ -> Singular(x, nf) * Pgg0_ -> Local(nf);
+           + 2 * Pgg0_->Singular(x, nf) * Pgg0_->Local(nf);
 }
 
 //==========================================================================================//
-//  Regular part of the analytical convolution between the splitting functions Pgg0 and Pgg0
+//  Regular part of the analytical convolution between the splitting functions
+//  Pgg0 and Pgg0
 //------------------------------------------------------------------------------------------//
 
-double ConvolutedSplittingFunctions::Pgg0_x_Pgg0_sing_integrated(double x, int nf) const {
+double ConvolutedSplittingFunctions::Pgg0_x_Pgg0_sing_integrated(
+    double x, int nf
+) const {
 
     return Pgg0sing_x_Pgg0sing_sing_integrated(x)
-           + 2 * Pgg0_ -> SingularIntegrated(x, nf) * Pgg0_ -> Local(nf);
+           + 2 * Pgg0_->SingularIntegrated(x, nf) * Pgg0_->Local(nf);
 }
 
 //==========================================================================================//
-//  Regular part of the analytical convolution between the splitting functions Pgg0 and Pgg0
+//  Regular part of the analytical convolution between the splitting functions
+//  Pgg0 and Pgg0
 //------------------------------------------------------------------------------------------//
 
 double ConvolutedSplittingFunctions::Pgg0_x_Pgg0_loc(int nf) const {
 
-    return Pgg0sing_x_Pgg0sing_loc() + Pgg0_ -> Local(nf);
+    double loc = Pgg0_->Local(nf);
+    return Pgg0sing_x_Pgg0sing_loc() + loc * loc;
 }
diff --git a/tests/test_approx_coeff_func.py b/tests/test_approx_coeff_func.py
index 7c83e92e..73ae9376 100644
--- a/tests/test_approx_coeff_func.py
+++ b/tests/test_approx_coeff_func.py
@@ -20,7 +20,7 @@ def test_mudependent_terms():
             for kind in ['2', 'L']:
                 for dim in [1000]:
                     for MCcalls in [20000]:
-                        for mf in [False, True]:
+                        for mf in ["analytical", "double_numerical", "monte_carlo"]:
                             for abserr in [1e-3]:
                                 relerr = abserr
                                 massive = ad.ExactCoefficientFunction(order, kind, channel, abserr, relerr, dim, mf, MCcalls)
diff --git a/tests/test_mc_integral.py b/tests/test_mc_integral.py
index 8b3e80d2..a88fdbe9 100644
--- a/tests/test_mc_integral.py
+++ b/tests/test_mc_integral.py
@@ -5,8 +5,9 @@
 
 def test_mc_integrals():
     for kind in ['2', 'L']:
-        massive_mc = ad.ExactCoefficientFunction(3, kind, 'g', 1e-3, 1e-3, 1000, True, 1000000)
-        massive_nomc = ad.ExactCoefficientFunction(3, kind, 'g', 1e-3, 1e-3, 1000, False, 25000)
+        massive_mc = ad.ExactCoefficientFunction(3, kind, 'g', 1e-3, 1e-3, 1000, "monte_carlo", 1000000)
+        massive_nomc = ad.ExactCoefficientFunction(3, kind, 'g', 1e-3, 1e-3, 1000, "double_numerical", 25000)
+        massive_analytical = ad.ExactCoefficientFunction(3, kind, 'g', 1e-3, 1e-3, 1000, "analytical")
         for xi in np.geomspace(1e-2, 1e4, 4, endpoint=True):
             m2Q2 = 1/xi
             xmax = 1. / (1. + 4 * m2Q2)
@@ -14,4 +15,6 @@ def test_mc_integrals():
                     for nf in range(1, 6 + 1):
                         res1 = massive_mc.MuDependentTerms(x, m2Q2, m2Q2, nf)
                         res2 = massive_nomc.MuDependentTerms(x, m2Q2, m2Q2, nf)
+                        res3 = massive_analytical.MuDependentTerms(x, m2Q2, m2Q2, nf)
                         np.testing.assert_allclose(res1, res2, rtol=5e-3)
+                        np.testing.assert_allclose(res2, res3, 1e-3)