1 #include "dimension.h"
2 #include "buchberger.h"
3 #include "log.h"
4 #include "printer.h"
5
6 const int bitsPerWord=64;
7 typedef int64 ExponentType;
8
allOnes(int n)9 ExponentType allOnes(int n)
10 {
11 if(n<64)return (((int64) 1)<<n)-1;
12 return -1;
13 }
14
sum(ExponentType v)15 int sum(ExponentType v)
16 {
17 int64 t0=v;
18 int64 t1=(t0&0x5555555555555555)+((t0>>1)&0x5555555555555555);
19 int64 t2=(t1&0x3333333333333333)+((t1>>2)&0x3333333333333333);
20 int64 t3=(t2&0x0f0f0f0f0f0f0f0f)+((t2>>4)&0x0f0f0f0f0f0f0f0f);
21 int64 t4=(t3&0x00ff00ff00ff00ff)+((t3>>8)&0x00ff00ff00ff00ff);
22 int64 t5=(t4&0x0000ffff0000ffff)+((t4>>16)&0x0000ffff0000ffff);
23 int64 t6=(t5&0x00000000ffffffff)+((t5>>32)&0x00000000ffffffff);
24 return t6;
25 }
26
rek64(ExponentType & ones,ExponentType & zeros,int nOnes,int nZeros,vector<ExponentType> vectors,int & best,int n)27 static void rek64(ExponentType &ones, ExponentType &zeros, int nOnes, int nZeros, vector<ExponentType> vectors, int &best, int n)
28 {
29 if(nOnes>best)best=nOnes;
30 if(n-nZeros<best)return;
31 if(nOnes+nZeros==n)return;
32
33 #if 1
34 int index=0;
35 for(int i=0;i<n;i++,index++)
36 // if((!ones[i])&&(!zeros[i]))break;
37 if((!((((int64)1)<<i)&ones))&&(!((((int64)1)<<i)&zeros)))break;
38 #else
39 int index=nOnes+nZeros;
40 #endif
41 assert(index<n);
42
43 ones|=((int64)1)<<index;
44 bool good=true;
45 for(vector<ExponentType>::const_iterator i=vectors.begin();i!=vectors.end();i++)
46 if(((*i)&~ones)==0)
47 {
48 good=false;
49 break;
50 }
51
52 if(good)
53 {
54 rek64(ones,zeros,nOnes+1,nZeros,vectors,best,n);
55 }
56 ones-=((int64)1)<<index;
57
58 vector<ExponentType> vectorsSubset;
59 vectorsSubset.reserve(vectors.size());
60
61 for(vector<ExponentType>::const_iterator i=vectors.begin();i!=vectors.end();i++)
62 {
63 // if((*i)[index]==0)vectorsSubset.push_back(*i);
64 if(((*i)&(((int64)1)<<index))==0)vectorsSubset.push_back(*i);
65 }
66
67 zeros|=((int64)1)<<index;
68 rek64(ones,zeros,nOnes,nZeros+1,vectorsSubset,best,n);
69 zeros-=((int64)1)<<index;
70 }
71 /*
72 * monomialGenerators should have zero-one exponent vectors with at most 64 variables in the ring.
73 */
krullDimensionOfMonomialIdeal64(PolynomialSet const & monomialGenerators)74 int krullDimensionOfMonomialIdeal64(PolynomialSet const &monomialGenerators)
75 {
76 int n=monomialGenerators.getRing().getNumberOfVariables();
77 assert(n<=64);
78 vector<ExponentType> generators;
79 for(PolynomialSet::const_iterator i=monomialGenerators.begin();i!=monomialGenerators.end();i++)
80 {
81 ExponentType v=0;
82 for(int j=0;j<n;j++)v=2*v+i->getMarked().m.exponent[j];
83 generators.push_back(v);
84 }
85
86 int best=0;
87
88 ExponentType zeros=0;
89 ExponentType ones=0;
90
91 ExponentType possiblyOne=0;
92 for(vector<ExponentType>::const_iterator i=generators.begin();i!=generators.end();i++)possiblyOne=possiblyOne|*i;
93 ones=allOnes(n)-possiblyOne;
94
95 rek64(ones,zeros,sum(ones),sum(zeros),generators,best,n);
96
97 return best;
98 }
99
100
101
radicalOfMonomialIdeal(PolynomialSet const & monomialGenerators)102 PolynomialSet radicalOfMonomialIdeal(PolynomialSet const &monomialGenerators)
103 {
104 PolynomialRing theRing=monomialGenerators.getRing();
105 PolynomialSet temp=monomialGenerators;
106
107 temp.markAndScale(LexicographicTermOrder()); //just to make sure that some term is marked
108
109 PolynomialSet ret(theRing);
110 for(PolynomialSet::const_iterator i=temp.begin();i!=temp.end();i++)
111 {
112 IntegerVector e=i->getMarked().m.exponent;
113 e=e.supportVector();
114 ret.push_back(Polynomial(Term(i->getMarked().c,Monomial(theRing,e))));
115 }
116 return ret;
117 }
118
increase(IntegerVector & v,int & numberOfOnes)119 static bool increase(IntegerVector &v, int &numberOfOnes)
120 {
121 int i=0;
122 while(i<v.size() && v[i]==1)
123 {
124 v[i]=0;
125 numberOfOnes--;
126 i++;
127 }
128 if(i==v.size())return false;
129 v[i]=1;
130 numberOfOnes++;
131 return true;
132 }
133
134 static int rekCallls;
rek(IntegerVector & ones,IntegerVector & zeros,int nOnes,int nZeros,IntegerVectorList const & vectors,int & best)135 static void rek(IntegerVector &ones, IntegerVector &zeros, int nOnes, int nZeros, IntegerVectorList const &vectors, int &best)
136 {
137 rekCallls++;
138 if(nOnes>best)best=nOnes;
139 if(ones.size()-nZeros<best)return;
140 if(nOnes+nZeros==ones.size())return;
141
142
143 log3
144 {
145 fprintf(Stderr,"Ones:\n");
146 AsciiPrinter(Stderr).printVector(ones);
147 fprintf(Stderr,"Zeros:\n");
148 AsciiPrinter(Stderr).printVector(zeros);
149
150 AsciiPrinter(Stderr).printVectorList(vectors);
151 }
152
153 int index=0;
154 for(int i=0;i<ones.size();i++,index++)
155 if((!ones[i])&&(!zeros[i]))break;
156
157 assert(index<ones.size());
158
159 ones[index]=1;
160 bool good=true;
161 for(IntegerVectorList::const_iterator i=vectors.begin();i!=vectors.end();i++)
162 if(i->divides(ones))
163 {
164 good=false;
165 break;
166 }
167
168 if(good)
169 {
170 rek(ones,zeros,nOnes+1,nZeros,vectors,best);
171 }
172 ones[index]=0;
173
174 IntegerVectorList vectorsSubset;
175
176 for(IntegerVectorList::const_iterator i=vectors.begin();i!=vectors.end();i++)
177 {
178 if((*i)[index]==0)vectorsSubset.push_back(*i);
179 }
180
181 if(0) // It is not clear that this is an improvement, even with an better implementation
182 {
183 IntegerVector newOnes(ones.size());
184 IntegerVector possiblyOne(ones.size());
185 for(IntegerVectorList::const_iterator i=vectors.begin();i!=vectors.end();i++)possiblyOne=max(possiblyOne,*i);
186
187 if((IntegerVector::allOnes(ones.size())-possiblyOne-ones-zeros).max()>0)
188 {
189 /* debug<<vectorsSubset;
190 debug<<"ones now\n"<<ones<<"\n";
191 debug<<"zeros now\n"<<zeros<<"\n";
192 debug<<"FORCED ONES:"<<IntegerVector::allOnes(ones.size())-possiblyOne<<"\n";
193 */ newOnes=max(IntegerVector::allOnes(ones.size())-possiblyOne-ones-zeros,IntegerVector(ones.size()));
194 // debug<<"NEW ONES:"<<newOnes<<"\n";
195 }
196
197 nOnes+=newOnes.sum();
198 ones+=newOnes;
199
200 zeros[index]=1;
201 rek(ones,zeros,nOnes,nZeros+1,vectorsSubset,best);
202
203 ones-=newOnes;
204 nOnes-=newOnes.sum();
205
206 zeros[index]=0;
207 }
208 else
209 {
210 zeros[index]=1;
211 rek(ones,zeros,nOnes,nZeros+1,vectorsSubset,best);
212 zeros[index]=0;
213 }
214
215 }
216
217
krullDimensionOfMonomialIdeal(PolynomialSet const & monomialGenerators)218 int krullDimensionOfMonomialIdeal(PolynomialSet const &monomialGenerators)
219 {
220 // debug<<"Taking radical\n";
221 PolynomialSet temp=radicalOfMonomialIdeal(monomialGenerators);
222 // debug<<"Minimizing\n";
223 minimize(&temp);
224 // debug<<"Done\n";
225 if(monomialGenerators.getRing().getNumberOfVariables()<=64)return krullDimensionOfMonomialIdeal64(temp);
226 IntegerVectorList vectors;
227 for(PolynomialSet::const_iterator i=temp.begin();i!=temp.end();i++)
228 vectors.push_back(i->getMarked().m.exponent);
229
230 int best=0;
231
232 int n=monomialGenerators.getRing().getNumberOfVariables();
233 IntegerVector zeros(n);
234 IntegerVector ones(n);
235
236 //Preprocessing step. This can be improved.
237 IntegerVector possiblyOne(n);
238 for(IntegerVectorList::const_iterator i=vectors.begin();i!=vectors.end();i++)possiblyOne=max(possiblyOne,*i);
239 ones=IntegerVector::allOnes(n)-possiblyOne;
240
241 rek(ones,zeros,ones.sum(),zeros.sum(),vectors,best);
242 //debug<<"REKCALLS"<<rekCallls<<"\n";
243
244 return best;
245 }
246
krullDimension(PolynomialSet const & groebnerBasis)247 int krullDimension(PolynomialSet const &groebnerBasis)
248 {
249 return krullDimensionOfMonomialIdeal(groebnerBasis.markedTermIdeal());
250 }
251