/* Vladimir Kotlyar  February 1996


 */

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "alpha_matrix.h"
#include "alpha_ineq.h"
#include "alpha_int.h"
#include "alpha_io.h"

#ifndef NDEBUG
#define XA_VERIFY_SORTED(i)	xa_verify_sorted(i)
#else
#define XA_VERIFY_SORTED(i)	((void)0)
#endif

#define XA_TAG_STEP		(32)

xa_ineq 
xa_alloc_ineq (arena_info *arena, int nterms)
{
  xa_ineq  ineq;

  assert(nterms >= 0);

  ineq = a_alloc0(arena, xa_ineq_t);
  xa_terms_ineq (ineq) = a_nalloc0(arena, xa_ineq_term_t, nterms+1);
  xa_nterm_ineq (ineq) = nterms;
  xa_set_sentinel_ineq (ineq);

  return ineq;
}

void 	
xa_verify_sorted (xa_ineq ineq)
{
  int k;
  int i1, i2;

  for (k = 2; k <= xa_nterm_ineq (ineq); k++) {
    i1 = xa_varind_ineq(ineq,k);
    i2 = xa_varind_ineq(ineq,k-1);
    assert(i1 > i2);
  }

  xa_verify_sentinel_ineq (ineq);
}

void
xa_compute_tag_ineq (xa_ineq ineq)
{
  int tag;
  int k;

  /* Just as in Omega the inequalities h^T x <= g and -h^T x <= g'
   * will have the tag with the same abs value but different signs */
  tag = 0;
  for (k = 1; k <= a_nterm_ineq (ineq); k++) {
    if ((tag < INT_MIN / 32) || (tag > INT_MAX / 32)) {
      /* this should deal gracefully with over/under-flow */
      tag = 0;
    }
    tag = 32 * tag + a_varind_ineq(ineq,k) * a_coeff_ineq(ineq,k);
  }
  
  a_tag_ineq (ineq) = tag;
}


AlphaIneq	
a_make_ineq_copy (arena_info *arena, AlphaIneq src)
{
  int		k;
  AlphaIneq 	ineq;

  ineq = xa_alloc_ineq (arena, a_nterm_ineq(src));

  a_tag_ineq (ineq) = a_tag_ineq (src);
  a_const_ineq (ineq) = a_const_ineq (src);
  
  for (k = 1; k <= a_nterm_ineq (ineq); k++) {
    a_varind_ineq(ineq,k) = a_varind_ineq(src,k);
    a_coeff_ineq(ineq,k) = a_coeff_ineq(src,k);
  }

  XA_VERIFY_SORTED(ineq);
  if (xa_is_normalized_ineq(src)) {
    xa_set_normalized_ineq (ineq);
  }
  else {
    xa_reset_normalized_ineq (ineq);
  }

  return ineq;
}

AlphaIneq	
a_make_ineq_saxpy (arena_info *arena, int a1, AlphaIneq i1, int a2, AlphaIneq i2)
{
  int		k1, k2, k, n, n1, n2, c;
  int		advance;
  AlphaIneq	ineq;

  /* first count the number of terms */
  n=0;
  k1=1;
  k2=1; 
  n1 = a_nterm_ineq(i1);
  n2 = a_nterm_ineq(i2);
  while((k1 <= n1+1) && (k2 <= n2+1)) {
    if (a_varind_ineq(i1,k1) > a_varind_ineq(i2,k2)) {
      advance = (a_coeff_ineq(i2,k2) != 0);
      k2++;
    }
    else if (a_varind_ineq(i1,k1) < a_varind_ineq(i2,k2)) {
      advance = (a_coeff_ineq(i1,k1) != 0);
      k1 ++;
    }
    else /* (a_varind_ineq(i1,k1) == a_varind_ineq(i2,k2)) */ {
      advance = ((a1*a_coeff_ineq(i1,k1) + a2*a_coeff_ineq(i2,k2)) != 0);
      k1++; k2++;
    }

    if (advance) {
      n++;
    }
  }

  /* this should be guaranteed by the sentinel:
   * k1 and k2 must have been bumped up at the same time 
   */
  assert ((k1 == n1+2) && (k2 == n2 + 2));

  ineq = xa_alloc_ineq (arena, n);

  /* now merge the terms */
  k=1;
  k1=1;
  k2=1; 
  while((k1 <= n1+1) && (k2 <= n2+1)) {
    advance = 0;
    if (a_varind_ineq(i1,k1) > a_varind_ineq(i2,k2)) {
      if (0 != a_coeff_ineq (i2, k2)) {
	a_coeff_ineq (ineq, k) = a1 * a_coeff_ineq (i2, k2);
	a_varind_ineq (ineq, k) = a_varind_ineq (i2, k2);
	advance = 1;
      }
      k2++;
    }
    else if (a_varind_ineq(i1,k1) < a_varind_ineq(i2,k2)) {
      if (0 != a_coeff_ineq (i1, k1)) {
	a_coeff_ineq (ineq, k) = a2 * a_coeff_ineq (i1, k1);
	a_varind_ineq (ineq, k) = a_varind_ineq (i1, k1);
	advance = 1;
      }
      k1 ++;
    }
    else if (a_varind_ineq(i1,k1) == a_varind_ineq(i2,k2)) {
      c = a1 * a_coeff_ineq (i1, k1) + a2 * a_coeff_ineq (i2, k2);
      if (0 != c) {
	a_coeff_ineq (ineq, k) = c;
	a_varind_ineq (ineq, k) = a_varind_ineq (i1, k1);
	advance = 1;
      }
      k1++; k2++;
    }
    if (advance) {
      k++;
    }
  }

  /* don't forget const coeffs */
  a_const_ineq (ineq) = a1 * a_const_ineq (i1) + a2 * a_const_ineq (i2);

  /* have to create a new tag */
  xa_compute_tag_ineq (ineq);
  xa_reset_normalized_ineq (ineq);

  XA_VERIFY_SORTED (ineq);
  
  return ineq;
}

AlphaIneq	
a_make_ineq_from_dense (arena_info *arena, AlphaVector a, int b)
{
  AlphaIneq 	ineq;
  int		n;
  int		j, k;

  for (n=0, j=1; j <= a_length_vec (a); j++) {
    if (0 != a_elt_vec(a,j)) {
      n++;
    }
  }

  ineq = xa_alloc_ineq (arena, n);
  a_const_ineq (ineq) = b;

  for (k=1, j=1; j <= a_length_vec (a); j++) {
    if (0 != a_elt_vec(a,j)) {
      a_varind_ineq (ineq, k) = j;
      a_coeff_ineq (ineq, k) = a_elt_vec (a, j);
      k++;
    }
  }
  
  xa_compute_tag_ineq (ineq);
  xa_reset_normalized_ineq (ineq);

  XA_VERIFY_SORTED (ineq);

  return ineq;
}

AlphaIneq
a_make_ineq_var_transform (arena_info *arena, AlphaIneq ineq, AlphaVector X0, AlphaMatrix H)
{
  /*
   * The inequality is a^T x <= b
   * The transform is x = x0 + H*t
   * So the new inequality is:
   * a^T * H * t <= b - a^T*x0
   */
  AlphaIneq	ineq_new;
  int		nt = a_cols_mx (H);
  int		mh = a_rows_mx (H);
  AlphaVector	Anew;
  int		sum;
  int		b_new;

  assert (nt > 0);
  assert (mh > 0);

  Anew = a_make_vector (arena, nt);
  A_FOR_COLS_MX(H,col,j) {
    sum = 0;
    A_FOR_TERMS_IN_INEQ (ineq, varind, coeff) {
      sum += coeff * a_elt_vec (col, varind);
    }
    A_END_TERMS_IN_INEQ;
    a_elt_vec (Anew,j) = sum;
  }
  A_END_COLS_MX;

  b_new = a_const_ineq (ineq);
  A_FOR_TERMS_IN_INEQ (ineq, varind, coeff) {
    b_new -= coeff * a_elt_vec (X0, varind);
  }
  A_END_TERMS_IN_INEQ;
  
  ineq_new = a_make_ineq_from_dense (arena, Anew, b_new);
  xa_reset_normalized_ineq (ineq_new);
  
  return ineq_new;
}

AlphaIneq
a_make_ineq_neg (arena_info *arena, AlphaIneq ineq)
{
  AlphaIneq 	ineq2 = a_make_ineq_copy (arena, ineq);
  int		n = xa_nterm_ineq (ineq);
  int		i;

  for (i = 1; i <= n; i++) {
    a_coeff_ineq (ineq2, i) = -a_coeff_ineq (ineq2, i);
  }
  a_const_ineq (ineq2) = -a_const_ineq (ineq2) - 1;
  
  xa_compute_tag_ineq (ineq2);
  xa_reset_normalized_ineq (ineq2);

  XA_VERIFY_SORTED (ineq2);

  return ineq2;
}

int		
a_coeff4var_ineq (AlphaIneq ineq, int varind)
{
  int k;

  /* simple linear search for now */
  for (k = 1; k <= a_nterm_ineq (ineq); k++) {
    if (a_varind_ineq (ineq, k) == varind) {
      return a_coeff_ineq (ineq, k);
    }
  }

  return 0;  
}

int
a_is_same_coeff_ineq (AlphaIneq i1, AlphaIneq i2)
{
  int	k;

  if (a_nterm_ineq(i1) != a_nterm_ineq(i2)) {
    return 0;
  }
  else {
    for (k = 1; k <= a_nterm_ineq(i1); k++) {
      if ((a_varind_ineq(i1,k) != a_varind_ineq(i2,k)) ||
	  (a_coeff_ineq(i1,k) != a_coeff_ineq(i2,k))) {
	return 0;
      }
    }
  }
  return 1;
}

int
a_is_opposite_coeff_ineq (AlphaIneq i1, AlphaIneq i2)
{
  int	k;

  if (a_nterm_ineq(i1) != a_nterm_ineq(i2)) {
    return 0;
  }
  else {
    for (k = 1; k <= a_nterm_ineq(i1); k++) {
      if ((a_varind_ineq(i1,k) != a_varind_ineq(i2,k)) ||
	  (a_coeff_ineq(i1,k) != -a_coeff_ineq(i2,k))) {
	return 0;
      }
    }
  }
  return 1;
}

int
a_normalize_ineq (AlphaIneq ineq)
{
  int g;
  int k;

  if (a_nterm_ineq(ineq) > 0) {
    g = a_abs(a_coeff_ineq(ineq,1));
    for (k = 2; k <= a_nterm_ineq(ineq); k++) {
      g = a_gcd (g, a_abs(a_coeff_ineq(ineq,k)));
    }
    for (k = 1; k <= a_nterm_ineq(ineq); k++) {
      a_coeff_ineq(ineq,k) = a_coeff_ineq(ineq,k) / g;
    }
    a_const_ineq(ineq) = a_floor_div(a_const_ineq(ineq), g);
  }
  else {
    g = 1;
  }

  xa_compute_tag_ineq (ineq);
  xa_set_normalized_ineq(ineq);

  return g;
}

void
a_dump_ineq (FILE *fp, AlphaIneq ineq)
{
  char 		plus;

  if (a_nterm_ineq (ineq) > 0) {
    plus = ' ';
    A_FOR_TERMS_IN_INEQ (ineq,v,c) {
      if (c == 1) {
	fprintf (fp, "%cx(%d) ", plus, v);
      }
      else if (c == -1) {
	fprintf (fp, "-x(%d) ", v);
      }
      else if (c > 0) {
	fprintf (fp, "%c%d*x(%d) ", plus, c, v);
      }
      else if (c < 0) {
	fprintf (fp, "-%d*x(%d) ", -c, v);
      }
      plus = '+';
    }
    A_END_TERMS_IN_INEQ;
  }
  else {
    fprintf (fp, "0 ");
  }

  fprintf (fp, " <= %d ", a_const_ineq (ineq));
}


/***********************************************************************
 *			LISTS OF INEQUALITIES
 ***********************************************************************/

xa_ineq_list
xa_make_ineq_list (arena_info *arena)
{
  xa_ineq_list l = a_alloc0 (arena, xa_ineq_list_t);

  return l;
}

void		
xa_clear_ineq_list (xa_ineq_list l)
{
  xa_ineq_list_elt	elt, next;

  for (elt = xa_head_ineq_list (l); elt; elt = next) {
    next = xa_next_ineq_list_elt (elt);
  }

  xa_length_ineq_list (l) = 0;
  xa_head_ineq_list (l) = 0;
}

void
xa_add_ineq_list(arena_info *arena, xa_ineq_list l, xa_ineq ineq)
{
  xa_ineq_list_elt elt, next;
  xa_ineq ineq_copy = a_make_ineq_copy (arena, ineq);

  next = xa_head_ineq_list(l);

  elt = a_alloc0(arena, xa_ineq_list_elt_t);
  xa_ineq_ineq_list_elt (elt) = ineq_copy;
  xa_next_ineq_list_elt (elt) = xa_head_ineq_list (l);
  xa_prev_ineq_list_elt (elt) = 0;
  xa_head_ineq_list (l) = elt;
  if (next) {
    xa_prev_ineq_list_elt (next) = elt;
  }

  xa_length_ineq_list (l) ++;
}

void
xa_delete_ineq_list_elt(xa_ineq_list l, xa_ineq_list_elt elt)
{ 
  xa_ineq_list_elt prev,next;
  prev = xa_prev_ineq_list_elt(elt);
  next = xa_next_ineq_list_elt(elt);
  if (prev) { xa_next_ineq_list_elt(prev) = next; }
  if (next) { xa_prev_ineq_list_elt(next) = prev; }
  if(xa_head_ineq_list(l) == elt) { 
    xa_head_ineq_list(l) = next; 
  }

  xa_length_ineq_list (l) --;
}


void 
a_normalize_ineq_list (AlphaIneqList l)
{
  A_FOR_INEQ_IN_LIST (l, ineq, elt) {
    a_normalize_ineq (ineq);
  }
  A_END_INEQ_IN_LIST;
}

AlphaIneqList 
a_make_ineq_list_copy (arena_info *arena, AlphaIneqList l)
{
  AlphaIneqList l2 = a_make_ineq_list(arena);
  AlphaIneq	ineq2;

  A_FOR_INEQ_IN_LIST (l, ineq, elt) {
    ineq2 = a_make_ineq_copy (arena, ineq);
    a_add_ineq_list (arena, l2, ineq2);
  }
  A_END_INEQ_IN_LIST;

  return l2;
}


void
a_dump_ineq_list (FILE *fp, AlphaIneqList l)
{
  fprintf (fp, "HEAD == %p LENGTH == %d\n", 
	   xa_head_ineq_list (l), 
	   xa_length_ineq_list (l));

  A_FOR_INEQ_IN_LIST(l,ineq,elt){
    fprintf (fp, "INEQ == %8p PREV == %8p NEXT == %8p ELT == %8p ",
	     ineq,
	     xa_prev_ineq_list_elt(elt),
	     xa_next_ineq_list_elt(elt),
	     elt);
    a_fprintf (fp, " %I \n", ineq);
  }
  A_END_INEQ_IN_LIST;
}

void
a_print_ineq_list (FILE *fp, AlphaIneqList l)
{
  fprintf (fp, "{ ");
  A_FOR_INEQ_IN_LIST(l,ineq,elt) {
    a_fprintf (fp, " %I,", ineq);
  }
  A_END_INEQ_IN_LIST;
  fprintf (fp, " }");
}

/*
 * the matrix and vector are assumed to be zero on entry 
 */
void
a_ineq_list2matrix_form (AlphaIneqList l, AlphaMatrix A, AlphaVector B)
{
  int  		m = a_rows_mx (A);
  int  		n = a_cols_mx (A);
  int		i;
  
  assert (m == a_length_vec (B));
  assert (m == a_length_ineq_list (l));

  i = 1;
  A_FOR_INEQ_IN_LIST (l, ineq, elt) {
    A_FOR_TERMS_IN_INEQ (ineq, varind, coeff) {
      a_elt_mx (A, i, varind) = coeff;
    }
    A_END_TERMS_IN_INEQ;
    a_elt_vec (B, i) = a_const_ineq (ineq);
    i++;
  }
  A_END_INEQ_IN_LIST;
}

void
a_matrix_form2ineq_list (arena_info *arena, AlphaMatrix A, AlphaVector B, AlphaIneqList l)
{
  int  		m = a_rows_mx (A);
  int  		n = a_cols_mx (A);
  int		i;
  AlphaIneq  	ineq;
  AlphaMatrix   At;
  
  assert (m == a_length_vec (B));
  assert (0 == a_length_ineq_list (l));

  At = a_make_matrix (arena, n, m);
  a_transpose (A, At);

  A_FOR_COLS_MX (At, col, j) {
    ineq = a_make_ineq_from_dense (arena, col, a_elt_vec (B, j));
    a_add_ineq_list (arena, l, ineq);
  }
  A_END_COLS_MX;

}


void 
a_transform_ineq_list (arena_info *arena, AlphaIneqList from, AlphaVector T0, 
		       AlphaMatrix H, AlphaIneqList to)
{
  AlphaIneq inew;

  A_FOR_INEQ_IN_LIST (from, ineq, elt) {
    inew = a_make_ineq_var_transform (arena, ineq, T0, H);
    a_normalize_ineq (ineq);
    a_add_ineq_list (arena, to, inew);
  }
  A_END_INEQ_IN_LIST;
}

int
a_is_satisfied_ineq_list (AlphaIneqList C, AlphaVector X)
{
  int sum;
  A_FOR_INEQ_IN_LIST (C, ineq, elt) {
    sum = 0;
    A_FOR_TERMS_IN_INEQ (ineq, varind, coeff) {
      sum += coeff * a_elt_vec (X, varind);
    }
    A_END_TERMS_IN_INEQ;
    
    if (sum > a_const_ineq (ineq)) {
      return 0;
    }
  }
  A_END_INEQ_IN_LIST;

  return 1;
}

