             for(long int out_y = min_y; out_y < max_y; out_y++)
             {
                 size_t mat_out = n*length_x*length_y*output_size_z + out_x*length_y*output_size_z + out_y*output_size_z + out_z;
-
-                const float_type prev_linf_coeff = expr_linf_coeff[mat_out];
-                const float_type prev_lsup_coeff = expr_lsup_coeff[mat_out];
-
-                elina_double_interval_mul_cst_coeff_const_expr(&tmp1, &tmp2, prev_linf_coeff, prev_lsup_coeff, aux_csts[out_z]);
+                const float_type aux_cst = aux_csts[out_z];
+
+                elina_double_interval_mul_cst_coeff_const_expr(&tmp1, &tmp2, expr_linf_coeff[mat_out], expr_lsup_coeff[mat_out], aux_cst);
 
                 maxRes = max(abs(linf_cst), abs(lsup_cst));
                 maxMul = max(abs(tmp1), abs(tmp2));

                linf_cst += tmp1 - (maxRes + maxMul)*ulp - min_denormal;
                lsup_cst += tmp2 + (maxRes + maxMul)*ulp + min_denormal;

                const float_type prev_uinf_coeff = expr_uinf_coeff[mat_out];
                const float_type prev_usup_coeff = expr_usup_coeff[mat_out];

                elina_double_interval_mul_cst_coeff_const_expr(&tmp1, &tmp2, prev_uinf_coeff, prev_usup_coeff, aux_csts[out_z]);

                maxRes = max(abs(uinf_cst), abs(usup_cst));
                maxMul = max(abs(tmp1), abs(tmp2));

        out_z += blockDim.x;
    }

    block_reduce_sum_csts(linf_cst, lsup_cst, blockDim.x);
    block_reduce_sum_csts(uinf_cst, usup_cst, blockDim.x);

    if(threadIdx.x == 0)
    {
