/*

 * Copyright 2006-2007 Universiteit Leiden

 * Copyright 2008-2009 Katholieke Universiteit Leuven

 * Copyright 2010      INRIA Saclay

 *

 * Use of this software is governed by the MIT license

 *

 * Written by Sven Verdoolaege, Leiden Institute of Advanced Computer Science,

 * Universiteit Leiden, Niels Bohrweg 1, 2333 CA Leiden, The Netherlands

 * and K.U.Leuven, Departement Computerwetenschappen, Celestijnenlaan 200A,

 * B-3001 Leuven, Belgium

 * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,

 * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France

 */

#include <isl_ctx_private.h>

#include <isl_map_private.h>

#include <isl/set.h>

#include <isl_seq.h>

#include <isl_morph.h>

#include <isl_factorization.h>

#include <isl_vertices_private.h>

#include <isl_polynomial_private.h>

#include <isl_options_private.h>

#include <isl_vec_private.h>

#include <isl_bernstein.h>

struct bernstein_data {

  enum isl_fold type;

  isl_qpolynomial *poly;

  int check_tight;

  isl_cell *cell;

  isl_qpolynomial_fold *fold;

  isl_qpolynomial_fold *fold_tight;

  isl_pw_qpolynomial_fold *pwf;

  isl_pw_qpolynomial_fold *pwf_tight;

};

static int vertex_is_integral(__isl_keep isl_basic_set *vertex)

{

  unsigned nvar;

  unsigned nparam;

  int i;

  nvar = isl_basic_set_dim(vertex, isl_dim_set);

  nparam = isl_basic_set_dim(vertex, isl_dim_param);

  for (i = 0; i < nvar; ++i) {

    int r = nvar - 1 - i;

    if (!isl_int_is_one(vertex->eq[r][1 + nparam + i]) &&

        !isl_int_is_negone(vertex->eq[r][1 + nparam + i]))

      return 0;

  }

  return 1;

}

static __isl_give isl_qpolynomial *vertex_coordinate(

  __isl_keep isl_basic_set *vertex, int i, __isl_take isl_space *dim)

{

  unsigned nvar;

  unsigned nparam;

  int r;

  isl_int denom;

  isl_qpolynomial *v;

  nvar = isl_basic_set_dim(vertex, isl_dim_set);

  nparam = isl_basic_set_dim(vertex, isl_dim_param);

  r = nvar - 1 - i;

  isl_int_init(denom);

  isl_int_set(denom, vertex->eq[r][1 + nparam + i]);

  isl_assert(vertex->ctx, !isl_int_is_zero(denom), goto error);

  if (isl_int_is_pos(denom))

    isl_seq_neg(vertex->eq[r], vertex->eq[r],

        1 + isl_basic_set_total_dim(vertex));

  else

    isl_int_neg(denom, denom);

  v = isl_qpolynomial_from_affine(dim, vertex->eq[r], denom);

  isl_int_clear(denom);

  return v;

error:

  isl_space_free(dim);

  isl_int_clear(denom);

  return NULL;

}

/* Check whether the bound associated to the selection "k" is tight,

 * which is the case if we select exactly one vertex and if that vertex

 * is integral for all values of the parameters.

 */

static int is_tight(int *k, int n, int d, isl_cell *cell)

{

  int i;

  for (i = 0; i < n; ++i) {

    int v;

    if (k[i] != d) {

      if (k[i])

        return 0;

      continue;

    }

    v = cell->ids[n - 1 - i];

    return vertex_is_integral(cell->vertices->v[v].vertex);

  }

  return 0;

}

static void add_fold(__isl_take isl_qpolynomial *b, __isl_keep isl_set *dom,

  int *k, int n, int d, struct bernstein_data *data)

{

  isl_qpolynomial_fold *fold;

  fold = isl_qpolynomial_fold_alloc(data->type, b);

  if (data->check_tight && is_tight(k, n, d, data->cell))

    data->fold_tight = isl_qpolynomial_fold_fold_on_domain(dom,

              data->fold_tight, fold);

  else

    data->fold = isl_qpolynomial_fold_fold_on_domain(dom,

              data->fold, fold);

}

/* Extract the coefficients of the Bernstein base polynomials and store

 * them in data->fold and data->fold_tight.

 *

 * In particular, the coefficient of each monomial

 * of multi-degree (k[0], k[1], ..., k[n-1]) is divided by the corresponding

 * multinomial coefficient d!/k[0]! k[1]! ... k[n-1]!

 *

 * c[i] contains the coefficient of the selected powers of the first i+1 vars.

 * multinom[i] contains the partial multinomial coefficient.

 */

static void extract_coefficients(isl_qpolynomial *poly,

  __isl_keep isl_set *dom, struct bernstein_data *data)

{

  int i;

  int d;

  int n;

  isl_ctx *ctx;

  isl_qpolynomial **c = NULL;

  int *k = NULL;

  int *left = NULL;

  isl_vec *multinom = NULL;

  if (!poly)

    return;

  ctx = isl_qpolynomial_get_ctx(poly);

  n = isl_qpolynomial_dim(poly, isl_dim_in);

  d = isl_qpolynomial_degree(poly);

  isl_assert(ctx, n >= 2, return);

  c = isl_calloc_array(ctx, isl_qpolynomial *, n);

  k = isl_alloc_array(ctx, int, n);

  left = isl_alloc_array(ctx, int, n);

  multinom = isl_vec_alloc(ctx, n);

  if (!c || !k || !left || !multinom)

    goto error;

  isl_int_set_si(multinom->el[0], 1);

  for (k[0] = d; k[0] >= 0; --k[0]) {

    int i = 1;

    isl_qpolynomial_free(c[0]);

    c[0] = isl_qpolynomial_coeff(poly, isl_dim_in, n - 1, k[0]);

    left[0] = d - k[0];

    k[1] = -1;

    isl_int_set(multinom->el[1], multinom->el[0]);

    while (i > 0) {

      if (i == n - 1) {

        int j;

        isl_space *dim;

        isl_qpolynomial *b;

        isl_qpolynomial *f;

        for (j = 2; j <= left[i - 1]; ++j)

          isl_int_divexact_ui(multinom->el[i],

            multinom->el[i], j);

        b = isl_qpolynomial_coeff(c[i - 1], isl_dim_in,

          n - 1 - i, left[i - 1]);

        b = isl_qpolynomial_project_domain_on_params(b);

        dim = isl_qpolynomial_get_domain_space(b);

        f = isl_qpolynomial_rat_cst_on_domain(dim, ctx->one,

          multinom->el[i]);

        b = isl_qpolynomial_mul(b, f);

        k[n - 1] = left[n - 2];

        add_fold(b, dom, k, n, d, data);

        --i;

        continue;

      }

      if (k[i] >= left[i - 1]) {

        --i;

        continue;

      }

      ++k[i];

      if (k[i])

        isl_int_divexact_ui(multinom->el[i],

          multinom->el[i], k[i]);

      isl_qpolynomial_free(c[i]);

      c[i] = isl_qpolynomial_coeff(c[i - 1], isl_dim_in,

          n - 1 - i, k[i]);

      left[i] = left[i - 1] - k[i];

      k[i + 1] = -1;

      isl_int_set(multinom->el[i + 1], multinom->el[i]);

      ++i;

    }

    isl_int_mul_ui(multinom->el[0], multinom->el[0], k[0]);

  }

  for (i = 0; i < n; ++i)

    isl_qpolynomial_free(c[i]);

  isl_vec_free(multinom);

  free(left);

  free(k);

  free(c);

  return;

error:

  isl_vec_free(multinom);

  free(left);

  free(k);

  if (c)

    for (i = 0; i < n; ++i)

      isl_qpolynomial_free(c[i]);

  free(c);

  return;

}

/* Perform bernstein expansion on the parametric vertices that are active

 * on "cell".

 *

 * data->poly has been homogenized in the calling function.

 *

 * We plug in the barycentric coordinates for the set variables

 *

 *    \vec x = \sum_i \alpha_i v_i(\vec p)

 *

 * and the constant "1 = \sum_i \alpha_i" for the homogeneous dimension.

 * Next, we extract the coefficients of the Bernstein base polynomials.

 */

static isl_stat bernstein_coefficients_cell(__isl_take isl_cell *cell,

  void *user)

{

  int i, j;

  struct bernstein_data *data = (struct bernstein_data *)user;

  isl_space *dim_param;

  isl_space *dim_dst;

  isl_qpolynomial *poly = data->poly;

  unsigned nvar;

  int n_vertices;

  isl_qpolynomial **subs;

  isl_pw_qpolynomial_fold *pwf;

  isl_set *dom;

  isl_ctx *ctx;

  if (!poly)

    goto error;

  nvar = isl_qpolynomial_dim(poly, isl_dim_in) - 1;

  n_vertices = cell->n_vertices;

  ctx = isl_qpolynomial_get_ctx(poly);

  if (n_vertices > nvar + 1 && ctx->opt->bernstein_triangulate)

    return isl_cell_foreach_simplex(cell,

              &bernstein_coefficients_cell, user);

  subs = isl_alloc_array(ctx, isl_qpolynomial *, 1 + nvar);

  if (!subs)

    goto error;

  dim_param = isl_basic_set_get_space(cell->dom);

  dim_dst = isl_qpolynomial_get_domain_space(poly);

  dim_dst = isl_space_add_dims(dim_dst, isl_dim_set, n_vertices);

  for (i = 0; i < 1 + nvar; ++i)

    subs[i] = isl_qpolynomial_zero_on_domain(isl_space_copy(dim_dst));

  for (i = 0; i < n_vertices; ++i) {

    isl_qpolynomial *c;

    c = isl_qpolynomial_var_on_domain(isl_space_copy(dim_dst), isl_dim_set,

          1 + nvar + i);

    for (j = 0; j < nvar; ++j) {

      int k = cell->ids[i];

      isl_qpolynomial *v;

      v = vertex_coordinate(cell->vertices->v[k].vertex, j,

            isl_space_copy(dim_param));

      v = isl_qpolynomial_add_dims(v, isl_dim_in,

              1 + nvar + n_vertices);

      v = isl_qpolynomial_mul(v, isl_qpolynomial_copy(c));

      subs[1 + j] = isl_qpolynomial_add(subs[1 + j], v);

    }

    subs[0] = isl_qpolynomial_add(subs[0], c);

  }

  isl_space_free(dim_dst);

  poly = isl_qpolynomial_copy(poly);

  poly = isl_qpolynomial_add_dims(poly, isl_dim_in, n_vertices);

  poly = isl_qpolynomial_substitute(poly, isl_dim_in, 0, 1 + nvar, subs);

  poly = isl_qpolynomial_drop_dims(poly, isl_dim_in, 0, 1 + nvar);

  data->cell = cell;

  dom = isl_set_from_basic_set(isl_basic_set_copy(cell->dom));

  data->fold = isl_qpolynomial_fold_empty(data->type, isl_space_copy(dim_param));

  data->fold_tight = isl_qpolynomial_fold_empty(data->type, dim_param);

  extract_coefficients(poly, dom, data);

  pwf = isl_pw_qpolynomial_fold_alloc(data->type, isl_set_copy(dom),

              data->fold);

  data->pwf = isl_pw_qpolynomial_fold_fold(data->pwf, pwf);

  pwf = isl_pw_qpolynomial_fold_alloc(data->type, dom, data->fold_tight);

  data->pwf_tight = isl_pw_qpolynomial_fold_fold(data->pwf_tight, pwf);

  isl_qpolynomial_free(poly);

  isl_cell_free(cell);

  for (i = 0; i < 1 + nvar; ++i)

    isl_qpolynomial_free(subs[i]);

  free(subs);

  return isl_stat_ok;

error:

  isl_cell_free(cell);

  return isl_stat_error;

}

/* Base case of applying bernstein expansion.

 *

 * We compute the chamber decomposition of the parametric polytope "bset"

 * and then perform bernstein expansion on the parametric vertices

 * that are active on each chamber.

 */

static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_base(

  __isl_take isl_basic_set *bset,

  __isl_take isl_qpolynomial *poly, struct bernstein_data *data, int *tight)

{

  unsigned nvar;

  isl_space *dim;

  isl_pw_qpolynomial_fold *pwf;

  isl_vertices *vertices;

  int covers;

  nvar = isl_basic_set_dim(bset, isl_dim_set);

  if (nvar == 0) {

    isl_set *dom;

    isl_qpolynomial_fold *fold;

    fold = isl_qpolynomial_fold_alloc(data->type, poly);

    dom = isl_set_from_basic_set(bset);

    if (tight)

      *tight = 1;

    pwf = isl_pw_qpolynomial_fold_alloc(data->type, dom, fold);

    return isl_pw_qpolynomial_fold_project_domain_on_params(pwf);

  }

  if (isl_qpolynomial_is_zero(poly)) {

    isl_set *dom;

    isl_qpolynomial_fold *fold;

    fold = isl_qpolynomial_fold_alloc(data->type, poly);

    dom = isl_set_from_basic_set(bset);

    pwf = isl_pw_qpolynomial_fold_alloc(data->type, dom, fold);

    if (tight)

      *tight = 1;

    return isl_pw_qpolynomial_fold_project_domain_on_params(pwf);

  }

  dim = isl_basic_set_get_space(bset);

  dim = isl_space_params(dim);

  dim = isl_space_from_domain(dim);

  dim = isl_space_add_dims(dim, isl_dim_set, 1);

  data->pwf = isl_pw_qpolynomial_fold_zero(isl_space_copy(dim), data->type);

  data->pwf_tight = isl_pw_qpolynomial_fold_zero(dim, data->type);

  data->poly = isl_qpolynomial_homogenize(isl_qpolynomial_copy(poly));

  vertices = isl_basic_set_compute_vertices(bset);

  if (isl_vertices_foreach_disjoint_cell(vertices,

          &bernstein_coefficients_cell, data) < 0)

    data->pwf = isl_pw_qpolynomial_fold_free(data->pwf);

  isl_vertices_free(vertices);

  isl_qpolynomial_free(data->poly);

  isl_basic_set_free(bset);

  isl_qpolynomial_free(poly);

  covers = isl_pw_qpolynomial_fold_covers(data->pwf_tight, data->pwf);

  if (covers < 0)

    goto error;

  if (tight)

    *tight = covers;

  if (covers) {

    isl_pw_qpolynomial_fold_free(data->pwf);

    return data->pwf_tight;

  }

  data->pwf = isl_pw_qpolynomial_fold_fold(data->pwf, data->pwf_tight);

  return data->pwf;

error:

  isl_pw_qpolynomial_fold_free(data->pwf_tight);

  isl_pw_qpolynomial_fold_free(data->pwf);

  return NULL;

}

/* Apply bernstein expansion recursively by working in on len[i]

 * set variables at a time, with i ranging from n_group - 1 to 0.

 */

static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_recursive(

  __isl_take isl_pw_qpolynomial *pwqp,

  int n_group, int *len, struct bernstein_data *data, int *tight)

{

  int i;

  unsigned nparam;

  unsigned nvar;

  isl_pw_qpolynomial_fold *pwf;

  if (!pwqp)

    return NULL;

  nparam = isl_pw_qpolynomial_dim(pwqp, isl_dim_param);

  nvar = isl_pw_qpolynomial_dim(pwqp, isl_dim_in);

  pwqp = isl_pw_qpolynomial_move_dims(pwqp, isl_dim_param, nparam,

          isl_dim_in, 0, nvar - len[n_group - 1]);

  pwf = isl_pw_qpolynomial_bound(pwqp, data->type, tight);

  for (i = n_group - 2; i >= 0; --i) {

    nparam = isl_pw_qpolynomial_fold_dim(pwf, isl_dim_param);

    pwf = isl_pw_qpolynomial_fold_move_dims(pwf, isl_dim_in, 0,

        isl_dim_param, nparam - len[i], len[i]);

    if (tight && !*tight)

      tight = NULL;

    pwf = isl_pw_qpolynomial_fold_bound(pwf, tight);

  }

  return pwf;

}

static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_factors(

  __isl_take isl_basic_set *bset,

  __isl_take isl_qpolynomial *poly, struct bernstein_data *data, int *tight)

{

  isl_factorizer *f;

  isl_set *set;

  isl_pw_qpolynomial *pwqp;

  isl_pw_qpolynomial_fold *pwf;

  f = isl_basic_set_factorizer(bset);

  if (!f)

    goto error;

  if (f->n_group == 0) {

    isl_factorizer_free(f);

    return  bernstein_coefficients_base(bset, poly, data, tight);

  }

  set = isl_set_from_basic_set(bset);

  pwqp = isl_pw_qpolynomial_alloc(set, poly);

  pwqp = isl_pw_qpolynomial_morph_domain(pwqp, isl_morph_copy(f->morph));

  pwf = bernstein_coefficients_recursive(pwqp, f->n_group, f->len, data,

            tight);

  isl_factorizer_free(f);

  return pwf;

error:

  isl_basic_set_free(bset);

  isl_qpolynomial_free(poly);

  return NULL;

}

static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_full_recursive(

  __isl_take isl_basic_set *bset,

  __isl_take isl_qpolynomial *poly, struct bernstein_data *data, int *tight)

{

  int i;

  int *len;

  unsigned nvar;

  isl_pw_qpolynomial_fold *pwf;

  isl_set *set;

  isl_pw_qpolynomial *pwqp;

  if (!bset || !poly)

    goto error;

  nvar = isl_basic_set_dim(bset, isl_dim_set);

  len = isl_alloc_array(bset->ctx, int, nvar);

  if (nvar && !len)

    goto error;

  for (i = 0; i < nvar; ++i)

    len[i] = 1;

  set = isl_set_from_basic_set(bset);

  pwqp = isl_pw_qpolynomial_alloc(set, poly);

  pwf = bernstein_coefficients_recursive(pwqp, nvar, len, data, tight);

  free(len);

  return pwf;

error:

  isl_basic_set_free(bset);

  isl_qpolynomial_free(poly);

  return NULL;

}

/* Compute a bound on the polynomial defined over the parametric polytope

 * using bernstein expansion and store the result

 * in bound->pwf and bound->pwf_tight.

 *

 * If bernstein_recurse is set to ISL_BERNSTEIN_FACTORS, we check if

 * the polytope can be factorized and apply bernstein expansion recursively

 * on the factors.

 * If bernstein_recurse is set to ISL_BERNSTEIN_INTERVALS, we apply

 * bernstein expansion recursively on each dimension.

 * Otherwise, we apply bernstein expansion on the entire polytope.

 */

isl_stat isl_qpolynomial_bound_on_domain_bernstein(

  __isl_take isl_basic_set *bset, __isl_take isl_qpolynomial *poly,

  struct isl_bound *bound)

{

  struct bernstein_data data;

  isl_pw_qpolynomial_fold *pwf;

  unsigned nvar;

  int tight = 0;

  int *tp = bound->check_tight ? 
&tight0
 : NULL;

  if (!bset || !poly)

    goto error;

  data.type = bound->type;

  data.check_tight = bound->check_tight;

  nvar = isl_basic_set_dim(bset, isl_dim_set);

  if (bset->ctx->opt->bernstein_recurse & ISL_BERNSTEIN_FACTORS)

    pwf = bernstein_coefficients_factors(bset, poly, &data, tp);

  else if (nvar > 1 &&

      (bset->ctx->opt->bernstein_recurse & ISL_BERNSTEIN_INTERVALS))

    pwf = bernstein_coefficients_full_recursive(bset, poly, &data, tp);

  else

    pwf = bernstein_coefficients_base(bset, poly, &data, tp);

  if (tight)

    bound->pwf_tight = isl_pw_qpolynomial_fold_fold(bound->pwf_tight, pwf);

  else

    bound->pwf = isl_pw_qpolynomial_fold_fold(bound->pwf, pwf);

  return isl_stat_ok;

error:

  isl_basic_set_free(bset);

  isl_qpolynomial_free(poly);

  return isl_stat_error;

}