         lsup_cst[n] = res_lsup_cst + lsup_cst[n];
         uinf_cst[n] = res_uinf_cst + uinf_cst[n];
         usup_cst[n] = res_usup_cst + usup_cst[n];
-    }
-}
-
-
-__global__
-void lexpr_replace_relu_bounds_conv_sparse(float_type* __restrict__ inf_coeff, float_type* __restrict__ sup_coeff, float_type* __restrict__ inf_cst, float_type* __restrict__ sup_cst, const float_type* __restrict__ lb_array, const float_type* __restrict__ ub_array, const size_t output_size_x, const size_t output_size_y, const size_t output_size_z, long int offset_x, long int offset_y, long int length_x, long int length_y, long int shift_x, long int shift_y, const bool use_area_heuristic)
-{
-    const size_t last_x = blockIdx.x;
-    const size_t last_y = blockIdx.y;
-    const size_t last_z = blockIdx.z;
-
-    size_t out_z = threadIdx.x;
-
-    __shared__ size_t n;
-
-    __shared__ long int min_out_x;
-    __shared__ long int min_out_y;
-
-    __shared__ long int min_x;
-    __shared__ long int min_y;
-
-    __shared__ long int max_x;
-    __shared__ long int max_y;
-
-    if(out_z == 0)
-    {
-        n = last_x*gridDim.y*gridDim.z + last_y*gridDim.z + last_z;
-
-        min_out_x = offset_x + last_x*shift_x;
-        min_out_y = offset_y + last_y*shift_y;
-
-        min_x = (min_out_x < 0) ? -min_out_x : 0;
-        min_y = (min_out_y < 0) ? -min_out_y : 0;
-
-        max_x = (length_x + min_out_x > output_size_x) ? output_size_x - min_out_x : length_x;
-        max_y = (length_y + min_out_y > output_size_y) ? output_size_y - min_out_y : length_y;
-    }
-
-    __syncthreads();
-
-    float_type res_inf_cst = 0;
-    float_type res_sup_cst = 0;
-
-    while(out_z < output_size_z)
-    {
-        for(long int out_x = min_x; out_x < max_x; out_x++)
-        {
-            for(long int out_y = min_y; out_y < max_y; out_y++)
-            {
-                size_t i = (out_x + min_out_x)*output_size_y*output_size_z + (out_y + min_out_y)*output_size_z + out_z;
-                size_t j = out_x*length_y*output_size_z + out_y*output_size_z + out_z;
-
-                size_t a = n*length_x*length_y*output_size_z + j;
-
-                const float_type lb = lb_array[i];
-                const float_type ub = ub_array[i];
-                const float_type width = ub - lb;
-                const float_type lambda_inf = ub/width;
-                const float_type lambda_sup = ub/width;
-
-                const float_type old_inf_coeff = inf_coeff[a];
-                const float_type old_sup_coeff = sup_coeff[a];
-
-                if((old_sup_coeff == 0) && (old_inf_coeff == 0))
-                {
-                    inf_coeff[a] = 0.0;
-                    sup_coeff[a] = 0.0;
-                }
-                else if(ub <= 0)
-                {
-                    inf_coeff[a] = 0.0;
-                    sup_coeff[a] = 0.0;
-                }
-                else if(lb > 0)
-                {
-                    inf_coeff[a] = old_inf_coeff;
-                    sup_coeff[a] = old_sup_coeff;
-                }
-                else if(old_sup_coeff < 0)
-                {
-                    const float_type mu_inf = -lambda_inf*lb;
-                    const float_type mu_sup = -lambda_sup*lb;
-                    elina_double_interval_mul_expr_coeff(&inf_coeff[a], &sup_coeff[a], lambda_inf, lambda_sup, old_inf_coeff, old_sup_coeff);
-                    float_type tmp1, tmp2;
-                    elina_double_interval_mul_cst_coeff(&tmp1, &tmp2, mu_inf, mu_sup, old_inf_coeff, old_sup_coeff);
-
-                    res_inf_cst += tmp1 - min_denormal;
-                    res_sup_cst += tmp2 + min_denormal;
-                }
-                else if (old_inf_coeff > 0)
-                {
-                    if(use_area_heuristic)
-                    {
-                        const float_type area1 = -lb*ub;
-                        const float_type area2 = 0.5*ub*width;
-                        const float_type area3 = -0.5*lb*width;
-
-                        if((area1 < area2) && (area1 < area3))
-                        {
-                            elina_double_interval_mul_expr_coeff(&inf_coeff[a], &sup_coeff[a], lambda_inf, lambda_sup, old_inf_coeff, old_sup_coeff);
-                        }
-                        else if((area2 < area1) && (area2 < area3))
-                        {
-                            inf_coeff[a] = 0.0;
-                            sup_coeff[a] = 0.0;
-                        }
-                        else
-                        {
-                            inf_coeff[a] = old_inf_coeff;
-                            sup_coeff[a] = old_sup_coeff;
-                        }
-                    }
-                    else
-                    {
-                        inf_coeff[a] = 0.0;
-                        sup_coeff[a] = 0.0;
-                    }
-                }
-                else
-                {
-                    inf_coeff[a] = 0.0;
-                    sup_coeff[a] = 0.0;
-                    float_type tmp1, tmp2;
-                    elina_double_interval_mul2(&tmp1, &tmp2, old_inf_coeff, old_sup_coeff, 0, ub);
-
-                    res_inf_cst += tmp1;
-                    res_sup_cst += tmp1;
-                }
-            }
-        }
-
-        out_z += blockDim.x;
-    }
-
-    block_reduce_sum(res_inf_cst, blockDim.x);
-    block_reduce_sum(res_sup_cst, blockDim.x);
-
-    if(threadIdx.x == 0)
-    {
-        inf_cst[n] = res_inf_cst + inf_cst[n];
-        sup_cst[n] = res_sup_cst + sup_cst[n];
-    }
-}
-
-
-__global__
-void uexpr_replace_relu_bounds_conv_sparse(float_type* __restrict__ inf_coeff, float_type* __restrict__ sup_coeff, float_type* __restrict__ inf_cst, float_type* __restrict__ sup_cst, const float_type* __restrict__ lb_array, const float_type* __restrict__ ub_array, const size_t output_size_x, const size_t output_size_y, const size_t output_size_z, long int offset_x, long int offset_y, long int length_x, long int length_y, long int shift_x, long int shift_y, const bool use_area_heuristic)
-{
-    const size_t last_x = blockIdx.x;
-    const size_t last_y = blockIdx.y;
-    const size_t last_z = blockIdx.z;
-
-    size_t out_z = threadIdx.x;
-
-    __shared__ size_t n;
-
-    __shared__ long int min_out_x;
-    __shared__ long int min_out_y;
-
-    __shared__ long int min_x;
-    __shared__ long int min_y;
-
-    __shared__ long int max_x;
-    __shared__ long int max_y;
-
-    if(out_z == 0)
-    {
-        n = last_x*gridDim.y*gridDim.z + last_y*gridDim.z + last_z;
-
-        min_out_x = offset_x + last_x*shift_x;
-        min_out_y = offset_y + last_y*shift_y;
-
-        min_x = (min_out_x < 0) ? -min_out_x : 0;
-        min_y = (min_out_y < 0) ? -min_out_y : 0;
-
-        max_x = (length_x + min_out_x > output_size_x) ? output_size_x - min_out_x : length_x;
-        max_y = (length_y + min_out_y > output_size_y) ? output_size_y - min_out_y : length_y;
-    }
-
-    __syncthreads();
-
-    float_type res_inf_cst = 0;
-    float_type res_sup_cst = 0;
-
-    while(out_z < output_size_z)
-    {
-        for(long int out_x = min_x; out_x < max_x; out_x++)
-        {
-            for(long int out_y = min_y; out_y < max_y; out_y++)
-            {
-                size_t i = (out_x + min_out_x)*output_size_y*output_size_z + (out_y + min_out_y)*output_size_z + out_z;
-                size_t j = out_x*length_y*output_size_z + out_y*output_size_z + out_z;
-
-                size_t a = n*length_x*length_y*output_size_z + j;
-
-                const float_type lb = lb_array[i];
-                const float_type ub = ub_array[i];
-                const float_type width = ub - lb;
-                const float_type lambda_inf = ub/width;
-                const float_type lambda_sup = ub/width;
-
-                const float_type old_inf_coeff = inf_coeff[a];
-                const float_type old_sup_coeff = sup_coeff[a];
-
-                if((old_sup_coeff == 0) && (old_inf_coeff == 0))
-                {
-                    inf_coeff[a] = 0.0;
-                    sup_coeff[a] = 0.0;
-                }
-                else if(ub <= 0)
-                {
-                    inf_coeff[a] = 0.0;
-                    sup_coeff[a] = 0.0;
-                }
-                else if(lb > 0)
-                {
-                    inf_coeff[a] = old_inf_coeff;
-                    sup_coeff[a] = old_sup_coeff;
-                }
-                else if(old_inf_coeff > 0)
-                {
-                    const float_type mu_inf = -lambda_inf*lb;
-                    const float_type mu_sup = -lambda_sup*lb;
-                    elina_double_interval_mul_expr_coeff(&inf_coeff[a], &sup_coeff[a], lambda_inf, lambda_sup, old_inf_coeff, old_sup_coeff);
-                    float_type tmp1, tmp2;
-                    elina_double_interval_mul_cst_coeff(&tmp1, &tmp2, mu_inf, mu_sup, old_inf_coeff, old_sup_coeff);
-
-                    res_inf_cst += tmp1 - min_denormal;
-                    res_sup_cst += tmp2 + min_denormal;
-                }
-                else if(old_sup_coeff < 0)
-                {
-                    if(use_area_heuristic)
-                    {
-                        const float_type area1 = -lb*ub;
-                        const float_type area2 = 0.5*ub*width;
-                        const float_type area3 = -0.5*lb*width;
-
-                        if((area1 < area2) && (area1 < area3))
-                        {
-                            elina_double_interval_mul_expr_coeff(&inf_coeff[a], &sup_coeff[a], lambda_inf, lambda_sup, old_inf_coeff, old_sup_coeff);
-                        }
-                        else if((area2 < area1) && (area2 < area3))
-                        {
-                            inf_coeff[a] = 0.0;
-                            sup_coeff[a] = 0.0;
-                        }
-                        else
-                        {
-                            inf_coeff[a] = old_inf_coeff;
-                            sup_coeff[a] = old_sup_coeff;
-                        }
-                    }
-                    else
-                    {
-                        inf_coeff[a] = 0.0;
-                        sup_coeff[a] = 0.0;
-                    }
-                }
-                else
-                {
-                    inf_coeff[a] = 0.0;
-                    sup_coeff[a] = 0.0;
-                    float_type tmp1, tmp2;
-                    elina_double_interval_mul2(&tmp1, &tmp2, old_inf_coeff, old_sup_coeff, 0, ub);
-
-                    res_inf_cst += tmp2;
-                    res_sup_cst += tmp2;
-                }
-            }
-        }
-
-        out_z += blockDim.x;
-    }
-
-    block_reduce_sum(res_inf_cst, blockDim.x);
-    block_reduce_sum(res_sup_cst, blockDim.x);
-
-    if(threadIdx.x == 0)
-    {
-        inf_cst[n] = res_inf_cst + inf_cst[n];
-        sup_cst[n] = res_sup_cst + sup_cst[n];
     }
 }
 

        res_lsup_cst[n] = lsup_cst + expr_lsup_cst[n];
        res_uinf_cst[n] = uinf_cst + expr_uinf_cst[n];
        res_usup_cst[n] = usup_cst + expr_usup_cst[n];
    }
}


__global__
void csts_from_previous_layer_conv(const float_type* __restrict__ expr_inf_coeff, const float_type* __restrict__ expr_sup_coeff, const float_type* __restrict__ expr_inf_cst, const float_type* __restrict__ expr_sup_cst, float_type* __restrict__ res_inf_cst, float_type* __restrict__ res_sup_cst, const float_type* __restrict__ aux_csts, const size_t current_layer_out_size_x, const size_t current_layer_out_size_y, const size_t current_layer_out_size_z)
{
    const size_t n = blockIdx.x;
    size_t j = threadIdx.x;

    float_type inf_cst = 0;
    float_type sup_cst = 0;

    float_type tmp1, tmp2;
    float_type maxRes, maxMul;

    while(j < current_layer_out_size_z)
    {
        for(size_t i = 0; i < current_layer_out_size_x*current_layer_out_size_y; i++)
        {
            size_t a = n*current_layer_out_size_x*current_layer_out_size_y*current_layer_out_size_z + i*current_layer_out_size_z + j;

            elina_double_interval_mul_cst_coeff_const_expr(&tmp1, &tmp2, expr_inf_coeff[a], expr_sup_coeff[a], aux_csts[j]);

            maxRes = max(abs(inf_cst), abs(sup_cst));
            maxMul = max(abs(tmp1), abs(tmp2));

            inf_cst += tmp1 - (maxRes + maxMul)*ulp - min_denormal;
            sup_cst += tmp2 + (maxRes + maxMul)*ulp + min_denormal;
        }

        j += blockDim.x;
    }

    block_reduce_sum_csts(inf_cst, sup_cst, blockDim.x);

    if(threadIdx.x == 0)
    {
        res_inf_cst[n] = inf_cst + expr_inf_cst[n];
        res_sup_cst[n] = sup_cst + expr_sup_cst[n];
    }
}


    if(fp->layers[k]->activation == RELU)
    {
        expr_replace_relu_bounds_conv_sparse<<<dim3(fp->layers[layerno]->output_size[0], fp->layers[layerno]->output_size[1], fp->layers[layerno]->output_size[2]/num_chunks), fp->layers[k]->output_size[2]>>>(*linf_coeff, *lsup_coeff, *uinf_coeff, *usup_coeff, *linf_cst, *lsup_cst, *uinf_cst, *usup_cst, aux_lb_array, aux_ub_array, fp->layers[k]->output_size[0], fp->layers[k]->output_size[1], fp->layers[k]->output_size[2], offset_x, offset_y, length_x, length_y, shift_x, shift_y, use_area_heuristic);
    }

    size_t missing_length = ((length_x - 1)*fp->layers[k]->strides[0] + fp->layers[k]->filter_size[0])*((length_y - 1)*fp->layers[k]->strides[1] + fp->layers[k]->filter_size[1])*fp->layers[k]->input_size[2];
