// 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; }