         uexpr[i] = copy_expr(out_neurons[i]->expr);
     }
 
-    for(size_t i = 0; i < num_out_neurons; i++)
-    {
-        expr_t** lexpr_ptr = &lexpr[i];
-        expr_t** uexpr_ptr = &uexpr[i];
-
-        for(int k = layerno - 1; k >= 0; k--)
-        {
-            neuron_t** aux_neurons = fp->layers[k]->neurons;
-            if(fp->layers[k]->type != MAXPOOL)
-            {
-                if(fp->layers[k]->activation == RELU)
-                {
-                    lexpr_replace_relu_bounds(lexpr_ptr, aux_neurons);
-                    uexpr_replace_relu_bounds(uexpr_ptr, aux_neurons);
-                }
-
-                expr_from_previous_layer(lexpr_ptr, aux_neurons);
-                expr_from_previous_layer(uexpr_ptr, aux_neurons);
-            }
-            else
-            {
-                lexpr_replace_maxpool_bounds(lexpr_ptr, aux_neurons);
-                uexpr_replace_maxpool_bounds(uexpr_ptr, aux_neurons);
-            }
-        }
-
+    for(int k = layerno - 1; k >= 0; k--)
+    {
+        neuron_t** aux_neurons = fp->layers[k]->neurons;
+        if(fp->layers[k]->type != MAXPOOL)
+        {
+            if(fp->layers[k]->activation == RELU)
+            {
+                lexpr_replace_relu_bounds(lexpr, aux_neurons, num_out_neurons);
+                uexpr_replace_relu_bounds(uexpr, aux_neurons, num_out_neurons);
+            }
+
+            expr_from_previous_layer(lexpr, aux_neurons, num_out_neurons);
+            expr_from_previous_layer(uexpr, aux_neurons, num_out_neurons);
+        }
+        else
+        {
+            lexpr_replace_maxpool_bounds(lexpr, aux_neurons, num_out_neurons);
+            uexpr_replace_maxpool_bounds(uexpr, aux_neurons, num_out_neurons);
+        }
     }
 
     for(size_t i = 0; i < num_out_neurons; i++)



__device__ __host__
void lexpr_replace_maxpool_bounds(expr_t** expr_ptr, neuron_t** neurons)
{
    //printf("begin\n");
    //fflush(stdout);
    expr_t* expr = *expr_ptr;
    const size_t num_neurons = expr->size;
    expr_t* res;

    size_t k;

    if(expr->type == DENSE)
    {
        k = 0;
    }
    else
    {
        k = expr->dim[0];
    }

    const neuron_t* neuron_k = neurons[k];

    if(expr->sup_coeff[0] < 0)
    {
        //expr_print(neuron_k->maxpool_uexpr);
        if(neuron_k->maxpool_uexpr == nullptr)
        {
            res = (expr_t*) malloc(sizeof(expr_t));
            res->inf_coeff = res->sup_coeff =  nullptr;
            res->dim = nullptr;
            res->size = 0;
            res->type = SPARSE;
            elina_double_interval_mul_cst_coeff(&res->inf_cst, &res->sup_cst, neuron_k->lb, neuron_k->ub, expr->inf_coeff[0], expr->sup_coeff[0]);
        }
        else
        {
            res = multiply_expr(neuron_k->maxpool_uexpr, expr->inf_coeff[0], expr->sup_coeff[0]);
        }
        //printf("multiply end %zu \n",k);
        //expr_print(res);
        //fflush(stdout);
    }
    else if(expr->inf_coeff[0] < 0)
    {
        //expr_print(neuron_k->maxpool_lexpr);
        if(neuron_k->maxpool_lexpr == nullptr)
        {
            res = (expr_t*) malloc(sizeof(expr_t));
            res->inf_coeff = nullptr;
            res->sup_coeff = nullptr;
            res->dim = nullptr;
            res->size = 0;
            res->type = SPARSE;
            elina_double_interval_mul_cst_coeff(&res->inf_cst, &res->sup_cst, neuron_k->lb, neuron_k->ub, expr->inf_coeff[0], expr->sup_coeff[0]);
        }
        else
        {
            res = multiply_expr(neuron_k->maxpool_lexpr, expr->inf_coeff[0], expr->sup_coeff[0]);
        }
        //printf("multiply end %zu \n",k);
        //expr_print(res);
        //fflush(stdout);
    }
    else
    {
        //printf("WTF1\n");
        //fflush(stdout);
        double tmp1, tmp2;
        elina_double_interval_mul_cst_coeff(&tmp1, &tmp2, neuron_k->lb, neuron_k->ub, expr->inf_coeff[0], expr->sup_coeff[0]);
        double coeff[1];
        size_t dim[1];
        coeff[0] = 0;
        dim[0] = 0;
        res = create_sparse_expr(coeff, -tmp1, dim, 1);
    }
    //printf("middle\n");
    //fflush(stdout);
    expr_t* mul_expr;

    for(size_t i = 1; i < num_neurons; i++)
    {
        if(expr->type == DENSE)
        {
            k = i;
        }
        else
        {
            k = expr->dim[i];
        }

        const neuron_t* neuron_k = neurons[k];

        if(expr->sup_coeff[i] < 0)
        {
            //expr_print(neuron_k->maxpool_uexpr);
            //printf("add start %zu %zu\n",k,i);

                //expr_print(res);

            if(neuron_k->maxpool_uexpr == nullptr)
            {
                mul_expr = (expr_t*) malloc(sizeof(expr_t));
                mul_expr->inf_coeff = nullptr;
                mul_expr->sup_coeff = nullptr;
                mul_expr->dim = nullptr;
                mul_expr->size = 0;
                mul_expr->type = SPARSE;
                //printf("lb: %g %g\n");
                elina_double_interval_mul_cst_coeff(&mul_expr->inf_cst, &mul_expr->sup_cst, neuron_k->lb, neuron_k->ub, expr->inf_coeff[i], expr->sup_coeff[i]);
                res->inf_cst += mul_expr->inf_cst;
                res->sup_cst += mul_expr->sup_cst;
            }
            else
            {
                mul_expr = multiply_expr(neuron_k->maxpool_uexpr, expr->inf_coeff[i], expr->sup_coeff[i]);
                add_expr(res, mul_expr);
            }
            //expr_print(mul_expr);
                //fflush(stdout);
            //printf("add finish\n");
                //expr_print(res);
                //fflush(stdout);
            free_expr(mul_expr);
        }
        else if (expr->inf_coeff[i] < 0)
        {
            //expr_print(neuron_k->maxpool_lexpr);
            //printf("add start %zu %zu\n",k,i);

                //expr_print(res);

            if(neuron_k->maxpool_lexpr == nullptr)
            {
                mul_expr = (expr_t*) malloc(sizeof(expr_t));
                mul_expr->inf_coeff = nullptr;
                mul_expr->sup_coeff = nullptr;
                mul_expr->dim = nullptr;
                mul_expr->size = 0;
                mul_expr->type = SPARSE;
                elina_double_interval_mul_cst_coeff(&mul_expr->inf_cst, &mul_expr->sup_cst, neuron_k->lb, neuron_k->ub, expr->inf_coeff[i], expr->sup_coeff[i]);
                res->inf_cst += mul_expr->inf_cst;
                res->sup_cst += mul_expr->sup_cst;
            }
            else
            {
                mul_expr = multiply_expr(neuron_k->maxpool_lexpr, expr->inf_coeff[i], expr->sup_coeff[i]);
                //printf("add start1 %zu %zu\n",k,i);
                //expr_print(res);
                //expr_print(mul_expr);
                //fflush(stdout);
                add_expr(res, mul_expr);
            }
            //expr_print(mul_expr);
            //    fflush(stdout);
            //printf("add finish1\n");
            //expr_print(res);
            //fflush(stdout);
            free_expr(mul_expr);
        }
        else
        {
            //printf("WTF2\n");
            //fflush(stdout);
            double tmp1, tmp2;
            elina_double_interval_mul_expr_coeff(&tmp1, &tmp2, neuron_k->lb, neuron_k->ub, expr->inf_coeff[i], expr->sup_coeff[i]);
            res->inf_cst = res->inf_cst + tmp1;
            res->sup_cst = res->sup_cst - tmp1;
        }
    }
    //printf("finish\n");
    //fflush(stdout);
    res->inf_cst = res->inf_cst + expr->inf_cst;
    res->sup_cst = res->sup_cst + expr->sup_cst;

    free_expr(expr);
    *expr_ptr = res;
}


__device__ __host__
void uexpr_replace_maxpool_bounds(expr_t** expr_ptr, neuron_t** neurons)
{
    expr_t* expr = *expr_ptr;
    const size_t num_neurons = expr->size;
    expr_t* res;

    size_t k;

    if(expr->type == DENSE)
    {
        k = 0;
    }
    else
    {
        k = expr->dim[0];
    }

    const neuron_t* neuron_k = neurons[k];

    if(expr->sup_coeff[0] < 0)
    {
        res = multiply_expr(neuron_k->maxpool_lexpr, expr->inf_coeff[0], expr->sup_coeff[0]);
    }
    else if(expr->inf_coeff[0] < 0)
    {
        res = multiply_expr(neuron_k->maxpool_uexpr, expr->inf_coeff[0], expr->sup_coeff[0]);
    }
    else
    {
        double tmp1, tmp2;
        elina_double_interval_mul_cst_coeff(&tmp1, &tmp2, neuron_k->lb, neuron_k->ub, expr->inf_coeff[0], expr->sup_coeff[0]);
        double coeff[1];
        size_t dim[1];
        coeff[0] = 0;
        dim[0] = 0;
        res = create_sparse_expr(coeff, tmp2, dim, 1);
    }

    for(size_t i = 1; i < num_neurons; i++)
    {
        if(expr->type == DENSE)
        {
            k = i;
        }
        else
        {
            k = expr->dim[i];
        }

        const neuron_t* neuron_k = neurons[k];

        if(expr->sup_coeff[i] < 0)
        {
            expr_t* mul_expr = multiply_expr(neuron_k->maxpool_lexpr, expr->inf_coeff[i], expr->sup_coeff[i]);
            add_expr(res, mul_expr);
            free_expr(mul_expr);
        }
        else if (expr->inf_coeff[i] < 0)
        {
            expr_t* mul_expr = multiply_expr(neuron_k->maxpool_uexpr, expr->inf_coeff[i], expr->sup_coeff[i]);
            add_expr(res, mul_expr);
            free_expr(mul_expr);
        }
        else
        {
            double tmp1, tmp2;
            elina_double_interval_mul_cst_coeff(&tmp1, &tmp2, neuron_k->lb, neuron_k->ub, expr->inf_coeff[i], expr->sup_coeff[i]);
            res->inf_cst = res->inf_cst - tmp2;
            res->sup_cst = res->sup_cst + tmp2;
        }
    }

    res->inf_cst = res->inf_cst + expr->inf_cst;
    res->sup_cst = res->sup_cst + expr->sup_cst;

    free_expr(expr);
    *expr_ptr = res;
}


__device__ __host__
void expr_from_previous_layer(expr_t** expr_ptr, neuron_t** neurons)
{
    expr_t* expr = *expr_ptr;

    if(expr->size == 0)
    {
        return;
    }

    if((expr->inf_coeff == nullptr) || (expr->sup_coeff == nullptr))
    {
        return;
    }

    const size_t in_num_neurons = expr->size;
    size_t k;
    expr_t* res;

    if(expr->type == DENSE)
    {
        k = 0;
    }
    else
    {
        k = expr->dim[0];
    }

    if(neurons[k]->expr->size == 0)
    {
        res = multiply_cst_expr(neurons[k]->expr, expr->inf_coeff[0], expr->sup_coeff[0]);
    }
    else
    {
        res = multiply_expr(neurons[k]->expr, expr->inf_coeff[0], expr->sup_coeff[0]);
    }

    for(size_t i = 1; i < in_num_neurons; i++)
    {
        if(expr->type == DENSE)
        {
            k = i;
        }
        else
        {
            k = expr->dim[i];
        }

        expr_t* mul_expr;

        if(neurons[k]->expr->size == 0)
        {
            mul_expr = multiply_cst_expr(neurons[k]->expr, expr->inf_coeff[i], expr->sup_coeff[i]);
            add_cst_expr(res, mul_expr);
            free_expr(mul_expr);
        }
        else if((expr->inf_coeff[i] != 0) || (expr->sup_coeff[i] != 0))
        {
            mul_expr = multiply_expr(neurons[k]->expr, expr->inf_coeff[i], expr->sup_coeff[i]);
            add_expr(res, mul_expr);
            free_expr(mul_expr);
        }
    }

    res->inf_cst = res->inf_cst + expr->inf_cst;
    res->sup_cst = res->sup_cst + expr->sup_cst;

    free_expr(expr);
    *expr_ptr = res;
}



#include <iostream>
#include <cuda.h>


// Global const variables can be used from CPU and from GPU.
constexpr const double min_denormal = 4.940656458412465441766e-324;
constexpr const double ulp =          2.220446049250313080848e-16;


__device__ __host__
void elina_double_interval_mul(double* const a_inf, double* const a_sup, const double b_inf, const double b_sup, const double c_inf, const double c_sup)
