binomialCD.sing

LIB "elim.lib";
proc binomialCD (ideal I)
{
  ideal sI = std(I);
  list rl = ringlist(basering);
  list Answer; // List that contains the decomposition
  list rvars = variables(maxideal(1)); // the ring variables
  ideal dummy = 1;
  list L = dummy,rvars,I; // Entry of the todolist
  ideal IntersectAnswer = (1); // Intersection of the answers up to now
  list todo = insert(list(),L);
  
  int compcount = 0;
  while (size(todo) > 0 ) {
    L = todo[1]; // Will decompose first entry of the todolist
    todo = delete(todo,1); // Shorten todolist
    if (size (reduce (IntersectAnswer , std(L[3]))) == 0){
      // This component is redundant
      print ("Redundant component");
      continue; // next on the todolist
    }
    else {
      // Test if there are variables remaining to check:
      if (size(L[2]) == 0) {
	  // We have an answer:
	  compcount++;
	  ideal newone = mstd(L[3])[2]; // Ideally compute minimal generating set
	  print ("component: ");
	  print (compcount);
	  Answer = insert(Answer, newone);
	  IntersectAnswer = intersect(IntersectAnswer, newone);
	  // Checks if we have everything:
	  if (matrix(std(IntersectAnswer)) == matrix (sI)){
	    todo = list();
	  }
	}
      else {
	// There is a remaining variable
	//	print (L);
	//	print (L[2]);
        def i = L[2][1];
	//	print ("Will now handle variable:");
	//	print (i);
	def newvars = delete(L[2],1);
	def result = sat (L[3],i);
	def J = result[1]; // The saturated ideal
	def k = result[2]; // The saturation exponent
	if (k>0) { // a saturation was needed:
	  ideal J2 = L[3] + i^k;
	  // Remove L[1] components from variables considiered
	  // earlier:
	  J2 = sat(J2,L[1])[1];
	  if (matrix(dummy) != matrix(std(J2))) {
	    // Something left ? New todo-item:
	    list newitem = L[1], newvars, J2;
	    todo = insert(todo, newitem );
	  }
	}
	// Continue with the next variable:
	if (matrix(std(J)) != matrix(dummy)) {
	  list newitem = L[1]*i,newvars,J;
	  todo = insert (todo, newitem);
	}
      }
    }
  }
  print ("Decomposition done");
  return (Answer);
}

int t=timer;
system("--ticks-per-sec", 100);         // give time in 1/100 sec

ring R = 0,(R00,U00,R01,D01,U01,R02,D02,R10,L10,U10,R11,L11,D11,U11,R12,L12,D12,L20,U20,L21,D21,U21,L22,D22),dp;
ideal I = U00*R01-R00*U10,R01*D11-D01*R00,D11*L10-L11*D01,L10*U00-U10*L11,U01*R02-R01*U11,R02*D12-D02*R01,D12*L11-L12*D02,L11*U01-U11*L12,U10*R11-R10*U20,R11*D21-D11*R10,D21*L20-L21*D11,L20*U10-U20*L21,U11*R12-R11*U21,R12*D22-D12*R11,D22*L21-L22*D12,L21*U11-U21*L22;

binomialCD(I);

print ("Time needed in 1/100 seconds");
print (timer-t);

quit();

// Questions:
// -) What is a fast and nice way to check I == J for ideals ?
// -) How can I compute minimal generating sets and avoid lenghty 
//    Groebner Bases?
// -) How do I convert the string "x" to the variable "x"?