| 1 | #ifndef POLYS_H |
|---|
| 2 | #define POLYS_H |
|---|
| 3 | /**************************************** |
|---|
| 4 | * Computer Algebra System SINGULAR * |
|---|
| 5 | ****************************************/ |
|---|
| 6 | /* $Id: polys.h,v 1.32 2000-08-14 12:56:46 obachman Exp $ */ |
|---|
| 7 | /* |
|---|
| 8 | * ABSTRACT - all basic methods to manipulate polynomials |
|---|
| 9 | */ |
|---|
| 10 | #include "polys-comp.h" |
|---|
| 11 | #include "polys-impl.h" |
|---|
| 12 | #define pFetchCopy(r,p) _pFetchCopy(r,p) |
|---|
| 13 | // Similar to pFetchCopy, except that poly p is deleted |
|---|
| 14 | #define pFetchCopyDelete(r, p) _pFetchCopyDelete(r, p) |
|---|
| 15 | |
|---|
| 16 | typedef poly* polyset; |
|---|
| 17 | extern int pVariables; |
|---|
| 18 | extern int pOrdSgn; |
|---|
| 19 | extern BOOLEAN pLexOrder; |
|---|
| 20 | extern BOOLEAN pMixedOrder; |
|---|
| 21 | extern poly ppNoether; |
|---|
| 22 | extern BOOLEAN pVectorOut; |
|---|
| 23 | // 1 for lex ordering (except ls), -1 otherwise |
|---|
| 24 | extern int pComponentOrder; |
|---|
| 25 | |
|---|
| 26 | /* function prototypes */ |
|---|
| 27 | |
|---|
| 28 | /* ---------------- General poly access and iteration macros -----------*/ |
|---|
| 29 | // IMPORTANT: Do _NOT_ make _ANY_ assumptions about their actual implementation |
|---|
| 30 | |
|---|
| 31 | #define pNext(p) _pNext(p) |
|---|
| 32 | #define pIter(p) _pIter(p) |
|---|
| 33 | |
|---|
| 34 | #define pGetCoeff(p) _pGetCoeff(p) |
|---|
| 35 | // deletes old coeff before setting the new |
|---|
| 36 | #define pSetCoeff(p,n) _pSetCoeff(p,n) |
|---|
| 37 | // sets new coeff without deleting the old one |
|---|
| 38 | #define pSetCoeff0(p,n) _pSetCoeff0(p,n) |
|---|
| 39 | |
|---|
| 40 | #define pGetOrder(p) _pGetOrder(p) |
|---|
| 41 | |
|---|
| 42 | // Gets/Sets Component field w.r.t. of polys of global ring |
|---|
| 43 | #define pSetComp(p,k) _pSetComp(p,k) |
|---|
| 44 | #define pGetComp(p) _pGetComp(p) |
|---|
| 45 | #define pIncrComp(p) _pIncrComp(p) |
|---|
| 46 | #define pDecrComp(p) _pDecrComp(p) |
|---|
| 47 | #define pAddComp(p,v) _pAddComp(p,v) |
|---|
| 48 | #define pSubComp(p,v) _pSubComp(p,v) |
|---|
| 49 | |
|---|
| 50 | // Gets/Sets Component field w.r.t. of polys of ring r |
|---|
| 51 | #define pRingSetComp(r,p,k) _pRingSetComp(r,p,k) |
|---|
| 52 | #define pRingGetComp(r,p) _pRingGetComp(r,p) |
|---|
| 53 | |
|---|
| 54 | // Gets/Sets ith exponent poly of global ring |
|---|
| 55 | #define pGetExp(p,i) _pGetExp(p,i) |
|---|
| 56 | #define pSetExp(p,i,v) _pSetExp(p,i,v) |
|---|
| 57 | |
|---|
| 58 | // Gets/Sets ith exponent of poly of ring r |
|---|
| 59 | #define pRingGetExp(r,p,i) _pRingGetExp(r,p,i) |
|---|
| 60 | #define pRingSetExp(r,p,i,v) _pRingSetExp(r,p,i,v) |
|---|
| 61 | |
|---|
| 62 | // Increments/decrements ith exponent by one |
|---|
| 63 | #define pIncrExp(p,i) _pIncrExp(p,i) |
|---|
| 64 | #define pDecrExp(p,i) _pDecrExp(p,i) |
|---|
| 65 | |
|---|
| 66 | // Gets Difference and sum of ith Exponent of m1, m2 |
|---|
| 67 | #define pGetExpSum(m1, m2, i) _pGetExpSum(p1, p2, i) |
|---|
| 68 | #define pGetExpDiff(p1, p2, i) _pGetExpDiff(p1, p2, i) |
|---|
| 69 | |
|---|
| 70 | // Adds/Substracts v to/from the ith Exponent |
|---|
| 71 | #define pAddExp(p,i,v) _pAddExp(p,i,v) |
|---|
| 72 | #define pSubExp(p,i,v) _pSubExp(p,i,v) |
|---|
| 73 | #define pMultExp(p,i,v) _pMultExp(p,i,v) |
|---|
| 74 | |
|---|
| 75 | // Gets a copy of (resp. sets) the exponent vector, where e is assumed |
|---|
| 76 | // to point to (pVariables+1)*sizeof(Exponent_t) memory. Exponents are |
|---|
| 77 | // filled in as follows: comp, e_1, .., e_n |
|---|
| 78 | #define pGetExpV(p, e) _pGetExpV(p, e) |
|---|
| 79 | #define pSetExpV(p, e) _pSetExpV(p, e) |
|---|
| 80 | |
|---|
| 81 | |
|---|
| 82 | /*-------------predicate on polys ----------------------*/ |
|---|
| 83 | BOOLEAN pIsConstant(const poly p); |
|---|
| 84 | BOOLEAN pIsConstantPoly(poly p); |
|---|
| 85 | BOOLEAN pIsConstantComp(const poly p); |
|---|
| 86 | int pIsPurePower(const poly p); |
|---|
| 87 | #define pIsVector(p) (pGetComp(p)!=0) |
|---|
| 88 | BOOLEAN pHasNotCF(poly p1, poly p2); /*has no common factor ?*/ |
|---|
| 89 | void pSplit(poly p, poly * r); /*p => IN(p), r => REST(p) */ |
|---|
| 90 | // Returns TRUE if m is monom of p, FALSE otherwise |
|---|
| 91 | BOOLEAN pIsMonomOf(poly p, poly m); |
|---|
| 92 | |
|---|
| 93 | |
|---|
| 94 | /*-------------ring management:----------------------*/ |
|---|
| 95 | //extern void pChangeRing(ring newRing); |
|---|
| 96 | extern void pSetGlobals(ring r, BOOLEAN complete = TRUE); |
|---|
| 97 | |
|---|
| 98 | /*-----------the ordering of monomials:-------------*/ |
|---|
| 99 | extern pSetmProc pSetm; |
|---|
| 100 | #define pComp0 rComp0 |
|---|
| 101 | #define pSetmComp pSetm |
|---|
| 102 | int pComp(poly p1,poly p2); |
|---|
| 103 | extern pLDegProc pLDeg; |
|---|
| 104 | extern pFDegProc pFDeg; |
|---|
| 105 | int pDeg(poly p); |
|---|
| 106 | int pTotaldegree(poly p); |
|---|
| 107 | int pWTotaldegree(poly p); |
|---|
| 108 | int pWDegree(poly p); |
|---|
| 109 | |
|---|
| 110 | /*-------------pComp for syzygies:-------------------*/ |
|---|
| 111 | |
|---|
| 112 | void pSetModDeg(intvec *w); |
|---|
| 113 | void pSetSchreyerOrdB(polyset nextorder, int length); |
|---|
| 114 | void pSetSchreyerOrdM(polyset nextorder, int length, int comps); |
|---|
| 115 | int pModuleOrder(); |
|---|
| 116 | |
|---|
| 117 | /*-------------storage management:-------------------*/ |
|---|
| 118 | // allocates the space for a new monomial -- no initialization !!! |
|---|
| 119 | #define pNew() _pNew(currPolyBin) |
|---|
| 120 | // allocates a new monomial and initializes everything to 0 |
|---|
| 121 | #define pInit() _pInit(currPolyBin) |
|---|
| 122 | // frees the space of the monomial m (coef is not freed) |
|---|
| 123 | #define pFree1(m) _pFree1(m, currPolyBin) |
|---|
| 124 | // frees the space of monomial and frees coefficient |
|---|
| 125 | #define pDelete1(m) _pDelete1(m, currPolyBin) |
|---|
| 126 | // deletes the whole polynomial p |
|---|
| 127 | #define pDelete(p) _pDelete(p, currPolyBin) |
|---|
| 128 | |
|---|
| 129 | // similar to routines above, except that monomials are assumed to be from heap h (resp. are allocated from heap h) |
|---|
| 130 | #define pHeapNew(h) _pNew(h) |
|---|
| 131 | #define pHeapInit(h) _pInit(h) |
|---|
| 132 | #define pHeapFree1(m, h) _pFree1(m, h) |
|---|
| 133 | #define pHeapDelete1(m, h) _pDelete1(m, h) |
|---|
| 134 | #define pHeapDelete(p, h) _pDelete(p, h) |
|---|
| 135 | |
|---|
| 136 | // makes a simple memcpy of m2 into m1 -- does _not_ allocate memory for m1 |
|---|
| 137 | #define pCopy2(m1, m2) _pCopy2(m1, m2) |
|---|
| 138 | // Returns a newly allocated copy of the monomial m, initializes coefficient |
|---|
| 139 | // and sets the next-field to NULL |
|---|
| 140 | #define pCopy1(m) _pCopy1(m) |
|---|
| 141 | // Returns a newly allocated copy of the monomial m, sets the coefficient to 0, |
|---|
| 142 | // and sets the next-field to NULL |
|---|
| 143 | #define pHead0(p) _pHead0(p) |
|---|
| 144 | // Returns a newly allocated copy of the monomial, copy includes copy of coeff |
|---|
| 145 | // sets the next-field to NULL |
|---|
| 146 | #define pHead(p) _pHead(currPolyBin, p) |
|---|
| 147 | // Similar to pHead, except that new monom is taken from heap dest_heap |
|---|
| 148 | #define pHeapHead(dest_heap, p) _pHead(dest_heap, p) |
|---|
| 149 | // Returns a shallow copy (i.e. only monomial is copied) of head of source_p |
|---|
| 150 | // Monom of copy is from dest_heap, assumes monoms of p are from source_heap |
|---|
| 151 | // Head of source_p is deleted along the way |
|---|
| 152 | #define pShallowCopyDeleteHead(dest_heap, source_p, source_heap) \ |
|---|
| 153 | _pShallowCopyDeleteHead(dest_heap, source_p, source_heap) |
|---|
| 154 | |
|---|
| 155 | extern poly pHeadProc(poly p); |
|---|
| 156 | // Returns copy of the whole polynomial |
|---|
| 157 | #define pCopy(p) _pCopy(currPolyBin, p) |
|---|
| 158 | // Returns copy of the whole poly, new monomials are taken from dest_heap |
|---|
| 159 | #define pHeapCopy(dest_heap, p) _pCopy(dest_heap, p) |
|---|
| 160 | |
|---|
| 161 | // Returns a shallow copy (i.e. only monomials are copied) of source_p |
|---|
| 162 | // Monoms of copy are from dest_heap, assumes monoms of p are from source_heap |
|---|
| 163 | // source_p is deleted along the way |
|---|
| 164 | #define pShallowCopyDelete(dest_heap, source_p, source_heap) \ |
|---|
| 165 | _pShallowCopyDelete(dest_heap, source_p, source_heap) |
|---|
| 166 | |
|---|
| 167 | |
|---|
| 168 | |
|---|
| 169 | // Adds exponents of p2 to exponents of p1 |
|---|
| 170 | // assumes that exponents >= 0 and and one Component != 0 |
|---|
| 171 | #define pMonAddOn(p1, p2) _pMonAddOn(p1, p2) |
|---|
| 172 | // Similar to pMonAddOn, excpet that new exponents are written in p1 |
|---|
| 173 | #define pMonAdd(p1, p2, p3) _prMonAdd(p1, p2, p3, currRing) |
|---|
| 174 | #define prMonAdd(p1, p2, p3, r) _prMonAdd(p1, p2, p3, r) |
|---|
| 175 | // Subtracts exponetns of p1 from p1, assumes both Components are Equal |
|---|
| 176 | #define pMonSubFrom(p1, p2) _pMonSubFrom(p1, p2) |
|---|
| 177 | |
|---|
| 178 | poly pmInit(char *s, BOOLEAN &ok); /* monom -> poly */ |
|---|
| 179 | void ppDelete(poly * a, ring r); |
|---|
| 180 | |
|---|
| 181 | /*-------------operations on polynomials:------------*/ |
|---|
| 182 | poly pAdd(poly p1, poly p2); |
|---|
| 183 | poly pNeg(poly p); |
|---|
| 184 | poly pSub(poly a, poly b); |
|---|
| 185 | poly pMult(poly a, poly b); |
|---|
| 186 | void pMultN(poly a, number c); |
|---|
| 187 | poly pMultCopyN(poly a, number c); |
|---|
| 188 | poly pPower(poly p, int i); |
|---|
| 189 | |
|---|
| 190 | // return TRUE, if exponent and component of Lm(p1) and Lm(p2) are equal, |
|---|
| 191 | // FALSE otherwise; |
|---|
| 192 | #define pEqual(p1, p2) _pEqual(p1, p2) // Assumes p1 != NULL & p2 != NULL |
|---|
| 193 | inline BOOLEAN pLmEqual(poly p1, poly p2) // no assumptions |
|---|
| 194 | { |
|---|
| 195 | if (p1 != NULL) |
|---|
| 196 | { |
|---|
| 197 | if (p2 != NULL) |
|---|
| 198 | return pEqual(p1, p2); |
|---|
| 199 | return FALSE; |
|---|
| 200 | } |
|---|
| 201 | if (p2 == NULL) |
|---|
| 202 | return TRUE; |
|---|
| 203 | return FALSE; |
|---|
| 204 | } |
|---|
| 205 | |
|---|
| 206 | // returns TRUE, if leading monom of a divides leading monom of b |
|---|
| 207 | // i.e., if there exists a expvector c > 0, s.t. b = a + c; assumes b != NULL |
|---|
| 208 | #define pDivisibleBy(a, b) _pDivisibleBy(a,b) |
|---|
| 209 | // like pDivisibleBy, except that it is assumed that a!=NULL |
|---|
| 210 | #define pDivisibleBy1(a,b) _pDivisibleBy1(a,b) |
|---|
| 211 | // like pDivisibleBy1, does not check components |
|---|
| 212 | #define pDivisibleBy2(a, b) _pDivisibleBy2(a,b) |
|---|
| 213 | |
|---|
| 214 | // Divisibility tests based on Short Exponent vectors |
|---|
| 215 | // sev_a == pGetShortExpVector(a) |
|---|
| 216 | // not_sev_b == ~ pGetShortExpVector(b) |
|---|
| 217 | #define pShortDivisibleBy(a, sev_a, b, not_sev_b) \ |
|---|
| 218 | _pShortDivisibleBy(a, sev_a, b, not_sev_b) |
|---|
| 219 | |
|---|
| 220 | |
|---|
| 221 | poly pDivide(poly a, poly b); |
|---|
| 222 | poly pDivideM(poly a, poly b); |
|---|
| 223 | void pLcm(poly a, poly b, poly m); |
|---|
| 224 | poly pDiff(poly a, int k); |
|---|
| 225 | poly pDiffOp(poly a, poly b,BOOLEAN multiply); |
|---|
| 226 | |
|---|
| 227 | int pLength(poly a); |
|---|
| 228 | int pMaxComp(poly p); |
|---|
| 229 | int pMinComp(poly p, ring r=currRing); |
|---|
| 230 | BOOLEAN pOneComp(poly p); |
|---|
| 231 | int pWeight(int c); |
|---|
| 232 | void pSetCompP(poly a, int i); |
|---|
| 233 | |
|---|
| 234 | // Quer sum (total degree) of all exponents of leading monomial |
|---|
| 235 | #define pExpQuerSum(p1) _pExpQuerSum(p1) |
|---|
| 236 | // sum from "from"th to "to"th exponent (including borders) |
|---|
| 237 | #define pExpQuerSum2(p1, from, to) _pExpQuerSum2(p1, from, to) |
|---|
| 238 | |
|---|
| 239 | char* pString(poly p); |
|---|
| 240 | char* pString0(poly p); |
|---|
| 241 | void pWrite(poly p); |
|---|
| 242 | void pWrite0(poly p); |
|---|
| 243 | void wrp(poly p); |
|---|
| 244 | |
|---|
| 245 | void pEnlargeSet(polyset *p, int length, int increment); |
|---|
| 246 | poly pISet(int i); |
|---|
| 247 | #define pOne() pISet(1) |
|---|
| 248 | |
|---|
| 249 | void pContent(poly p); |
|---|
| 250 | void pCleardenom(poly p); |
|---|
| 251 | void pNormalize(poly p); |
|---|
| 252 | |
|---|
| 253 | // homogenizes p by multiplying certain powers of the varnum-th variable |
|---|
| 254 | poly pHomogen (poly p, int varnum); |
|---|
| 255 | |
|---|
| 256 | // replaces the maximal powers of the leading monomial of p2 in p1 by |
|---|
| 257 | // the same powers of n, utility for dehomogenization |
|---|
| 258 | poly pDehomogen (poly p1,poly p2,number n); |
|---|
| 259 | BOOLEAN pIsHomogeneous (poly p); |
|---|
| 260 | |
|---|
| 261 | // returns the leading monomial of p1 divided by the maximal power of |
|---|
| 262 | // that of p2 |
|---|
| 263 | poly pDivByMonom (poly p1,poly p2); |
|---|
| 264 | |
|---|
| 265 | // Returns as i-th entry of P the coefficient of the (i-1) power of |
|---|
| 266 | // the leading monomial of p2 in p1 |
|---|
| 267 | void pCancelPolyByMonom (poly p1,poly p2,polyset * P,int * SizeOfSet); |
|---|
| 268 | |
|---|
| 269 | // orders monoms of poly using insertion sort, performs pSetm on each |
|---|
| 270 | // monom (i.e. sets Order field) |
|---|
| 271 | poly pOrdPolyInsertSetm(poly p); |
|---|
| 272 | |
|---|
| 273 | // orders monoms of poly using merge sort (ususally faster than |
|---|
| 274 | // insertion sort). ASSUMES that pSetm was performed on monoms |
|---|
| 275 | // (i.e. that Order field is set correctly) |
|---|
| 276 | poly pOrdPolyMerge(poly p); |
|---|
| 277 | |
|---|
| 278 | poly pPermPoly (poly p, int * perm, ring OldRing, |
|---|
| 279 | int *par_perm=NULL, int OldPar=0); |
|---|
| 280 | |
|---|
| 281 | /*BOOLEAN pVectorHasUnitM(poly p, int * k);*/ |
|---|
| 282 | BOOLEAN pVectorHasUnitB(poly p, int * k); |
|---|
| 283 | void pVectorHasUnit(poly p, int * k, int * len); |
|---|
| 284 | poly pTakeOutComp1(poly * p, int k); |
|---|
| 285 | // Splits *p into two polys: *q which consists of all monoms with |
|---|
| 286 | // component == comp and *p of all other monoms *lq == pLength(*q) |
|---|
| 287 | // On rreturn all components pf *q == 0 |
|---|
| 288 | void pTakeOutComp(poly *p, Exponent_t comp, poly *q, int *lq); |
|---|
| 289 | // Similar to pTakeOutComp, except that only those components are |
|---|
| 290 | // taken out whose Order == order |
|---|
| 291 | // ASSUME: monomial ordering is Order compatible, i.e., if m1, m2 Monoms then |
|---|
| 292 | // m1 >= m2 ==> pGetOrder(m1) >= pGetOrder(m2) |
|---|
| 293 | void pDecrOrdTakeOutComp(poly *p, Exponent_t comp, Order_t order, |
|---|
| 294 | poly *q, int *lq); |
|---|
| 295 | // This is something weird -- Don't use it, unless you know what you are doing |
|---|
| 296 | poly pTakeOutComp(poly * p, int k); |
|---|
| 297 | void pSetPolyComp(poly p, int comp); |
|---|
| 298 | void pDeleteComp(poly * p,int k); |
|---|
| 299 | void pNorm(poly p); |
|---|
| 300 | poly pSubst(poly p, int n, poly e); |
|---|
| 301 | poly pJet(poly p, int m); |
|---|
| 302 | poly pJetW(poly p, int m, short * iv); |
|---|
| 303 | // maximum weigthed degree of all monomials of p, w is indexed from |
|---|
| 304 | // 1..pVariables |
|---|
| 305 | int pDegW(poly p, short *w); |
|---|
| 306 | |
|---|
| 307 | /*-----------type conversions ----------------------------*/ |
|---|
| 308 | poly pPolys2Vec(polyset p, int len); |
|---|
| 309 | void pVec2Polys(poly v, polyset *p, int *len); |
|---|
| 310 | int pVar(poly m); |
|---|
| 311 | |
|---|
| 312 | /*-----------specials for spoly-computations--------------*/ |
|---|
| 313 | poly pMultT(poly a, poly e); |
|---|
| 314 | int pDivComp(poly p, poly q); |
|---|
| 315 | BOOLEAN pCompareChain (poly p,poly p1,poly p2,poly lcm); |
|---|
| 316 | BOOLEAN pEqualPolys(poly p1,poly p2); |
|---|
| 317 | BOOLEAN pComparePolys(poly p1,poly p2); |
|---|
| 318 | |
|---|
| 319 | /*-----------Misc stuff-----------------------------------*/ |
|---|
| 320 | // returns the "Short Exponent Vector" -- used to speed up divisibility |
|---|
| 321 | // tests (see polys-impl.cc ) |
|---|
| 322 | unsigned long pGetShortExpVector(poly p); |
|---|
| 323 | |
|---|
| 324 | |
|---|
| 325 | |
|---|
| 326 | |
|---|
| 327 | #ifdef PDEBUG |
|---|
| 328 | #define pHeapTest(A,B) pDBTest(A, B, __FILE__,__LINE__) |
|---|
| 329 | #define pTest(A) pDBTest(A, currPolyBin, __FILE__,__LINE__) |
|---|
| 330 | #define prTest(p, r) prDBTest(p, r, __FILE__, __LINE__) |
|---|
| 331 | |
|---|
| 332 | BOOLEAN pDBTest(poly p, char *f, int l); |
|---|
| 333 | BOOLEAN pDBTest(poly p, omBin tail_heap, char *f, int l); |
|---|
| 334 | BOOLEAN prDBTest(poly p, ring r, char *f, int l); |
|---|
| 335 | BOOLEAN pDBTest(poly p, omBin tail_heap, omBin lm_heap, char *f, int l); |
|---|
| 336 | #else |
|---|
| 337 | #define prTest(p, r) (TRUE) |
|---|
| 338 | #define pHeapTest(A,B) (TRUE) |
|---|
| 339 | #define pTest(A) (TRUE) |
|---|
| 340 | #define pDBTest(A,B,C) (TRUE) |
|---|
| 341 | #endif |
|---|
| 342 | #endif |
|---|
