1 | // emacs edit mode for this file is -*- C++ -*- |
---|
2 | // $Id: fglmhom.cc,v 1.7 1998-01-05 16:39:18 Singular Exp $ |
---|
3 | |
---|
4 | /**************************************** |
---|
5 | * Computer Algebra System SINGULAR * |
---|
6 | ****************************************/ |
---|
7 | /* |
---|
8 | * ABSTRACT - The FGLM-Algorithm extended for homogeneous ideals. |
---|
9 | * Calculates via the hilbert-function a groebner basis. |
---|
10 | */ |
---|
11 | |
---|
12 | #ifndef NOSTREAMIO |
---|
13 | #include <iostream.h> |
---|
14 | #endif |
---|
15 | |
---|
16 | #include "mod2.h" |
---|
17 | #ifdef HAVE_FGLM |
---|
18 | #include "tok.h" |
---|
19 | #include "structs.h" |
---|
20 | #include "subexpr.h" |
---|
21 | #include "polys.h" |
---|
22 | #include "ideals.h" |
---|
23 | #include "ring.h" |
---|
24 | #include "ipid.h" |
---|
25 | #include "ipshell.h" |
---|
26 | #include "febase.h" |
---|
27 | #include "maps.h" |
---|
28 | #include "mmemory.h" |
---|
29 | #include "fglm.h" |
---|
30 | #include "fglmvec.h" |
---|
31 | #include "fglmgauss.h" |
---|
32 | #include "intvec.h" |
---|
33 | #include "kstd1.h" |
---|
34 | #include "stairc.h" // -> hHStdSeries, hFirstSeries usw. |
---|
35 | #include <templates/ftmpl_list.h> |
---|
36 | |
---|
37 | // obachman: Got rid off those "redefiende messages by includeing fglm.h |
---|
38 | #include "fglm.h" |
---|
39 | #if 0 |
---|
40 | #define PROT(msg) if (BTEST1(OPT_PROT)) Print(msg) |
---|
41 | #define STICKYPROT(msg) if (BTEST1(OPT_PROT)) Print(msg) |
---|
42 | #define PROT2(msg,arg) if (BTEST1(OPT_PROT)) Print(msg,arg) |
---|
43 | #define STICKYPROT2(msg,arg) if (BTEST1(OPT_PROT)) Print(msg,arg) |
---|
44 | #define fglmASSERT(ignore1,ignore2) |
---|
45 | #endif |
---|
46 | |
---|
47 | struct doublepoly |
---|
48 | { |
---|
49 | poly sm; |
---|
50 | poly dm; |
---|
51 | }; |
---|
52 | |
---|
53 | class homogElem |
---|
54 | { |
---|
55 | public: |
---|
56 | doublepoly mon; |
---|
57 | fglmVector v; |
---|
58 | fglmVector dv; |
---|
59 | int basis; |
---|
60 | int destbasis; |
---|
61 | BOOLEAN inDest; |
---|
62 | homogElem() : v(), dv(), basis(0), destbasis(0), inDest(FALSE) {} |
---|
63 | homogElem( poly m, int b, BOOLEAN ind ) : |
---|
64 | basis(b), inDest(ind) |
---|
65 | { |
---|
66 | mon.dm= m; |
---|
67 | mon.sm= NULL; |
---|
68 | } |
---|
69 | }; |
---|
70 | |
---|
71 | struct homogData |
---|
72 | { |
---|
73 | ideal sourceIdeal; |
---|
74 | doublepoly * sourceHeads; |
---|
75 | int numSourceHeads; |
---|
76 | ideal destIdeal; |
---|
77 | int numDestPolys; |
---|
78 | homogElem * monlist; |
---|
79 | int monlistmax; |
---|
80 | int monlistblock; |
---|
81 | int numMonoms; |
---|
82 | int basisSize; |
---|
83 | int overall; // nur zum testen. |
---|
84 | int numberofdestbasismonoms; |
---|
85 | // homogData() : sourceHeads(NULL), numSourceHeads(0), monlist(NULL), |
---|
86 | // numMonoms(0), basisSize(0) {} |
---|
87 | }; |
---|
88 | |
---|
89 | int |
---|
90 | hfglmNextdegree( intvec * source, ideal current, int & deg ) |
---|
91 | { |
---|
92 | int numelems; |
---|
93 | intvec * newhilb = hHstdSeries( current, NULL, currQuotient ); |
---|
94 | |
---|
95 | loop |
---|
96 | { |
---|
97 | if ( deg < newhilb->length() ) |
---|
98 | { |
---|
99 | if ( deg < source->length() ) |
---|
100 | numelems= (*newhilb)[deg]-(*source)[deg]; |
---|
101 | else |
---|
102 | numelems= (*newhilb)[deg]; |
---|
103 | } |
---|
104 | else |
---|
105 | { |
---|
106 | if (deg < source->length()) |
---|
107 | numelems= -(*source)[deg]; |
---|
108 | else |
---|
109 | { |
---|
110 | deg= 0; |
---|
111 | return 0; |
---|
112 | } |
---|
113 | } |
---|
114 | if (numelems != 0) |
---|
115 | return numelems; |
---|
116 | deg++; |
---|
117 | } |
---|
118 | delete newhilb; |
---|
119 | } |
---|
120 | |
---|
121 | void |
---|
122 | generateMonoms( poly m, int var, int deg, homogData * dat ) |
---|
123 | { |
---|
124 | if ( var == pVariables ) { |
---|
125 | BOOLEAN inSource = FALSE; |
---|
126 | BOOLEAN inDest = FALSE; |
---|
127 | poly mon = pCopy( m ); |
---|
128 | pSetExp( mon, var, deg ); |
---|
129 | pSetm( mon ); |
---|
130 | ++dat->overall; |
---|
131 | int i; |
---|
132 | for ( i= dat->numSourceHeads - 1; (i >= 0) && (inSource==FALSE); i-- ) { |
---|
133 | if ( pDivisibleBy( dat->sourceHeads[i].dm, mon ) ) { |
---|
134 | inSource= TRUE; |
---|
135 | } |
---|
136 | } |
---|
137 | for ( i= dat->numDestPolys - 1; (i >= 0) && (inDest==FALSE); i-- ) { |
---|
138 | if ( pDivisibleBy( (dat->destIdeal->m)[i], mon ) ) { |
---|
139 | inDest= TRUE; |
---|
140 | } |
---|
141 | } |
---|
142 | if ( (!inSource) || (!inDest) ) { |
---|
143 | int basis = 0; |
---|
144 | if ( !inSource ) |
---|
145 | basis= ++(dat->basisSize); |
---|
146 | if ( !inDest ) |
---|
147 | ++dat->numberofdestbasismonoms; |
---|
148 | if ( dat->numMonoms == dat->monlistmax ) { |
---|
149 | dat->monlist= (homogElem * )ReAlloc( dat->monlist, (dat->monlistmax)*sizeof( homogElem ), (dat->monlistmax+dat->monlistblock) * sizeof( homogElem ) ); |
---|
150 | int k; |
---|
151 | for ( k= dat->monlistmax; k < (dat->monlistmax+dat->monlistblock); k++ ) |
---|
152 | dat->monlist[k].homogElem(); |
---|
153 | dat->monlistmax+= dat->monlistblock; |
---|
154 | } |
---|
155 | dat->monlist[dat->numMonoms]= homogElem( mon, basis, inDest ); |
---|
156 | dat->numMonoms++; |
---|
157 | if ( inSource && ! inDest ) PROT( "\\" ); |
---|
158 | if ( ! inSource && inDest ) PROT( "/" ); |
---|
159 | if ( ! inSource && ! inDest ) PROT( "." ); |
---|
160 | } |
---|
161 | else { |
---|
162 | pDelete( & mon ); |
---|
163 | } |
---|
164 | return; |
---|
165 | } |
---|
166 | else { |
---|
167 | poly newm = pCopy( m ); |
---|
168 | while ( deg >= 0 ) { |
---|
169 | generateMonoms( newm, var+1, deg, dat ); |
---|
170 | pIncrExp( newm, var ); |
---|
171 | pSetm( newm ); |
---|
172 | deg--; |
---|
173 | } |
---|
174 | pDelete( & newm ); |
---|
175 | } |
---|
176 | return; |
---|
177 | } |
---|
178 | |
---|
179 | void |
---|
180 | mapMonoms( ring oldRing, homogData & dat ) |
---|
181 | { |
---|
182 | int * vperm = (int *)Alloc( (currRing->N + 1)*sizeof(int) ); |
---|
183 | maFindPerm( oldRing->names, oldRing->N, NULL, 0, currRing->names, currRing->N, NULL, 0, vperm, NULL ); |
---|
184 | nSetMap( oldRing->ch, oldRing->parameter, oldRing->P, oldRing->minpoly ); |
---|
185 | int s; |
---|
186 | for ( s= dat.numMonoms - 1; s >= 0; s-- ) { |
---|
187 | // dat.monlist[s].mon.sm= pPermPoly( dat.monlist[s].mon.dm, vperm, currRing->N, NULL, 0 ); |
---|
188 | // obachman: changed the folowing to reflect the new calling interface of |
---|
189 | // pPermPoly -- Tim please check whether this is correct! |
---|
190 | dat.monlist[s].mon.sm= pPermPoly( dat.monlist[s].mon.dm, vperm, oldRing, NULL, 0 ); |
---|
191 | } |
---|
192 | } |
---|
193 | |
---|
194 | void |
---|
195 | getVectorRep( homogData & dat ) |
---|
196 | { |
---|
197 | // Calculate the NormalForms |
---|
198 | int s; |
---|
199 | for ( s= 0; s < dat.numMonoms; s++ ) { |
---|
200 | if ( dat.monlist[s].inDest == FALSE ) { |
---|
201 | fglmVector v; |
---|
202 | if ( dat.monlist[s].basis == 0 ) { |
---|
203 | v= fglmVector( dat.basisSize ); |
---|
204 | // now the monom is in L(source) |
---|
205 | PROT( "(" ); |
---|
206 | poly nf = kNF( dat.sourceIdeal, NULL, dat.monlist[s].mon.sm ); |
---|
207 | PROT( ")" ); |
---|
208 | poly temp = nf; |
---|
209 | while (temp != NULL ) { |
---|
210 | int t; |
---|
211 | for ( t= dat.numMonoms - 1; t >= 0; t-- ) { |
---|
212 | if ( dat.monlist[t].basis > 0 ) { |
---|
213 | if ( pEqual( dat.monlist[t].mon.sm, temp ) ) { |
---|
214 | number coeff= nCopy( pGetCoeff( temp ) ); |
---|
215 | v.setelem( dat.monlist[t].basis, coeff ); |
---|
216 | } |
---|
217 | } |
---|
218 | } |
---|
219 | temp= pIter(temp); |
---|
220 | } |
---|
221 | pDelete( & nf ); |
---|
222 | } |
---|
223 | else { |
---|
224 | PROT( "." ); |
---|
225 | v= fglmVector( dat.basisSize, dat.monlist[s].basis ); |
---|
226 | } |
---|
227 | dat.monlist[s].v= v; |
---|
228 | } |
---|
229 | } |
---|
230 | } |
---|
231 | |
---|
232 | void |
---|
233 | remapVectors( ring oldring, homogData & dat ) |
---|
234 | { |
---|
235 | nSetMap( oldring->ch, oldring->parameter, oldring->P, oldring->minpoly ); |
---|
236 | int s; |
---|
237 | for ( s= dat.numMonoms - 1; s >= 0; s-- ) { |
---|
238 | if ( dat.monlist[s].inDest == FALSE ) { |
---|
239 | int k; |
---|
240 | fglmVector newv( dat.basisSize ); |
---|
241 | for ( k= dat.basisSize; k > 0; k-- ){ |
---|
242 | number newnum= nMap( dat.monlist[s].v.getelem( k ) ); |
---|
243 | newv.setelem( k, newnum ); |
---|
244 | } |
---|
245 | dat.monlist[s].dv= newv; |
---|
246 | } |
---|
247 | } |
---|
248 | } |
---|
249 | |
---|
250 | void |
---|
251 | gaussreduce( homogData & dat, int maxnum, int BS ) |
---|
252 | { |
---|
253 | int s; |
---|
254 | int found= 0; |
---|
255 | |
---|
256 | int destbasisSize = 0; |
---|
257 | gaussReducer gauss( dat.basisSize ); |
---|
258 | |
---|
259 | for ( s= 0; (s < dat.numMonoms) && (found < maxnum); s++ ) { |
---|
260 | if ( dat.monlist[s].inDest == FALSE ) { |
---|
261 | if ( gauss.reduce( dat.monlist[s].dv ) == FALSE ) { |
---|
262 | destbasisSize++; |
---|
263 | dat.monlist[s].destbasis= destbasisSize; |
---|
264 | gauss.store(); |
---|
265 | PROT( "." ); |
---|
266 | } |
---|
267 | else { |
---|
268 | fglmVector p= gauss.getDependence(); |
---|
269 | poly result = pCopy( dat.monlist[s].mon.dm ); |
---|
270 | pSetCoeff( result, nCopy( p.getconstelem( p.size() ) ) ); |
---|
271 | int l = 0; |
---|
272 | int k; |
---|
273 | for ( k= 1; k < p.size(); k++ ) { |
---|
274 | if ( ! p.elemIsZero( k ) ) { |
---|
275 | while ( dat.monlist[l].destbasis != k ) |
---|
276 | l++; |
---|
277 | poly temp = pCopy( dat.monlist[l].mon.dm ); |
---|
278 | pSetCoeff( temp, nCopy( p.getconstelem( k ) ) ); |
---|
279 | result= pAdd( result, temp ); |
---|
280 | } |
---|
281 | } |
---|
282 | if ( ! nGreaterZero( pGetCoeff( result ) ) ) result= pNeg( result ); |
---|
283 | // PROT2( "(%s)", pString( result ) ); |
---|
284 | PROT( "+" ); |
---|
285 | found++; |
---|
286 | (dat.destIdeal->m)[dat.numDestPolys]= result; |
---|
287 | dat.numDestPolys++; |
---|
288 | if ( IDELEMS(dat.destIdeal) == dat.numDestPolys ) { |
---|
289 | pEnlargeSet( & dat.destIdeal->m, IDELEMS( dat.destIdeal ), BS ); |
---|
290 | IDELEMS( dat.destIdeal )+= BS; |
---|
291 | } |
---|
292 | |
---|
293 | } |
---|
294 | |
---|
295 | } |
---|
296 | } |
---|
297 | PROT2( "(%i", s ); |
---|
298 | PROT2( "/%i)", dat.numberofdestbasismonoms ); |
---|
299 | } |
---|
300 | |
---|
301 | |
---|
302 | BOOLEAN |
---|
303 | fglmhomog( idhdl sourceRingHdl, ideal sourceIdeal, idhdl destRingHdl, ideal & destIdeal ) |
---|
304 | { |
---|
305 | #define groebnerBS 16 |
---|
306 | int numGBelems; |
---|
307 | int deg = 0; |
---|
308 | |
---|
309 | homogData dat; |
---|
310 | |
---|
311 | // get the hilbert series and the leading monomials of the sourceIdeal: |
---|
312 | rSetHdl( sourceRingHdl, TRUE ); |
---|
313 | ring sourceRing = currRing; |
---|
314 | |
---|
315 | intvec * hilb = hHstdSeries( sourceIdeal, NULL, currQuotient ); |
---|
316 | int s; |
---|
317 | dat.sourceIdeal= sourceIdeal; |
---|
318 | dat.sourceHeads= (doublepoly *)Alloc( IDELEMS( sourceIdeal ) * sizeof( doublepoly ) ); |
---|
319 | for ( s= IDELEMS( sourceIdeal ) - 1; s >= 0; s-- ) { |
---|
320 | dat.sourceHeads[s].sm= pHead( (sourceIdeal->m)[s] ); |
---|
321 | } |
---|
322 | dat.numSourceHeads= IDELEMS( sourceIdeal ); |
---|
323 | rSetHdl( destRingHdl, TRUE ); |
---|
324 | ring destRing = currRing; |
---|
325 | |
---|
326 | // Map the sourceHeads to the destRing |
---|
327 | int * vperm = (int *)Alloc( (sourceRing->N + 1)*sizeof(int) ); |
---|
328 | maFindPerm( sourceRing->names, sourceRing->N, NULL, 0, currRing->names, currRing->N, NULL, 0, vperm, NULL ); |
---|
329 | nSetMap( sourceRing->ch, sourceRing->parameter, sourceRing->P, sourceRing->minpoly ); |
---|
330 | for ( s= IDELEMS( sourceIdeal ) - 1; s >= 0; s-- ) { |
---|
331 | dat.sourceHeads[s].dm= pPermPoly( dat.sourceHeads[s].sm, vperm, sourceRing, NULL, 0 ); |
---|
332 | } |
---|
333 | |
---|
334 | dat.destIdeal= idInit( groebnerBS, 1 ); |
---|
335 | dat.numDestPolys= 0; |
---|
336 | |
---|
337 | while ( (numGBelems= hfglmNextdegree( hilb, dat.destIdeal, deg )) != 0 ) { |
---|
338 | int num = 0; // the number of monoms of degree deg |
---|
339 | PROT2( "deg= %i ", deg ); |
---|
340 | PROT2( "num= %i\ngen>", numGBelems ); |
---|
341 | dat.monlistblock= 512; |
---|
342 | dat.monlistmax= dat.monlistblock; |
---|
343 | dat.monlist= (homogElem *)Alloc( dat.monlistmax*sizeof( homogElem ) ); |
---|
344 | int j; |
---|
345 | for ( j= dat.monlistmax - 1; j >= 0; j-- ) dat.monlist[j].homogElem(); |
---|
346 | dat.numMonoms= 0; |
---|
347 | dat.basisSize= 0; |
---|
348 | dat.overall= 0; |
---|
349 | dat.numberofdestbasismonoms= 0; |
---|
350 | |
---|
351 | poly start= pOne(); |
---|
352 | generateMonoms( start, 1, deg, &dat ); |
---|
353 | pDelete( & start ); |
---|
354 | |
---|
355 | PROT2( "(%i/", dat.basisSize ); |
---|
356 | PROT2( "%i)\nvec>", dat.overall ); |
---|
357 | // switch to sourceRing and map monoms |
---|
358 | rSetHdl( sourceRingHdl, TRUE ); |
---|
359 | mapMonoms( destRing, dat ); |
---|
360 | getVectorRep( dat ); |
---|
361 | |
---|
362 | // switch to destination Ring and remap the vectors |
---|
363 | rSetHdl( destRingHdl, TRUE ); |
---|
364 | remapVectors( sourceRing, dat ); |
---|
365 | |
---|
366 | PROT( "<\nred>" ); |
---|
367 | // now do gaussian reduction |
---|
368 | gaussreduce( dat, numGBelems, groebnerBS ); |
---|
369 | |
---|
370 | Free( (ADDRESS)dat.monlist, dat.monlistmax*sizeof( homogElem ) ); |
---|
371 | PROT( "<\n" ); |
---|
372 | } |
---|
373 | PROT( "\n" ); |
---|
374 | destIdeal= dat.destIdeal; |
---|
375 | idSkipZeroes( destIdeal ); |
---|
376 | return TRUE; |
---|
377 | } |
---|
378 | |
---|
379 | ideal |
---|
380 | fglmhomProc(leftv first, leftv second) |
---|
381 | { |
---|
382 | idhdl dest= currRingHdl; |
---|
383 | ideal result; |
---|
384 | // in den durch das erste Argument angegeben Ring schalten: |
---|
385 | rSetHdl( (idhdl)first->data, TRUE ); |
---|
386 | idhdl ih= currRing->idroot->get( second->name, myynest ); |
---|
387 | ASSERT( ih!=NULL, "Error: Can't find ideal in ring"); |
---|
388 | rSetHdl( dest, TRUE ); |
---|
389 | |
---|
390 | ideal i= IDIDEAL(ih); |
---|
391 | fglmhomog( (idhdl)first->data, i, dest, result ); |
---|
392 | |
---|
393 | return( result ); |
---|
394 | } |
---|
395 | |
---|
396 | #endif |
---|
397 | |
---|
398 | // Questions: |
---|
399 | // Muss ich einen intvec freigeben? |
---|
400 | // ---------------------------------------------------------------------------- |
---|
401 | // Local Variables: *** |
---|
402 | // compile-command: "make Singular" *** |
---|
403 | // page-delimiter: "^\\(\\|//!\\)" *** |
---|
404 | // fold-internal-margins: nil *** |
---|
405 | // End: *** |
---|