source: git/Singular/LIB/cisimplicial.lib @ 33b509

// IB/IG, last modified: 29.07.2010
//////////////////////////////////////////////////////////////////////

category= "Commutative Algebra";
version = "$Id$";
info="
LIBRARY: cisimplicial.lib. Determines if the toric ideal of
a simplicial toric variety is a complete intersection

AUTHORS: I.Bermejo,               ibermejo@ull.es
@*       I.Garcia-Marco,          iggarcia@ull.es

OVERVIEW:
  A library for determining if a simplicial toric ideal is a complete
  intersection with NO NEED of computing explicitly a system of generators
  of such ideal. The procedures are based on two papers:
  I. Bermejo, I. Garcia-Marco and J.J. Salazar-Gonzalez: 'An algorithm for
  checking whether the toric ideal of an affine monomial curve is a complete
  intersection', J. Symbolic Computation 42 (2007) pags: 971--991 and
  I.Bermejo and I. Garcia-Marco: 'Complete  intersections in simplicial toric
  varieties', Preprint (2010)

PROCEDURES:
 minMult(a,b);        computes the minimum multiple of a that belongs to the
                semigroup generated by b
 belongSemigroup(v,A[,n]);     checks whether A*x = v has a nonnegative
                               integral solution
 oneDimBelongSemigroup(n,v[,m]);      checks whether v*x = n has a
                                      nonnegative integral solution
 cardGroup(A);           computes the cardinal of Z^m / ZA
 isCI(A);     checks wether I(A) is a complete intersection
";

LIB "general.lib";


/////////////////////////////////////////////////////////////////////


static proc Multiple(intvec v, intmat A, list #)
"
USAGE:   Multiple(v,A[,n]); v is an integral vector, A is an integral matrix, n is an integer.
RETURN:  0 if none of the [n first] columns of A divides the vector v or an intvec otherwise.
         In this case v = answer[2] * (answer[1]-th column of A).
ASSUME:  nrows(v) = nrows(A) [and n <= nrows(A)].
"
{

   intvec answer;
   if (v == 0)
   {
      answer[1] = 1;
      answer[2] = 0;
      return (answer);
   }

   int last;
   int e = size(#);
   if  (e > 0)
   {
      last = #[1];
   }
   else
   {
      last = ncols(A);
   }

   int i,j,s;
   for (j = 1; j <= last; j++)
   {
      s = 0;
      for (i = 1; i <= nrows(A); i++)
      {
          if ((v[i] == 0)&&(A[i,j] != 0))
          {
             // it is not multiple of A_j
            break;
          }
          if (v[i] != 0)
          {
            if (A[i,j] == 0)
            {
              // it is not multiple of A_j
              break;
            }
            if (s == 0)
            {
              s = v[i] div A[i,j];
            }
            if (v[i] != s * A[i,j])
            {
              break;
            }
      }
      if (i == nrows(A))
      {
        answer[1] = j;
        answer[2] = s;
        // v = s * A_j
        return (answer);
      }
    }
   }

   // None of the columns of A divides v
   return (0);

}

/////////////////////////////////////////////////////////////////////


proc oneDimBelongSemigroup(int n, intvec v, list #)
"
USAGE:   oneDimBelongSemigroup(n,v[,m]); v is an integral vector,
         n is a positive integer[, m is a positive integer].
RETURN:  counters, a vector with nonnegative entries such that
         v*counters = n. If it does not exist such a vector, it returns 0.
         If a third parameter m is introduced, it will only consider the
         first m entries of v.
ASSUME:  v is an integral vector with positive entries.
EXAMPLE: example oneDimBelongSemigroup; shows some examples.
"
{
//--------------------------- initialisation ---------------------------------


   int i, j;
   int PartialSum;
   int num;
   int e = size(#);

   if  (e > 0)
   {
      num = #[1];
   }
   else
   {
      num = nrows(v);
   }

   intvec v2 = v[1..num];
   intvec counter, belong;
   belong[num] = 0;
   counter[num] = 0;
   for (i = 1; i <= num; i++)
   {
      if (n % v[i] == 0)
      {
         // ---- n is multiple of v[i]
         belong[i] = n div v[i];
         return (belong);
      }
   }

   if (num == 1)
   {
      // ---- num = 1 and n is not multiple of v[1] --> FALSE
      return(0);
   }

   PartialSum = 0;

   intvec w = sort(v2)[1];
   intvec cambio = sort(v2)[2];

   // ---- Iterative procedure to determine if n is in the semigroup generated by v
   while (1)
   {
      if (n  >= PartialSum)
      {
         if (((n - PartialSum) % w[1]) == 0)
         {
            // ---- n belongs to the semigroup generated by v,
            belong[cambio[1]] = (n - PartialSum) div w[1];
            for (j = 2; j <= num; j++)
            {
               belong[cambio[j]] = counter[j];
            }
            return(belong);
         }
      }
      i = num;
      while (!defined(end))
      {
         if (i == 1)
         {
            // ---- Stop, n is not in the semigroup
            return(0);
         }
         if (i > 1)
         {
            // counters control
            if (w[i] > n - PartialSum)
            {
               PartialSum = PartialSum - (counter[i]*w[i]);
               counter[i] = 0;
               i--;
            }
            else
            {
               counter[i] = counter[i] + 1;
               PartialSum = PartialSum + w[i];
               int end;
            }
         }
      }
      kill end;
   }
}
example
{ "EXAMPLE:";echo=2;
   int a = 95;
   intvec v = 18,51,13;
   oneDimBelongSemigroup(a,v);
   "// 95 = 1*18 + 1*25 + 2*13";
   oneDimBelongSemigroup(a,v,2);
   "// 95 is not a combination of 18 and 52;";
}



/////////////////////////////////////////////////////////////////////


proc SBelongSemigroup (intvec v, intmat A, list #)
"
USAGE:  SBelongSemigroup(v,A[,k]); v is an integral vector, A is an
        integral matrix, n is a positive integer.
RETURN: counters, a vector with nonnegative entries such that
        A*counters = v. If it does not exist such vector, it returns 0.
        If a third parameter k is introduced, it will only consider the
        first k columns of A.
ASSUME: A is an m x n matrix with nonnegative entries, n >= m, for every
        i,j <= m, i != j then A[i,j] = 0, v has nonnegative entries and
        nrows(v) = nrows(A);
"
{

   int last;
   int e = size(#);
   if  (e > 0)
   {
      last = #[1];
   }
   else
   {
      last = ncols(A);
   }

   intvec counters;
   counters[last] = 0;

   int i, j, k, l;
   intvec d;
   for (i = 1; i <= nrows(A); i++)
   {
      d[i] = A[i,i];
   }

   i = 1;
   while ((i < nrows(v)) && (v[i] % d[i] == 0))
   {
      i++;
   }
   if (v[i] % d[i] == 0)
   {
      // v is a combination of the first nrows(A) columns
      for (j = 1; j <= nrows(v); j++)
      {
         counters[j] = v[j] div d[j];
      }
      return(counters);
   }

   int gcdrow;
   for (i = 1; i <= nrows(A); i++)
   {
      gcdrow = d[i];
      for (j = nrows(A)+1; j <= last; j++)
      {
         if (A[i,j] != 0)
         {
            gcdrow = gcd(gcdrow,A[i,j]);
         }
      }
      if (v[i] % gcdrow != 0)
      {
         return (0);
      }
   }

   intvec swap;
   for (i = 1; i <= last; i++)
   {
      swap[i] = i;
   }
   for (i = nrows(A) + 1; i <= last; i++)
   {
      for (j = 1; j <= nrows(v); j++)
      {
         if (A[j,i] > v[j])
         {
            swap[i] = 0;
            for (k = 1; k <= nrows(A); k++)
            {
               A[k,i] = A[k,last];
            }
            swap[i] = swap[last];
            last--;
            i--;
            break;
         }
      }
   }
   if (nrows(A) == last)
   {
      return (0);
   }


   intvec order;
   order[last] = 0;
   for (i = nrows(A) + 1; i <= last; i++)
   {
      order[i] = 1;
      for (j = 1; j <= nrows(A); j++)
      {
         if (A[j,i] > 0)
         {
            order[i] = lcm(order[i], d[j] div gcd(A[j,i],d[j]));
         }
      }
   }

   intvec counters2;
   counters2[last] = 0;

   // A full enumeration is performed
   while(1)
   {
      i = nrows(counters2);
      while (1)
      {
         j = 1;
         if (counters2[i] < order[i] - 1)
         {
            while ((j < nrows(A)) and (v[j] >= A[j,i]))
            {
               j++;
            }
         }
         if ((v[j] < A[j,i])||(counters2[i] == order[i] - 1))
         {
            // A_j is not < v componentwise or counters2 = order[i]
            // we cannot increase counters2[i]
            if (counters2[i] != 0)
            {
               for (k = 1; k <= nrows(v); k++)
               {
                  v[k] = v[k] + counters2[i] * A[k,i];
               }
               counters2[i] = 0;
            }
            i--;
            if (i <= nrows(A))
            {
               // A*x = v has not nonnegative integral solution
               return(0);
            }
         }
         else
         {
            // j = nrows(A), then A_j < v (componentwise)
            // we add one unit to counters2[i]
            for (k = 1; k <= nrows(v); k++)
            {
               v[k] = v[k] - A[k,i];
            }
            counters2[i] = counters2[i] + 1;

            l = 1;
            while ((l < nrows(v)) and (v[l] % d[l] == 0))
            {
               l++;
            }
            if (v[l] % d[l] == 0)
            {
               // v is a combination of the first nrows(A) columns
               for (k = 1; k <= nrows(v); k++)
               {
                  counters[k] = v[k] div d[k];
               }
               for (k = nrows(v) + 1; k <= nrows(counters2); k++)
               {
                  counters[swap[k]] = counters2[k];
               }
               // A*counters = v
               return(counters);
            }

            break;
         }
      }
   }
}


/////////////////////////////////////////////////////////////////////



proc belongSemigroup (intvec v, intmat A, list #)
"
USAGE:   belongSemigroup(v,A[,k]); v is an integral vector, A is an
         integral matrix, n is a positive integer.
RETURN:  counters, a vector with nonnegative entries such that
         A*counters = v. If it does not exist such a vector, it returns 0.
         If a third parameter k is introduced, it will only consider the
         first k columns of A.
ASSUME:  A is a matrix with nonnegative entries, nonzero colums,
         v is a nonnegative vector and nrows(v) = nrows(A).
EXAMPLE: example belongSemigroup; shows some examples.
"
{
   int inputlast;
   int e = size(#);
   if  (e > 0)
   {
      inputlast = #[1];
   }
   else
   {
      inputlast = ncols(A);
   }

   int i, j, k;
   intvec counters;
   int last = inputlast;

   if (last == 0)
   {
      return (counters);
   }

   intvec swap;
   for (i = 1; i <= last; i++)
   {
      swap[i] = i;
   }
   for (i = 1; i <= last; i++)
   {
      for (j = 1; j <= nrows(v); j++)
      {
         if (A[j,i] > v[j])
         {
            swap[i] = 0;
            for (k = 1; k <= nrows(A); k++)
            {
               A[k,i] = A[k,last];
            }
            swap[i] = swap[last];
            last--;
            i--;
            break;
         }
      }
   }

   if (last == 0)
   {
      return (0);
   }

   intvec multip = Multiple(v, A, last);

   if (multip != 0)
   {
      intvec counters2;
      counters2[inputlast] = 0;
      // v = mult.value * ((mult.column)-th column of A)
      counters2[swap[multip[1]]] = multip[2];
      return (counters2);
   }

   counters[last] = 0;
   // A full enumeration is performed
   // Example: si v = (3,2,6), a1 = (1,0,0), a2 = (0,1,1), a3 = (0,0,1), then counters will take the
   // following values:
   // 000, 001, 002, 003, 004, 005, 006, 010, 011, 012, 013, 014, 015, 020, 021, 022, 023, 024 ---> 324
   while(1)
   {
      i = nrows(counters);
      while (1)
      {
         j = 1;
         while (j <= nrows(A))
         {
            if (v[j] < A[j,i])
            {
               break;
            }
            j++;
         }
         if (j <= nrows(A))
         {
            if (counters[i] != 0)
            {
               for (k = 1; k <= nrows(A); k++)
               {
                  v[k] = v[k] + counters[i] * A[k,i];
               }
               counters[i] = 0;
            }
            i--;
            if (i < 2)
            {
               // Does not belong
               return (0);
            }
         }
         else
         {
            for (k = 1; k <= nrows(A); k++)
            {
               v[k] = v[k] - A[k,i];
            }
            counters[i] = counters[i] + 1;
            multip = Multiple(v, A, i);
            if (multip != 0)
            {
               // v belongs, we return the solution counters so that A * counters = v
               counters[multip[1]] = multip[2];
               intvec counters2;
               counters2[inputlast] = 0;
               for (i = 1; i <= last; i++)
               {
                  counters2[swap[i]] = counters[i];
               }
               return (counters2);
            }

            break;
         }
      }
   }
}
example
{ "EXAMPLE:"; echo=2;
   intmat A[3][4] = 10,3,2,1,2,1,1,3,5,0,1,2;
   print(A);
   intvec v = 23,12,10;
   belongSemigroup(v,A);
   "// A * (1,3,1,2) = v";
   belongSemigroup(v,A,3);
   "// v is not a combination of the first 3 columns of A";
   intvec w = 12,4,1;
   belongSemigroup(w,A);
   "// w is not a combination of the columns of A";
}

/////////////////////////////////////////////////////////////////////



proc cardGroup(intmat A, list #)
"
USAGE:   cardGroup(A[,n]); A is a matrix with integral coefficients.
RETURN:         It returns a bigint. If we denote by ZA the group generated
         by the columns of the matrix A, then it returns the number of
         elements of the group of Z^m / ZA, where m = number of rows of A.
         If a second parameter n is introduced, it will
         only consider the first n columns of A. It returns 0 if Z^m / ZA
         is infinite; this is, when rank ZA < m.
EXAMPLE: example cardGroup; shows an example.
"
{
   int i, j, k, l;
   int coef1, coef2, aux;
   bigint aux2, aux3;
   list gcdiv;

   int last;
   int e = size(#);
   if  (e > 0)
   {
      last = #[1];
   }
   else
   {
      last = ncols(A);
   }

   // First we put the matrix A in diagonal form.
   for (i = 1; i <= nrows(A); i++)
   {
      j = i;
      while (A[i,j] == 0)
      {
         j++;
         if (j > last)
         {
            return (0);
            // Group is infinite
         }
      }
      if (j > i)
      {
         for (k = i; k <= nrows(A); k++)
         {
            // swap columns to have a nonzero pivot
            aux = A[k,i];
            A[k,i] = A[k,j];
            A[k,j] = aux;
         }
      }
      for (k = j+1; k <= last; k++)
      {
         if (A[i,k] != 0)
         {
            gcdiv = extgcd(A[i,i],A[i,k]);
            coef1 = A[i,k] div gcdiv[1];
            coef2 = A[i,i] div gcdiv[1];
            for (l = i; l <= nrows(A); l++)
            {
               // Perform elemental operations in the matrix
               // to put A in diagonal form
               aux2 = bigint(gcdiv[2]) * bigint(A[l,i]) + bigint(gcdiv[3]) * bigint(A[l,k]);
               aux3 = bigint(coef1) * bigint(A[l,i]) - bigint(coef2) * bigint(A[l,k]);
               A[l,k] = int(aux3);
               A[l,i] = int(aux2);
            }
         }
      }
   }

   // Once the matrix is in diagonal form, we only have to multiply
   // the diagonal elements

   bigint determinant = bigint(A[1,1]);
   bigint entry;
   for (i = 2; i <= nrows(A); i++)
   {
      entry = bigint(A[i,i]);
      determinant = determinant * entry;
   }
   determinant = absValue(determinant);
   return(determinant);

}
example
{ "EXAMPLE:"; echo=2;
   intmat A[3][5] = 24,  0,  0,  8, 3,
                     0, 24,  0, 10, 6,
                     0,  0, 24,  5, 9;
   cardGroup(A);
}

///////////////////////////////////////////////////////////////////////////////////////////////////////////

proc minMult(int a, intvec b)
"
USAGE:   minMult (a, b); a integer, b integer vector.
RETURN:  an integer k, the minimum positive integer such that k*a belongs to the
         semigroup generated by the integers in b.
ASSUME:  a is a positive integer, b is a vector of positive integers.
EXAMPLE: example minMult; shows an example.
"
{
//--------------------------- initialisation ---------------------------------
   int i, j, min, max;
   int n = nrows(b);

   if (n == 1)
   {
      // ---- trivial case
      return(b[1]/gcd(a,b[1]));
   }

   max = b[1];
   for (i = 2; i <= n; i++)
   {
      if (b[i] > max)
      {
         max = b[i];
      }
   }
   int NumNodes = a + max; //----Number of nodes in the graph

   int dist = 1;
   // ---- Auxiliary structures to obtain the shortest path between the nodes 1 and a+1 of this graph
   intvec queue = 1;
   intvec queue2;

   // ---- Control vector:
   //      control[i] = 0 -> node not reached yet
   //      control[i] = 1 -> node in queue1
   //      control[i] = 2 -> node in queue2
   //      control[i] = 3 -> node already processed
   intvec control;
   control[1] = 3;         // Starting node
   control[a + max] = 0;   // Ending node
   int current = 1;        // Current node
   int next;               // Node connected to corrent by arc (current, next)

   int ElemQueue, ElemQueue2;
   int PosQueue = 1;

   // Algoritmo de Dijkstra
   while (1)
   {
      if (current <= a)
      {
         // ---- current <= a, arcs are (current, current + b[i])
         for (i = 1; i <= n; i++)
         {
            next = current + b[i];
            if (next == a+1)
            {
               kill control;
               kill queue;
               kill queue2;
               return (dist);
            }
            if ((control[next] == 0)||(control[next] == 2))
            {
               control[next] = 1;
               queue = queue, next;
            }
         }
      }
      if (current > a)
      {
         // ---- current > a, the only possible ars is (current, current - a)
         next = current - a;
         if (control[next] == 0)
         {
            control[next] = 2;
            queue2[nrows(queue2) + 1] = next;
         }
      }
      PosQueue++;
      if (PosQueue <= nrows(queue))
      {
         current = queue[PosQueue];
      }
      else
      {
         dist++;
         if (control[a+1] == 2)
         {
            return(dist);
         }
         queue = queue2[2..nrows(queue2)];
         current = queue[1];
         PosQueue = 1;
         queue2 = 0;
      }
      control[current] = 3;
   }
}
example
{ "EXAMPLE:";echo=2;
   int a = 46;
   intvec b = 13,17,59;
   minMult(a,b);
   "// 3*a = 8*b[1] + 2*b[2]"
}


/////////////////////////////////////////////////////////////////////


static proc CheckMin (int posiblemin, intmat A, int column, list #)
"
USAGE:  CheckMin(posiblemin,A,column[,n]); posiblemin is an integer,
        A is an integral matrix and column and last are integers.
RETURN: 1 if posiblemin is the minimum value x such that x * (column-th colum of A)
        belongs to the semigroup generated by all the columns of A except
        A_column. It returns 0 otherwise. If an extra parameter n is
        introduced then it will only consider the first n columns of A.
ASSUME: 1 <= column <= ncols(A), A does not have negative entries or zero columns
"
{

   // If one can write (posiblemin-1)*A_column as a non-trivial combination of the
   // colums of A, then posiblemin is > to the real minimum
   intvec counters, multip;
   counters[ncols(A)] = 0;

   int i,j,k;
   int last;

   int e = size(#);
   if  (e > 0)
   {
      last = #[1];
   }
   else
   {
      last = ncols(A);
   }

   intvec v, aux;
   for (i = 1; i <= nrows(A); i++)
   {
      v[i] = (posiblemin-1)*A[i,column];
      // We swap A_column with A_1
      aux[i] = A[i,1];
      A[i,1] = A[i,column];
      A[i,column] = aux[i];
   }

   for (i = 2; i <= last; i++)
   {
      for (j = 1; j <= nrows(v); j++)
      {
         if (A[j,i] > v[j])
         {
            for (k = 1; k <= nrows(A); k++)
            {
               A[k,i] = A[k,last];
            }
            last--;
            i--;
            break;
         }
      }
   }

   // A full enumeration is performed
   while(1)
   {
      i = last;
      while (1)
      {
         j = 1;
         while (j <= nrows(A))
         {
            if (v[j] < A[j,i])
            {
               break;
            }
            j++;
         }
         if (j <= nrows(A))
         {
            if (counters[i] != 0)
            {
               for (k = 1; k <= nrows(A); k++)
               {
                  v[k] = v[k] + counters[i] * A[k,i];
               }
               counters[i] = 0;
            }
            i--;
            if (i == 1)
            {
               // The only solution is that v = (posiblemin-1)*A.col[1]
               return (1);
            }
         }
         else
         {
            for (k = 1; k <= nrows(A); k++)
            {
               v[k] = v[k] - A[k,i];
            }
            counters[i] = counters[i]+1;
            multip = Multiple(v, A, i);
            if (multip != 0)
            {
               return (0);
            }
            break;
         }
      }
   }
}

/////////////////////////////////////////////////////////////////////


static proc SimplicialCheckMin (int posiblemin, intmat A, int column, list #)
"
USAGE:  SimplicialCheckMin(posiblemin,A,column[,last]); posiblemin is an integer,
        A is an integral matrix and column and last are integers.
RETURN: 1 if posiblemin is less or equal to the minimum value x such that
        x * A_column belongs to the semigroup generated by all the columns of A
        except A_column. It returns 0 otherwise. If an extra parameter last is
        introduced then it will only consider the first n columns of A.
ASSUME: 1 <= column <= ncols(A), A does not have negative entries or zero columns
        A[i,j] = 0 for all 1 <= i,j <= nrows(A) where i != j
"
{
   // If one can write (posiblemin-1)*A_column as a non-trivial combination of the
   // colums of A, then posiblemin is > than the real minimum

   int last;
   int e = size(#);
   if  (e > 0)
   {
      last = #[1];
   }
   else
   {
      last = ncols(A);
   }

   int i, j, k, l;
   intvec d, v;
   for (i = 1; i <= nrows(A); i++)
   {
      d[i] = A[i,i];
   }

   for (i = 1; i <= nrows(A); i++)
   {
      v[i] = A[i,column] * (posiblemin-1);
   }

   i = 1;
   while ((i < nrows(v)) && (v[i] % d[i] == 0))
   {
      i++;
   }
   if (v[i] % d[i] == 0)
   {
      // v is a combination of the first nrows(A) columns
      return(0);
   }

   int aux;
   for (i = 1; i <= nrows(A); i++)
   {
      aux = A[i,nrows(A)+1];
      A[i,nrows(A)+1] = A[i,column];
      A[i,column] = aux;
   }

   for (i = nrows(A) + 2; i <= last; i++)
   {
       for (j = 1; j <= nrows(v); j++)
       {
          if (A[j,i] > v[j])
          {
             for (k = 1; k <= nrows(A); k++)
             {
                A[k,i] = A[k,last];
             }
             last--;
             i--;
             break;
          }
       }
    }

    intvec order;
    order[last] = 0;
    for (i = nrows(A) + 1; i <= last; i++)
    {
       order[i] = 1;
       for (j = 1; j <= nrows(A); j++)
       {
          if (A[j,i] > 0)
          {
             order[i] = lcm(order[i], d[j] div gcd(A[j,i],d[j]));
          }
       }
    }


    if (order[nrows(A)+1] < posiblemin-1)
    {
       return (0);
    }

    order[nrows(A)+1] = posiblemin-1;

    intvec counters;
    counters[last] = 0;

    // A full enumeration is performed
    while(1)
    {
       i = last;
       while (1)
       {
          j = 1;
          if (counters[i] < order[i] - 1)
          {
             while ((j < nrows(A)) and (v[j] >= A[j,i]))
             {
                j++;
             }
          }
          if ((v[j] < A[j,i])||(counters[i] == order[i] - 1))
          {
             // A_j is not < v componentwise or counters = order[i]-1
             // we cannot increase counters[i]
             if (counters[i] != 0)
             {
                for (k = 1; k <= nrows(v); k++)
                {
                   v[k] = v[k] + counters[i] * A[k,i];
                }
                counters[i] = 0;
             }
             i--;
            if (i <= nrows(A))
            {
               // A*x = v has not nonnegative integral solution different
               // from the trivial one
               return(1);
            }
         }
         else
         {
            // j = nrows(A), then A_j < v (componentwise)
            // we add one unit to counters[i]
            for (k = 1; k <= nrows(v); k++)
            {
               v[k] = v[k] - A[k,i];
            }
            counters[i] = counters[i] + 1;

            l = 1;
            while ((l < nrows(v)) and (v[l] % d[l] == 0))
            {
               l++;
            }
            if (v[l] % d[l] == 0)
            {
               // v is a combination of the first nrows(A) columns
               return(0);
            }

            break;
         }
      }
   }
}



/////////////////////////////////////////////////////////////////////


static proc Proportional(intvec a, intvec b)
"
USAGE:  Proportional(a,b); a,b integral vectors
RETURN: 1 if nrows(a) = nrows(b) and the vectors a and b are proportional;
        this is, there exist a rational number k such that k*a = b, and
        0 otherwise
ASSUME: a, b are nonzero vectors
"
{

   if (nrows(a) != nrows(b))
   {
      return (0)
   }

   int i, pivot;
   pivot = 1;
   while (a[pivot] == 0)
   {
      if (b[pivot] != 0)
      {
         // Not proportional
         return (0);
      }
      pivot++;
   }

   if (b[pivot] == 0)
   {
      // Not proportional
      return (0);
   }

   for (i = pivot + 1; i <= nrows(a); i++)
   {
      if (a[i] * b[pivot] != a[pivot] * b[i])
      {
         // Not proportional
         return (0);
      }
   }
   return (1);

}

/////////////////////////////////////////////////////////////////////


static proc StimatesMin(intmat A, int column, intvec line, list #)
"
USAGE:  StimatesMin(A,column,line[,n]); A is an integral matrix, column is
        an integer, line is an integral vector and n is an integer.
RETURN: The minimum integer k such that k * A_column belongs to the
        semigroup generated by all the columns of A except A_column. It
        returns 0 if it is not necessary to compute for the main program.
        If an extra parameter n is introuduced it considers the first n
        colums of A.
ASSUME: 1 <= column [<= n] <= ncols(A), A has nonnegative entries, line is
        a vector such that line[i] = line[j] if and only if the i-th and
        j-th columns of A are proportional
"
{

   int last;
   int e = size(#);
   if  (e > 0)
   {
      last = #[1];
   }
   else
   {
      last = nrows(line);
   }

   intvec current;
   int i,j,k;
   for (i = 1; i <= nrows(A); i++)
   {
      current[i] = A[i,column];
   }

   int nonzero = 1;
   while (current[nonzero] == 0)
   {
      nonzero++;
   }

   // We will only consider those colums A_j such that line[j] = line[column]
   intvec prop, jthcolumn;
   for (j = 1; j <= last; j++)
   {
      if (j != column)
      {
         if (line[column] == line[j])
         {
            prop[nrows(prop)+1] = A[nonzero,j];
         }
      }
   }
   int posiblemin = 0;
   if (prop[nrows(prop)] > 0)
   {
      // 1-dim minimum
      posiblemin = minMult(current[nonzero],prop);
   }

   if (line[column] <= nrows(A))
   {
      // It is not necessary to do CheckMin
      return(posiblemin);
   }


   if (posiblemin > 0)
   {
      if (SimplicialCheckMin(posiblemin, A, column, last))
      {
         // It is the real minimum, otherwise minimum < posiblemin
         return (posiblemin);
      }
   }
   // Not necessary to compute the minimum explicitly
   return (0);
}

/////////////////////////////////////////////////////////////////////


proc isCI(intmat A)
"
USAGE:   isCI(A); A is an integral matrix
RETURN:  1 if the simplicial toric ideal I(A) is a complete intersection
         and 0 otherwise. If printlevel > 0 and I(A) is a complete
         intersection it also shows a minimal set of generators of I(A)
ASSUME:  A is an m x n integral matrix with nonnegative entries and for
         every 1 <= i <= m, there exist a column in A whose i-th coordinate
         is not null and the rest are 0.
EXAMPLE: example isCI; shows some examples
"
{
//--------------------------- initialisation ---------------------------------

   intvec d;
   intvec minimum;
   intvec swap;                // swap[i] = j if and only if the i-th column of B equals the j-th column of A
   intvec line;                // line[i] = line[j] if and only if the i-th and the j-th column of B are proportional


   int n, m;           // Size of the input
   int i,j,k,l,t;

   n = ncols(A);
   m = nrows(A);
   intmat B[m][n];                                // auxiliary matrix
   intmat support[2*n-m][n];        // support[i,j] = 1 if and only if a_j belongs to V_i
   if (printlevel > 0)
   {
      ring r = 0,x(1..n),dp;
      ideal toric;                        // In case I(A) is a complete intersection, we obtain a
                                        // minimal set of generators
   }

   for (i = 1; i <= n; i++)
   {
      int zero = 0;
      swap[i] = i;
      for (j = 1; j <= m; j++)
      {
         B[j,i] = A[j,i];
         if (B[j,i] > 0)
         {
            zero++;
         }
         if (B[j,i] < 0)
         {
            print("// There are negative entries in the matrix");
            return (0);
         }
      }
      if (zero == 0)
      {
         print ("// There is a zero column in the matrix");
         return (0);
      }
   kill zero;
   }

//--------------------------- preprocessing the input ---------------------------------


   int aux, found;

   // We write B in a standard form; this is, for 1 <= i,j <= m, j != i then B[i,j] = 0
   for (i = 1; i <= m; i++)
   {
      j = i;
      found = 0;
      while (j <= n)
      {
         if (B[i,j] != 0)
         {
            k = 1;
            while ((k == i)||(B[k,j] == 0))
            {
               k++;
               if (k == m+1)
               {
                  for (l = 1; l <= m; l++)
                  {
                     aux = B[l,j];
                     B[l,j] = B[l,i];
                     B[l,i] = aux;
                  }
                  aux = swap[j];
                  swap[j] = swap[i];
                  swap[i] = aux;
                  found = 1;
                  break;
               }
            }
         }
         if (found == 1)
         {
            break;
         }
         j++;
      }
      if (j == n+1)
      {
         print("// There exists an i such that no column in A has the i-th coordinate positive and the rest are 0.");
            // It is not simplicial
         return (0);
      }
   }

   // Initialisation of variables
   int numgens = 0;                // number of generators built

   for (i = 1; i <= n; i++)
   {
      support[i,i] = 1;
   }

   intvec ithcolumn, jthcolumn, belong;
   int numcols = ncols(B);
   line[numcols] = 0;

   // line[i] = line[j] if and only if B_i and B_j are proportional.
   // Moreover for i = 1,...,nrows(B) we have that line[i]= i
   for (i = 1; i <= numcols; i++)
   {
      if (line[i] == 0)
      {
         for (j = 1; j <= nrows(B); j++)
         {
            ithcolumn[j] = B[j,i];
         }
         line[i] = i;
         for (j = i+1; j <= numcols; j++)
         {
            for (k = 1; k <= nrows(B); k++)
            {
               jthcolumn[k] = B[k,j];
            }
            if (Proportional(ithcolumn, jthcolumn))
            {
               line[j] = i;
            }
         }
      }
   }

//----------------- We apply reduction ---------------

   bigint det1, det2;
   int minim, swapiold, lineiold;
   int change = 1;

   det1 = cardGroup(B,numcols);
   while (change == 1)
   {
      change = 0;
      for (i = 1; i <= numcols; i++)
      {
         for (j = 1; j <= m; j++)
         {
            ithcolumn[j] = B[j,i];
            B[j,i] = B[j,numcols];
         }
         swapiold = swap[i];
         swap[i] = swap[numcols];
         lineiold = line[i];
         line[i] = line[numcols];
         det2 = cardGroup(B,numcols-1);
         minim = int(det2/det1);
         if (lineiold > m)
         {
            belong = SBelongSemigroup(minim*ithcolumn,B,numcols-1);
            if (belong != 0)
            {
               // It belongs, we remove the ith column
               if (printlevel > 0)
               {
                  // Create a generator
                  poly mon1 = x(swapiold)^(minim);
                  poly mon2 = 1;
                  for (j = 1; j <= nrows(belong); j++)
                  {
                     mon2 = mon2 * x(swap[j])^(belong[j]);
                  }
                  toric = toric, mon1-mon2;
                  kill mon1;
                  kill mon2;
               }
               det1 = det2;
               change = 1;
               numgens++;
               numcols--;
               i--;
            }
         }
         else
         {
            // line[i] <= m
            intvec semigroup, supportmon;
            int position;
            for (j = 1; j <= numcols-1; j++)
            {
               if (line[j] == lineiold)
               {
                  position = j;
                  semigroup = semigroup, B[lineiold,j];
                  supportmon = supportmon, swap[j];
               }
            }
            if (semigroup == 0)
            {
               belong = 0;
            }
            else
            {
               semigroup = semigroup[2..nrows(semigroup)];
               supportmon = supportmon[2..nrows(supportmon)];
               belong = oneDimBelongSemigroup(minim*ithcolumn[lineiold],semigroup);
            }
            if (belong != 0)
            {
               // It belongs, we remove the ith column
               if (printlevel > 0)
               {
                  // We create a generator
                  poly mon1,mon2;
                  mon1 = x(swapiold)^(minim);
                  mon2 = 1;
                  for (j = 1; j <= nrows(supportmon); j++)
                  {
                     mon2 = mon2 * x(supportmon[j])^(belong[j]);
                  }
                  toric = toric, mon1-mon2;
                  kill mon1, mon2;
               }
               det1 = det2;
               numcols--;
               numgens++;
               change = 1;
               if (i <= m)
               {
                  // We put again B in standard form
                  if (position != i)
                  {
                     for (j = 1; j <= m; j++)
                     {
                        aux = B[j,position];
                        B[j,position] = B[j,i];
                        B[j,i] = aux;
                     }
                     aux = swap[i];
                     swap[i] = swap[position];
                     swap[position] = aux;
                     aux = line[i];
                     line[i] = line[position];
                     line[position] = aux;
                  }
               }
               i--;
            }
            kill position;
            kill semigroup;
            kill supportmon;
         }
         if (belong == 0)
         {
            for (j = 1; j <= m; j++)
            {
               B[j,i] = ithcolumn[j];
            }
            swap[i] = swapiold;
            line[i] = lineiold;
         }
      }
   }

    // Initialisation of variables
    minimum[numcols] = 0;

//----------------- We run the first part of the algorithm ---------------

   // Estimation of m_i with i = 1,...,n
   for (i = 1; i <= numcols; i++)
   {
      minimum[i] = StimatesMin(B,i,line,numcols);
   }

   int nonzero;
   int stagenumber = 0;
   change = 1;
   while (change == 1)
   {
      // If for every i,j we have that m_i*B_i != m_j*B_j we leave the while loop
      change = 0;
      for (i = 1; (i <= numcols) && (change == 0); i++)
      {
         if (minimum[i] != 0)
         {
            for (j = i+1; (j <= numcols) && (change == 0); j++)
            {
               if ((minimum[j] != 0)&&(line[i] == line[j]))
               {
                  // We look for a nonzero entry in B_i
                  nonzero = 1;
                  while (B[nonzero,i] == 0)
                  {
                     nonzero++;
                  }
                  if (minimum[i]*B[nonzero,i] == minimum[j]*B[nonzero,j])
                  {
                     // m_i b_i = m_j b_j
                     numgens++;
                     stagenumber++;
                     if (swap[i] <= n)
                     {
                        // For k = swap[i], we have that V_k = {b_i}, so m_i b_i belongs to V_k
                              if (printlevel > 0)
                        {
                           poly mon1 = x(swap[i])^(minimum[i]);
                        }
                     }
                     else
                     {
                        // We check wether m_i b_i belongs to the semigroup generated by V_k
                        // where k = swap[i]. All vectors in V_k are proportional to b_i
                        intvec checkbelong;
                        int miai;
                        intvec supporti;
                        miai = minimum[i]*B[nonzero,i];
                        for (k = 1; k <= n; k++)
                        {
                           if (support[swap[i],k])
                           {
                              supporti = supporti, k;
                              checkbelong[nrows(supporti)-1] = A[nonzero,k];
                           }
                        }
                        // 1-dim belong semigroup
                        belong = oneDimBelongSemigroup(miai,checkbelong);
                        if (belong == 0)
                        {
                           // It does not belong
                           print ("// It is NOT a complete intersection");
                           return (0);
                        }
                        if (printlevel > 0)
                        {
                           poly mon1 = 1;
                           // It belongs, we construct a monomial of the new generator
                           for (k = 1; k < nrows(supporti); k++)
                           {
                              mon1 = mon1 * x(supporti[k+1])^(belong[k]);
                           }
                        }
                        kill miai;
                        kill checkbelong;
                        kill supporti;
                     }
                     if (swap[j] <= n)
                     {
                        // For k = swap[j], we have that V_k = {b_j}, so m_j b_j belongs to V_k
                        if (printlevel > 0)
                        {
                           poly mon2 = x(swap[j])^(minimum[j]);
                           toric = toric, mon1-mon2;
                           kill mon1;
                           kill mon2;
                        }
                     }
                     else
                     {
                        // We check wether m_j b_j belongs to the semigroup generated by V_k
                        // where k = swap[j]. All vectors in V_k are proportional to b_j
                        intvec checkbelong;
                        int mjaj;
                        intvec supportj;
                        nonzero = 1;
                        while (B[nonzero,j] == 0)
                        {
                           nonzero++;
                        }
                        mjaj = minimum[j]*B[nonzero,j];
                        for (k = 1; k <= n; k++)
                        {
                           if (support[swap[j],k])
                           {
                              supportj = supportj, k;
                              checkbelong[nrows(supportj)-1] = A[nonzero,k];
                           }
                        }
                        // 1-dim belong semigroup
                        belong = oneDimBelongSemigroup(mjaj,checkbelong);
                        if (belong == 0)
                        {
                           // It does not belong
                           print ("// It is NOT a complete intersection");
                           return (0);
                        }
                        if (printlevel > 0)
                        {
                           poly mon2 = 1;
                           // It belongs, we construct the second monomial of the generator
                           for (k = 1; k < nrows(supportj); k++)
                           {
                              mon2 = mon2 * x(supportj[k+1])^(belong[k]);
                           }
                           toric = toric,mon1-mon2;
                           kill mon1;
                           kill mon2;
                        }
                        kill checkbelong;
                        kill mjaj;
                        kill supportj;
                     }

                     // Now we remove b_i, b_j from B and we add gcd(b_i,b_j)
                     change = 1;
                     for (k = 1; k <= n; k++)
                     {
                        // V_{support.nrows} = V_i + V_j
                        support[n+stagenumber,k] = support[swap[i],k] + support[swap[j],k];
                     }
                     // line[i] does not change
                     line[j] = line[numcols];
                     swap[i] = n+stagenumber;
                     swap[j] = swap[numcols];
                     k = 1;
                     while (B[k,i] == 0)
                     {
                        k++;
                     }
                     int dp;
                     dp = gcd(B[k,i], B[k,j]);
                     int factor = B[k,i] div dp;
                     B[k,i] = dp;
                     k++;
                     kill dp;
                     while (k <= nrows(B))
                     {
                        B[k,i] = B[k,i] / factor;
                        k++;
                     }
                     kill factor;
                     for (k = 1; k <= nrows(B); k++)
                     {
                        B[k,j] = B[k,numcols];
                     }
                     minimum[j] = minimum[numcols];
                     numcols--;
                     // We compute a new m_i
                     minimum[i] = StimatesMin(B,i,line,numcols);
                  }
               }
            }
         }
      }
   }

//--------------------------- We apply reduction ---------------------------------

   intvec minZ;
   for (i = 1; i <= n+stagenumber; i++)
   {
      minZ[i] = 1;
   }
   det1 = cardGroup(B,numcols);
   int posaux;
   intvec counters1,counters2;

   while (numgens < n - m)
   {
      change = 0;
      i = nrows(B) + 1;
      while ((i <= numcols)&&(numgens < n-m))
      {
         // The vector we are going to study
         for (j = 1; j <= nrows(B); j++)
         {
            ithcolumn[j] = B[j,i];
            B[j,i] = B[j,numcols];
         }
         posaux = swap[i];
         swap[i] = swap[numcols];

         numcols--;

         det2 = cardGroup(B,numcols);
         minim = int(det2/det1);

         if (minim != 1)
         {
            // If minimum = 1, we have previously checked that current does not belong to
            // the semigroup generated by the columns of B.

            for (j = 1; j <= nrows(B); j++)
            {
               ithcolumn[j] = ithcolumn[j]*minim;
            }
            minZ[posaux] = minim * minZ[posaux];
            det1 = det2;  // det1 *= minimum

            intvec support1, support2;
            nonzero = 1;
            while (ithcolumn[nonzero] == 0)
            {
               nonzero++;
               if (nonzero > nrows(ithcolumn))
               {
                  print("// ERROR");
                  return (0);
               }
            }
            int currentnonzero = ithcolumn[nonzero];
            intvec B1;
            for (j = 1; j <= n; j++)
            {
               if (support[posaux,j])
               {
                  // a_j is proportional to b_i = a_posaux
                  support1 = support1, j;
                  B1[nrows(support1)-1] = A[nonzero,j];
               }
            }
            // 1-dim belongsemigroup
            counters1 = oneDimBelongSemigroup(currentnonzero, B1);
            kill currentnonzero;
            kill B1;
            if (counters1 != 0)
            {
               intmat B2[m][n];
               // It belongs, now we have to check if current belongs to the semigroup
               // generated by the columns of B2
               for (j = 1; j <= numcols; j++)
               {
                  for (k = 1; k <= n; k++)
                  {
                     if (support[swap[j],k])
                     {
                        // a_k may not be proportional to b_i = a_posaux
                        support2 = support2, k;
                        for (l = 1; l <= m; l++)
                        {
                           B2[l,nrows(support2)-1] = A[l,k];
                        }
                     }
                  }
               }
               // We write B2 in standard form
               for (l = 1; l <= m; l++)
               {
                  j = l;
                  found = 0;
                  while (found == 0)
                  {
                     if (B2[l,j] != 0)
                     {
                        k = 1;
                        while ((k == l)||(B2[k,j] == 0))
                        {
                           k++;
                           if (k == m + 1)
                           {
                              for (t = 1; t <= m; t++)
                              {
                                 jthcolumn[t] = B2[t,j];
                                 B2[t,j] = B2[t,l];
                                 B2[t,l] = jthcolumn[t];
                              }
                              aux = support2[j+1];
                              support2[j+1] = support2[l+1];
                              support2[l+1] = aux;
                              found = 1;
                              break;
                           }
                        }
                     }
                     j++;
                  }
               }

               // m-dim belong semigroup
               counters2 = SBelongSemigroup(ithcolumn, B2, nrows(support2)-1);
               kill B2;

               if (counters2 != 0)
               {
                  // current belongs, we construct the new generator
                  numgens++;
                  if (printlevel > 0)
                  {
                     poly mon1, mon2;
                     mon1 = 1;
                     mon2 = 1;
                     for (j = 1; j < nrows(support1); j++)
                     {
                        mon1 = mon1 * x(support1[j+1])^(counters1[j]);
                     }
                     for (j = 1; j < nrows(support2); j++)
                     {
                        mon2 = mon2 * x(support2[j+1])^(counters2[j]);
                     }
                     toric = toric, mon1-mon2;
                     kill mon1;
                     kill mon2;
                  }

                  if (i != numcols)
                  {
                     i--;
                  }
                  change = 1; // We have removed a column of B

               }
            }
            kill support1,support2;
         }

         if ((counters1 == 0)||(counters2 == 0)||(minim == 1))
         {
            // We swap again the columns in B
            for (j = 1; j <= m; j++)
            {
               B[j,i] = ithcolumn[j];
            }
            numcols++;
            swap[numcols] = swap[i];
            swap[i] = posaux;
         }

         if ((i == numcols)&&(change == 0))
         {
            print(" // It is NOT a Complete Intersection.");
            return (0);
         }
         i++;
      }
   }

   // We have removed all possible columns
   if (printlevel > 0)
   {
      print("// Generators of the toric ideal");
      toric = simplify(toric,2);
      toric;
   }
   print("// It is a complete intersection");
   return(1);
}
example
{ "EXAMPLE:"; echo=2;
   intmat A[2][5] = 60,0,140,150,21,0,60,140,150,21;
   print(A);
   printlevel = 0;
   isCI(A);
   printlevel = 1;
   isCI(A);
   intmat B[3][5] = 12,0,0,1,2,0,10,0,3,2,0,0,8,3,3;
   print(B);
   isCI(B);
   printlevel=0;
}