Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I would like to compute, using Groebner bases, the dimension of the variety defined by a set of polynomials in several variables. In the wikipedia page the method is described, and an implementation is possible, but I was hoping there was a predefined function that would do that.

share|improve this question
    
Is this about math or Mathematica? –  Yves Klett Nov 14 '13 at 12:52
    
Dear Yves, this was a question about Mathematica, I wanted to know if there is any predefined function that computes the dimension by implementing the algorithm cited above. –  Ion Nechita Nov 15 '13 at 17:05
    
Thank you for the info - in that case you are in good hands with Daniel :) –  Yves Klett Nov 16 '13 at 6:45

2 Answers 2

up vote 6 down vote accepted

I show some aged code cribbed from this MathGroup post (See also this library item). There is some amount of explanation in each. But the most salient quote, from the post, is this. "It is unlikely that I can do justice to the explanation. So far as correctness goes, let's just say the code is on the honor system."

firstContainsSecond[l1_, l2_] := (Union[l1, l2] === l1)

isIndependentSet[set_, sets_] := 
 Map[! firstContainsSecond[set, #] &, sets]

getMaxIndependentSets[vars_, inset_, heds_, maxlen_, indx_, sets_] := 
 Module[{currentset, vlen = Length[vars], ilen = Length[inset], 
   enlarged = False, newmax = maxlen, maxsets = sets}, 
  Do[If[ilen + vlen - i + 1 <= maxlen, Break[]];
   currentset = Append[inset, vars[[i]]];
   If[isIndependentSet[currentset, 
     heds], {maxsets, enlarged, newmax} = 
     getMaxIndependentSets[vars, currentset, heds, newmax, i + 1, 
      maxsets];
    If[! enlarged, maxsets = {maxsets, currentset};
     newmax = Max[newmax, Length[currentset]];
     enlarged = True;];], {i, indx, vlen}];
  {maxsets, enlarged, newmax}]

idealDimension[ideal_, vars_] := 
 Module[{ord = MonomialOrder -> DegreeReverseLexicographic, gb, 
   pheads, heds, maxsets, el, mlen}, 
  gb = GroebnerBasis[ideal, vars, ord];
  pheads = 
   Map[First[First[#]] &, 
    First[Internal`DistributedTermsList[gb, vars, ord]]];
  heds = Apply[
    And, (Map[(vars*#) &, (pheads /. _?Positive -> 1)]/.0 -> 
      Sequence[])];
  {maxsets, el, mlen} = 
   getMaxIndependentSets[vars, {}, heds, 0, 1, {}];
  mlen]

This comes with no warranties. If ever there were any, they must have expired by now.

share|improve this answer
    
Thanks a lot Daniel! I will try the code and let you know how it worked ;) –  Ion Nechita Nov 15 '13 at 10:25
    
Daniel, I managed to get the code working in Mathematica 9, I will post it here in a couple of days after some debugging. Thank you very much for your help! –  Ion Nechita Nov 16 '13 at 8:27

Here's the code kindly provided by Daniel Lichtblau in the accepted answer, modified to work with Mathematica 9. All the credit goes to him, and, of course, this comes with no warranties ;) The code is inspired from T.Becker, V.Weisphenning - Groebner bases, (Table 9.6, pp 449).

firstContainsSecond[l1_, l2_] := (Union[l1, l2] === l1);
isIndependentSet[set_, sets_] := 
  Map[! firstContainsSecond[set, #] &, sets];
getMaxIndependentSets[vars_, inset_, heds_, maxlen_, indx_, sets_] := 
  Module[{currentset, vlen = Length[vars], ilen = Length[inset], 
    enlarged = False, newmax = maxlen, maxsets = sets},
   Do[
    If[ilen + vlen - i + 1 <= maxlen, Break[]];
    currentset = Append[inset, vars[[i]]];
    If[And @@ 
      isIndependentSet[currentset, heds], {maxsets, enlarged, 
       newmax} = 
      getMaxIndependentSets[vars, currentset, heds, newmax, i + 1, 
       maxsets];
     If[! enlarged, maxsets = Append[maxsets, currentset];
      newmax = Max[newmax, Length[currentset]];
      enlarged = True;];], {i, indx, vlen}];
   {maxsets, enlarged, newmax}
   ];
idealDimension[ideal_, vars_] := 
  Module[{ord = DegreeReverseLexicographic, gb, pheads, heds, maxsets,
     el, mlen}, 
   gb = GroebnerBasis[ideal, vars, MonomialOrder -> ord];
   pheads = Map[First[MonomialList[#, vars, ord]] &, gb];
   heds = Map[Variables[#] &, pheads];
   {maxsets, el, mlen} = 
    getMaxIndependentSets[vars, {}, heds, 0, 1, {}];
   mlen
   ];
share|improve this answer
    
Thanks. I had forgotten about the context change to DistributedTermsList. And the "new" MonomialList did not exist when I wrote that. Ironically, version 3 also had a MonomialList that was related to the current one. I no longer recall what internal misgivings caused it to go away, only to reappear a few versions later. –  Daniel Lichtblau Nov 16 '13 at 23:16

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.