/**************************************** * Computer Algebra System SINGULAR * ****************************************/ /*************************************************************** * File: p_polys.cc * Purpose: implementation of ring independent poly procedures? * Author: obachman (Olaf Bachmann) * Created: 8/00 *******************************************************************/ #include "config.h" #include #include #include #include #include #include // ??? #include #include #define TRANSEXT_PRIVATES #include #include #include #include #include "ring.h" #include "p_polys.h" #include #include #include // #include // #include #ifndef NDEBUG // #include #endif #ifdef HAVE_PLURAL #include "nc/nc.h" #include "nc/sca.h" #endif #include "coeffrings.h" #ifdef HAVE_FACTORY #include "clapsing.h" #endif /* * lift ideal with coeffs over Z (mod N) to Q via Farey */ poly p_Farey(poly p, number N, const ring r) { poly h=p_Copy(p,r); poly hh=h; while(h!=NULL) { number c=pGetCoeff(h); pSetCoeff0(h,n_Farey(c,N,r->cf)); n_Delete(&c,r->cf); pIter(h); } while((hh!=NULL)&&(n_IsZero(pGetCoeff(hh),r->cf))) { p_LmDelete(&hh,r); } h=hh; while((h!=NULL) && (pNext(h)!=NULL)) { if(n_IsZero(pGetCoeff(pNext(h)),r->cf)) { p_LmDelete(&pNext(h),r); } else pIter(h); } return hh; } /*2 * xx,q: arrays of length 0..rl-1 * xx[i]: SB mod q[i] * assume: char=0 * assume: q[i]!=0 * destroys xx */ poly p_ChineseRemainder(poly *xx, number *x,number *q, int rl, const ring R) { poly r,h,hh; int j; poly res_p=NULL; loop { /* search the lead term */ r=NULL; for(j=rl-1;j>=0;j--) { h=xx[j]; if ((h!=NULL) &&((r==NULL)||(p_LmCmp(r,h,R)==-1))) r=h; } /* nothing found -> return */ if (r==NULL) break; /* create the monomial in h */ h=p_Head(r,R); /* collect the coeffs in x[..]*/ for(j=rl-1;j>=0;j--) { hh=xx[j]; if ((hh!=NULL) && (p_LmCmp(r,hh,R)==0)) { x[j]=pGetCoeff(hh); hh=p_LmFreeAndNext(hh,R); xx[j]=hh; } else x[j]=n_Init(0, R); } number n=n_ChineseRemainderSym(x,q,rl,TRUE,R->cf); for(j=rl-1;j>=0;j--) { x[j]=NULL; // nlInit(0...) takes no memory } if (n_IsZero(n,R)) p_Delete(&h,R); else { //Print("new mon:");pWrite(h); p_SetCoeff(h,n,R); pNext(h)=res_p; res_p=h; // building res_p in reverse order! } } res_p=pReverse(res_p); p_Test(res_p, R); return res_p; } /*************************************************************** * * Completing what needs to be set for the monomial * ***************************************************************/ // this is special for the syz stuff static int* _components = NULL; static long* _componentsShifted = NULL; static int _componentsExternal = 0; BOOLEAN pSetm_error=0; #ifndef NDEBUG # define MYTEST 0 #else /* ifndef NDEBUG */ # define MYTEST 0 #endif /* ifndef NDEBUG */ void p_Setm_General(poly p, const ring r) { p_LmCheckPolyRing(p, r); int pos=0; if (r->typ!=NULL) { loop { long ord=0; sro_ord* o=&(r->typ[pos]); switch(o->ord_typ) { case ro_dp: { int a,e; a=o->data.dp.start; e=o->data.dp.end; for(int i=a;i<=e;i++) ord+=p_GetExp(p,i,r); p->exp[o->data.dp.place]=ord; break; } case ro_wp_neg: ord=POLY_NEGWEIGHT_OFFSET; // no break; case ro_wp: { int a,e; a=o->data.wp.start; e=o->data.wp.end; int *w=o->data.wp.weights; #if 1 for(int i=a;i<=e;i++) ord+=p_GetExp(p,i,r)*w[i-a]; #else long ai; int ei,wi; for(int i=a;i<=e;i++) { ei=p_GetExp(p,i,r); wi=w[i-a]; ai=ei*wi; if (ai/ei!=wi) pSetm_error=TRUE; ord+=ai; if (ordexp[o->data.wp.place]=ord; break; } case ro_am: { ord = POLY_NEGWEIGHT_OFFSET; const short a=o->data.am.start; const short e=o->data.am.end; const int * w=o->data.am.weights; #if 1 for(short i=a; i<=e; i++, w++) ord += ((*w) * p_GetExp(p,i,r)); #else long ai; int ei,wi; for(short i=a;i<=e;i++) { ei=p_GetExp(p,i,r); wi=w[i-a]; ai=ei*wi; if (ai/ei!=wi) pSetm_error=TRUE; ord += ai; if (orddata.am.len_gen; if ((c > 0) && (c <= len_gen)) { assume( w == o->data.am.weights_m ); assume( w[0] == len_gen ); ord += w[c]; } p->exp[o->data.am.place] = ord; break; } case ro_wp64: { int64 ord=0; int a,e; a=o->data.wp64.start; e=o->data.wp64.end; int64 *w=o->data.wp64.weights64; int64 ei,wi,ai; for(int i=a;i<=e;i++) { //Print("exp %d w %d \n",p_GetExp(p,i,r),(int)w[i-a]); //ord+=((int64)p_GetExp(p,i,r))*w[i-a]; ei=(int64)p_GetExp(p,i,r); wi=w[i-a]; ai=ei*wi; if(ei!=0 && ai/ei!=wi) { pSetm_error=TRUE; #if SIZEOF_LONG == 4 Print("ai %lld, wi %lld\n",ai,wi); #else Print("ai %ld, wi %ld\n",ai,wi); #endif } ord+=ai; if (ord>31 ); /*(ord/(mask+1));*/ //Print("mask: %x, ord: %d, a_0: %d, a_1: %d\n" //,(int)mask,(int)ord,(int)a_0,(int)a_1); //Print("mask: %d",mask); p->exp[o->data.wp64.place]=a_1; p->exp[o->data.wp64.place+1]=a_0; // if(p_Setm_error) Print("***************************\n // ***************************\n // **WARNING: overflow error**\n // ***************************\n // ***************************\n"); break; } case ro_cp: { int a,e; a=o->data.cp.start; e=o->data.cp.end; int pl=o->data.cp.place; for(int i=a;i<=e;i++) { p->exp[pl]=p_GetExp(p,i,r); pl++; } break; } case ro_syzcomp: { int c=p_GetComp(p,r); long sc = c; int* Components = (_componentsExternal ? _components : o->data.syzcomp.Components); long* ShiftedComponents = (_componentsExternal ? _componentsShifted: o->data.syzcomp.ShiftedComponents); if (ShiftedComponents != NULL) { assume(Components != NULL); assume(c == 0 || Components[c] != 0); sc = ShiftedComponents[Components[c]]; assume(c == 0 || sc != 0); } p->exp[o->data.syzcomp.place]=sc; break; } case ro_syz: { const unsigned long c = p_GetComp(p, r); const short place = o->data.syz.place; const int limit = o->data.syz.limit; if (c > (unsigned long)limit) p->exp[place] = o->data.syz.curr_index; else if (c > 0) { assume( (1 <= c) && (c <= (unsigned long)limit) ); p->exp[place]= o->data.syz.syz_index[c]; } else { assume(c == 0); p->exp[place]= 0; } break; } // Prefix for Induced Schreyer ordering case ro_isTemp: // Do nothing?? (to be removed into suffix later on...?) { assume(p != NULL); #ifndef NDEBUG #if MYTEST Print("p_Setm_General: ro_isTemp ord: pos: %d, p: ", pos); p_DebugPrint(p, r, r, 1); #endif #endif int c = p_GetComp(p, r); assume( c >= 0 ); // Let's simulate case ro_syz above.... // Should accumulate (by Suffix) and be a level indicator const int* const pVarOffset = o->data.isTemp.pVarOffset; assume( pVarOffset != NULL ); // TODO: Can this be done in the suffix??? for( int i = 1; i <= r->N; i++ ) // No v[0] here!!! { const int vo = pVarOffset[i]; if( vo != -1) // TODO: optimize: can be done once! { // Hans! Please don't break it again! p_SetExp(p, ..., r, vo) is correct: p_SetExp(p, p_GetExp(p, i, r), r, vo); // copy put them verbatim // Hans! Please don't break it again! p_GetExp(p, r, vo) is correct: assume( p_GetExp(p, r, vo) == p_GetExp(p, i, r) ); // copy put them verbatim } } #ifndef NDEBUG for( int i = 1; i <= r->N; i++ ) // No v[0] here!!! { const int vo = pVarOffset[i]; if( vo != -1) // TODO: optimize: can be done once! { // Hans! Please don't break it again! p_GetExp(p, r, vo) is correct: assume( p_GetExp(p, r, vo) == p_GetExp(p, i, r) ); // copy put them verbatim } } #if MYTEST // if( p->exp[o->data.isTemp.start] > 0 ) PrintS("after Values: "); p_DebugPrint(p, r, r, 1); #endif #endif break; } // Suffix for Induced Schreyer ordering case ro_is: { #ifndef NDEBUG #if MYTEST Print("p_Setm_General: ro_is ord: pos: %d, p: ", pos); p_DebugPrint(p, r, r, 1); #endif #endif assume(p != NULL); int c = p_GetComp(p, r); assume( c >= 0 ); const ideal F = o->data.is.F; const int limit = o->data.is.limit; assume( limit >= 0 ); const int start = o->data.is.start; if( F != NULL && c > limit ) { #ifndef NDEBUG #if MYTEST Print("p_Setm_General: ro_is : in rSetm: pos: %d, c: %d > limit: %d\n", c, pos, limit); // p_DebugPrint(p, r, r, 1); PrintS("preComputed Values: "); p_DebugPrint(p, r, r, 1); #endif #endif // if( c > limit ) // BUG??? p->exp[start] = 1; // else // p->exp[start] = 0; c -= limit; assume( c > 0 ); c--; if( c >= IDELEMS(F) ) break; assume( c < IDELEMS(F) ); // What about others??? const poly pp = F->m[c]; // get reference monomial!!! if(pp == NULL) break; assume(pp != NULL); #ifndef NDEBUG #if MYTEST Print("Respective F[c - %d: %d] pp: ", limit, c); p_DebugPrint(pp, r, r, 1); #endif #endif const int end = o->data.is.end; assume(start <= end); // const int st = o->data.isTemp.start; #ifndef NDEBUG Print("p_Setm_General: is(-Temp-) :: c: %d, limit: %d, [st:%d] ===>>> %ld\n", c, limit, start, p->exp[start]); #endif // p_ExpVectorAdd(p, pp, r); for( int i = start; i <= end; i++) // v[0] may be here... p->exp[i] += pp->exp[i]; // !!!!!!!! ADD corresponding LT(F) // p_MemAddAdjust(p, ri); if (r->NegWeightL_Offset != NULL) { for (int i=r->NegWeightL_Size-1; i>=0; i--) { const int _i = r->NegWeightL_Offset[i]; if( start <= _i && _i <= end ) p->exp[_i] -= POLY_NEGWEIGHT_OFFSET; } } #ifndef NDEBUG const int* const pVarOffset = o->data.is.pVarOffset; assume( pVarOffset != NULL ); for( int i = 1; i <= r->N; i++ ) // No v[0] here!!! { const int vo = pVarOffset[i]; if( vo != -1) // TODO: optimize: can be done once! // Hans! Please don't break it again! p_GetExp(p/pp, r, vo) is correct: assume( p_GetExp(p, r, vo) == (p_GetExp(p, i, r) + p_GetExp(pp, r, vo)) ); } // TODO: how to check this for computed values??? #if MYTEST PrintS("Computed Values: "); p_DebugPrint(p, r, r, 1); #endif #endif } else { p->exp[start] = 0; //!!!!????? where????? const int* const pVarOffset = o->data.is.pVarOffset; // What about v[0] - component: it will be added later by // suffix!!! // TODO: Test it! const int vo = pVarOffset[0]; if( vo != -1 ) p->exp[vo] = c; // initial component v[0]! #ifndef NDEBUG #if MYTEST Print("ELSE p_Setm_General: ro_is :: c: %d <= limit: %d, vo: %d, exp: %d\n", c, limit, vo, p->exp[vo]); p_DebugPrint(p, r, r, 1); #endif #endif } break; } default: dReportError("wrong ord in rSetm:%d\n",o->ord_typ); return; } pos++; if (pos == r->OrdSize) return; } } } void p_Setm_Syz(poly p, ring r, int* Components, long* ShiftedComponents) { _components = Components; _componentsShifted = ShiftedComponents; _componentsExternal = 1; p_Setm_General(p, r); _componentsExternal = 0; } // dummy for lp, ls, etc void p_Setm_Dummy(poly p, const ring r) { p_LmCheckPolyRing(p, r); } // for dp, Dp, ds, etc void p_Setm_TotalDegree(poly p, const ring r) { p_LmCheckPolyRing(p, r); p->exp[r->pOrdIndex] = p_Totaldegree(p, r); } // for wp, Wp, ws, etc void p_Setm_WFirstTotalDegree(poly p, const ring r) { p_LmCheckPolyRing(p, r); p->exp[r->pOrdIndex] = p_WFirstTotalDegree(p, r); } p_SetmProc p_GetSetmProc(ring r) { // covers lp, rp, ls, if (r->typ == NULL) return p_Setm_Dummy; if (r->OrdSize == 1) { if (r->typ[0].ord_typ == ro_dp && r->typ[0].data.dp.start == 1 && r->typ[0].data.dp.end == r->N && r->typ[0].data.dp.place == r->pOrdIndex) return p_Setm_TotalDegree; if (r->typ[0].ord_typ == ro_wp && r->typ[0].data.wp.start == 1 && r->typ[0].data.wp.end == r->N && r->typ[0].data.wp.place == r->pOrdIndex && r->typ[0].data.wp.weights == r->firstwv) return p_Setm_WFirstTotalDegree; } return p_Setm_General; } /* -------------------------------------------------------------------*/ /* several possibilities for pFDeg: the degree of the head term */ /* comptible with ordering */ long p_Deg(poly a, const ring r) { p_LmCheckPolyRing(a, r); // assume(p_GetOrder(a, r) == p_WTotaldegree(a, r)); // WRONG assume! return p_GetOrder(a, r); } // p_WTotalDegree for weighted orderings // whose first block covers all variables long p_WFirstTotalDegree(poly p, const ring r) { int i; long sum = 0; for (i=1; i<= r->firstBlockEnds; i++) { sum += p_GetExp(p, i, r)*r->firstwv[i-1]; } return sum; } /*2 * compute the degree of the leading monomial of p * with respect to weigths from the ordering * the ordering is not compatible with degree so do not use p->Order */ long p_WTotaldegree(poly p, const ring r) { p_LmCheckPolyRing(p, r); int i, k; long j =0; // iterate through each block: for (i=0;r->order[i]!=0;i++) { int b0=r->block0[i]; int b1=r->block1[i]; switch(r->order[i]) { case ringorder_M: for (k=b0 /*r->block0[i]*/;k<=b1 /*r->block1[i]*/;k++) { // in jedem block: j+= p_GetExp(p,k,r)*r->wvhdl[i][k - b0 /*r->block0[i]*/]*r->OrdSgn; } break; case ringorder_wp: case ringorder_ws: case ringorder_Wp: case ringorder_Ws: for (k=b0 /*r->block0[i]*/;k<=b1 /*r->block1[i]*/;k++) { // in jedem block: j+= p_GetExp(p,k,r)*r->wvhdl[i][k - b0 /*r->block0[i]*/]; } break; case ringorder_lp: case ringorder_ls: case ringorder_rs: case ringorder_dp: case ringorder_ds: case ringorder_Dp: case ringorder_Ds: case ringorder_rp: for (k=b0 /*r->block0[i]*/;k<=b1 /*r->block1[i]*/;k++) { j+= p_GetExp(p,k,r); } break; case ringorder_a64: { int64* w=(int64*)r->wvhdl[i]; for (k=0;k<=(b1 /*r->block1[i]*/ - b0 /*r->block0[i]*/);k++) { //there should be added a line which checks if w[k]>2^31 j+= p_GetExp(p,k+1, r)*(long)w[k]; } //break; return j; } case ringorder_c: case ringorder_C: case ringorder_S: case ringorder_s: case ringorder_aa: case ringorder_IS: break; case ringorder_a: for (k=b0 /*r->block0[i]*/;k<=b1 /*r->block1[i]*/;k++) { // only one line j+= p_GetExp(p,k, r)*r->wvhdl[i][ k- b0 /*r->block0[i]*/]; } //break; return j; #ifndef NDEBUG default: Print("missing order %d in p_WTotaldegree\n",r->order[i]); break; #endif } } return j; } long p_DegW(poly p, const short *w, const ring R) { assume( p_Test(p, R) ); assume( w != NULL ); if( p == NULL ) // TODO: ??? return -1; long r = totaldegreeWecart_IV(p, R, w); pIter(p); for ( ; p!=NULL; pIter(p) ) { const long t = totaldegreeWecart_IV(p, R, w); if (t > r) r = t; } return r; } int p_Weight(int i, const ring r) { if ((r->firstwv==NULL) || (i>r->firstBlockEnds)) { return 1; } return r->firstwv[i-1]; } long p_WDegree(poly p, const ring r) { if (r->firstwv==NULL) return p_Totaldegree(p, r); p_LmCheckPolyRing(p, r); int i; long j =0; for(i=1;i<=r->firstBlockEnds;i++) j+=p_GetExp(p, i, r)*r->firstwv[i-1]; for (;i<=r->N;i++) j+=p_GetExp(p,i, r)*p_Weight(i, r); return j; } /* ---------------------------------------------------------------------*/ /* several possibilities for pLDeg: the maximal degree of a monomial in p*/ /* compute in l also the pLength of p */ /*2 * compute the length of a polynomial (in l) * and the degree of the monomial with maximal degree: the last one */ long pLDeg0(poly p,int *l, const ring r) { p_CheckPolyRing(p, r); long k= p_GetComp(p, r); int ll=1; if (k > 0) { while ((pNext(p)!=NULL) && (p_GetComp(pNext(p), r)==k)) { pIter(p); ll++; } } else { while (pNext(p)!=NULL) { pIter(p); ll++; } } *l=ll; return r->pFDeg(p, r); } /*2 * compute the length of a polynomial (in l) * and the degree of the monomial with maximal degree: the last one * but search in all components before syzcomp */ long pLDeg0c(poly p,int *l, const ring r) { assume(p!=NULL); #ifdef PDEBUG _p_Test(p,r,PDEBUG); #endif p_CheckPolyRing(p, r); long o; int ll=1; if (! rIsSyzIndexRing(r)) { while (pNext(p) != NULL) { pIter(p); ll++; } o = r->pFDeg(p, r); } else { int curr_limit = rGetCurrSyzLimit(r); poly pp = p; while ((p=pNext(p))!=NULL) { if (p_GetComp(p, r)<=curr_limit/*syzComp*/) ll++; else break; pp = p; } #ifdef PDEBUG _p_Test(pp,r,PDEBUG); #endif o = r->pFDeg(pp, r); } *l=ll; return o; } /*2 * compute the length of a polynomial (in l) * and the degree of the monomial with maximal degree: the first one * this works for the polynomial case with degree orderings * (both c,dp and dp,c) */ long pLDegb(poly p,int *l, const ring r) { p_CheckPolyRing(p, r); long k= p_GetComp(p, r); long o = r->pFDeg(p, r); int ll=1; if (k != 0) { while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k)) { ll++; } } else { while ((p=pNext(p)) !=NULL) { ll++; } } *l=ll; return o; } /*2 * compute the length of a polynomial (in l) * and the degree of the monomial with maximal degree: * this is NOT the last one, we have to look for it */ long pLDeg1(poly p,int *l, const ring r) { p_CheckPolyRing(p, r); long k= p_GetComp(p, r); int ll=1; long t,max; max=r->pFDeg(p, r); if (k > 0) { while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k)) { t=r->pFDeg(p, r); if (t>max) max=t; ll++; } } else { while ((p=pNext(p))!=NULL) { t=r->pFDeg(p, r); if (t>max) max=t; ll++; } } *l=ll; return max; } /*2 * compute the length of a polynomial (in l) * and the degree of the monomial with maximal degree: * this is NOT the last one, we have to look for it * in all components */ long pLDeg1c(poly p,int *l, const ring r) { p_CheckPolyRing(p, r); int ll=1; long t,max; max=r->pFDeg(p, r); if (rIsSyzIndexRing(r)) { long limit = rGetCurrSyzLimit(r); while ((p=pNext(p))!=NULL) { if (p_GetComp(p, r)<=limit) { if ((t=r->pFDeg(p, r))>max) max=t; ll++; } else break; } } else { while ((p=pNext(p))!=NULL) { if ((t=r->pFDeg(p, r))>max) max=t; ll++; } } *l=ll; return max; } // like pLDeg1, only pFDeg == pDeg long pLDeg1_Deg(poly p,int *l, const ring r) { assume(r->pFDeg == p_Deg); p_CheckPolyRing(p, r); long k= p_GetComp(p, r); int ll=1; long t,max; max=p_GetOrder(p, r); if (k > 0) { while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k)) { t=p_GetOrder(p, r); if (t>max) max=t; ll++; } } else { while ((p=pNext(p))!=NULL) { t=p_GetOrder(p, r); if (t>max) max=t; ll++; } } *l=ll; return max; } long pLDeg1c_Deg(poly p,int *l, const ring r) { assume(r->pFDeg == p_Deg); p_CheckPolyRing(p, r); int ll=1; long t,max; max=p_GetOrder(p, r); if (rIsSyzIndexRing(r)) { long limit = rGetCurrSyzLimit(r); while ((p=pNext(p))!=NULL) { if (p_GetComp(p, r)<=limit) { if ((t=p_GetOrder(p, r))>max) max=t; ll++; } else break; } } else { while ((p=pNext(p))!=NULL) { if ((t=p_GetOrder(p, r))>max) max=t; ll++; } } *l=ll; return max; } // like pLDeg1, only pFDeg == pTotoalDegree long pLDeg1_Totaldegree(poly p,int *l, const ring r) { p_CheckPolyRing(p, r); long k= p_GetComp(p, r); int ll=1; long t,max; max=p_Totaldegree(p, r); if (k > 0) { while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k)) { t=p_Totaldegree(p, r); if (t>max) max=t; ll++; } } else { while ((p=pNext(p))!=NULL) { t=p_Totaldegree(p, r); if (t>max) max=t; ll++; } } *l=ll; return max; } long pLDeg1c_Totaldegree(poly p,int *l, const ring r) { p_CheckPolyRing(p, r); int ll=1; long t,max; max=p_Totaldegree(p, r); if (rIsSyzIndexRing(r)) { long limit = rGetCurrSyzLimit(r); while ((p=pNext(p))!=NULL) { if (p_GetComp(p, r)<=limit) { if ((t=p_Totaldegree(p, r))>max) max=t; ll++; } else break; } } else { while ((p=pNext(p))!=NULL) { if ((t=p_Totaldegree(p, r))>max) max=t; ll++; } } *l=ll; return max; } // like pLDeg1, only pFDeg == p_WFirstTotalDegree long pLDeg1_WFirstTotalDegree(poly p,int *l, const ring r) { p_CheckPolyRing(p, r); long k= p_GetComp(p, r); int ll=1; long t,max; max=p_WFirstTotalDegree(p, r); if (k > 0) { while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k)) { t=p_WFirstTotalDegree(p, r); if (t>max) max=t; ll++; } } else { while ((p=pNext(p))!=NULL) { t=p_WFirstTotalDegree(p, r); if (t>max) max=t; ll++; } } *l=ll; return max; } long pLDeg1c_WFirstTotalDegree(poly p,int *l, const ring r) { p_CheckPolyRing(p, r); int ll=1; long t,max; max=p_WFirstTotalDegree(p, r); if (rIsSyzIndexRing(r)) { long limit = rGetCurrSyzLimit(r); while ((p=pNext(p))!=NULL) { if (p_GetComp(p, r)<=limit) { if ((t=p_Totaldegree(p, r))>max) max=t; ll++; } else break; } } else { while ((p=pNext(p))!=NULL) { if ((t=p_Totaldegree(p, r))>max) max=t; ll++; } } *l=ll; return max; } /*************************************************************** * * Maximal Exponent business * ***************************************************************/ static inline unsigned long p_GetMaxExpL2(unsigned long l1, unsigned long l2, const ring r, unsigned long number_of_exp) { const unsigned long bitmask = r->bitmask; unsigned long ml1 = l1 & bitmask; unsigned long ml2 = l2 & bitmask; unsigned long max = (ml1 > ml2 ? ml1 : ml2); unsigned long j = number_of_exp - 1; if (j > 0) { unsigned long mask = bitmask << r->BitsPerExp; while (1) { ml1 = l1 & mask; ml2 = l2 & mask; max |= ((ml1 > ml2 ? ml1 : ml2) & mask); j--; if (j == 0) break; mask = mask << r->BitsPerExp; } } return max; } static inline unsigned long p_GetMaxExpL2(unsigned long l1, unsigned long l2, const ring r) { return p_GetMaxExpL2(l1, l2, r, r->ExpPerLong); } poly p_GetMaxExpP(poly p, const ring r) { p_CheckPolyRing(p, r); if (p == NULL) return p_Init(r); poly max = p_LmInit(p, r); pIter(p); if (p == NULL) return max; int i, offset; unsigned long l_p, l_max; unsigned long divmask = r->divmask; do { offset = r->VarL_Offset[0]; l_p = p->exp[offset]; l_max = max->exp[offset]; // do the divisibility trick to find out whether l has an exponent if (l_p > l_max || (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask))) max->exp[offset] = p_GetMaxExpL2(l_max, l_p, r); for (i=1; iVarL_Size; i++) { offset = r->VarL_Offset[i]; l_p = p->exp[offset]; l_max = max->exp[offset]; // do the divisibility trick to find out whether l has an exponent if (l_p > l_max || (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask))) max->exp[offset] = p_GetMaxExpL2(l_max, l_p, r); } pIter(p); } while (p != NULL); return max; } unsigned long p_GetMaxExpL(poly p, const ring r, unsigned long l_max) { unsigned long l_p, divmask = r->divmask; int i; while (p != NULL) { l_p = p->exp[r->VarL_Offset[0]]; if (l_p > l_max || (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask))) l_max = p_GetMaxExpL2(l_max, l_p, r); for (i=1; iVarL_Size; i++) { l_p = p->exp[r->VarL_Offset[i]]; // do the divisibility trick to find out whether l has an exponent if (l_p > l_max || (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask))) l_max = p_GetMaxExpL2(l_max, l_p, r); } pIter(p); } return l_max; } /*************************************************************** * * Misc things * ***************************************************************/ // returns TRUE, if all monoms have the same component BOOLEAN p_OneComp(poly p, const ring r) { if(p!=NULL) { long i = p_GetComp(p, r); while (pNext(p)!=NULL) { pIter(p); if(i != p_GetComp(p, r)) return FALSE; } } return TRUE; } /*2 *test if a monomial /head term is a pure power */ int p_IsPurePower(const poly p, const ring r) { int i,k=0; for (i=r->N;i;i--) { if (p_GetExp(p,i, r)!=0) { if(k!=0) return 0; k=i; } } return k; } /*2 *test if a polynomial is univariate * return -1 for constant, * 0 for not univariate,s * i if dep. on var(i) */ int p_IsUnivariate(poly p, const ring r) { int i,k=-1; while (p!=NULL) { for (i=r->N;i;i--) { if (p_GetExp(p,i, r)!=0) { if((k!=-1)&&(k!=i)) return 0; k=i; } } pIter(p); } return k; } // set entry e[i] to 1 if var(i) occurs in p, ignore var(j) if e[j]>0 int p_GetVariables(poly p, int * e, const ring r) { int i; int n=0; while(p!=NULL) { n=0; for(i=r->N; i>0; i--) { if(e[i]==0) { if (p_GetExp(p,i,r)>0) { e[i]=1; n++; } } else n++; } if (n==r->N) break; pIter(p); } return n; } /*2 * returns a polynomial representing the integer i */ poly p_ISet(long i, const ring r) { poly rc = NULL; if (i!=0) { rc = p_Init(r); pSetCoeff0(rc,n_Init(i,r->cf)); if (n_IsZero(pGetCoeff(rc),r->cf)) p_LmDelete(&rc,r); } return rc; } /*2 * an optimized version of p_ISet for the special case 1 */ poly p_One(const ring r) { poly rc = p_Init(r); pSetCoeff0(rc,n_Init(1,r->cf)); return rc; } void p_Split(poly p, poly *h) { *h=pNext(p); pNext(p)=NULL; } /*2 * pair has no common factor ? or is no polynomial */ BOOLEAN p_HasNotCF(poly p1, poly p2, const ring r) { if (p_GetComp(p1,r) > 0 || p_GetComp(p2,r) > 0) return FALSE; int i = rVar(r); loop { if ((p_GetExp(p1, i, r) > 0) && (p_GetExp(p2, i, r) > 0)) return FALSE; i--; if (i == 0) return TRUE; } } /*2 * convert monomial given as string to poly, e.g. 1x3y5z */ const char * p_Read(const char *st, poly &rc, const ring r) { if (r==NULL) { rc=NULL;return st;} int i,j; rc = p_Init(r); const char *s = r->cf->cfRead(st,&(rc->coef),r->cf); if (s==st) /* i.e. it does not start with a coeff: test if it is a ringvar*/ { j = r_IsRingVar(s,r); if (j >= 0) { p_IncrExp(rc,1+j,r); while (*s!='\0') s++; goto done; } } while (*s!='\0') { char ss[2]; ss[0] = *s++; ss[1] = '\0'; j = r_IsRingVar(ss,r); if (j >= 0) { const char *s_save=s; s = eati(s,&i); if (((unsigned long)i) > r->bitmask) { // exponent to large: it is not a monomial p_LmDelete(&rc,r); return s_save; } p_AddExp(rc,1+j, (long)i, r); } else { // 1st char of is not a varname // We return the parsed polynomial nevertheless. This is needed when // we are parsing coefficients in a rational function field. s--; break; } } done: if (n_IsZero(pGetCoeff(rc),r->cf)) p_LmDelete(&rc,r); else { #ifdef HAVE_PLURAL // in super-commutative ring // squares of anti-commutative variables are zeroes! if(rIsSCA(r)) { const unsigned int iFirstAltVar = scaFirstAltVar(r); const unsigned int iLastAltVar = scaLastAltVar(r); assume(rc != NULL); for(unsigned int k = iFirstAltVar; k <= iLastAltVar; k++) if( p_GetExp(rc, k, r) > 1 ) { p_LmDelete(&rc, r); goto finish; } } #endif p_Setm(rc,r); } finish: return s; } poly p_mInit(const char *st, BOOLEAN &ok, const ring r) { poly p; const char *s=p_Read(st,p,r); if (*s!='\0') { if ((s!=st)&&isdigit(st[0])) { errorreported=TRUE; } ok=FALSE; p_Delete(&p,r); return NULL; } #ifdef PDEBUG _p_Test(p,r,PDEBUG); #endif ok=!errorreported; return p; } /*2 * returns a polynomial representing the number n * destroys n */ poly p_NSet(number n, const ring r) { if (n_IsZero(n,r->cf)) { n_Delete(&n, r->cf); return NULL; } else { poly rc = p_Init(r); pSetCoeff0(rc,n); return rc; } } /*2 * assumes that LM(a) = LM(b)*m, for some monomial m, * returns the multiplicant m, * leaves a and b unmodified */ poly p_Divide(poly a, poly b, const ring r) { assume((p_GetComp(a,r)==p_GetComp(b,r)) || (p_GetComp(b,r)==0)); int i; poly result = p_Init(r); for(i=(int)r->N; i; i--) p_SetExp(result,i, p_GetExp(a,i,r)- p_GetExp(b,i,r),r); p_SetComp(result, p_GetComp(a,r) - p_GetComp(b,r),r); p_Setm(result,r); return result; } poly p_Div_nn(poly p, const number n, const ring r) { pAssume(!n_IsZero(n,r->cf)); p_Test(p, r); poly q = p; while (p != NULL) { number nc = pGetCoeff(p); pSetCoeff0(p, n_Div(nc, n, r->cf)); n_Delete(&nc, r->cf); pIter(p); } p_Test(q, r); return q; } /*2 * divides a by the monomial b, ignores monomials which are not divisible * assumes that b is not NULL, destroyes b */ poly p_DivideM(poly a, poly b, const ring r) { if (a==NULL) { p_Delete(&b,r); return NULL; } poly result=a; poly prev=NULL; int i; #ifdef HAVE_RINGS number inv=pGetCoeff(b); #else number inv=n_Invers(pGetCoeff(b),r->cf); #endif while (a!=NULL) { if (p_DivisibleBy(b,a,r)) { for(i=(int)r->N; i; i--) p_SubExp(a,i, p_GetExp(b,i,r),r); p_SubComp(a, p_GetComp(b,r),r); p_Setm(a,r); prev=a; pIter(a); } else { if (prev==NULL) { p_LmDelete(&result,r); a=result; } else { p_LmDelete(&pNext(prev),r); a=pNext(prev); } } } #ifdef HAVE_RINGS if (n_IsUnit(inv,r->cf)) { inv = n_Invers(inv,r->cf); p_Mult_nn(result,inv,r); n_Delete(&inv, r->cf); } else { p_Div_nn(result,inv,r); } #else p_Mult_nn(result,inv,r); n_Delete(&inv, r->cf); #endif p_Delete(&b, r); return result; } #ifdef HAVE_RINGS /* TRUE iff LT(f) | LT(g) */ BOOLEAN p_DivisibleByRingCase(poly f, poly g, const ring r) { int exponent; for(int i = (int)rVar(r); i>0; i--) { exponent = p_GetExp(g, i, r) - p_GetExp(f, i, r); if (exponent < 0) return FALSE; } return n_DivBy(pGetCoeff(g), pGetCoeff(f), r->cf); } #endif // returns the LCM of the head terms of a and b in *m void p_Lcm(const poly a, const poly b, poly m, const ring r) { for (int i=rVar(r); i; --i) p_SetExp(m,i, si_max( p_GetExp(a,i,r), p_GetExp(b,i,r)),r); p_SetComp(m, si_max(p_GetComp(a,r), p_GetComp(b,r)),r); /* Don't do a pSetm here, otherwise hres/lres chockes */ } /* assumes that p and divisor are univariate polynomials in r, mentioning the same variable; assumes divisor != NULL; p may be NULL; assumes a global monomial ordering in r; performs polynomial division of p by divisor: - afterwards p contains the remainder of the division, i.e., p_before = result * divisor + p_afterwards; - if needResult == TRUE, then the method computes and returns 'result', otherwise NULL is returned (This parametrization can be used when one is only interested in the remainder of the division. In this case, the method will be slightly faster.) leaves divisor unmodified */ poly p_PolyDiv(poly &p, const poly divisor, const BOOLEAN needResult, const ring r) { assume(divisor != NULL); if (p == NULL) return NULL; poly result = NULL; number divisorLC = p_GetCoeff(divisor, r); int divisorLE = p_GetExp(divisor, 1, r); while ((p != NULL) && (p_Deg(p, r) >= p_Deg(divisor, r))) { /* determine t = LT(p) / LT(divisor) */ poly t = p_ISet(1, r); number c = n_Div(p_GetCoeff(p, r), divisorLC, r->cf); n_Normalize(c,r->cf); p_SetCoeff(t, c, r); int e = p_GetExp(p, 1, r) - divisorLE; p_SetExp(t, 1, e, r); p_Setm(t, r); if (needResult) result = p_Add_q(result, p_Copy(t, r), r); p = p_Add_q(p, p_Neg(p_Mult_q(t, p_Copy(divisor, r), r), r), r); } return result; } /*2 * returns the partial differentiate of a by the k-th variable * does not destroy the input */ poly p_Diff(poly a, int k, const ring r) { poly res, f, last; number t; last = res = NULL; while (a!=NULL) { if (p_GetExp(a,k,r)!=0) { f = p_LmInit(a,r); t = n_Init(p_GetExp(a,k,r),r->cf); pSetCoeff0(f,n_Mult(t,pGetCoeff(a),r->cf)); n_Delete(&t,r->cf); if (n_IsZero(pGetCoeff(f),r->cf)) p_LmDelete(&f,r); else { p_DecrExp(f,k,r); p_Setm(f,r); if (res==NULL) { res=last=f; } else { pNext(last)=f; last=f; } } } pIter(a); } return res; } static poly p_DiffOpM(poly a, poly b,BOOLEAN multiply, const ring r) { int i,j,s; number n,h,hh; poly p=p_One(r); n=n_Mult(pGetCoeff(a),pGetCoeff(b),r->cf); for(i=rVar(r);i>0;i--) { s=p_GetExp(b,i,r); if (scf); p_LmDelete(&p,r); return NULL; } if (multiply) { for(j=p_GetExp(a,i,r); j>0;j--) { h = n_Init(s,r->cf); hh=n_Mult(n,h,r->cf); n_Delete(&h,r->cf); n_Delete(&n,r->cf); n=hh; s--; } p_SetExp(p,i,s,r); } else { p_SetExp(p,i,s-p_GetExp(a,i,r),r); } } p_Setm(p,r); /*if (multiply)*/ p_SetCoeff(p,n,r); if (n_IsZero(n,r->cf)) p=p_LmDeleteAndNext(p,r); // return NULL as p is a monomial return p; } poly p_DiffOp(poly a, poly b,BOOLEAN multiply, const ring r) { poly result=NULL; poly h; for(;a!=NULL;pIter(a)) { for(h=b;h!=NULL;pIter(h)) { result=p_Add_q(result,p_DiffOpM(a,h,multiply,r),r); } } return result; } /*2 * subtract p2 from p1, p1 and p2 are destroyed * do not put attention on speed: the procedure is only used in the interpreter */ poly p_Sub(poly p1, poly p2, const ring r) { return p_Add_q(p1, p_Neg(p2,r),r); } /*3 * compute for a monomial m * the power m^exp, exp > 1 * destroys p */ static poly p_MonPower(poly p, int exp, const ring r) { int i; if(!n_IsOne(pGetCoeff(p),r->cf)) { number x, y; y = pGetCoeff(p); n_Power(y,exp,&x,r->cf); n_Delete(&y,r->cf); pSetCoeff0(p,x); } for (i=rVar(r); i!=0; i--) { p_MultExp(p,i, exp,r); } p_Setm(p,r); return p; } /*3 * compute for monomials p*q * destroys p, keeps q */ static void p_MonMult(poly p, poly q, const ring r) { number x, y; y = pGetCoeff(p); x = n_Mult(y,pGetCoeff(q),r->cf); n_Delete(&y,r->cf); pSetCoeff0(p,x); //for (int i=pVariables; i!=0; i--) //{ // pAddExp(p,i, pGetExp(q,i)); //} //p->Order += q->Order; p_ExpVectorAdd(p,q,r); } /*3 * compute for monomials p*q * keeps p, q */ static poly p_MonMultC(poly p, poly q, const ring rr) { number x; poly r = p_Init(rr); x = n_Mult(pGetCoeff(p),pGetCoeff(q),rr->cf); pSetCoeff0(r,x); p_ExpVectorSum(r,p, q, rr); return r; } /*3 * create binomial coef. */ static number* pnBin(int exp, const ring r) { int e, i, h; number x, y, *bin=NULL; x = n_Init(exp,r->cf); if (n_IsZero(x,r->cf)) { n_Delete(&x,r->cf); return bin; } h = (exp >> 1) + 1; bin = (number *)omAlloc0(h*sizeof(number)); bin[1] = x; if (exp < 4) return bin; i = exp - 1; for (e=2; ecf); i--; y = n_Mult(x,bin[e-1],r->cf); n_Delete(&x,r->cf); x = n_Init(e,r->cf); bin[e] = n_IntDiv(y,x,r->cf); n_Delete(&x,r->cf); n_Delete(&y,r->cf); } return bin; } static void pnFreeBin(number *bin, int exp,const coeffs r) { int e, h = (exp >> 1) + 1; if (bin[1] != NULL) { for (e=1; e 1 * with binomial coef. */ static poly p_TwoMonPower(poly p, int exp, const ring r) { int eh, e; long al; poly *a; poly tail, b, res, h; number x; number *bin = pnBin(exp,r); tail = pNext(p); if (bin == NULL) { p_MonPower(p,exp,r); p_MonPower(tail,exp,r); #ifdef PDEBUG p_Test(p,r); #endif return p; } eh = exp >> 1; al = (exp + 1) * sizeof(poly); a = (poly *)omAlloc(al); a[1] = p; for (e=1; eeh; e--) { h = a[e]; x = n_Mult(bin[exp-e],pGetCoeff(h),r->cf); p_SetCoeff(h,x,r); p_MonMult(h,b,r); res = pNext(res) = h; p_MonMult(b,tail,r); } for (e=eh; e!=0; e--) { h = a[e]; x = n_Mult(bin[e],pGetCoeff(h),r->cf); p_SetCoeff(h,x,r); p_MonMult(h,b,r); res = pNext(res) = h; p_MonMult(b,tail,r); } p_LmDelete(&tail,r); pNext(res) = b; pNext(b) = NULL; res = a[exp]; omFreeSize((ADDRESS)a, al); pnFreeBin(bin, exp, r->cf); // tail=res; // while((tail!=NULL)&&(pNext(tail)!=NULL)) // { // if(nIsZero(pGetCoeff(pNext(tail)))) // { // pLmDelete(&pNext(tail)); // } // else // pIter(tail); // } #ifdef PDEBUG p_Test(res,r); #endif return res; } static poly p_Pow(poly p, int i, const ring r) { poly rc = p_Copy(p,r); i -= 2; do { rc = p_Mult_q(rc,p_Copy(p,r),r); p_Normalize(rc,r); i--; } while (i != 0); return p_Mult_q(rc,p,r); } /*2 * returns the i-th power of p * p will be destroyed */ poly p_Power(poly p, int i, const ring r) { poly rc=NULL; if (i==0) { p_Delete(&p,r); return p_One(r); } if(p!=NULL) { if ( (i > 0) && ((unsigned long ) i > (r->bitmask))) { Werror("exponent %d is too large, max. is %ld",i,r->bitmask); return NULL; } switch (i) { // cannot happen, see above // case 0: // { // rc=pOne(); // pDelete(&p); // break; // } case 1: rc=p; break; case 2: rc=p_Mult_q(p_Copy(p,r),p,r); break; default: if (i < 0) { p_Delete(&p,r); return NULL; } else { #ifdef HAVE_PLURAL if (rIsPluralRing(r)) /* in the NC case nothing helps :-( */ { int j=i; rc = p_Copy(p,r); while (j>1) { rc = p_Mult_q(p_Copy(p,r),rc,r); j--; } p_Delete(&p,r); return rc; } #endif rc = pNext(p); if (rc == NULL) return p_MonPower(p,i,r); /* else: binom ?*/ int char_p=rChar(r); if ((pNext(rc) != NULL) #ifdef HAVE_RINGS || rField_is_Ring(r) #endif ) return p_Pow(p,i,r); if ((char_p==0) || (i<=char_p)) return p_TwoMonPower(p,i,r); return p_Pow(p,i,r); } /*end default:*/ } } return rc; } /* --------------------------------------------------------------------------------*/ /* content suff */ static number p_InitContent(poly ph, const ring r); #define CLEARENUMERATORS 1 void p_Content(poly ph, const ring r) { assume( ph != NULL ); assume( r != NULL ); assume( r->cf != NULL ); #if CLEARENUMERATORS if( 0 ) { const coeffs C = r->cf; // experimentall (recursive enumerator treatment) of alg. Ext! CPolyCoeffsEnumerator itr(ph); n_ClearContent(itr, r->cf); p_Test(ph, r); n_Test(pGetCoeff(ph), C); assume(n_GreaterZero(pGetCoeff(ph), C)); // ?? // if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r); return; } #endif #ifdef HAVE_RINGS if (rField_is_Ring(r)) { if (rField_has_Units(r)) { number k = n_GetUnit(pGetCoeff(ph),r->cf); if (!n_IsOne(k,r->cf)) { number tmpGMP = k; k = n_Invers(k,r->cf); n_Delete(&tmpGMP,r->cf); poly h = pNext(ph); p_SetCoeff(ph, n_Mult(pGetCoeff(ph), k,r->cf),r); while (h != NULL) { p_SetCoeff(h, n_Mult(pGetCoeff(h), k,r->cf),r); pIter(h); } // assume( n_GreaterZero(pGetCoeff(ph),r->cf) ); // if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r); } n_Delete(&k,r->cf); } return; } #endif number h,d; poly p; if(TEST_OPT_CONTENTSB) return; if(pNext(ph)==NULL) { p_SetCoeff(ph,n_Init(1,r->cf),r); } else { assume( pNext(ph) != NULL ); #if CLEARENUMERATORS if( nCoeff_is_Q(r->cf) || nCoeff_is_Q_a(r->cf) ) { // experimentall (recursive enumerator treatment) of alg. Ext! CPolyCoeffsEnumerator itr(ph); n_ClearContent(itr, r->cf); p_Test(ph, r); n_Test(pGetCoeff(ph), r->cf); assume(n_GreaterZero(pGetCoeff(ph), r->cf)); // ?? // if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r); return; } #endif n_Normalize(pGetCoeff(ph),r->cf); if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r); if (rField_is_Q(r)) // should not be used anymore if CLEARENUMERATORS is 1 { h=p_InitContent(ph,r); p=ph; } else { h=n_Copy(pGetCoeff(ph),r->cf); p = pNext(ph); } while (p!=NULL) { n_Normalize(pGetCoeff(p),r->cf); d=n_Gcd(h,pGetCoeff(p),r->cf); n_Delete(&h,r->cf); h = d; if(n_IsOne(h,r->cf)) { break; } pIter(p); } p = ph; //number tmp; if(!n_IsOne(h,r->cf)) { while (p!=NULL) { //d = nDiv(pGetCoeff(p),h); //tmp = nIntDiv(pGetCoeff(p),h); //if (!nEqual(d,tmp)) //{ // StringSetS("** div0:");nWrite(pGetCoeff(p));StringAppendS("/"); // nWrite(h);StringAppendS("=");nWrite(d);StringAppendS(" int:"); // nWrite(tmp);Print(StringEndS("\n")); // NOTE/TODO: use StringAppendS("\n"); omFree(s); //} //nDelete(&tmp); d = n_IntDiv(pGetCoeff(p),h,r->cf); p_SetCoeff(p,d,r); pIter(p); } } n_Delete(&h,r->cf); #ifdef HAVE_FACTORY // if ( (n_GetChar(r) == 1) || (n_GetChar(r) < 0) ) /* Q[a],Q(a),Zp[a],Z/p(a) */ // { // singclap_divide_content(ph, r); // if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r); // } #endif if (rField_is_Q_a(r)) { // we only need special handling for alg. ext. if (getCoeffType(r->cf)==n_algExt) { number hzz = n_Init(1, r->cf->extRing->cf); p=ph; while (p!=NULL) { // each monom: coeff in Q_a poly c_n_n=(poly)pGetCoeff(p); poly c_n=c_n_n; while (c_n!=NULL) { // each monom: coeff in Q d=n_Lcm(hzz,pGetCoeff(c_n),r->cf->extRing->cf); n_Delete(&hzz,r->cf->extRing->cf); hzz=d; pIter(c_n); } pIter(p); } /* hzz contains the 1/lcm of all denominators in c_n_n*/ h=n_Invers(hzz,r->cf->extRing->cf); n_Delete(&hzz,r->cf->extRing->cf); n_Normalize(h,r->cf->extRing->cf); if(!n_IsOne(h,r->cf->extRing->cf)) { p=ph; while (p!=NULL) { // each monom: coeff in Q_a poly c_n=(poly)pGetCoeff(p); while (c_n!=NULL) { // each monom: coeff in Q d=n_Mult(h,pGetCoeff(c_n),r->cf->extRing->cf); n_Normalize(d,r->cf->extRing->cf); n_Delete(&pGetCoeff(c_n),r->cf->extRing->cf); pGetCoeff(c_n)=d; pIter(c_n); } pIter(p); } } n_Delete(&h,r->cf->extRing->cf); } } } if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r); } // Not yet? #if 1 // currently only used by Singular/janet void p_SimpleContent(poly ph, int smax, const ring r) { if(TEST_OPT_CONTENTSB) return; if (ph==NULL) return; if (pNext(ph)==NULL) { p_SetCoeff(ph,n_Init(1,r->cf),r); return; } if ((pNext(pNext(ph))==NULL)||(!rField_is_Q(r))) { return; } number d=p_InitContent(ph,r); if (n_Size(d,r->cf)<=smax) { //if (TEST_OPT_PROT) PrintS("G"); return; } poly p=ph; number h=d; if (smax==1) smax=2; while (p!=NULL) { #if 0 d=nlGcd(h,pGetCoeff(p),r->cf); nlDelete(&h,r->cf); h = d; #else nlInpGcd(h,pGetCoeff(p),r->cf); #endif if(nlSize(h,r->cf)cf)) h=nlNeg(h,r->cf); if(nlIsOne(h,r->cf)) return; //if (TEST_OPT_PROT) PrintS("c"); while (p!=NULL) { #if 1 d = nlIntDiv(pGetCoeff(p),h,r->cf); p_SetCoeff(p,d,r); #else nlInpIntDiv(pGetCoeff(p),h,r->cf); #endif pIter(p); } nlDelete(&h,r->cf); } #endif static number p_InitContent(poly ph, const ring r) // only for coefficients in Q #if 0 { assume(!TEST_OPT_CONTENTSB); assume(ph!=NULL); assume(pNext(ph)!=NULL); assume(rField_is_Q(r)); if (pNext(pNext(ph))==NULL) { return nlGetNom(pGetCoeff(pNext(ph)),r->cf); } poly p=ph; number n1=nlGetNom(pGetCoeff(p),r->cf); pIter(p); number n2=nlGetNom(pGetCoeff(p),r->cf); pIter(p); number d; number t; loop { nlNormalize(pGetCoeff(p),r->cf); t=nlGetNom(pGetCoeff(p),r->cf); if (nlGreaterZero(t,r->cf)) d=nlAdd(n1,t,r->cf); else d=nlSub(n1,t,r->cf); nlDelete(&t,r->cf); nlDelete(&n1,r->cf); n1=d; pIter(p); if (p==NULL) break; nlNormalize(pGetCoeff(p),r->cf); t=nlGetNom(pGetCoeff(p),r->cf); if (nlGreaterZero(t,r->cf)) d=nlAdd(n2,t,r->cf); else d=nlSub(n2,t,r->cf); nlDelete(&t,r->cf); nlDelete(&n2,r->cf); n2=d; pIter(p); if (p==NULL) break; } d=nlGcd(n1,n2,r->cf); nlDelete(&n1,r->cf); nlDelete(&n2,r->cf); return d; } #else { number d=pGetCoeff(ph); if(SR_HDL(d)&SR_INT) return d; int s=mpz_size1(d->z); int s2=-1; number d2; loop { pIter(ph); if(ph==NULL) { if (s2==-1) return nlCopy(d,r->cf); break; } if (SR_HDL(pGetCoeff(ph))&SR_INT) { s2=s; d2=d; s=0; d=pGetCoeff(ph); if (s2==0) break; } else if (mpz_size1((pGetCoeff(ph)->z))<=s) { s2=s; d2=d; d=pGetCoeff(ph); s=mpz_size1(d->z); } } return nlGcd(d,d2,r->cf); } #endif //void pContent(poly ph) //{ // number h,d; // poly p; // // p = ph; // if(pNext(p)==NULL) // { // pSetCoeff(p,nInit(1)); // } // else // { //#ifdef PDEBUG // if (!pTest(p)) return; //#endif // nNormalize(pGetCoeff(p)); // if(!nGreaterZero(pGetCoeff(ph))) // { // ph = pNeg(ph); // nNormalize(pGetCoeff(p)); // } // h=pGetCoeff(p); // pIter(p); // while (p!=NULL) // { // nNormalize(pGetCoeff(p)); // if (nGreater(h,pGetCoeff(p))) h=pGetCoeff(p); // pIter(p); // } // h=nCopy(h); // p=ph; // while (p!=NULL) // { // d=n_Gcd(h,pGetCoeff(p)); // nDelete(&h); // h = d; // if(nIsOne(h)) // { // break; // } // pIter(p); // } // p = ph; // //number tmp; // if(!nIsOne(h)) // { // while (p!=NULL) // { // d = nIntDiv(pGetCoeff(p),h); // pSetCoeff(p,d); // pIter(p); // } // } // nDelete(&h); //#ifdef HAVE_FACTORY // if ( (nGetChar() == 1) || (nGetChar() < 0) ) /* Q[a],Q(a),Zp[a],Z/p(a) */ // { // pTest(ph); // singclap_divide_content(ph); // pTest(ph); // } //#endif // } //} #if 0 void p_Content(poly ph, const ring r) { number h,d; poly p; if(pNext(ph)==NULL) { p_SetCoeff(ph,n_Init(1,r->cf),r); } else { n_Normalize(pGetCoeff(ph),r->cf); if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r); h=n_Copy(pGetCoeff(ph),r->cf); p = pNext(ph); while (p!=NULL) { n_Normalize(pGetCoeff(p),r->cf); d=n_Gcd(h,pGetCoeff(p),r->cf); n_Delete(&h,r->cf); h = d; if(n_IsOne(h,r->cf)) { break; } pIter(p); } p = ph; //number tmp; if(!n_IsOne(h,r->cf)) { while (p!=NULL) { //d = nDiv(pGetCoeff(p),h); //tmp = nIntDiv(pGetCoeff(p),h); //if (!nEqual(d,tmp)) //{ // StringSetS("** div0:");nWrite(pGetCoeff(p));StringAppendS("/"); // nWrite(h);StringAppendS("=");nWrite(d);StringAppendS(" int:"); // nWrite(tmp);Print(StringEndS("\n")); // NOTE/TODO: use StringAppendS("\n"); omFree(s); //} //nDelete(&tmp); d = n_IntDiv(pGetCoeff(p),h,r->cf); p_SetCoeff(p,d,r->cf); pIter(p); } } n_Delete(&h,r->cf); #ifdef HAVE_FACTORY //if ( (n_GetChar(r) == 1) || (n_GetChar(r) < 0) ) /* Q[a],Q(a),Zp[a],Z/p(a) */ //{ // singclap_divide_content(ph); // if(!n_GreaterZero(pGetCoeff(ph),r)) ph = p_Neg(ph,r); //} #endif } } #endif /* ---------------------------------------------------------------------------*/ /* cleardenom suff */ poly p_Cleardenom(poly ph, const ring r) { if( ph == NULL ) return NULL; assume( r != NULL ); assume( r->cf != NULL ); const coeffs C = r->cf; #if CLEARENUMERATORS if( 0 ) { CPolyCoeffsEnumerator itr(ph); n_ClearDenominators(itr, C); n_ClearContent(itr, C); // divide out the content p_Test(ph, r); n_Test(pGetCoeff(ph), C); assume(n_GreaterZero(pGetCoeff(ph), C)); // ?? // if(!n_GreaterZero(pGetCoeff(ph),C)) ph = p_Neg(ph,r); return ph; } #endif poly start=ph; number d, h; poly p; #ifdef HAVE_RINGS if (rField_is_Ring(r)) { p_Content(ph,r); assume( n_GreaterZero(pGetCoeff(ph),C) ); if(!n_GreaterZero(pGetCoeff(ph),C)) ph = p_Neg(ph,r); return start; } #endif if (rField_is_Zp(r) && TEST_OPT_INTSTRATEGY) { if(!n_GreaterZero(pGetCoeff(ph),C)) ph = p_Neg(ph,r); return start; } p = ph; assume(p != NULL); if(pNext(p)==NULL) { if (TEST_OPT_CONTENTSB) { number n=n_GetDenom(pGetCoeff(p),r->cf); if (!n_IsOne(n,r->cf)) { number nn=n_Mult(pGetCoeff(p),n,r->cf); n_Normalize(nn,r->cf); p_SetCoeff(p,nn,r); } n_Delete(&n,r->cf); } else p_SetCoeff(p,n_Init(1,r->cf),r); assume( n_GreaterZero(pGetCoeff(ph),C) ); if(!n_GreaterZero(pGetCoeff(ph),C)) ph = p_Neg(ph,r); return start; } assume(pNext(p)!=NULL); #if CLEARENUMERATORS if( nCoeff_is_Q(C) || nCoeff_is_Q_a(C) ) { CPolyCoeffsEnumerator itr(ph); n_ClearDenominators(itr, C); n_ClearContent(itr, C); // divide out the content p_Test(ph, r); n_Test(pGetCoeff(ph), C); assume(n_GreaterZero(pGetCoeff(ph), C)); // ?? // if(!n_GreaterZero(pGetCoeff(ph),C)) ph = p_Neg(ph,r); return start; } #endif if(1) { h = n_Init(1,r->cf); while (p!=NULL) { n_Normalize(pGetCoeff(p),r->cf); d=n_Lcm(h,pGetCoeff(p),r->cf); n_Delete(&h,r->cf); h=d; pIter(p); } /* contains the 1/lcm of all denominators */ if(!n_IsOne(h,r->cf)) { p = ph; while (p!=NULL) { /* should be: * number hh; * nGetDenom(p->coef,&hh); * nMult(&h,&hh,&d); * nNormalize(d); * nDelete(&hh); * nMult(d,p->coef,&hh); * nDelete(&d); * nDelete(&(p->coef)); * p->coef =hh; */ d=n_Mult(h,pGetCoeff(p),r->cf); n_Normalize(d,r->cf); p_SetCoeff(p,d,r); pIter(p); } n_Delete(&h,r->cf); if (n_GetChar(r->cf)==1) { loop { h = n_Init(1,r->cf); p=ph; while (p!=NULL) { d=n_Lcm(h,pGetCoeff(p),r->cf); n_Delete(&h,r->cf); h=d; pIter(p); } /* contains the 1/lcm of all denominators */ if(!n_IsOne(h,r->cf)) { p = ph; while (p!=NULL) { /* should be: * number hh; * nGetDenom(p->coef,&hh); * nMult(&h,&hh,&d); * nNormalize(d); * nDelete(&hh); * nMult(d,p->coef,&hh); * nDelete(&d); * nDelete(&(p->coef)); * p->coef =hh; */ d=n_Mult(h,pGetCoeff(p),r->cf); n_Normalize(d,r->cf); p_SetCoeff(p,d,r); pIter(p); } n_Delete(&h,r->cf); } else { n_Delete(&h,r->cf); break; } } } } if (h!=NULL) n_Delete(&h,r->cf); p_Content(ph,r); #ifdef HAVE_RATGRING if (rIsRatGRing(r)) { /* quick unit detection in the rational case is done in gr_nc_bba */ pContentRat(ph); start=ph; } #endif } assume( n_GreaterZero(pGetCoeff(ph),C) ); if(!n_GreaterZero(pGetCoeff(ph),C)) ph = p_Neg(ph,r); return start; } void p_Cleardenom_n(poly ph,const ring r,number &c) { const coeffs C = r->cf; number d, h; assume( ph != NULL ); poly p = ph; #if CLEARENUMERATORS if( 0 ) { CPolyCoeffsEnumerator itr(ph); n_ClearDenominators(itr, d, C); // multiply with common denom. d n_ClearContent(itr, h, C); // divide by the content h c = n_Div(d, h, C); // d/h n_Delete(&d, C); n_Delete(&h, C); n_Test(c, C); p_Test(ph, r); n_Test(pGetCoeff(ph), C); assume(n_GreaterZero(pGetCoeff(ph), C)); // ?? /* if(!n_GreaterZero(pGetCoeff(ph),C)) { ph = p_Neg(ph,r); c = n_Neg(c, C); } */ return; } #endif if( pNext(p) == NULL ) { c=n_Invers(pGetCoeff(p), C); p_SetCoeff(p, n_Init(1, C), r); assume( n_GreaterZero(pGetCoeff(ph),C) ); if(!n_GreaterZero(pGetCoeff(ph),C)) { ph = p_Neg(ph,r); c = n_Neg(c, C); } return; } assume( pNext(p) != NULL ); #if CLEARENUMERATORS if( nCoeff_is_Q(C) || nCoeff_is_Q_a(C) ) { CPolyCoeffsEnumerator itr(ph); n_ClearDenominators(itr, d, C); // multiply with common denom. d n_ClearContent(itr, h, C); // divide by the content h c = n_Div(d, h, C); // d/h n_Delete(&d, C); n_Delete(&h, C); n_Test(c, C); p_Test(ph, r); n_Test(pGetCoeff(ph), C); assume(n_GreaterZero(pGetCoeff(ph), C)); // ?? /* if(!n_GreaterZero(pGetCoeff(ph),C)) { ph = p_Neg(ph,r); c = n_Neg(c, C); } */ return; } #endif if(1) { h = n_Init(1,r->cf); while (p!=NULL) { n_Normalize(pGetCoeff(p),r->cf); d=n_Lcm(h,pGetCoeff(p),r->cf); n_Delete(&h,r->cf); h=d; pIter(p); } c=h; /* contains the 1/lcm of all denominators */ if(!n_IsOne(h,r->cf)) { p = ph; while (p!=NULL) { /* should be: * number hh; * nGetDenom(p->coef,&hh); * nMult(&h,&hh,&d); * nNormalize(d); * nDelete(&hh); * nMult(d,p->coef,&hh); * nDelete(&d); * nDelete(&(p->coef)); * p->coef =hh; */ d=n_Mult(h,pGetCoeff(p),r->cf); n_Normalize(d,r->cf); p_SetCoeff(p,d,r); pIter(p); } if (rField_is_Q_a(r)) { loop { h = n_Init(1,r->cf); p=ph; while (p!=NULL) { d=n_Lcm(h,pGetCoeff(p),r->cf); n_Delete(&h,r->cf); h=d; pIter(p); } /* contains the 1/lcm of all denominators */ if(!n_IsOne(h,r->cf)) { p = ph; while (p!=NULL) { /* should be: * number hh; * nGetDenom(p->coef,&hh); * nMult(&h,&hh,&d); * nNormalize(d); * nDelete(&hh); * nMult(d,p->coef,&hh); * nDelete(&d); * nDelete(&(p->coef)); * p->coef =hh; */ d=n_Mult(h,pGetCoeff(p),r->cf); n_Normalize(d,r->cf); p_SetCoeff(p,d,r); pIter(p); } number t=n_Mult(c,h,r->cf); n_Delete(&c,r->cf); c=t; } else { break; } n_Delete(&h,r->cf); } } } } if(!n_GreaterZero(pGetCoeff(ph),C)) { ph = p_Neg(ph,r); c = n_Neg(c, C); } } number p_GetAllDenom(poly ph, const ring r) { number d=n_Init(1,r->cf); poly p = ph; while (p!=NULL) { number h=n_GetDenom(pGetCoeff(p),r->cf); if (!n_IsOne(h,r->cf)) { number dd=n_Mult(d,h,r->cf); n_Delete(&d,r->cf); d=dd; } n_Delete(&h,r->cf); pIter(p); } return d; } int p_Size(poly p, const ring r) { int count = 0; while ( p != NULL ) { count+= n_Size( pGetCoeff( p ), r->cf ); pIter( p ); } return count; } /*2 *make p homogeneous by multiplying the monomials by powers of x_varnum *assume: deg(var(varnum))==1 */ poly p_Homogen (poly p, int varnum, const ring r) { pFDegProc deg; if (r->pLexOrder && (r->order[0]==ringorder_lp)) deg=p_Totaldegree; else deg=r->pFDeg; poly q=NULL, qn; int o,ii; sBucket_pt bp; if (p!=NULL) { if ((varnum < 1) || (varnum > rVar(r))) { return NULL; } o=deg(p,r); q=pNext(p); while (q != NULL) { ii=deg(q,r); if (ii>o) o=ii; pIter(q); } q = p_Copy(p,r); bp = sBucketCreate(r); while (q != NULL) { ii = o-deg(q,r); if (ii!=0) { p_AddExp(q,varnum, (long)ii,r); p_Setm(q,r); } qn = pNext(q); pNext(q) = NULL; sBucket_Add_p(bp, q, 1); q = qn; } sBucketDestroyAdd(bp, &q, &ii); } return q; } /*2 *tests if p is homogeneous with respect to the actual weigths */ BOOLEAN p_IsHomogeneous (poly p, const ring r) { poly qp=p; int o; if ((p == NULL) || (pNext(p) == NULL)) return TRUE; pFDegProc d; if (r->pLexOrder && (r->order[0]==ringorder_lp)) d=p_Totaldegree; else d=r->pFDeg; o = d(p,r); do { if (d(qp,r) != o) return FALSE; pIter(qp); } while (qp != NULL); return TRUE; } /*----------utilities for syzygies--------------*/ BOOLEAN p_VectorHasUnitB(poly p, int * k, const ring r) { poly q=p,qq; int i; while (q!=NULL) { if (p_LmIsConstantComp(q,r)) { i = p_GetComp(q,r); qq = p; while ((qq != q) && (p_GetComp(qq,r) != i)) pIter(qq); if (qq == q) { *k = i; return TRUE; } else pIter(q); } else pIter(q); } return FALSE; } void p_VectorHasUnit(poly p, int * k, int * len, const ring r) { poly q=p,qq; int i,j=0; *len = 0; while (q!=NULL) { if (p_LmIsConstantComp(q,r)) { i = p_GetComp(q,r); qq = p; while ((qq != q) && (p_GetComp(qq,r) != i)) pIter(qq); if (qq == q) { j = 0; while (qq!=NULL) { if (p_GetComp(qq,r)==i) j++; pIter(qq); } if ((*len == 0) || (j<*len)) { *len = j; *k = i; } } } pIter(q); } } poly p_TakeOutComp1(poly * p, int k, const ring r) { poly q = *p; if (q==NULL) return NULL; poly qq=NULL,result = NULL; if (p_GetComp(q,r)==k) { result = q; /* *p */ while ((q!=NULL) && (p_GetComp(q,r)==k)) { p_SetComp(q,0,r); p_SetmComp(q,r); qq = q; pIter(q); } *p = q; pNext(qq) = NULL; } if (q==NULL) return result; // if (pGetComp(q) > k) pGetComp(q)--; while (pNext(q)!=NULL) { if (p_GetComp(pNext(q),r)==k) { if (result==NULL) { result = pNext(q); qq = result; } else { pNext(qq) = pNext(q); pIter(qq); } pNext(q) = pNext(pNext(q)); pNext(qq) =NULL; p_SetComp(qq,0,r); p_SetmComp(qq,r); } else { pIter(q); // if (pGetComp(q) > k) pGetComp(q)--; } } return result; } poly p_TakeOutComp(poly * p, int k, const ring r) { poly q = *p,qq=NULL,result = NULL; if (q==NULL) return NULL; BOOLEAN use_setmcomp=rOrd_SetCompRequiresSetm(r); if (p_GetComp(q,r)==k) { result = q; do { p_SetComp(q,0,r); if (use_setmcomp) p_SetmComp(q,r); qq = q; pIter(q); } while ((q!=NULL) && (p_GetComp(q,r)==k)); *p = q; pNext(qq) = NULL; } if (q==NULL) return result; if (p_GetComp(q,r) > k) { p_SubComp(q,1,r); if (use_setmcomp) p_SetmComp(q,r); } poly pNext_q; while ((pNext_q=pNext(q))!=NULL) { if (p_GetComp(pNext_q,r)==k) { if (result==NULL) { result = pNext_q; qq = result; } else { pNext(qq) = pNext_q; pIter(qq); } pNext(q) = pNext(pNext_q); pNext(qq) =NULL; p_SetComp(qq,0,r); if (use_setmcomp) p_SetmComp(qq,r); } else { /*pIter(q);*/ q=pNext_q; if (p_GetComp(q,r) > k) { p_SubComp(q,1,r); if (use_setmcomp) p_SetmComp(q,r); } } } return result; } // Splits *p into two polys: *q which consists of all monoms with // component == comp and *p of all other monoms *lq == pLength(*q) void p_TakeOutComp(poly *r_p, long comp, poly *r_q, int *lq, const ring r) { spolyrec pp, qq; poly p, q, p_prev; int l = 0; #ifdef HAVE_ASSUME int lp = pLength(*r_p); #endif pNext(&pp) = *r_p; p = *r_p; p_prev = &pp; q = &qq; while(p != NULL) { while (p_GetComp(p,r) == comp) { pNext(q) = p; pIter(q); p_SetComp(p, 0,r); p_SetmComp(p,r); pIter(p); l++; if (p == NULL) { pNext(p_prev) = NULL; goto Finish; } } pNext(p_prev) = p; p_prev = p; pIter(p); } Finish: pNext(q) = NULL; *r_p = pNext(&pp); *r_q = pNext(&qq); *lq = l; #ifdef HAVE_ASSUME assume(pLength(*r_p) + pLength(*r_q) == lp); #endif p_Test(*r_p,r); p_Test(*r_q,r); } void p_DeleteComp(poly * p,int k, const ring r) { poly q; while ((*p!=NULL) && (p_GetComp(*p,r)==k)) p_LmDelete(p,r); if (*p==NULL) return; q = *p; if (p_GetComp(q,r)>k) { p_SubComp(q,1,r); p_SetmComp(q,r); } while (pNext(q)!=NULL) { if (p_GetComp(pNext(q),r)==k) p_LmDelete(&(pNext(q)),r); else { pIter(q); if (p_GetComp(q,r)>k) { p_SubComp(q,1,r); p_SetmComp(q,r); } } } } /*2 * convert a vector to a set of polys, * allocates the polyset, (entries 0..(*len)-1) * the vector will not be changed */ void p_Vec2Polys(poly v, poly* *p, int *len, const ring r) { poly h; int k; *len=p_MaxComp(v,r); if (*len==0) *len=1; *p=(poly*)omAlloc0((*len)*sizeof(poly)); while (v!=NULL) { h=p_Head(v,r); k=p_GetComp(h,r); p_SetComp(h,0,r); (*p)[k-1]=p_Add_q((*p)[k-1],h,r); pIter(v); } } // // resets the pFDeg and pLDeg: if pLDeg is not given, it is // set to currRing->pLDegOrig, i.e. to the respective LDegProc which // only uses pFDeg (and not p_Deg, or pTotalDegree, etc) void pSetDegProcs(ring r, pFDegProc new_FDeg, pLDegProc new_lDeg) { assume(new_FDeg != NULL); r->pFDeg = new_FDeg; if (new_lDeg == NULL) new_lDeg = r->pLDegOrig; r->pLDeg = new_lDeg; } // restores pFDeg and pLDeg: void pRestoreDegProcs(ring r, pFDegProc old_FDeg, pLDegProc old_lDeg) { assume(old_FDeg != NULL && old_lDeg != NULL); r->pFDeg = old_FDeg; r->pLDeg = old_lDeg; } /*-------- several access procedures to monomials -------------------- */ /* * the module weights for std */ static pFDegProc pOldFDeg; static pLDegProc pOldLDeg; static BOOLEAN pOldLexOrder; static long pModDeg(poly p, ring r) { long d=pOldFDeg(p, r); int c=p_GetComp(p, r); if ((c>0) && ((r->pModW)->range(c-1))) d+= (*(r->pModW))[c-1]; return d; //return pOldFDeg(p, r)+(*pModW)[p_GetComp(p, r)-1]; } void p_SetModDeg(intvec *w, ring r) { if (w!=NULL) { r->pModW = w; pOldFDeg = r->pFDeg; pOldLDeg = r->pLDeg; pOldLexOrder = r->pLexOrder; pSetDegProcs(r,pModDeg); r->pLexOrder = TRUE; } else { r->pModW = NULL; pRestoreDegProcs(r,pOldFDeg, pOldLDeg); r->pLexOrder = pOldLexOrder; } } /*2 * handle memory request for sets of polynomials (ideals) * l is the length of *p, increment is the difference (may be negative) */ void pEnlargeSet(poly* *p, int l, int increment) { poly* h; h=(poly*)omReallocSize((poly*)*p,l*sizeof(poly),(l+increment)*sizeof(poly)); if (increment>0) { //for (i=l; icf)) return; // Werror("p_Norm not possible in the case of coefficient rings."); } else #endif if (p1!=NULL) { if (pNext(p1)==NULL) { p_SetCoeff(p1,n_Init(1,r->cf),r); return; } poly h; if (!n_IsOne(pGetCoeff(p1),r->cf)) { number k, c; n_Normalize(pGetCoeff(p1),r->cf); k = pGetCoeff(p1); c = n_Init(1,r->cf); pSetCoeff0(p1,c); h = pNext(p1); while (h!=NULL) { c=n_Div(pGetCoeff(h),k,r->cf); // no need to normalize: Z/p, R // normalize already in nDiv: Q_a, Z/p_a // remains: Q if (rField_is_Q(r) && (!n_IsOne(c,r->cf))) n_Normalize(c,r->cf); p_SetCoeff(h,c,r); pIter(h); } n_Delete(&k,r->cf); } else { //if (r->cf->cfNormalize != nDummy2) //TODO: OPTIMIZE { h = pNext(p1); while (h!=NULL) { n_Normalize(pGetCoeff(h),r->cf); pIter(h); } } } } } /*2 *normalize all coefficients */ void p_Normalize(poly p,const ring r) { if (rField_has_simple_inverse(r)) return; /* Z/p, GF(p,n), R, long R/C */ while (p!=NULL) { #ifdef LDEBUG n_Test(pGetCoeff(p), r->cf); #endif n_Normalize(pGetCoeff(p),r->cf); pIter(p); } } // splits p into polys with Exp(n) == 0 and Exp(n) != 0 // Poly with Exp(n) != 0 is reversed static void p_SplitAndReversePoly(poly p, int n, poly *non_zero, poly *zero, const ring r) { if (p == NULL) { *non_zero = NULL; *zero = NULL; return; } spolyrec sz; poly z, n_z, next; z = &sz; n_z = NULL; while(p != NULL) { next = pNext(p); if (p_GetExp(p, n,r) == 0) { pNext(z) = p; pIter(z); } else { pNext(p) = n_z; n_z = p; } p = next; } pNext(z) = NULL; *zero = pNext(&sz); *non_zero = n_z; } /*3 * substitute the n-th variable by 1 in p * destroy p */ static poly p_Subst1 (poly p,int n, const ring r) { poly qq=NULL, result = NULL; poly zero=NULL, non_zero=NULL; // reverse, so that add is likely to be linear p_SplitAndReversePoly(p, n, &non_zero, &zero,r); while (non_zero != NULL) { assume(p_GetExp(non_zero, n,r) != 0); qq = non_zero; pIter(non_zero); qq->next = NULL; p_SetExp(qq,n,0,r); p_Setm(qq,r); result = p_Add_q(result,qq,r); } p = p_Add_q(result, zero,r); p_Test(p,r); return p; } /*3 * substitute the n-th variable by number e in p * destroy p */ static poly p_Subst2 (poly p,int n, number e, const ring r) { assume( ! n_IsZero(e,r->cf) ); poly qq,result = NULL; number nn, nm; poly zero, non_zero; // reverse, so that add is likely to be linear p_SplitAndReversePoly(p, n, &non_zero, &zero,r); while (non_zero != NULL) { assume(p_GetExp(non_zero, n, r) != 0); qq = non_zero; pIter(non_zero); qq->next = NULL; n_Power(e, p_GetExp(qq, n, r), &nn,r->cf); nm = n_Mult(nn, pGetCoeff(qq),r->cf); #ifdef HAVE_RINGS if (n_IsZero(nm,r->cf)) { p_LmFree(&qq,r); n_Delete(&nm,r->cf); } else #endif { p_SetCoeff(qq, nm,r); p_SetExp(qq, n, 0,r); p_Setm(qq,r); result = p_Add_q(result,qq,r); } n_Delete(&nn,r->cf); } p = p_Add_q(result, zero,r); p_Test(p,r); return p; } /* delete monoms whose n-th exponent is different from zero */ static poly p_Subst0(poly p, int n, const ring r) { spolyrec res; poly h = &res; pNext(h) = p; while (pNext(h)!=NULL) { if (p_GetExp(pNext(h),n,r)!=0) { p_LmDelete(&pNext(h),r); } else { pIter(h); } } p_Test(pNext(&res),r); return pNext(&res); } /*2 * substitute the n-th variable by e in p * destroy p */ poly p_Subst(poly p, int n, poly e, const ring r) { if (e == NULL) return p_Subst0(p, n,r); if (p_IsConstant(e,r)) { if (n_IsOne(pGetCoeff(e),r->cf)) return p_Subst1(p,n,r); else return p_Subst2(p, n, pGetCoeff(e),r); } #ifdef HAVE_PLURAL if (rIsPluralRing(r)) { return nc_pSubst(p,n,e,r); } #endif int exponent,i; poly h, res, m; int *me,*ee; number nu,nu1; me=(int *)omAlloc((rVar(r)+1)*sizeof(int)); ee=(int *)omAlloc((rVar(r)+1)*sizeof(int)); if (e!=NULL) p_GetExpV(e,ee,r); res=NULL; h=p; while (h!=NULL) { if ((e!=NULL) || (p_GetExp(h,n,r)==0)) { m=p_Head(h,r); p_GetExpV(m,me,r); exponent=me[n]; me[n]=0; for(i=rVar(r);i>0;i--) me[i]+=exponent*ee[i]; p_SetExpV(m,me,r); if (e!=NULL) { n_Power(pGetCoeff(e),exponent,&nu,r->cf); nu1=n_Mult(pGetCoeff(m),nu,r->cf); n_Delete(&nu,r->cf); p_SetCoeff(m,nu1,r); } res=p_Add_q(res,m,r); } p_LmDelete(&h,r); } omFreeSize((ADDRESS)me,(rVar(r)+1)*sizeof(int)); omFreeSize((ADDRESS)ee,(rVar(r)+1)*sizeof(int)); return res; } /*2 * returns a re-ordered convertion of a number as a polynomial, * with permutation of parameters * NOTE: this only works for Frank's alg. & trans. fields */ poly n_PermNumber(const number z, const int *par_perm, const int , const ring src, const ring dst) { #if 0 PrintS("\nSource Ring: \n"); rWrite(src); if(0) { number zz = n_Copy(z, src->cf); PrintS("z: "); n_Write(zz, src); n_Delete(&zz, src->cf); } PrintS("\nDestination Ring: \n"); rWrite(dst); /*Print("\nOldPar: %d\n", OldPar); for( int i = 1; i <= OldPar; i++ ) { Print("par(%d) -> par/var (%d)\n", i, par_perm[i-1]); }*/ #endif if( z == NULL ) return NULL; const coeffs srcCf = src->cf; assume( srcCf != NULL ); assume( !nCoeff_is_GF(srcCf) ); assume( rField_is_Extension(src) ); poly zz = NULL; const ring srcExtRing = srcCf->extRing; assume( srcExtRing != NULL ); const coeffs dstCf = dst->cf; assume( dstCf != NULL ); if( nCoeff_is_algExt(srcCf) ) // nCoeff_is_GF(srcCf)? { zz = (poly) z; if( zz == NULL ) return NULL; } else if (nCoeff_is_transExt(srcCf)) { assume( !IS0(z) ); zz = NUM(z); p_Test (zz, srcExtRing); if( zz == NULL ) return NULL; //if( !DENIS1(z) ) //WarnS("Not implemented yet: Cannot permute a rational fraction and make a polynomial out of it! Ignorring the denumerator."); } else { assume (FALSE); Werror("Number permutation is not implemented for this data yet!"); return NULL; } assume( zz != NULL ); p_Test (zz, srcExtRing); nMapFunc nMap = n_SetMap(srcExtRing->cf, dstCf); assume( nMap != NULL ); poly qq; if ((par_perm == NULL) && (rPar(dst) != 0 && rVar (srcExtRing) > 0)) { int* perm; perm=(int *)omAlloc0((rVar(srcExtRing)+1)*sizeof(int)); perm[0]= 0; for(int i=si_min(rVar(srcExtRing),rPar(dst));i>0;i--) perm[i]=-i; qq = p_PermPoly(zz, perm, srcExtRing, dst, nMap, NULL, rVar(srcExtRing)-1); omFreeSize ((ADDRESS)perm, (rVar(srcExtRing)+1)*sizeof(int)); } else qq = p_PermPoly(zz, par_perm-1, srcExtRing, dst, nMap, NULL, rVar (srcExtRing)-1); assume (p_Test (qq, dst)); // poly p_PermPoly (poly p, int * perm, const ring oldRing, const ring dst, nMapFunc nMap, int *par_perm, int OldPar) // assume( FALSE ); WarnS("longalg missing 2"); return qq; } /*2 *returns a re-ordered copy of a polynomial, with permutation of the variables */ poly p_PermPoly (poly p, const int * perm, const ring oldRing, const ring dst, nMapFunc nMap, const int *par_perm, int OldPar) { #if 0 p_Test(p, oldRing); PrintS("\np_PermPoly::p: "); p_Write(p, oldRing, oldRing); PrintLn(); #endif const int OldpVariables = rVar(oldRing); poly result = NULL; poly result_last = NULL; poly aq = NULL; /* the map coefficient */ poly qq; /* the mapped monomial */ assume(dst != NULL); assume(dst->cf != NULL); while (p != NULL) { // map the coefficient if ( ((OldPar == 0) || (par_perm == NULL) || rField_is_GF(oldRing)) && (nMap != NULL) ) { qq = p_Init(dst); assume( nMap != NULL ); number n = nMap(p_GetCoeff(p, oldRing), oldRing->cf, dst->cf); assume (n_Test (n,dst->cf)); if ( nCoeff_is_algExt(dst->cf) ) n_Normalize(n, dst->cf); p_GetCoeff(qq, dst) = n;// Note: n can be a ZERO!!! // coef may be zero: // p_Test(qq, dst); } else { qq = p_One(dst); // aq = naPermNumber(p_GetCoeff(p, oldRing), par_perm, OldPar, oldRing); // no dst??? // poly n_PermNumber(const number z, const int *par_perm, const int P, const ring src, const ring dst) aq = n_PermNumber(p_GetCoeff(p, oldRing), par_perm, OldPar, oldRing, dst); p_Test(aq, dst); if ( nCoeff_is_algExt(dst->cf) ) p_Normalize(aq,dst); if (aq == NULL) p_SetCoeff(qq, n_Init(0, dst->cf),dst); // Very dirty trick!!! p_Test(aq, dst); } if (rRing_has_Comp(dst)) p_SetComp(qq, p_GetComp(p, oldRing), dst); if ( n_IsZero(pGetCoeff(qq), dst->cf) ) { p_LmDelete(&qq,dst); qq = NULL; } else { // map pars: int mapped_to_par = 0; for(int i = 1; i <= OldpVariables; i++) { int e = p_GetExp(p, i, oldRing); if (e != 0) { if (perm==NULL) p_SetExp(qq, i, e, dst); else if (perm[i]>0) p_AddExp(qq,perm[i], e/*p_GetExp( p,i,oldRing)*/, dst); else if (perm[i]<0) { number c = p_GetCoeff(qq, dst); if (rField_is_GF(dst)) { assume( dst->cf->extRing == NULL ); number ee = n_Param(1, dst); number eee; n_Power(ee, e, &eee, dst->cf); //nfDelete(ee,dst); ee = n_Mult(c, eee, dst->cf); //nfDelete(c,dst);nfDelete(eee,dst); pSetCoeff0(qq,ee); } else if (nCoeff_is_Extension(dst->cf)) { const int par = -perm[i]; assume( par > 0 ); // WarnS("longalg missing 3"); #if 1 const coeffs C = dst->cf; assume( C != NULL ); const ring R = C->extRing; assume( R != NULL ); assume( par <= rVar(R) ); poly pcn; // = (number)c assume( !n_IsZero(c, C) ); if( nCoeff_is_algExt(C) ) pcn = (poly) c; else // nCoeff_is_transExt(C) pcn = NUM(c); if (pNext(pcn) == NULL) // c->z p_AddExp(pcn, -perm[i], e, R); else /* more difficult: we have really to multiply: */ { poly mmc = p_ISet(1, R); p_SetExp(mmc, -perm[i], e, R); p_Setm(mmc, R); number nnc; // convert back to a number: number nnc = mmc; if( nCoeff_is_algExt(C) ) nnc = (number) mmc; else // nCoeff_is_transExt(C) nnc = ntInit(mmc, C); p_GetCoeff(qq, dst) = n_Mult((number)c, nnc, C); n_Delete((number *)&c, C); n_Delete((number *)&nnc, C); } mapped_to_par=1; #endif } } else { /* this variable maps to 0 !*/ p_LmDelete(&qq, dst); break; } } } if ( mapped_to_par && qq!= NULL && nCoeff_is_algExt(dst->cf) ) { number n = p_GetCoeff(qq, dst); n_Normalize(n, dst->cf); p_GetCoeff(qq, dst) = n; } } pIter(p); #if 0 p_Test(aq,dst); PrintS("\naq: "); p_Write(aq, dst, dst); PrintLn(); #endif #if 1 if (qq!=NULL) { p_Setm(qq,dst); p_Test(aq,dst); p_Test(qq,dst); #if 0 p_Test(qq,dst); PrintS("\nqq: "); p_Write(qq, dst, dst); PrintLn(); #endif if (aq!=NULL) qq=p_Mult_q(aq,qq,dst); aq = qq; while (pNext(aq) != NULL) pIter(aq); if (result_last==NULL) { result=qq; } else { pNext(result_last)=qq; } result_last=aq; aq = NULL; } else if (aq!=NULL) { p_Delete(&aq,dst); } } result=p_SortAdd(result,dst); #else // if (qq!=NULL) // { // pSetm(qq); // pTest(qq); // pTest(aq); // if (aq!=NULL) qq=pMult(aq,qq); // aq = qq; // while (pNext(aq) != NULL) pIter(aq); // pNext(aq) = result; // aq = NULL; // result = qq; // } // else if (aq!=NULL) // { // pDelete(&aq); // } //} //p = result; //result = NULL; //while (p != NULL) //{ // qq = p; // pIter(p); // qq->next = NULL; // result = pAdd(result, qq); //} #endif p_Test(result,dst); #if 0 p_Test(result,dst); PrintS("\nresult: "); p_Write(result,dst,dst); PrintLn(); #endif return result; } /************************************************************** * * Jet * **************************************************************/ poly pp_Jet(poly p, int m, const ring R) { poly r=NULL; poly t=NULL; while (p!=NULL) { if (p_Totaldegree(p,R)<=m) { if (r==NULL) r=p_Head(p,R); else if (t==NULL) { pNext(r)=p_Head(p,R); t=pNext(r); } else { pNext(t)=p_Head(p,R); pIter(t); } } pIter(p); } return r; } poly p_Jet(poly p, int m,const ring R) { while((p!=NULL) && (p_Totaldegree(p,R)>m)) p_LmDelete(&p,R); if (p==NULL) return NULL; poly r=p; while (pNext(p)!=NULL) { if (p_Totaldegree(pNext(p),R)>m) { p_LmDelete(&pNext(p),R); } else pIter(p); } return r; } poly pp_JetW(poly p, int m, short *w, const ring R) { poly r=NULL; poly t=NULL; while (p!=NULL) { if (totaldegreeWecart_IV(p,R,w)<=m) { if (r==NULL) r=p_Head(p,R); else if (t==NULL) { pNext(r)=p_Head(p,R); t=pNext(r); } else { pNext(t)=p_Head(p,R); pIter(t); } } pIter(p); } return r; } poly p_JetW(poly p, int m, short *w, const ring R) { while((p!=NULL) && (totaldegreeWecart_IV(p,R,w)>m)) p_LmDelete(&p,R); if (p==NULL) return NULL; poly r=p; while (pNext(p)!=NULL) { if (totaldegreeWecart_IV(pNext(p),R,w)>m) { p_LmDelete(&pNext(p),R); } else pIter(p); } return r; } /*************************************************************/ int p_MinDeg(poly p,intvec *w, const ring R) { if(p==NULL) return -1; int d=-1; while(p!=NULL) { int d0=0; for(int j=0;j=w->length()) d0+=p_GetExp(p,j+1,R); else d0+=(*w)[j]*p_GetExp(p,j+1,R); if(d0cf); poly v=p_NSet(u0,R); if(n==0) return v; short *ww=iv2array(w,R); poly u1=p_JetW(p_Sub(p_One(R),p_Mult_nn(u,u0,R),R),n,ww,R); if(u1==NULL) { omFreeSize((ADDRESS)ww,(rVar(R)+1)*sizeof(short)); return v; } poly v1=p_Mult_nn(p_Copy(u1,R),u0,R); v=p_Add_q(v,p_Copy(v1,R),R); for(int i=n/p_MinDeg(u1,w,R);i>1;i--) { v1=p_JetW(p_Mult_q(v1,p_Copy(u1,R),R),n,ww,R); v=p_Add_q(v,p_Copy(v1,R),R); } p_Delete(&u1,R); p_Delete(&v1,R); omFreeSize((ADDRESS)ww,(rVar(R)+1)*sizeof(short)); return v; } BOOLEAN p_EqualPolys(poly p1,poly p2, const ring r) { while ((p1 != NULL) && (p2 != NULL)) { if (! p_LmEqual(p1, p2,r)) return FALSE; if (! n_Equal(p_GetCoeff(p1,r), p_GetCoeff(p2,r),r->cf )) return FALSE; pIter(p1); pIter(p2); } return (p1==p2); } static inline BOOLEAN p_ExpVectorEqual(poly p1, poly p2, const ring r1, const ring r2) { assume( r1 == r2 || rSamePolyRep(r1, r2) ); p_LmCheckPolyRing1(p1, r1); p_LmCheckPolyRing1(p2, r2); int i = r1->ExpL_Size; assume( r1->ExpL_Size == r2->ExpL_Size ); unsigned long *ep = p1->exp; unsigned long *eq = p2->exp; do { i--; if (ep[i] != eq[i]) return FALSE; } while (i); return TRUE; } BOOLEAN p_EqualPolys(poly p1,poly p2, const ring r1, const ring r2) { assume( r1 == r2 || rSamePolyRep(r1, r2) ); // will be used in rEqual! assume( r1->cf == r2->cf ); while ((p1 != NULL) && (p2 != NULL)) { // returns 1 if ExpVector(p)==ExpVector(q): does not compare numbers !! // #define p_LmEqual(p1, p2, r) p_ExpVectorEqual(p1, p2, r) if (! p_ExpVectorEqual(p1, p2, r1, r2)) return FALSE; if (! n_Equal(p_GetCoeff(p1,r1), p_GetCoeff(p2,r2), r1->cf )) return FALSE; pIter(p1); pIter(p2); } return (p1==p2); } /*2 *returns TRUE if p1 is a skalar multiple of p2 *assume p1 != NULL and p2 != NULL */ BOOLEAN p_ComparePolys(poly p1,poly p2, const ring r) { number n,nn; pAssume(p1 != NULL && p2 != NULL); if (!p_LmEqual(p1,p2,r)) //compare leading mons return FALSE; if ((pNext(p1)==NULL) && (pNext(p2)!=NULL)) return FALSE; if ((pNext(p2)==NULL) && (pNext(p1)!=NULL)) return FALSE; if (pLength(p1) != pLength(p2)) return FALSE; #ifdef HAVE_RINGS if (rField_is_Ring(r)) { if (!n_DivBy(p_GetCoeff(p1, r), p_GetCoeff(p2, r), r)) return FALSE; } #endif n=n_Div(p_GetCoeff(p1,r),p_GetCoeff(p2,r),r); while ((p1 != NULL) /*&& (p2 != NULL)*/) { if ( ! p_LmEqual(p1, p2,r)) { n_Delete(&n, r); return FALSE; } if (!n_Equal(p_GetCoeff(p1, r), nn = n_Mult(p_GetCoeff(p2, r),n, r->cf), r->cf)) { n_Delete(&n, r); n_Delete(&nn, r); return FALSE; } n_Delete(&nn, r); pIter(p1); pIter(p2); } n_Delete(&n, r); return TRUE; } /*2 * returns the length of a (numbers of monomials) * respect syzComp */ poly p_Last(const poly p, int &l, const ring r) { if (p == NULL) { l = 0; return NULL; } l = 1; poly a = p; if (! rIsSyzIndexRing(r)) { poly next = pNext(a); while (next!=NULL) { a = next; next = pNext(a); l++; } } else { int curr_limit = rGetCurrSyzLimit(r); poly pp = a; while ((a=pNext(a))!=NULL) { if (p_GetComp(a,r)<=curr_limit/*syzComp*/) l++; else break; pp = a; } a=pp; } return a; } int p_Var(poly m,const ring r) { if (m==NULL) return 0; if (pNext(m)!=NULL) return 0; int i,e=0; for (i=rVar(r); i>0; i--) { int exp=p_GetExp(m,i,r); if (exp==1) { if (e==0) e=i; else return 0; } else if (exp!=0) { return 0; } } return e; } /*2 *the minimal index of used variables - 1 */ int p_LowVar (poly p, const ring r) { int k,l,lex; if (p == NULL) return -1; k = 32000;/*a very large dummy value*/ while (p != NULL) { l = 1; lex = p_GetExp(p,l,r); while ((l < (rVar(r))) && (lex == 0)) { l++; lex = p_GetExp(p,l,r); } l--; if (l < k) k = l; pIter(p); } return k; } /*2 * verschiebt die Indizees der Modulerzeugenden um i */ void p_Shift (poly * p,int i, const ring r) { poly qp1 = *p,qp2 = *p;/*working pointers*/ int j = p_MaxComp(*p,r),k = p_MinComp(*p,r); if (j+i < 0) return ; while (qp1 != NULL) { if ((p_GetComp(qp1,r)+i > 0) || ((j == -i) && (j == k))) { p_AddComp(qp1,i,r); p_SetmComp(qp1,r); qp2 = qp1; pIter(qp1); } else { if (qp2 == *p) { pIter(*p); p_LmDelete(&qp2,r); qp2 = *p; qp1 = *p; } else { qp2->next = qp1->next; if (qp1!=NULL) p_LmDelete(&qp1,r); qp1 = qp2->next; } } } } /*************************************************************** * * Storage Managament Routines * ***************************************************************/ static inline unsigned long GetBitFields(long e, unsigned int s, unsigned int n) { #define Sy_bit_L(x) (((unsigned long)1L)<<(x)) unsigned int i = 0; unsigned long ev = 0L; assume(n > 0 && s < BIT_SIZEOF_LONG); do { assume(s+i < BIT_SIZEOF_LONG); if (e > (long) i) ev |= Sy_bit_L(s+i); else break; i++; } while (i < n); return ev; } // Short Exponent Vectors are used for fast divisibility tests // ShortExpVectors "squeeze" an exponent vector into one word as follows: // Let n = BIT_SIZEOF_LONG / pVariables. // If n == 0 (i.e. pVariables > BIT_SIZE_OF_LONG), let m == the number // of non-zero exponents. If (m>BIT_SIZEOF_LONG), then sev = ~0, else // first m bits of sev are set to 1. // Otherwise (i.e. pVariables <= BIT_SIZE_OF_LONG) // represented by a bit-field of length n (resp. n+1 for some // exponents). If the value of an exponent is greater or equal to n, then // all of its respective n bits are set to 1. If the value of an exponent // is smaller than n, say m, then only the first m bits of the respective // n bits are set to 1, the others are set to 0. // This way, we have: // exp1 / exp2 ==> (ev1 & ~ev2) == 0, i.e., // if (ev1 & ~ev2) then exp1 does not divide exp2 unsigned long p_GetShortExpVector(poly p, const ring r) { assume(p != NULL); if (p == NULL) return 0; unsigned long ev = 0; // short exponent vector unsigned int n = BIT_SIZEOF_LONG / r->N; // number of bits per exp unsigned int m1; // highest bit which is filled with (n+1) unsigned int i = 0, j=1; if (n == 0) { if (r->N <2*BIT_SIZEOF_LONG) { n=1; m1=0; } else { for (; j<=(unsigned long) r->N; j++) { if (p_GetExp(p,j,r) > 0) i++; if (i == BIT_SIZEOF_LONG) break; } if (i>0) ev = ~((unsigned long)0) >> ((unsigned long) (BIT_SIZEOF_LONG - i)); return ev; } } else { m1 = (n+1)*(BIT_SIZEOF_LONG - n*r->N); } n++; while (i