diff --git a/poincare/src/integral.cpp b/poincare/src/integral.cpp
index 299eccf44..77e59259f 100644
--- a/poincare/src/integral.cpp
+++ b/poincare/src/integral.cpp
@@ -113,72 +113,81 @@ IntegralNode::DetailedResult<T> IntegralNode::kronrodGaussQuadrature(T a, T b, C
   static T max = sizeof(T) == sizeof(double) ? DBL_MAX : FLT_MAX;
   /* We here use Kronrod-Legendre quadrature with n = 21
    * The abscissa and weights are taken from QUADPACK library. */
-  const static T wg[5]= {0.066671344308688137593568809893332, 0.149451349150580593145776339657697,
-    0.219086362515982043995534934228163, 0.269266719309996355091226921569469, 0.295524224714752870173892994651338};
-  const static T xgk[11]= {0.995657163025808080735527280689003, 0.973906528517171720077964012084452,
+
+  // Abscissae for the gauss (odd weights) and kronrod rules (all weights)
+  const static T x[11]= {0.995657163025808080735527280689003, 0.973906528517171720077964012084452,
     0.930157491355708226001207180059508, 0.865063366688984510732096688423493, 0.780817726586416897063717578345042,
     0.679409568299024406234327365114874, 0.562757134668604683339000099272694, 0.433395394129247190799265943165784,
     0.294392862701460198131126603103866, 0.148874338981631210884826001129720, 0.000000000000000000000000000000000};
-  const static T wgk[11]= {0.011694638867371874278064396062192, 0.032558162307964727478818972459390,
+
+  // Weights for the gauss integral
+  const static T wGauss[5]= {0.066671344308688137593568809893332, 0.149451349150580593145776339657697,
+    0.219086362515982043995534934228163, 0.269266719309996355091226921569469, 0.295524224714752870173892994651338};
+
+  // Weights for the kronrod rule
+  const static T wKronrod[11]= {0.011694638867371874278064396062192, 0.032558162307964727478818972459390,
     0.054755896574351996031381300244580, 0.075039674810919952767043140916190, 0.093125454583697605535065465083366,
     0.109387158802297641899210590325805, 0.123491976262065851077958109831074, 0.134709217311473325928054001771707,
     0.142775938577060080797094273138717, 0.147739104901338491374841515972068, 0.149445554002916905664936468389821};
+
   T fv1[10];
   T fv2[10];
 
-  T centr = 0.5*(a+b);
-  T hlgth = 0.5*(b-a);
-  T dhlgth = std::fabs(hlgth);
+  T center = 0.5 * (a+b);
+  T halfLength = 0.5 * (b-a);
+  T absHalfLength = std::fabs(halfLength);
 
   DetailedResult<T> errorResult;
   errorResult.integral = NAN;
   errorResult.absoluteError = 0;
 
-  T resg = 0;
-  T fc = functionValueAtAbscissa(centr, context, complexFormat, angleUnit);
-  if (std::isnan(fc)) {
+  T gaussIntegral = 0;
+  T fCenter = functionValueAtAbscissa(center, context, complexFormat, angleUnit);
+  if (std::isnan(fCenter)) {
     return errorResult;
   }
-  T resk = wgk[10]*fc;
-  T resabs = std::fabs(resk);
+  T kronrodIntegral = wKronrod[10] * fCenter;
+  T absKronrodIntegral = std::fabs(kronrodIntegral);
   for (int j = 0; j < 10; j++) {
-    T absc = hlgth*xgk[j];
-    T fval1 = functionValueAtAbscissa(centr-absc, context, complexFormat, angleUnit);
+    T xDelta = halfLength * x[j];
+    T fval1 = functionValueAtAbscissa(center - xDelta, context, complexFormat, angleUnit);
     if (std::isnan(fval1)) {
       return errorResult;
     }
-    T fval2 = functionValueAtAbscissa(centr+absc, context, complexFormat, angleUnit);
+    T fval2 = functionValueAtAbscissa(center + xDelta, context, complexFormat, angleUnit);
     if (std::isnan(fval2)) {
       return errorResult;
     }
     fv1[j] = fval1;
     fv2[j] = fval2;
-    T fsum = fval1+fval2;
-    if (j%2 == 1) {
-      resg += wg[j/2]*fsum;
+    T fsum = fval1 + fval2;
+    if (j % 2 == 1) {
+      gaussIntegral += wGauss[j/2] * fsum;
     }
-    resk += wgk[j]*fsum;
-    resabs += wgk[j]*(std::fabs(fval1)+std::fabs(fval2));
+    kronrodIntegral += wKronrod[j] * fsum;
+    absKronrodIntegral += wKronrod[j] * (std::fabs(fval1) + std::fabs(fval2));
   }
 
-  T reskh = resk*0.5;
-  T resasc = wgk[10]*std::fabs(fc-reskh);
+  T halfKronrodIntegral = 0.5 * kronrodIntegral;
+  T kronrodIntegralDifference = wKronrod[10] * std::fabs(fCenter - halfKronrodIntegral);
   for (int j = 0; j < 10; j++) {
-    resasc += wgk[j]*(std::fabs(fv1[j]-reskh)+std::fabs(fv2[j]-reskh));
+    kronrodIntegralDifference += wKronrod[j] * (std::fabs(fv1[j] - halfKronrodIntegral) + std::fabs(fv2[j] - halfKronrodIntegral));
   }
-  T integral = resk*hlgth;
-  resabs = resabs*dhlgth;
-  resasc = resasc*dhlgth;
-  T abserr = std::fabs((resk-resg)*hlgth);
-  if (resasc != 0 && abserr != 0) {
-    abserr = 1 > std::pow((T)(200*abserr/resasc), (T)1.5)? resasc*std::pow((T)(200*abserr/resasc), (T)1.5) : resasc;
+  T integral = kronrodIntegral * halfLength;
+  absKronrodIntegral = absKronrodIntegral * absHalfLength;
+  kronrodIntegralDifference = kronrodIntegralDifference * absHalfLength;
+  T absError = std::fabs((kronrodIntegral - gaussIntegral) * halfLength);
+  if (kronrodIntegralDifference != 0 && absError != 0) {
+    T errorCoefficient = std::pow((T)(200*absError/kronrodIntegralDifference), (T)1.5);
+    absError = 1 > errorCoefficient ? kronrodIntegralDifference * errorCoefficient : kronrodIntegralDifference;
   }
-  if (resabs > max/(50.0*epsilon)) {
-    abserr = abserr > epsilon*50*resabs ? abserr : epsilon*50*resabs;
+  if (absKronrodIntegral > max/(50.0 * epsilon)) {
+    T minError = epsilon * 50 * absKronrodIntegral;
+    absError = absError > minError ? absError : minError;
   }
   DetailedResult<T> result;
   result.integral = integral;
-  result.absoluteError = abserr;
+  result.absoluteError = absError;
   return result;
 }