17 #include <factoryconf.h>
23 #include <kernel/subexpr.h>
27 #include <kernel/ipid.h>
28 #include <kernel/ipshell.h>
42 #define PROT(msg) if (BTEST1(OPT_PROT)) Print(msg)
43 #define STICKYPROT(msg) if (BTEST1(OPT_PROT)) Print(msg)
44 #define PROT2(msg,arg) if (BTEST1(OPT_PROT)) Print(msg,arg)
45 #define STICKYPROT2(msg,arg) if (BTEST1(OPT_PROT)) Print(msg,arg)
46 #define fglmASSERT(ignore1,ignore2)
64 homogElem() : v(), dv(), basis(0), destbasis(0), inDest(
FALSE) {}
71 #ifndef HAVE_EXPLICIT_CONSTR
92 doublepoly * sourceHeads;
102 int numberofdestbasismonoms;
108 hfglmNextdegree(
intvec * source, ideal current,
int & deg )
115 if ( deg < newhilb->length() )
117 if ( deg < source->length() )
118 numelems= (*newhilb)[deg]-(*source)[deg];
120 numelems= (*newhilb)[deg];
124 if (deg < source->length())
125 numelems= -(*source)[deg];
140 generateMonoms(
poly m,
int var,
int deg, homogData * dat )
150 for ( i= dat->numSourceHeads - 1; (i >= 0) && (inSource==
FALSE); i-- ) {
155 for ( i= dat->numDestPolys - 1; (i >= 0) && (inDest==
FALSE); i-- ) {
160 if ( (!inSource) || (!inDest) ) {
163 basis= ++(dat->basisSize);
165 ++dat->numberofdestbasismonoms;
166 if ( dat->numMonoms == dat->monlistmax ) {
168 #ifdef HAVE_EXPLICIT_CONSTR
172 (dat->monlistmax)*
sizeof( homogElem ),
173 (dat->monlistmax+dat->monlistblock) *
sizeof( homogElem ) );
174 for ( k= dat->monlistmax; k < (dat->monlistmax+dat->monlistblock); k++ )
175 dat->monlist[k].homogElem();
178 int newsize = dat->monlistmax + dat->monlistblock;
179 homogElem * tempelem =
new homogElem[ newsize ];
181 for ( k= dat->monlistmax - 1; k >= 0; k-- )
182 tempelem[
k].initialize( dat->monlist[k] );
184 homogElem = tempelem;
186 dat->monlistmax+= dat->monlistblock;
188 #ifdef HAVE_EXPLICIT_CONSTR
189 dat->monlist[dat->numMonoms]= homogElem( mon, basis, inDest );
191 dat->monlist[dat->numMonoms].initialize( mon, basis, inDest );
194 if ( inSource && ! inDest )
PROT(
"\\" );
195 if ( ! inSource && inDest )
PROT(
"/" );
196 if ( ! inSource && ! inDest )
PROT(
"." );
206 generateMonoms( newm, var+1, deg, dat );
217 mapMonoms( ring oldRing, homogData & dat )
224 for ( s= dat.numMonoms - 1; s >= 0; s-- ) {
228 dat.monlist[
s].mon.sm= pPermPoly( dat.monlist[s].mon.dm, vperm, oldRing,
NULL, 0 );
233 getVectorRep( homogData & dat )
237 for ( s= 0; s < dat.numMonoms; s++ ) {
238 if ( dat.monlist[s].inDest ==
FALSE ) {
240 if ( dat.monlist[s].basis == 0 ) {
244 poly nf =
kNF( dat.sourceIdeal,
NULL, dat.monlist[s].mon.sm );
247 while (temp !=
NULL ) {
249 for ( t= dat.numMonoms - 1; t >= 0; t-- ) {
250 if ( dat.monlist[t].basis > 0 ) {
251 if (
pLmEqual( dat.monlist[t].mon.sm, temp ) ) {
253 v.
setelem( dat.monlist[t].basis, coeff );
263 v=
fglmVector( dat.basisSize, dat.monlist[s].basis );
271 remapVectors( ring oldring, homogData & dat )
276 for ( s= dat.numMonoms - 1; s >= 0; s-- ) {
277 if ( dat.monlist[s].inDest ==
FALSE ) {
280 for ( k= dat.basisSize; k > 0; k-- ){
281 number newnum= nMap( dat.monlist[s].v.getelem( k ) );
282 newv.setelem( k, newnum );
284 dat.monlist[
s].dv= newv;
290 gaussreduce( homogData & dat,
int maxnum,
int BS )
295 int destbasisSize = 0;
298 for ( s= 0; (s < dat.numMonoms) && (found < maxnum); s++ ) {
299 if ( dat.monlist[s].inDest ==
FALSE ) {
300 if ( gauss.reduce( dat.monlist[s].dv ) ==
FALSE ) {
302 dat.monlist[
s].destbasis= destbasisSize;
312 for ( k= 1; k < p.
size(); k++ ) {
314 while ( dat.monlist[l].destbasis != k )
316 poly temp =
pCopy( dat.monlist[l].mon.dm );
318 result=
pAdd( result, temp );
325 (dat.destIdeal->m)[dat.numDestPolys]= result;
327 if (
IDELEMS(dat.destIdeal) == dat.numDestPolys ) {
337 PROT2(
"/%i)", dat.numberofdestbasismonoms );
342 fglmhomog( ring sourceRing, ideal sourceIdeal, ring destRing, ideal & destIdeal )
344 #define groebnerBS 16
355 dat.sourceIdeal= sourceIdeal;
356 dat.sourceHeads= (doublepoly *)
omAlloc(
IDELEMS( sourceIdeal ) *
sizeof( doublepoly ) );
357 for ( s=
IDELEMS( sourceIdeal ) - 1; s >= 0; s-- )
359 dat.sourceHeads[
s].sm=
pHead( (sourceIdeal->m)[s] );
361 dat.numSourceHeads=
IDELEMS( sourceIdeal );
365 int * vperm = (
int *)
omAlloc( (sourceRing->N + 1)*
sizeof(int) );
370 for ( s=
IDELEMS( sourceIdeal ) - 1; s >= 0; s-- )
372 dat.sourceHeads[
s].dm= pPermPoly( dat.sourceHeads[s].sm, vperm, sourceRing,
NULL, 0 );
375 dat.destIdeal=
idInit( groebnerBS, 1 );
378 while ( (numGBelems= hfglmNextdegree( hilb, dat.destIdeal, deg )) != 0 )
381 PROT2(
"deg= %i ", deg );
382 PROT2(
"num= %i\ngen>", numGBelems );
383 dat.monlistblock= 512;
384 dat.monlistmax= dat.monlistblock;
385 #ifdef HAVE_EXPLICIT_CONSTR
386 dat.monlist= (homogElem *)
omAlloc( dat.monlistmax*
sizeof( homogElem ) );
388 for ( j= dat.monlistmax - 1; j >= 0; j-- ) dat.monlist[
j].homogElem();
390 dat.monlist =
new homogElem[ dat.monlistmax ];
395 dat.numberofdestbasismonoms= 0;
398 generateMonoms( start, 1, deg, &dat );
401 PROT2(
"(%i/", dat.basisSize );
402 PROT2(
"%i)\nvec>", dat.overall );
405 mapMonoms( destRing, dat );
410 remapVectors( sourceRing, dat );
414 gaussreduce( dat, numGBelems, groebnerBS );
416 #ifdef HAVE_EXPLICIT_CONSTR
419 delete [] dat.monlist;
424 destIdeal= dat.destIdeal;
const CanonicalForm int s
poly kNF(ideal F, ideal Q, poly p, int syzComp, int lazyReduce)
Compatiblity layer for legacy polynomial operations (over currRing)
#define omFreeSize(addr, size)
number getconstelem(int i) const
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy ...
intvec * hHstdSeries(ideal S, intvec *modulweight, intvec *wdegree, ideal Q, ring tailRing)
int int kStrategy strat if(h==NULL) return NULL
#define omReallocSize(addr, o_size, size)
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
#define pHead(p)
returns newly allocated copy of Lm(p), coef is copied, next=NULL, p might be NULL ...
void idSkipZeroes(ideal ide)
gives an ideal/module the minimal possible size
void rChangeCurrRing(ring r)
ideal idInit(int idsize, int rank)
initialise an ideal / module
const Variable & v
< [in] a sqrfree bivariate poly
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)
#define pDivisibleBy(a, b)
returns TRUE, if leading monom of a divides leading monom of b i.e., if there exists a expvector c > ...
void pEnlargeSet(poly **p, int l, int increment)
void setelem(int i, number &n)
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
#define pCopy(p)
return a copy of the poly