My Project
Loading...
Searching...
No Matches
fglmhom.cc
Go to the documentation of this file.
1// emacs edit mode for this file is -*- C++ -*-
2
3/****************************************
4* Computer Algebra System SINGULAR *
5****************************************/
6/*
7* ABSTRACT - The FGLM-Algorithm extended for homogeneous ideals.
8* Calculates via the hilbert-function a groebner basis.
9*/
10
11
12
13
14
15#include "kernel/mod2.h"
16#if 0
17#include "factoryconf.h"
18#ifndef NOSTREAMIO
19#include "iostream.h"
20#endif
21// #include "Singular/tok.h"
22#include "kernel/structs.h"
23#include "kernel/subexpr.h"
24#include "kernel/polys.h"
25#include "kernel/ideals.h"
27#include "kernel/ipid.h"
28#include "kernel/ipshell.h"
30#include "fglm.h"
31#include "fglmvec.h"
32#include "fglmgauss.h"
33#include "misc/intvec.h"
37
38// obachman: Got rid of those "redefined" messages by including fglm.h
39#include "fglm.h"
40#if 0
41#define PROT(msg) if (BTEST1(OPT_PROT)) Print(msg)
42#define STICKYPROT(msg) if (BTEST1(OPT_PROT)) Print(msg)
43#define PROT2(msg,arg) if (BTEST1(OPT_PROT)) Print(msg,arg)
44#define STICKYPROT2(msg,arg) if (BTEST1(OPT_PROT)) Print(msg,arg)
45#define fglmASSERT(ignore1,ignore2)
46#endif
47
48struct doublepoly
49{
50 poly sm;
51 poly dm;
52};
53
54class homogElem
55{
56public:
57 doublepoly mon;
59 fglmVector dv;
60 int basis;
61 int destbasis;
62 BOOLEAN inDest;
63 homogElem() : v(), dv(), basis(0), destbasis(0), inDest(FALSE) {}
64 homogElem( poly m, int b, BOOLEAN ind ) :
65 basis(b), inDest(ind)
66 {
67 mon.dm= m;
68 mon.sm= NULL;
69 }
70#ifndef HAVE_EXPLICIT_CONSTR
71 void initialize( poly m, int b, BOOLEAN ind )
72 {
73 basis = b;
74 inDest = ind;
75 mon.dm = m;
76 mon.sm = NULL;
77 }
78 void initialize( const homogElem h )
79 {
80 basis = h.basis;
81 inDest = h.inDest;
82 mon.dm = h.mon.dm;
83 mon.sm = h.mon.sm;
84 }
85#endif
86};
87
88struct homogData
89{
90 ideal sourceIdeal;
91 doublepoly * sourceHeads;
92 int numSourceHeads;
93 ideal destIdeal;
94 int numDestPolys;
95 homogElem * monlist;
96 int monlistmax;
97 int monlistblock;
98 int numMonoms;
99 int basisSize;
100 int overall; // nur zum testen.
101 int numberofdestbasismonoms;
102// homogData() : sourceHeads(NULL), numSourceHeads(0), monlist(NULL),
103// numMonoms(0), basisSize(0) {}
104};
105
106int
107hfglmNextdegree( intvec * source, ideal current, int & deg )
108{
109 int numelems;
110 intvec * newhilb = hFirstSeries( current, NULL, currRing->qideal );
111
112 loop
113 {
114 if ( deg < newhilb->length() )
115 {
116 if ( deg < source->length() )
117 numelems= (*newhilb)[deg]-(*source)[deg];
118 else
119 numelems= (*newhilb)[deg];
120 }
121 else
122 {
123 if (deg < source->length())
124 numelems= -(*source)[deg];
125 else
126 {
127 deg= 0;
128 return 0;
129 }
130 }
131 if (numelems != 0)
132 return numelems;
133 deg++;
134 }
135 delete newhilb;
136}
137
138void
139generateMonoms( poly m, int var, int deg, homogData * dat )
140{
141 if ( var == (currRing->N) ) {
142 BOOLEAN inSource = FALSE;
143 BOOLEAN inDest = FALSE;
144 poly mon = pCopy( m );
145 pSetExp( mon, var, deg );
146 pSetm( mon );
147 ++dat->overall;
148 int i;
149 for ( i= dat->numSourceHeads - 1; (i >= 0) && (inSource==FALSE); i-- ) {
150 if ( pDivisibleBy( dat->sourceHeads[i].dm, mon ) ) {
151 inSource= TRUE;
152 }
153 }
154 for ( i= dat->numDestPolys - 1; (i >= 0) && (inDest==FALSE); i-- ) {
155 if ( pDivisibleBy( (dat->destIdeal->m)[i], mon ) ) {
156 inDest= TRUE;
157 }
158 }
159 if ( (!inSource) || (!inDest) ) {
160 int basis = 0;
161 if ( !inSource )
162 basis= ++(dat->basisSize);
163 if ( !inDest )
164 ++dat->numberofdestbasismonoms;
165 if ( dat->numMonoms == dat->monlistmax ) {
166 int k;
167#ifdef HAVE_EXPLICIT_CONSTR
168 // Expand array using Singulars ReAlloc function
169 dat->monlist=
170 (homogElem * )omReallocSize( dat->monlist,
171 (dat->monlistmax)*sizeof( homogElem ),
172 (dat->monlistmax+dat->monlistblock) * sizeof( homogElem ) );
173 for ( k= dat->monlistmax; k < (dat->monlistmax+dat->monlistblock); k++ )
174 dat->monlist[k].homogElem();
175#else
176 // Expand array by generating new one and copying
177 int newsize = dat->monlistmax + dat->monlistblock;
178 homogElem * tempelem = new homogElem[ newsize ];
179 // Copy old elements
180 for ( k= dat->monlistmax - 1; k >= 0; k-- )
181 tempelem[k].initialize( dat->monlist[k] );
182 delete [] homogElem;
183 homogElem = tempelem;
184#endif
185 dat->monlistmax+= dat->monlistblock;
186 }
187#ifdef HAVE_EXPLICIT_CONSTR
188 dat->monlist[dat->numMonoms]= homogElem( mon, basis, inDest );
189#else
190 dat->monlist[dat->numMonoms].initialize( mon, basis, inDest );
191#endif
192 dat->numMonoms++;
193 if ( inSource && ! inDest ) PROT( "\\" );
194 if ( ! inSource && inDest ) PROT( "/" );
195 if ( ! inSource && ! inDest ) PROT( "." );
196 }
197 else {
198 pDelete( & mon );
199 }
200 return;
201 }
202 else {
203 poly newm = pCopy( m );
204 while ( deg >= 0 ) {
205 generateMonoms( newm, var+1, deg, dat );
206 pIncrExp( newm, var );
207 pSetm( newm );
208 deg--;
209 }
210 pDelete( & newm );
211 }
212 return;
213}
214
215void
216mapMonoms( ring oldRing, homogData & dat )
217{
218 int * vperm = (int *)omAlloc( (currRing->N + 1)*sizeof(int) );
219 maFindPerm( oldRing->names, oldRing->N, NULL, 0, currRing->names, currRing->N, NULL, 0, vperm, NULL );
220 //nSetMap( oldRing->ch, oldRing->parameter, oldRing->P, oldRing->minpoly );
221 nSetMap( oldRing );
222 int s;
223 for ( s= dat.numMonoms - 1; s >= 0; s-- ) {
224// dat.monlist[s].mon.sm= pPermPoly( dat.monlist[s].mon.dm, vperm, currRing->N, NULL, 0 );
225 // obachman: changed the following to reflect the new calling interface of
226 // pPermPoly -- Tim please check whether this is correct!
227 dat.monlist[s].mon.sm= pPermPoly( dat.monlist[s].mon.dm, vperm, oldRing, NULL, 0 );
228 }
229}
230
231void
232getVectorRep( homogData & dat )
233{
234 // Calculate the NormalForms
235 int s;
236 for ( s= 0; s < dat.numMonoms; s++ ) {
237 if ( dat.monlist[s].inDest == FALSE ) {
239 if ( dat.monlist[s].basis == 0 ) {
240 v= fglmVector( dat.basisSize );
241 // now the monom is in L(source)
242 PROT( "(" );
243 poly nf = kNF( dat.sourceIdeal, NULL, dat.monlist[s].mon.sm );
244 PROT( ")" );
245 poly temp = nf;
246 while (temp != NULL ) {
247 int t;
248 for ( t= dat.numMonoms - 1; t >= 0; t-- ) {
249 if ( dat.monlist[t].basis > 0 ) {
250 if ( pLmEqual( dat.monlist[t].mon.sm, temp ) ) {
251 number coeff= nCopy( pGetCoeff( temp ) );
252 v.setelem( dat.monlist[t].basis, coeff );
253 }
254 }
255 }
256 pIter(temp);
257 }
258 pDelete( & nf );
259 }
260 else {
261 PROT( "." );
262 v= fglmVector( dat.basisSize, dat.monlist[s].basis );
263 }
264 dat.monlist[s].v= v;
265 }
266 }
267}
268
269void
270remapVectors( ring oldring, homogData & dat )
271{
272 //nSetMap( oldring->ch, oldring->parameter, oldring->P, oldring->minpoly );
273 nSetMap( oldring );
274 int s;
275 for ( s= dat.numMonoms - 1; s >= 0; s-- ) {
276 if ( dat.monlist[s].inDest == FALSE ) {
277 int k;
278 fglmVector newv( dat.basisSize );
279 for ( k= dat.basisSize; k > 0; k-- ){
280 number newnum= nMap( dat.monlist[s].v.getelem( k ) );
281 newv.setelem( k, newnum );
282 }
283 dat.monlist[s].dv= newv;
284 }
285 }
286}
287
288void
289gaussreduce( homogData & dat, int maxnum, int BS )
290{
291 int s;
292 int found= 0;
293
294 int destbasisSize = 0;
295 gaussReducer gauss( dat.basisSize );
296
297 for ( s= 0; (s < dat.numMonoms) && (found < maxnum); s++ ) {
298 if ( dat.monlist[s].inDest == FALSE ) {
299 if ( gauss.reduce( dat.monlist[s].dv ) == FALSE ) {
300 destbasisSize++;
301 dat.monlist[s].destbasis= destbasisSize;
302 gauss.store();
303 PROT( "." );
304 }
305 else {
306 fglmVector p= gauss.getDependence();
307 poly result = pCopy( dat.monlist[s].mon.dm );
308 pSetCoeff( result, nCopy( p.getconstelem( p.size() ) ) );
309 int l = 0;
310 int k;
311 for ( k= 1; k < p.size(); k++ ) {
312 if ( ! p.elemIsZero( k ) ) {
313 while ( dat.monlist[l].destbasis != k )
314 l++;
315 poly temp = pCopy( dat.monlist[l].mon.dm );
316 pSetCoeff( temp, nCopy( p.getconstelem( k ) ) );
317 result= pAdd( result, temp );
318 }
319 }
320 if ( ! nGreaterZero( pGetCoeff( result ) ) ) result= pNeg( result );
321// PROT2( "(%s)", pString( result ) );
322 PROT( "+" );
323 found++;
324 (dat.destIdeal->m)[dat.numDestPolys]= result;
325 dat.numDestPolys++;
326 if ( IDELEMS(dat.destIdeal) == dat.numDestPolys ) {
327 pEnlargeSet( & dat.destIdeal->m, IDELEMS( dat.destIdeal ), BS );
328 IDELEMS( dat.destIdeal )+= BS;
329 }
330
331 }
332
333 }
334 }
335 PROT2( "(%i", s );
336 PROT2( "/%i)", dat.numberofdestbasismonoms );
337}
338
339
341fglmhomog( ring sourceRing, ideal sourceIdeal, ring destRing, ideal & destIdeal )
342{
343#define groebnerBS 16
344 int numGBelems;
345 int deg = 0;
346
347 homogData dat;
348
349 // get the hilbert series and the leading monomials of the sourceIdeal:
350 rChangeCurrRing( sourceRing );
351
352 intvec * hilb = hFirstSeries( sourceIdeal, NULL, currRing->qideal );
353 int s;
354 dat.sourceIdeal= sourceIdeal;
355 dat.sourceHeads= (doublepoly *)omAlloc( IDELEMS( sourceIdeal ) * sizeof( doublepoly ) );
356 for ( s= IDELEMS( sourceIdeal ) - 1; s >= 0; s-- )
357 {
358 dat.sourceHeads[s].sm= pHead( (sourceIdeal->m)[s] );
359 }
360 dat.numSourceHeads= IDELEMS( sourceIdeal );
361 rChangeCurrRing( destRing );
362
363 // Map the sourceHeads to the destRing
364 int * vperm = (int *)omAlloc( (sourceRing->N + 1)*sizeof(int) );
365 maFindPerm( sourceRing->names, sourceRing->N, NULL, 0, currRing->names,
366 currRing->N, NULL, 0, vperm, NULL, currRing->ch);
367 //nSetMap( sourceRing->ch, sourceRing->parameter, sourceRing->P, sourceRing->minpoly );
368 nSetMap( sourceRing );
369 for ( s= IDELEMS( sourceIdeal ) - 1; s >= 0; s-- )
370 {
371 dat.sourceHeads[s].dm= pPermPoly( dat.sourceHeads[s].sm, vperm, sourceRing, NULL, 0 );
372 }
373
374 dat.destIdeal= idInit( groebnerBS, 1 );
375 dat.numDestPolys= 0;
376
377 while ( (numGBelems= hfglmNextdegree( hilb, dat.destIdeal, deg )) != 0 )
378 {
379 int num = 0; // the number of monoms of degree deg
380 PROT2( "deg= %i ", deg );
381 PROT2( "num= %i\ngen>", numGBelems );
382 dat.monlistblock= 512;
383 dat.monlistmax= dat.monlistblock;
384#ifdef HAVE_EXPLICIT_CONSTR
385 dat.monlist= (homogElem *)omAlloc( dat.monlistmax*sizeof( homogElem ) );
386 int j;
387 for ( j= dat.monlistmax - 1; j >= 0; j-- ) dat.monlist[j].homogElem();
388#else
389 dat.monlist = new homogElem[ dat.monlistmax ];
390#endif
391 dat.numMonoms= 0;
392 dat.basisSize= 0;
393 dat.overall= 0;
394 dat.numberofdestbasismonoms= 0;
395
396 poly start= pOne();
397 generateMonoms( start, 1, deg, &dat );
398 pDelete( & start );
399
400 PROT2( "(%i/", dat.basisSize );
401 PROT2( "%i)\nvec>", dat.overall );
402 // switch to sourceRing and map monoms
403 rChangeCurrRing( sourceRing );
404 mapMonoms( destRing, dat );
405 getVectorRep( dat );
406
407 // switch to destination Ring and remap the vectors
408 rChangeCurrRing( destRing );
409 remapVectors( sourceRing, dat );
410
411 PROT( "<\nred>" );
412 // now do gaussian reduction
413 gaussreduce( dat, numGBelems, groebnerBS );
414
415#ifdef HAVE_EXPLICIT_CONSTR
416 omFreeSize( (ADDRESS)dat.monlist, dat.monlistmax*sizeof( homogElem ) );
417#else
418 delete [] dat.monlist;
419#endif
420 PROT( "<\n" );
421 }
422 PROT( "\n" );
423 destIdeal= dat.destIdeal;
424 idSkipZeroes( destIdeal );
425 return TRUE;
426}
427
428/* ideal fglmhomProc(leftv first, leftv second) // TODO: Move to Singular/
429{
430 idhdl dest= currRingHdl;
431 ideal result;
432 // in den durch das erste Argument angegeben Ring schalten:
433 rSetHdl( (idhdl)first->data, TRUE );
434 idhdl ih= currRing->idroot->get( second->name, myynest );
435 ASSERT( ih!=NULL, "Error: Can't find ideal in ring");
436 rSetHdl( dest, TRUE );
437
438 ideal i= IDIDEAL(ih);
439 fglmhomog( IDRING((idhdl)first->data), i, IDRING(dest), result );
440
441 return( result );
442} */
443
444#endif
445
446// Questions:
447// Muss ich einen intvec freigeben?
448// ----------------------------------------------------------------------------
449// Local Variables: ***
450// compile-command: "make Singular" ***
451// page-delimiter: "^\\‍(␌\\|//!\\‍)" ***
452// fold-internal-margins: nil ***
453// End: ***
int BOOLEAN
Definition: auxiliary.h:87
#define TRUE
Definition: auxiliary.h:100
#define FALSE
Definition: auxiliary.h:96
void * ADDRESS
Definition: auxiliary.h:119
CanonicalForm num(const CanonicalForm &f)
int l
Definition: cfEzgcd.cc:100
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
int p
Definition: cfModGcd.cc:4078
CanonicalForm b
Definition: cfModGcd.cc:4103
Definition: intvec.h:23
return result
Definition: facAbsBiFact.cc:75
const CanonicalForm int s
Definition: facAbsFact.cc:51
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
bool found
Definition: facFactorize.cc:55
int j
Definition: facHensel.cc:110
if(!FE_OPT_NO_SHELL_FLAG)(void) system(sys)
intvec * hFirstSeries(ideal A, intvec *module_w, ideal Q, intvec *wdegree)
Definition: hilb.cc:2036
static BOOLEAN length(leftv result, leftv arg)
Definition: interval.cc:257
STATIC_VAR Poly * h
Definition: janet.cc:971
#define PROT(msg)
Definition: fglm.h:19
#define PROT2(msg, arg)
Definition: fglm.h:21
poly kNF(ideal F, ideal Q, poly p, int syzComp, int lazyReduce)
Definition: kstd1.cc:3182
void maFindPerm(char const *const *const preim_names, int preim_n, char const *const *const preim_par, int preim_p, char const *const *const names, int n, char const *const *const par, int nop, int *perm, int *par_perm, n_coeffType ch)
Definition: maps.cc:163
#define pIter(p)
Definition: monomials.h:37
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:44
void initialize()
#define nSetMap(R)
Definition: numbers.h:43
#define nCopy(n)
Definition: numbers.h:15
#define nGreaterZero(n)
Definition: numbers.h:27
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
#define omAlloc(size)
Definition: omAllocDecl.h:210
#define omReallocSize(addr, o_size, size)
Definition: omAllocDecl.h:220
#define NULL
Definition: omList.c:12
void pEnlargeSet(poly **p, int l, int increment)
Definition: p_polys.cc:3696
void rChangeCurrRing(ring r)
Definition: polys.cc:15
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
Compatibility layer for legacy polynomial operations (over currRing)
#define pAdd(p, q)
Definition: polys.h:203
#define pDelete(p_ptr)
Definition: polys.h:186
#define pHead(p)
returns newly allocated copy of Lm(p), coef is copied, next=NULL, p might be NULL
Definition: polys.h:67
#define pSetm(p)
Definition: polys.h:271
#define pNeg(p)
Definition: polys.h:198
#define pLmEqual(p1, p2)
Definition: polys.h:111
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition: polys.h:31
#define pIncrExp(p, i)
Definition: polys.h:43
#define pDivisibleBy(a, b)
returns TRUE, if leading monom of a divides leading monom of b i.e., if there exists a expvector c > ...
Definition: polys.h:138
#define pSetExp(p, i, v)
Definition: polys.h:42
#define pCopy(p)
return a copy of the poly
Definition: polys.h:185
#define pOne()
Definition: polys.h:315
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
void idSkipZeroes(ideal ide)
gives an ideal/module the minimal possible size
#define IDELEMS(i)
Definition: simpleideals.h:23
#define loop
Definition: structs.h:75
Definition: gnumpfl.cc:25