/*

 * Copyright 2008-2009 Katholieke Universiteit Leuven

 *

 * Use of this software is governed by the MIT license

 *

 * Written by Sven Verdoolaege, K.U.Leuven, Departement

 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium

 */

#include <assert.h>

#include <isl_map_private.h>

#include <isl_seq.h>

#include "isl_tab.h"

#include <isl_int.h>

#include <isl_config.h>

struct tab_lp {

  struct isl_ctx  *ctx;

  struct isl_vec  *row;

  struct isl_tab  *tab;

  struct isl_tab_undo **stack;

  isl_int   *obj;

  isl_int    opt;

  isl_int    opt_denom;

  isl_int    tmp;

  isl_int    tmp2;

  int          neq;

  unsigned   dim;

  /* number of constraints in initial product tableau */

  int    con_offset;

  /* objective function has fixed or no integer value */

  int    is_fixed;

};

#ifdef USE_GMP_FOR_MP

#define GBR_type              mpq_t

#define GBR_init(v)             mpq_init(v)

#define GBR_clear(v)              mpq_clear(v)

#define GBR_set(a,b)          mpq_set(a,b)

#define GBR_set_ui(a,b)         mpq_set_ui(a,b,1)

#define GBR_mul(a,b,c)          mpq_mul(a,b,c)

#define GBR_lt(a,b)         (mpq_cmp(a,b) < 0)

#define GBR_is_zero(a)          (mpq_sgn(a) == 0)

#define GBR_numref(a)         mpq_numref(a)

#define GBR_denref(a)         mpq_denref(a)

#define GBR_floor(a,b)          mpz_fdiv_q(a,GBR_numref(b),GBR_denref(b))

#define GBR_ceil(a,b)         mpz_cdiv_q(a,GBR_numref(b),GBR_denref(b))

#define GBR_set_num_neg(a, b)       mpz_neg(GBR_numref(*a), b);

#define GBR_set_den(a, b)       mpz_set(GBR_denref(*a), b);

#endif /* USE_GMP_FOR_MP */

#ifdef USE_IMATH_FOR_MP

#include <imrat.h>

#define GBR_type              mp_rat

#define GBR_init(v)             v = mp_rat_alloc()

#define GBR_clear(v)              mp_rat_free(v)

#define GBR_set(a,b)          mp_rat_copy(b,a)

#define GBR_set_ui(a,b)         mp_rat_set_uvalue(a,b,1)

#define GBR_mul(a,b,c)          mp_rat_mul(b,c,a)

#define GBR_lt(a,b)         (mp_rat_compare(a,b) < 0)

#define GBR_is_zero(a)          (mp_rat_compare_zero(a) == 0)

#ifdef USE_SMALL_INT_OPT

#define GBR_numref(a) isl_sioimath_encode_big(mp_rat_numer_ref(a))

#define GBR_denref(a) isl_sioimath_encode_big(mp_rat_denom_ref(a))

#define GBR_floor(a, b) isl_sioimath_fdiv_q((a), GBR_numref(b), GBR_denref(b))

#define GBR_ceil(a, b)  isl_sioimath_cdiv_q((a), GBR_numref(b), GBR_denref(b))

#define GBR_set_num_neg(a, b)                              \

  do {                                               \

    isl_sioimath_scratchspace_t scratch;       \

    impz_neg(mp_rat_numer_ref(*a),             \

        isl_sioimath_bigarg_src(*b, &scratch));\

  } while (0)

#define GBR_set_den(a, b)                                  \

  do {                                               \

    isl_sioimath_scratchspace_t scratch;       \

    impz_set(mp_rat_denom_ref(*a),             \

        isl_sioimath_bigarg_src(*b, &scratch));\

  } while (0)

#else /* USE_SMALL_INT_OPT */

#define GBR_numref(a)   mp_rat_numer_ref(a)

#define GBR_denref(a)   mp_rat_denom_ref(a)

#define GBR_floor(a,b)    impz_fdiv_q(a,GBR_numref(b),GBR_denref(b))

#define GBR_ceil(a,b)   impz_cdiv_q(a,GBR_numref(b),GBR_denref(b))

#define GBR_set_num_neg(a, b) impz_neg(GBR_numref(*a), b)

#define GBR_set_den(a, b) impz_set(GBR_denref(*a), b)

#endif /* USE_SMALL_INT_OPT */

#endif /* USE_IMATH_FOR_MP */

static struct tab_lp *init_lp(struct isl_tab *tab);

static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim);

static int solve_lp(struct tab_lp *lp);

static void get_obj_val(struct tab_lp* lp, GBR_type *F);

static void delete_lp(struct tab_lp *lp);

static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim);

static void get_alpha(struct tab_lp* lp, int row, GBR_type *alpha);

static int del_lp_row(struct tab_lp *lp) WARN_UNUSED;

static int cut_lp_to_hyperplane(struct tab_lp *lp, isl_int *row);

#define GBR_LP                struct tab_lp

#define GBR_lp_init(P)              init_lp(P)

#define GBR_lp_set_obj(lp, obj, dim)      set_lp_obj(lp, obj, dim)

#define GBR_lp_solve(lp)        solve_lp(lp)

#define GBR_lp_get_obj_val(lp, F)     get_obj_val(lp, F)

#define GBR_lp_delete(lp)       delete_lp(lp)

#define GBR_lp_next_row(lp)       lp->neq

#define GBR_lp_add_row(lp, row, dim)      add_lp_row(lp, row, dim)

#define GBR_lp_get_alpha(lp, row, alpha)    get_alpha(lp, row, alpha)

#define GBR_lp_del_row(lp)        del_lp_row(lp)

#define GBR_lp_is_fixed(lp)       (lp)->is_fixed

#define GBR_lp_cut(lp, obj)           cut_lp_to_hyperplane(lp, obj)

#include "basis_reduction_templ.c"

/* Set up a tableau for the Cartesian product of bset with itself.

 * This could be optimized by first setting up a tableau for bset

 * and then performing the Cartesian product on the tableau.

 */

static struct isl_tab *gbr_tab(struct isl_tab *tab, struct isl_vec *row)

{

  unsigned dim;

  struct isl_tab *prod;

  if (!tab || !row)

    return NULL;

  dim = tab->n_var;

  prod = isl_tab_product(tab, tab);

  if (isl_tab_extend_cons(prod, 3 * dim + 1) < 0) {

    isl_tab_free(prod);

    return NULL;

  }

  return prod;

}

static struct tab_lp *init_lp(struct isl_tab *tab)

{

  struct tab_lp *lp = NULL;

  if (!tab)

    return NULL;

  lp = isl_calloc_type(tab->mat->ctx, struct tab_lp);

  if (!lp)

    return NULL;

  isl_int_init(lp->opt);

  isl_int_init(lp->opt_denom);

  isl_int_init(lp->tmp);

  isl_int_init(lp->tmp2);

  lp->dim = tab->n_var;

  lp->ctx = tab->mat->ctx;

  isl_ctx_ref(lp->ctx);

  lp->stack = isl_alloc_array(lp->ctx, struct isl_tab_undo *, lp->dim);

  lp->row = isl_vec_alloc(lp->ctx, 1 + 2 * lp->dim);

  if (!lp->row)

    goto error;

  lp->tab = gbr_tab(tab, lp->row);

  if (!lp->tab)

    goto error;

  lp->con_offset = lp->tab->n_con;

  lp->obj = NULL;

  lp->neq = 0;

  return lp;

error:

  delete_lp(lp);

  return NULL;

}

static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim)

{

  lp->obj = row;

}

static int solve_lp(struct tab_lp *lp)

{

  enum isl_lp_result res;

  unsigned flags = 0;

  lp->is_fixed = 0;

  isl_int_set_si(lp->row->el[0], 0);

  isl_seq_cpy(lp->row->el + 1, lp->obj, lp->dim);

  isl_seq_neg(lp->row->el + 1 + lp->dim, lp->obj, lp->dim);

  if (lp->neq)

    flags = ISL_TAB_SAVE_DUAL;

  res = isl_tab_min(lp->tab, lp->row->el, lp->ctx->one,

        &lp->opt, &lp->opt_denom, flags);

  isl_int_mul_ui(lp->opt_denom, lp->opt_denom, 2);

  if (isl_int_abs_lt(lp->opt, lp->opt_denom)) {

    struct isl_vec *sample = isl_tab_get_sample_value(lp->tab);

    if (!sample)

      return -1;

    isl_seq_inner_product(lp->obj, sample->el + 1, lp->dim, &lp->tmp);

    isl_seq_inner_product(lp->obj, sample->el + 1 + lp->dim, lp->dim, &lp->tmp2);

    isl_int_cdiv_q(lp->tmp, lp->tmp, sample->el[0]);

    isl_int_fdiv_q(lp->tmp2, lp->tmp2, sample->el[0]);

    if (isl_int_ge(lp->tmp, lp->tmp2))

      
lp->is_fixed = 110.6k
;

    isl_vec_free(sample);

  }

  isl_int_divexact_ui(lp->opt_denom, lp->opt_denom, 2);

  if (res < 0)

    return -1;

  if (res != isl_lp_ok)

    
isl_die0
(lp->ctx, isl_error_internal,

      "unexpected missing (bounded) solution", return -1);

  return 0;

}

/* The current objective function has a fixed (or no) integer value.

 * Cut the tableau to the hyperplane that fixes this value in

 * both halves of the tableau.

 * Return 1 if the resulting tableau is empty.

 */

static int cut_lp_to_hyperplane(struct tab_lp *lp, isl_int *row)

{

  enum isl_lp_result res;

  isl_int_set_si(lp->row->el[0], 0);

  isl_seq_cpy(lp->row->el + 1, row, lp->dim);

  isl_seq_clr(lp->row->el + 1 + lp->dim, lp->dim);

  res = isl_tab_min(lp->tab, lp->row->el, lp->ctx->one,

        &lp->tmp, NULL, 0);

  if (res != isl_lp_ok)

    return -1;

  isl_int_neg(lp->row->el[0], lp->tmp);

  if (isl_tab_add_eq(lp->tab, lp->row->el) < 0)

    return -1;

  isl_seq_cpy(lp->row->el + 1 + lp->dim, row, lp->dim);

  isl_seq_clr(lp->row->el + 1, lp->dim);

  if (isl_tab_add_eq(lp->tab, lp->row->el) < 0)

    return -1;

  lp->con_offset += 2;

  return lp->tab->empty;

}

static void get_obj_val(struct tab_lp* lp, GBR_type *F)

{

  GBR_set_num_neg(F, lp->opt);

  GBR_set_den(F, lp->opt_denom);

}

static void delete_lp(struct tab_lp *lp)

{

  if (!lp)

    return;

  isl_int_clear(lp->opt);

  isl_int_clear(lp->opt_denom);

  isl_int_clear(lp->tmp);

  isl_int_clear(lp->tmp2);

  isl_vec_free(lp->row);

  free(lp->stack);

  isl_tab_free(lp->tab);

  isl_ctx_deref(lp->ctx);

  free(lp);

}

static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim)

{

  lp->stack[lp->neq] = isl_tab_snap(lp->tab);

  isl_int_set_si(lp->row->el[0], 0);

  isl_seq_cpy(lp->row->el + 1, row, lp->dim);

  isl_seq_neg(lp->row->el + 1 + lp->dim, row, lp->dim);

  if (isl_tab_add_valid_eq(lp->tab, lp->row->el) < 0)

    return -1;

  return lp->neq++;

}

static void get_alpha(struct tab_lp* lp, int row, GBR_type *alpha)

{

  row += lp->con_offset;

  GBR_set_num_neg(alpha, lp->tab->dual->el[1 + row]);

  GBR_set_den(alpha, lp->tab->dual->el[0]);

}

static int del_lp_row(struct tab_lp *lp)

{

  lp->neq--;

  return isl_tab_rollback(lp->tab, lp->stack[lp->neq]);

}