[51c163] | 1 | /**************************************** |
---|
| 2 | * Computer Algebra System SINGULAR * |
---|
| 3 | ****************************************/ |
---|
[4301a7] | 4 | /* $Id: polys-impl.cc,v 1.38 2000-03-22 17:07:25 Singular Exp $ */ |
---|
[51c163] | 5 | |
---|
| 6 | /*************************************************************** |
---|
| 7 | * |
---|
| 8 | * File: polys-impl.cc |
---|
| 9 | * Purpose: low-level procuders for polys. |
---|
| 10 | * |
---|
| 11 | * If you touch anything here, you better know what you are doing. |
---|
| 12 | * What is here should not be used directly from other routines -- the |
---|
| 13 | * encapsulations in polys.[h,cc] should be used, instead. |
---|
| 14 | * |
---|
| 15 | ***************************************************************/ |
---|
| 16 | #ifndef POLYS_IMPL_CC |
---|
| 17 | #define POLYS_IMPL_CC |
---|
| 18 | |
---|
| 19 | #include <stdio.h> |
---|
| 20 | #include <string.h> |
---|
| 21 | #include "mod2.h" |
---|
[ca7a56] | 22 | #include "tok.h" |
---|
[51c163] | 23 | #include "structs.h" |
---|
[52a5e8] | 24 | #include "mmprivate.h" |
---|
[51c163] | 25 | #include "mmemory.h" |
---|
| 26 | #include "febase.h" |
---|
| 27 | #include "numbers.h" |
---|
| 28 | #include "polys.h" |
---|
| 29 | #include "ring.h" |
---|
[e38489] | 30 | #include "polys-impl.h" |
---|
[eb17bd3] | 31 | |
---|
[e95eaa7] | 32 | #ifdef HAVE_SHIFTED_EXPONENTS |
---|
| 33 | #ifdef PDEBUG |
---|
| 34 | int pDBsyzComp=0; |
---|
| 35 | #endif |
---|
| 36 | #endif |
---|
| 37 | |
---|
[51c163] | 38 | /*************************************************************** |
---|
| 39 | * |
---|
| 40 | * Storage Managament Routines |
---|
| 41 | * |
---|
| 42 | ***************************************************************/ |
---|
| 43 | |
---|
| 44 | |
---|
| 45 | /*2 |
---|
| 46 | * delete a poly, resets pointer |
---|
| 47 | * put the monomials in the freelist |
---|
| 48 | */ |
---|
[52a5e8] | 49 | #ifdef MDEBUG |
---|
[a9a7be] | 50 | void pDBDelete(poly * p, memHeap heap, char * f, int l) |
---|
[51c163] | 51 | { |
---|
| 52 | poly h = *p; |
---|
| 53 | |
---|
| 54 | while (h!=NULL) |
---|
| 55 | { |
---|
| 56 | #ifdef LDEBUG |
---|
| 57 | nDBDelete(&(h->coef),f,l); |
---|
| 58 | #else |
---|
| 59 | nDelete(&(h->coef)); |
---|
| 60 | #endif |
---|
| 61 | pIter(h); |
---|
[a9a7be] | 62 | pDBFree1((ADDRESS)*p, heap, f,l); |
---|
[51c163] | 63 | *p=h; |
---|
| 64 | if (l>0) l= -l; |
---|
| 65 | } |
---|
| 66 | *p = NULL; |
---|
| 67 | } |
---|
| 68 | #else |
---|
[a9a7be] | 69 | void _pDelete(poly* p, memHeap heap) |
---|
[51c163] | 70 | { |
---|
| 71 | poly h = *p; |
---|
| 72 | poly pp; |
---|
| 73 | |
---|
| 74 | while (h!=NULL) |
---|
| 75 | { |
---|
| 76 | nDelete(&(h->coef)); |
---|
| 77 | pp=h; |
---|
| 78 | pIter(h); |
---|
[a9a7be] | 79 | _pFree1((ADDRESS)pp, heap); |
---|
[51c163] | 80 | } |
---|
| 81 | *p = NULL; |
---|
| 82 | } |
---|
| 83 | #endif |
---|
| 84 | |
---|
| 85 | /*2 |
---|
| 86 | * remove first monom |
---|
| 87 | */ |
---|
[52a5e8] | 88 | #ifdef MDEBUG |
---|
[a9a7be] | 89 | void pDBDelete1(poly * p, memHeap heap, char * f, int l) |
---|
[51c163] | 90 | { |
---|
| 91 | poly h = *p; |
---|
| 92 | |
---|
| 93 | if (h==NULL) return; |
---|
| 94 | nDelete(&(h->coef)); |
---|
| 95 | *p = pNext(h); |
---|
[a9a7be] | 96 | pDBFree1((ADDRESS)h, heap, f,l); |
---|
[51c163] | 97 | } |
---|
| 98 | #else |
---|
[a9a7be] | 99 | void _pDelete1(poly* p, memHeap heap) |
---|
[51c163] | 100 | { |
---|
| 101 | poly h = *p; |
---|
| 102 | |
---|
| 103 | if (h==NULL) return; |
---|
| 104 | nDelete(&(h->coef)); |
---|
| 105 | *p = pNext(h); |
---|
[a9a7be] | 106 | _pFree1((ADDRESS)h, heap); |
---|
[51c163] | 107 | } |
---|
| 108 | #endif |
---|
| 109 | |
---|
[a9a7be] | 110 | |
---|
[51c163] | 111 | void ppDelete(poly* p, ring rg) |
---|
| 112 | { |
---|
[47faf56] | 113 | ring origRing = currRing; |
---|
| 114 | rChangeCurrRing(rg, FALSE); |
---|
[51c163] | 115 | pDelete(p); |
---|
[47faf56] | 116 | rChangeCurrRing(origRing, FALSE); |
---|
[51c163] | 117 | } |
---|
| 118 | |
---|
| 119 | /*2 |
---|
| 120 | * creates a copy of p |
---|
| 121 | */ |
---|
[52a5e8] | 122 | #ifdef MDEBUG |
---|
[a9a7be] | 123 | poly pDBCopy(memHeap d_h, poly s_p, char *f,int l) |
---|
[51c163] | 124 | #else |
---|
[a9a7be] | 125 | poly _pCopy(memHeap d_h, poly s_p) |
---|
[51c163] | 126 | #endif |
---|
| 127 | { |
---|
[a9a7be] | 128 | spolyrec dp; |
---|
| 129 | poly d_p = &dp; |
---|
[51c163] | 130 | |
---|
[a9a7be] | 131 | assume(d_h != NULL && (d_h == mm_specHeap) || |
---|
[b7b08c] | 132 | d_h->size == mm_specHeap->size); |
---|
[9d06971] | 133 | pTest(s_p); |
---|
[a9a7be] | 134 | while (s_p != NULL) |
---|
| 135 | { |
---|
[52a5e8] | 136 | #ifdef MDEBUG |
---|
[a9a7be] | 137 | d_p->next = (poly) mmDBAllocHeap(d_h, f, l); |
---|
[51c163] | 138 | #else |
---|
[a9a7be] | 139 | mmAllocHeap(d_p->next, d_h); |
---|
[51c163] | 140 | #endif |
---|
[a9a7be] | 141 | d_p = d_p->next; |
---|
| 142 | pSetCoeff0(d_p, nCopy(pGetCoeff(s_p))); |
---|
| 143 | memcpyW(&(d_p->exp.l[0]), &(s_p->exp.l[0]), currRing->ExpLSize); |
---|
| 144 | pIter(s_p); |
---|
| 145 | } |
---|
| 146 | pNext(d_p) = NULL; |
---|
[87bef42] | 147 | #if defined(MDEBUG) && defined(PDEBUG) |
---|
| 148 | pDBTest(dp.next, d_h, f,l); |
---|
| 149 | #else |
---|
[a9a7be] | 150 | pHeapTest(dp.next, d_h); |
---|
[87bef42] | 151 | #endif |
---|
[a9a7be] | 152 | return dp.next; |
---|
| 153 | } |
---|
| 154 | |
---|
| 155 | #ifdef MDEBUG |
---|
| 156 | poly pDBCopy(poly s_p, char *f,int l) |
---|
| 157 | { |
---|
| 158 | return pDBCopy(mm_specHeap, s_p, f, l); |
---|
| 159 | } |
---|
| 160 | #else |
---|
| 161 | poly _pCopy(poly s_p) |
---|
| 162 | { |
---|
| 163 | return _pCopy(mm_specHeap, s_p); |
---|
| 164 | } |
---|
| 165 | #endif |
---|
| 166 | |
---|
| 167 | |
---|
| 168 | #ifdef MDEBUG |
---|
| 169 | poly pDBShallowCopyDelete(memHeap d_h,poly *p,memHeap s_h, char *f,int l) |
---|
| 170 | #else |
---|
| 171 | poly _pShallowCopyDelete(memHeap d_h, poly *p, memHeap s_h) |
---|
| 172 | #endif |
---|
| 173 | { |
---|
| 174 | spolyrec dp; |
---|
| 175 | poly d_p = &dp, tmp; |
---|
| 176 | poly s_p = *p; |
---|
| 177 | |
---|
| 178 | assume(d_h != NULL && s_h != NULL && |
---|
[b7b08c] | 179 | d_h->size == s_h->size); |
---|
[a9a7be] | 180 | |
---|
| 181 | if (currRing->ExpLSize <= 2) |
---|
[51c163] | 182 | { |
---|
[a9a7be] | 183 | if (currRing->ExpLSize == 1) |
---|
[51c163] | 184 | { |
---|
[a9a7be] | 185 | while (s_p != NULL) |
---|
| 186 | { |
---|
[52a5e8] | 187 | #ifdef MDEBUG |
---|
[a9a7be] | 188 | d_p->next = (poly) mmDBAllocHeap(d_h, f, l); |
---|
[51c163] | 189 | #else |
---|
[a9a7be] | 190 | mmAllocHeap(d_p->next, d_h); |
---|
[51c163] | 191 | #endif |
---|
[a9a7be] | 192 | d_p = d_p->next; |
---|
| 193 | |
---|
| 194 | d_p->coef = s_p->coef; |
---|
| 195 | d_p->exp.l[0] = s_p->exp.l[0]; |
---|
| 196 | |
---|
| 197 | tmp = pNext(s_p); |
---|
| 198 | #ifdef MDEBUG |
---|
| 199 | mmDBFreeHeap(s_p, s_h, f, l); |
---|
| 200 | #else |
---|
| 201 | mmFreeHeap(s_p, s_h); |
---|
| 202 | #endif |
---|
| 203 | s_p = tmp; |
---|
| 204 | } |
---|
| 205 | } |
---|
| 206 | else |
---|
| 207 | { |
---|
| 208 | while (s_p != NULL) |
---|
| 209 | { |
---|
| 210 | #ifdef MDEBUG |
---|
| 211 | d_p->next = (poly) mmDBAllocHeap(d_h, f, l); |
---|
| 212 | #else |
---|
| 213 | mmAllocHeap(d_p->next, d_h); |
---|
| 214 | #endif |
---|
| 215 | d_p = d_p->next; |
---|
| 216 | |
---|
| 217 | d_p->coef = s_p->coef; |
---|
| 218 | d_p->exp.l[0] = s_p->exp.l[0]; |
---|
| 219 | d_p->exp.l[1] = s_p->exp.l[1]; |
---|
| 220 | |
---|
| 221 | tmp = pNext(s_p); |
---|
| 222 | #ifdef MDEBUG |
---|
| 223 | mmDBFreeHeap(s_p, s_h, f, l); |
---|
| 224 | #else |
---|
| 225 | mmFreeHeap(s_p, s_h); |
---|
| 226 | #endif |
---|
| 227 | s_p = tmp; |
---|
| 228 | } |
---|
[51c163] | 229 | } |
---|
| 230 | } |
---|
[a9a7be] | 231 | else |
---|
| 232 | { |
---|
| 233 | if (currRing->ExpLSize & 1) |
---|
| 234 | { |
---|
| 235 | while (s_p != NULL) |
---|
| 236 | { |
---|
| 237 | |
---|
| 238 | #ifdef MDEBUG |
---|
| 239 | d_p->next = (poly) mmDBAllocHeap(d_h, f, l); |
---|
| 240 | #else |
---|
| 241 | mmAllocHeap(d_p->next, d_h); |
---|
| 242 | #endif |
---|
| 243 | d_p = d_p->next; |
---|
| 244 | |
---|
| 245 | d_p->coef = s_p->coef; |
---|
| 246 | memcpy_nwODD(&(d_p->exp.l[0]), &(s_p->exp.l[1]), currRing->ExpLSize); |
---|
| 247 | |
---|
| 248 | tmp = pNext(s_p); |
---|
| 249 | #ifdef MDEBUG |
---|
| 250 | mmDBFreeHeap(s_p, s_h, f, l); |
---|
| 251 | #else |
---|
| 252 | mmFreeHeap(s_p, s_h); |
---|
| 253 | #endif |
---|
| 254 | s_p = tmp; |
---|
| 255 | } |
---|
| 256 | } |
---|
| 257 | else |
---|
| 258 | { |
---|
| 259 | while (s_p != NULL) |
---|
| 260 | { |
---|
| 261 | |
---|
| 262 | #ifdef MDEBUG |
---|
| 263 | d_p->next = (poly) mmDBAllocHeap(d_h, f, l); |
---|
| 264 | #else |
---|
| 265 | mmAllocHeap(d_p->next, d_h); |
---|
| 266 | #endif |
---|
| 267 | d_p = d_p->next; |
---|
| 268 | |
---|
| 269 | d_p->coef = s_p->coef; |
---|
| 270 | memcpy_nwEVEN(&(d_p->exp.l[0]), &(s_p->exp.l[1]), currRing->ExpLSize); |
---|
| 271 | |
---|
| 272 | tmp = pNext(s_p); |
---|
| 273 | #ifdef MDEBUG |
---|
| 274 | mmDBFreeHeap(s_p, s_h, f, l); |
---|
| 275 | #else |
---|
| 276 | mmFreeHeap(s_p, s_h); |
---|
| 277 | #endif |
---|
| 278 | s_p = tmp; |
---|
| 279 | } |
---|
| 280 | } |
---|
| 281 | } |
---|
| 282 | pNext(d_p) = NULL; |
---|
| 283 | pHeapTest(dp.next, d_h); |
---|
| 284 | *p = NULL; |
---|
| 285 | return pNext(dp.next); |
---|
[51c163] | 286 | } |
---|
| 287 | |
---|
| 288 | |
---|
| 289 | /*2 |
---|
| 290 | * creates a copy of the initial monomial of p |
---|
| 291 | * sets the coeff of the copy to a defined value |
---|
| 292 | */ |
---|
[52a5e8] | 293 | #ifdef MDEBUG |
---|
[51c163] | 294 | poly pDBCopy1(poly p,char *f,int l) |
---|
[eb17bd3] | 295 | #else |
---|
[51c163] | 296 | poly _pCopy1(poly p) |
---|
[eb17bd3] | 297 | #endif |
---|
[51c163] | 298 | { |
---|
| 299 | poly w; |
---|
[52a5e8] | 300 | #ifdef MDEBUG |
---|
[a9a7be] | 301 | w = pDBNew(mm_specHeap, f,l); |
---|
[eb17bd3] | 302 | #else |
---|
[51c163] | 303 | w = pNew(); |
---|
[eb17bd3] | 304 | #endif |
---|
[51c163] | 305 | pCopy2(w,p); |
---|
| 306 | nNew(&(w->coef)); |
---|
| 307 | pNext(w) = NULL; |
---|
| 308 | return w; |
---|
| 309 | } |
---|
| 310 | |
---|
| 311 | /*2 |
---|
| 312 | * returns (a copy of) the head term of a |
---|
| 313 | */ |
---|
[52a5e8] | 314 | #ifdef MDEBUG |
---|
[a9a7be] | 315 | poly pDBHead(memHeap heap, poly p,char *f, int l) |
---|
[51c163] | 316 | #else |
---|
[a9a7be] | 317 | poly _pHead(memHeap heap, poly p) |
---|
[51c163] | 318 | #endif |
---|
| 319 | { |
---|
| 320 | poly w=NULL; |
---|
| 321 | |
---|
| 322 | if (p!=NULL) |
---|
| 323 | { |
---|
[a9a7be] | 324 | assume(heap != NULL && (heap == mm_specHeap) || |
---|
[b7b08c] | 325 | heap->size == mm_specHeap->size); |
---|
[a9a7be] | 326 | |
---|
[52a5e8] | 327 | #ifdef MDEBUG |
---|
[a9a7be] | 328 | w = (poly) mmDBAllocHeap(heap, f, l); |
---|
[51c163] | 329 | #else |
---|
[a9a7be] | 330 | mmAllocHeap(w, heap); |
---|
[51c163] | 331 | #endif |
---|
[a9a7be] | 332 | memcpyW(&(w->exp.l[0]), &(p->exp.l[0]), currRing->ExpLSize); |
---|
[51c163] | 333 | pSetCoeff0(w,nCopy(pGetCoeff(p))); |
---|
| 334 | pNext(w) = NULL; |
---|
| 335 | } |
---|
| 336 | return w; |
---|
| 337 | } |
---|
| 338 | |
---|
[a9a7be] | 339 | #ifdef MDEBUG |
---|
| 340 | poly pDBShallowCopyDeleteHead(memHeap d_h,poly *s_p,memHeap s_h, char *f,int l) |
---|
| 341 | #else |
---|
| 342 | poly _pShallowCopyDeleteHead(memHeap d_h, poly *s_p, memHeap s_h) |
---|
| 343 | #endif |
---|
| 344 | { |
---|
| 345 | poly w = NULL; |
---|
| 346 | poly p = *s_p; |
---|
| 347 | |
---|
| 348 | if (p!=NULL) |
---|
| 349 | { |
---|
| 350 | assume(d_h != NULL && s_h != NULL && |
---|
[b7b08c] | 351 | d_h->size == s_h->size); |
---|
[a9a7be] | 352 | |
---|
| 353 | #ifdef MDEBUG |
---|
| 354 | w = (poly) mmDBAllocHeap(d_h, f, l); |
---|
| 355 | #else |
---|
| 356 | mmAllocHeap(w, d_h); |
---|
| 357 | #endif |
---|
| 358 | memcpyW(&(w->exp.l[0]), &(p->exp.l[0]), currRing->ExpLSize); |
---|
| 359 | pSetCoeff0(w,pGetCoeff(p)); |
---|
| 360 | pNext(w) = NULL; |
---|
| 361 | |
---|
| 362 | *s_p = pNext(p); |
---|
| 363 | #ifdef MDEBUG |
---|
| 364 | mmDBFreeHeap(p, s_h, f, l); |
---|
| 365 | #else |
---|
| 366 | mmFreeHeap(p, s_h); |
---|
| 367 | #endif |
---|
| 368 | } |
---|
| 369 | return w; |
---|
| 370 | } |
---|
| 371 | |
---|
| 372 | |
---|
| 373 | |
---|
[51c163] | 374 | poly pHeadProc(poly p) |
---|
| 375 | { |
---|
| 376 | return pHead(p); |
---|
| 377 | } |
---|
| 378 | |
---|
| 379 | /*2 |
---|
| 380 | * returns (a copy of) the head term of a without the coef |
---|
| 381 | */ |
---|
[52a5e8] | 382 | #ifdef MDEBUG |
---|
[51c163] | 383 | poly pDBHead0(poly p,char *f, int l) |
---|
| 384 | #else |
---|
| 385 | poly _pHead0(poly p) |
---|
| 386 | #endif |
---|
| 387 | { |
---|
| 388 | poly w=NULL; |
---|
| 389 | |
---|
| 390 | if (p!=NULL) |
---|
| 391 | { |
---|
[a9a7be] | 392 | #if defined(PDEBUG) && defined(MDEBUG) |
---|
| 393 | w = pDBNew(mm_specHeap, f,l); |
---|
[51c163] | 394 | #else |
---|
| 395 | w = pNew(); |
---|
| 396 | #endif |
---|
| 397 | pCopy2(w,p); |
---|
| 398 | pSetCoeff0(w,NULL); |
---|
| 399 | pNext(w) = NULL; |
---|
| 400 | } |
---|
| 401 | return w; |
---|
| 402 | } |
---|
| 403 | |
---|
| 404 | |
---|
| 405 | /*************************************************************** |
---|
| 406 | * |
---|
| 407 | * Routines for turned on debugging |
---|
| 408 | * |
---|
| 409 | ***************************************************************/ |
---|
| 410 | |
---|
[e38489] | 411 | #if defined(PDEBUG) && PDEBUG > 1 |
---|
[51c163] | 412 | Exponent_t pPDSetExp(poly p, int v, Exponent_t e, char* f, int l) |
---|
| 413 | { |
---|
| 414 | if (v == 0) |
---|
| 415 | { |
---|
[4f011a] | 416 | Warn("zero index to exponent in %s:%d\n", f, l); |
---|
[51c163] | 417 | } |
---|
| 418 | if (v > pVariables) |
---|
| 419 | { |
---|
[4f011a] | 420 | Warn("index %d to exponent too large in %s:%d\n", v, f, l); |
---|
[51c163] | 421 | } |
---|
[a9a7be] | 422 | return (p)->exp.e[_pExpIndex(v)]=(e); |
---|
[51c163] | 423 | } |
---|
| 424 | |
---|
| 425 | Exponent_t pPDGetExp(poly p, int v, char* f, int l) |
---|
| 426 | { |
---|
| 427 | if (v == 0) |
---|
| 428 | { |
---|
[4f011a] | 429 | Warn("zero index to exponent in %s:%d\n", f, l); |
---|
[51c163] | 430 | } |
---|
| 431 | if (v > pVariables) |
---|
| 432 | { |
---|
[4f011a] | 433 | Warn("index %d to exponent too large in %s:%d\n", v, f, l); |
---|
[51c163] | 434 | } |
---|
[a9a7be] | 435 | return (p)->exp.e[_pExpIndex(v)]; |
---|
[51c163] | 436 | } |
---|
| 437 | |
---|
[e78cce] | 438 | Exponent_t pPDRingSetExp(ring r, poly p, int v, Exponent_t e, char* f, int l) |
---|
| 439 | { |
---|
| 440 | if (v == 0) |
---|
| 441 | { |
---|
[4f011a] | 442 | Warn("zero index to exponent in %s:%d\n", f, l); |
---|
[e78cce] | 443 | } |
---|
| 444 | if (v > r->N) |
---|
| 445 | { |
---|
[4f011a] | 446 | Warn("index %d to exponent too large in %s:%d\n", v, f, l); |
---|
[e78cce] | 447 | } |
---|
[a9a7be] | 448 | return (p)->exp.e[_pRingExpIndex(r, v)]=(e); |
---|
[e78cce] | 449 | } |
---|
| 450 | |
---|
| 451 | Exponent_t pPDRingGetExp(ring r, poly p, int v, char* f, int l) |
---|
| 452 | { |
---|
| 453 | if (v == 0) |
---|
| 454 | { |
---|
[4f011a] | 455 | Warn("zero index to exponent in %s:%d\n", f, l); |
---|
[e78cce] | 456 | } |
---|
| 457 | if (v > r->N) |
---|
| 458 | { |
---|
[4f011a] | 459 | Warn("index %d to exponent too large in %s:%d\n", v, f, l); |
---|
[e78cce] | 460 | } |
---|
[a9a7be] | 461 | return (p)->exp.e[_pRingExpIndex(r,v)]; |
---|
[e78cce] | 462 | } |
---|
| 463 | |
---|
[51c163] | 464 | Exponent_t pPDIncrExp(poly p, int v, char* f, int l) |
---|
| 465 | { |
---|
| 466 | if (v == 0) |
---|
| 467 | { |
---|
[4f011a] | 468 | Warn("zero index to exponent in %s:%d\n", f, l); |
---|
[51c163] | 469 | } |
---|
| 470 | if (v > pVariables) |
---|
| 471 | { |
---|
[4f011a] | 472 | Warn("index %d to exponent too large in %s:%d\n", v, f, l); |
---|
[51c163] | 473 | } |
---|
[a9a7be] | 474 | return ((p)->exp.e[_pExpIndex(v)])++; |
---|
[51c163] | 475 | } |
---|
| 476 | |
---|
| 477 | Exponent_t pPDDecrExp(poly p, int v, char* f, int l) |
---|
| 478 | { |
---|
| 479 | if (v == 0) |
---|
| 480 | { |
---|
[4f011a] | 481 | Warn("zero index to exponent in %s:%d\n", f, l); |
---|
[51c163] | 482 | } |
---|
| 483 | if (v > pVariables) |
---|
| 484 | { |
---|
[4f011a] | 485 | Warn("index %d to exponent too large in %s:%d\n", v, f, l); |
---|
[51c163] | 486 | } |
---|
[a9a7be] | 487 | return ((p)->exp.e[_pExpIndex(v)])--; |
---|
[51c163] | 488 | } |
---|
| 489 | |
---|
| 490 | Exponent_t pPDAddExp(poly p, int v, Exponent_t e, char* f, int l) |
---|
| 491 | { |
---|
| 492 | if (v == 0) |
---|
| 493 | { |
---|
[4f011a] | 494 | Warn("zero index to exponent in %s:%d\n", f, l); |
---|
[51c163] | 495 | } |
---|
| 496 | if (v > pVariables) |
---|
| 497 | { |
---|
[4f011a] | 498 | Warn("index %d to exponent too large in %s:%d\n", v, f, l); |
---|
[51c163] | 499 | } |
---|
[a9a7be] | 500 | return ((p)->exp.e[_pExpIndex(v)]) += (e); |
---|
[51c163] | 501 | } |
---|
| 502 | |
---|
| 503 | Exponent_t pPDSubExp(poly p, int v, Exponent_t e, char* f, int l) |
---|
| 504 | { |
---|
| 505 | if (v == 0) |
---|
| 506 | { |
---|
[4f011a] | 507 | Warn("zero index to exponent in %s:%d\n", f, l); |
---|
[51c163] | 508 | } |
---|
| 509 | if (v > pVariables) |
---|
| 510 | { |
---|
[4f011a] | 511 | Warn("index %d to exponent too large in %s:%d\n", v, f, l); |
---|
[51c163] | 512 | } |
---|
[a9a7be] | 513 | return ((p)->exp.e[_pExpIndex(v)]) -= (e); |
---|
[51c163] | 514 | } |
---|
| 515 | |
---|
| 516 | Exponent_t pPDMultExp(poly p, int v, Exponent_t e, char* f, int l) |
---|
| 517 | { |
---|
| 518 | if (v == 0) |
---|
| 519 | { |
---|
[4f011a] | 520 | Warn("zero index to exponent in %s:%d\n", f, l); |
---|
[51c163] | 521 | } |
---|
| 522 | if (v > pVariables) |
---|
| 523 | { |
---|
[4f011a] | 524 | Warn("index %d to exponent too large in %s:%d\n", v, f, l); |
---|
[51c163] | 525 | } |
---|
[a9a7be] | 526 | return ((p)->exp.e[_pExpIndex(v)]) *= (e); |
---|
| 527 | } |
---|
| 528 | |
---|
| 529 | // routines on components |
---|
| 530 | Exponent_t pDBSetComp(poly p, Exponent_t k, int le, char* f, int l) |
---|
| 531 | { |
---|
| 532 | if (k < 0) |
---|
| 533 | { |
---|
| 534 | Warn("set negative component %d in %s:%d", k, f, l); |
---|
| 535 | } |
---|
| 536 | if (currRing->order[1] == ringorder_S) |
---|
| 537 | { |
---|
| 538 | if (le <= 0) le = currRing->typ[1].data.syzcomp.length; |
---|
| 539 | if (k > l) |
---|
| 540 | { |
---|
| 541 | Warn("component %d larger then max %d in %s:%d", |
---|
| 542 | k, le, f, l); |
---|
| 543 | } |
---|
| 544 | } |
---|
| 545 | return _pGetComp(p) = (k); |
---|
| 546 | } |
---|
| 547 | |
---|
| 548 | Exponent_t pDBDecrComp(poly p, char* f, int l) |
---|
| 549 | { |
---|
| 550 | if (_pGetComp(p) < 1) |
---|
| 551 | { |
---|
[4f011a] | 552 | Warn("decrement to negative component %d in %s:%d\n", _pGetComp(p), f, l); |
---|
[a9a7be] | 553 | } |
---|
| 554 | return _pGetComp(p)--; |
---|
| 555 | } |
---|
| 556 | |
---|
| 557 | Exponent_t pDBAddComp(poly p, Exponent_t k, int le, char* f, int l) |
---|
| 558 | { |
---|
| 559 | if (_pGetComp(p) + k < 0) |
---|
| 560 | { |
---|
[4f011a] | 561 | Warn("add to negative component %d + %d = %d in %s:%d\n", _pGetComp(p), |
---|
[a9a7be] | 562 | k, _pGetComp(p) + k, f, l); |
---|
| 563 | } |
---|
| 564 | _pGetComp(p) += (k); |
---|
| 565 | |
---|
| 566 | if (currRing->order[1] == ringorder_S) |
---|
| 567 | { |
---|
| 568 | if (le <= 0) le = currRing->typ[1].data.syzcomp.length; |
---|
| 569 | if (_pGetComp(p) > le) |
---|
| 570 | { |
---|
[4f011a] | 571 | Warn("sum of components %d larger then max %d in %s:%d\n", |
---|
[a9a7be] | 572 | _pGetComp(p), le, f, l); |
---|
| 573 | assume(0); |
---|
| 574 | } |
---|
| 575 | } |
---|
| 576 | return _pGetComp(p); |
---|
| 577 | } |
---|
| 578 | |
---|
| 579 | Exponent_t pDBSubComp(poly p, Exponent_t k, char* f, int l) |
---|
| 580 | { |
---|
| 581 | if (_pGetComp(p) - k < 0) |
---|
| 582 | { |
---|
[4f011a] | 583 | Warn("sub to negative component %d - %d = %d in %s:%d\n", _pGetComp(p), |
---|
[a9a7be] | 584 | k, _pGetComp(p) - k, f, l); |
---|
| 585 | } |
---|
| 586 | return _pGetComp(p) -= (k); |
---|
| 587 | } |
---|
| 588 | |
---|
| 589 | Exponent_t pDBRingSetComp(ring r, poly p, Exponent_t k, char* f, int l) |
---|
| 590 | { |
---|
| 591 | if (k < 0) |
---|
| 592 | { |
---|
[4f011a] | 593 | Warn("set negative component %d in %s:%d\n", k, f, l); |
---|
[a9a7be] | 594 | } |
---|
| 595 | return _pRingGetComp(r, p) = (k); |
---|
[51c163] | 596 | } |
---|
| 597 | |
---|
| 598 | // checks whether fast monom add did not overflow |
---|
[e38489] | 599 | void pDBMonAddOn(poly p1, poly p2, char* f, int l) |
---|
[51c163] | 600 | { |
---|
[4f011a] | 601 | poly ptemp = pHead0(p1); |
---|
[c2b529] | 602 | |
---|
[a9a7be] | 603 | if (pGetComp(p1) != 0 && pGetComp(p2) != 0) |
---|
| 604 | { |
---|
[e38489] | 605 | Warn("Error in pMonAddOn: both components %d:%d !=0 in %s:%d", |
---|
[a9a7be] | 606 | pGetComp(p1), pGetComp(p2), f, l); |
---|
| 607 | } |
---|
| 608 | |
---|
[e38489] | 609 | __pMonAddOn(p1, p2); |
---|
[51c163] | 610 | |
---|
| 611 | for (int i=1; i<=pVariables; i++) |
---|
[eb17bd3] | 612 | { |
---|
[51c163] | 613 | pAddExp(ptemp, i, pGetExp(p2, i)); |
---|
[eb17bd3] | 614 | } |
---|
[4f011a] | 615 | pAddComp(ptemp, pGetComp(p2)); |
---|
[a9a7be] | 616 | pSetm(ptemp); |
---|
[51c163] | 617 | |
---|
| 618 | if (! pEqual(ptemp, p1)) |
---|
[196a7b] | 619 | { |
---|
[4f011a] | 620 | Warn("Error in pMonAddOn in %s:%d\n", f, l); |
---|
[196a7b] | 621 | } |
---|
[c2b529] | 622 | |
---|
[51c163] | 623 | pFree1(ptemp); |
---|
| 624 | } |
---|
| 625 | |
---|
[e38489] | 626 | void pDBMonSubFrom(poly p1, poly p2, char* f, int l) |
---|
[a9a7be] | 627 | { |
---|
| 628 | poly ptemp = pNew(); |
---|
| 629 | pCopy2(ptemp, p1); |
---|
| 630 | |
---|
| 631 | if ((pGetComp(p1) != pGetComp(p2)) && (pGetComp(p2)!=0)) |
---|
| 632 | { |
---|
[e38489] | 633 | Warn("Error in pMonSubFrom: components are different %d:%d in %s:%d", |
---|
[a9a7be] | 634 | pGetComp(p1), pGetComp(p2), f, l); |
---|
| 635 | } |
---|
| 636 | |
---|
[e38489] | 637 | __pMonSubFrom(p1, p2); |
---|
[a9a7be] | 638 | |
---|
| 639 | for (int i=1; i<=pVariables; i++) |
---|
| 640 | { |
---|
[e38489] | 641 | if (pGetExp(ptemp, i) < pGetExp(p2, i)) |
---|
| 642 | { |
---|
| 643 | Warn("Error in pMonSubFrom: %dth exponent %d of p1 smaller than %d of p2", i, pGetExp(ptemp, i), pGetExp(p2, i)); |
---|
| 644 | } |
---|
[a9a7be] | 645 | pSubExp(ptemp, i, pGetExp(p2, i)); |
---|
| 646 | } |
---|
| 647 | pSetComp(ptemp, pGetComp(ptemp)-pGetComp(p2)); |
---|
| 648 | pSetm(ptemp); |
---|
| 649 | |
---|
| 650 | if (! pEqual(ptemp, p1)) |
---|
| 651 | { |
---|
[e38489] | 652 | Warn("Error in pMonSubFrom in %s:%d", f, l); |
---|
[a9a7be] | 653 | } |
---|
| 654 | |
---|
| 655 | pFree1(ptemp); |
---|
| 656 | } |
---|
| 657 | |
---|
[416465] | 658 | void pDBMonAdd(poly p1, poly p2, poly p3, ring r, char* f, int l) |
---|
[51c163] | 659 | { |
---|
[416465] | 660 | if (r == currRing) |
---|
[51c163] | 661 | { |
---|
[416465] | 662 | if (pGetComp(p3) != 0 && pGetComp(p2) != 0) |
---|
| 663 | { |
---|
| 664 | Warn("Error in pMonAdd: both components %d:%d !=0 in %s:%d", |
---|
| 665 | pGetComp(p3), pGetComp(p2), f, l); |
---|
| 666 | } |
---|
| 667 | if (p2 == p1 || p3 == p1) |
---|
| 668 | { |
---|
| 669 | Warn("Error in pMonAdd: Destination equals source in %s:%d", f, l); |
---|
| 670 | } |
---|
[51c163] | 671 | } |
---|
[416465] | 672 | __pMonAdd(p1, p2, p3, r); |
---|
[51c163] | 673 | |
---|
[416465] | 674 | if (r == currRing) |
---|
[e56c23] | 675 | { |
---|
[416465] | 676 | poly ptemp = pInit(); |
---|
| 677 | for (int i=1; i<=pVariables; i++) |
---|
[e56c23] | 678 | { |
---|
[416465] | 679 | pSetExp(ptemp, i, pGetExp(p2, i) + pGetExp(p3, i)); |
---|
| 680 | if (pGetExp(ptemp, i) != pGetExp(p1, i)) |
---|
| 681 | { |
---|
| 682 | Warn("Error in pMonAdd: %th exponent: %d != (%d == %d + %d)", |
---|
| 683 | i, pGetExp(p1, i), pGetExp(ptemp, i), pGetExp(p2, i), |
---|
| 684 | pGetExp(p3, i)); |
---|
| 685 | } |
---|
[e56c23] | 686 | } |
---|
[416465] | 687 | pSetComp(ptemp, pGetComp(p2) + pGetComp(p3)); |
---|
| 688 | pSetm(ptemp); |
---|
[51c163] | 689 | |
---|
[416465] | 690 | if (! pEqual(ptemp, p1)) |
---|
| 691 | Warn("Error in pMonAdd in %s:%d", f, l); |
---|
| 692 | pFree1(ptemp); |
---|
| 693 | } |
---|
[51c163] | 694 | } |
---|
[eb17bd3] | 695 | |
---|
[e78cce] | 696 | static BOOLEAN OldpDivisibleBy(poly a, poly b) |
---|
| 697 | { |
---|
[51c163] | 698 | if ((a!=NULL)&&((pGetComp(a)==0) || (pGetComp(a) == pGetComp(b)))) |
---|
| 699 | { |
---|
| 700 | for (int i = 1; i<=pVariables; i++) |
---|
[e78cce] | 701 | if (pGetExp(a, i) > pGetExp(b,i)) return FALSE; |
---|
| 702 | return TRUE; |
---|
[51c163] | 703 | } |
---|
[e78cce] | 704 | return FALSE; |
---|
| 705 | } |
---|
| 706 | |
---|
| 707 | |
---|
| 708 | BOOLEAN pDBDivisibleBy(poly a, poly b, char* f, int l) |
---|
| 709 | { |
---|
| 710 | BOOLEAN istrue = OldpDivisibleBy(a,b); |
---|
| 711 | BOOLEAN f_istrue = _pDivisibleBy_orig(a, b); |
---|
| 712 | |
---|
[51c163] | 713 | if (istrue != f_istrue) |
---|
| 714 | { |
---|
[4f011a] | 715 | Warn("Error in pDivisibleBy in %s:%d\n", f, l); |
---|
[51c163] | 716 | _pDivisibleBy_orig(a, b); |
---|
| 717 | } |
---|
[a9a7be] | 718 | return istrue; |
---|
[51c163] | 719 | } |
---|
| 720 | |
---|
[e78cce] | 721 | BOOLEAN pDBDivisibleBy1(poly a, poly b, char* f, int l) |
---|
| 722 | { |
---|
| 723 | BOOLEAN istrue = OldpDivisibleBy(a,b); |
---|
| 724 | BOOLEAN f_istrue = _pDivisibleBy1_orig(a, b); |
---|
| 725 | |
---|
| 726 | if (istrue != f_istrue) |
---|
| 727 | { |
---|
[4f011a] | 728 | Warn("Error in pDivisibleBy1 in %s:%d\n", f, l); |
---|
[e78cce] | 729 | _pDivisibleBy1_orig(a, b); |
---|
| 730 | } |
---|
[a9a7be] | 731 | return istrue; |
---|
[e78cce] | 732 | } |
---|
| 733 | |
---|
| 734 | BOOLEAN pDBDivisibleBy2(poly a, poly b, char* f, int l) |
---|
| 735 | { |
---|
| 736 | BOOLEAN istrue = OldpDivisibleBy(a,b); |
---|
| 737 | BOOLEAN f_istrue = __pDivisibleBy(a, b); |
---|
[eb17bd3] | 738 | |
---|
[e78cce] | 739 | if (istrue != f_istrue) |
---|
| 740 | { |
---|
[4f011a] | 741 | Warn("Error in pDivisibleBy2 in %s:%d\n", f, l); |
---|
[e78cce] | 742 | __pDivisibleBy(a, b); |
---|
| 743 | } |
---|
| 744 | return f_istrue; |
---|
| 745 | } |
---|
| 746 | |
---|
[e38489] | 747 | #endif // defined(PDEBUG) && PDEBUG > 1 |
---|
[eb17bd3] | 748 | |
---|
[e38489] | 749 | #ifdef PDEBUG |
---|
[c54075] | 750 | BOOLEAN mmDBEqual(poly p, poly q, char *f, int l) |
---|
| 751 | { |
---|
| 752 | int i; |
---|
[e95eaa7] | 753 | |
---|
[c54075] | 754 | for (i = 1; i<=pVariables; i++) |
---|
| 755 | { |
---|
| 756 | if (pGetExp(p,i) != pGetExp(q, i)) return FALSE; |
---|
| 757 | } |
---|
| 758 | if (pGetComp(p) != pGetComp(q)) return FALSE; |
---|
| 759 | if (__pEqual(p, q) != TRUE) |
---|
| 760 | { |
---|
| 761 | Warn("Error in pEqual: exp/comp same, bug monoms different in %s:%d", |
---|
| 762 | f, l); |
---|
| 763 | } |
---|
| 764 | return TRUE; |
---|
| 765 | } |
---|
| 766 | |
---|
[416465] | 767 | BOOLEAN prDBTest(poly p, ring r, char* f, int l) |
---|
| 768 | { |
---|
| 769 | ring old_ring = NULL; |
---|
| 770 | BOOLEAN res; |
---|
[e95eaa7] | 771 | |
---|
[416465] | 772 | if (r != currRing) |
---|
| 773 | { |
---|
| 774 | old_ring = currRing; |
---|
| 775 | rChangeCurrRing(r); |
---|
| 776 | } |
---|
| 777 | res = pDBTest(p, currRing->mm_specHeap, f, l); |
---|
| 778 | if (old_ring != NULL) |
---|
| 779 | { |
---|
| 780 | rChangeCurrRing(old_ring); |
---|
| 781 | } |
---|
| 782 | return res; |
---|
| 783 | } |
---|
| 784 | |
---|
| 785 | |
---|
[51c163] | 786 | BOOLEAN pDBTest(poly p, char *f, int l) |
---|
[52a5e8] | 787 | { |
---|
| 788 | return pDBTest(p, mm_specHeap, f,l); |
---|
| 789 | } |
---|
| 790 | |
---|
| 791 | BOOLEAN pDBTest(poly p, memHeap heap, char *f, int l) |
---|
[51c163] | 792 | { |
---|
| 793 | poly old=NULL; |
---|
| 794 | BOOLEAN ismod=FALSE; |
---|
[a9a7be] | 795 | if (heap == NULL) heap = mm_specHeap; |
---|
| 796 | |
---|
[51c163] | 797 | while (p!=NULL) |
---|
| 798 | { |
---|
[e38489] | 799 | if (heap != MM_UNKNOWN_HEAP) |
---|
| 800 | { |
---|
[51c163] | 801 | #ifdef MDEBUG |
---|
[52a5e8] | 802 | if (!mmDBTestHeapBlock(p, heap, f,l)) |
---|
[51c163] | 803 | return FALSE; |
---|
[52a5e8] | 804 | #elif defined(HEAP_DEBUG) |
---|
| 805 | if (! mmDebugCheckHeapAddr(p, heap, MM_HEAP_ADDR_USED_FLAG, f, l)) |
---|
| 806 | return FALSE; |
---|
[51c163] | 807 | #endif |
---|
[e38489] | 808 | } |
---|
| 809 | |
---|
[51c163] | 810 | #ifdef LDEBUG |
---|
| 811 | if (!nDBTest(p->coef,f,l)) |
---|
[a38f5ea] | 812 | return FALSE; |
---|
[51c163] | 813 | #endif |
---|
| 814 | if ((p->coef==NULL)&&(nGetChar()<2)) |
---|
| 815 | { |
---|
[4f011a] | 816 | Warn("NULL coef in poly in %s:%d\n",f,l); |
---|
[51c163] | 817 | return FALSE; |
---|
| 818 | } |
---|
| 819 | if (nIsZero(p->coef)) |
---|
| 820 | { |
---|
[4f011a] | 821 | Warn("zero coef in poly in %s:%d\n",f,l); |
---|
[51c163] | 822 | return FALSE; |
---|
| 823 | } |
---|
| 824 | int i=pVariables; |
---|
[69618d] | 825 | #ifndef HAVE_SHIFTED_EXPONENTS |
---|
| 826 | // can not hapen for SHIFTED_EXPONENTS |
---|
[51c163] | 827 | for(;i;i--) |
---|
| 828 | { |
---|
| 829 | if (pGetExp(p,i)<0) |
---|
| 830 | { |
---|
[4f011a] | 831 | Warn("neg. Exponent %d of x(%d) in %s:%d\n",pGetExp(p,i),i,f,l); |
---|
[51c163] | 832 | return FALSE; |
---|
| 833 | } |
---|
| 834 | } |
---|
[69618d] | 835 | #endif |
---|
[51c163] | 836 | if (pGetComp(p)<0) |
---|
| 837 | { |
---|
[4f011a] | 838 | Warn("neg Component in %s:%d\n",f,l); |
---|
[51c163] | 839 | return FALSE; |
---|
| 840 | } |
---|
| 841 | if (ismod==0) |
---|
| 842 | { |
---|
| 843 | if (pGetComp(p)>0) ismod=1; |
---|
| 844 | else ismod=2; |
---|
| 845 | } |
---|
| 846 | else if (ismod==1) |
---|
| 847 | { |
---|
| 848 | if (pGetComp(p)==0) |
---|
| 849 | { |
---|
[4f011a] | 850 | Warn("mix vec./poly in %s:%d\n",f,l); |
---|
[51c163] | 851 | return FALSE; |
---|
| 852 | } |
---|
| 853 | } |
---|
| 854 | else if (ismod==2) |
---|
| 855 | { |
---|
| 856 | if (pGetComp(p)!=0) |
---|
| 857 | { |
---|
[4f011a] | 858 | Warn("mix poly/vec. in %s:%d\n",f,l); |
---|
[51c163] | 859 | return FALSE; |
---|
| 860 | } |
---|
| 861 | } |
---|
[a9a7be] | 862 | if (currRing->order[1] == ringorder_S) |
---|
[51c163] | 863 | { |
---|
[a9a7be] | 864 | long c1, cc1, ccc1, ec1; |
---|
| 865 | sro_ord* o = &(currRing->typ[1]); |
---|
| 866 | |
---|
| 867 | c1 = pGetComp(p); |
---|
| 868 | cc1 = o->data.syzcomp.Components[c1]; |
---|
| 869 | ccc1 = o->data.syzcomp.ShiftedComponents[cc1]; |
---|
| 870 | if (! (c1 == 0 || cc1 != 0)) |
---|
| 871 | { |
---|
[4f011a] | 872 | Warn("Component <-> TrueComponent zero mismatch\n", f, l); |
---|
[a9a7be] | 873 | return FALSE; |
---|
| 874 | } |
---|
| 875 | if (! (c1 == 0 || ccc1 != 0)) |
---|
| 876 | { |
---|
[4f011a] | 877 | Warn("Component <-> ShiftedComponent zero mismatch\n", f, l); |
---|
[a9a7be] | 878 | return FALSE; |
---|
| 879 | } |
---|
| 880 | ec1 = p->exp.l[currRing->typ[1].data.syzcomp.place]; |
---|
| 881 | if (ec1 != ccc1) |
---|
| 882 | { |
---|
[4f011a] | 883 | Warn("Shifted comp out of sync. should %d, is %d", ccc1, ec1); |
---|
[a9a7be] | 884 | return FALSE; |
---|
| 885 | } |
---|
[51c163] | 886 | } |
---|
[9d06971] | 887 | if (currRing->order[0] == ringorder_s) |
---|
| 888 | { |
---|
[416465] | 889 | unsigned long syzindex = p->exp.l[currRing->typ[0].data.syz.place]; |
---|
| 890 | pSetm(p); |
---|
| 891 | if (p->exp.l[currRing->typ[0].data.syz.place] != syzindex) |
---|
[9d06971] | 892 | { |
---|
[416465] | 893 | Warn("Syzindex wrong: Was %dl but should be %d in %s:%d\n", |
---|
| 894 | syzindex, p->exp.l[currRing->typ[0].data.syz.place], f, l); |
---|
[9d06971] | 895 | } |
---|
| 896 | } |
---|
[51c163] | 897 | old=p; |
---|
| 898 | pIter(p); |
---|
| 899 | if (pComp(old,p)!=1) |
---|
| 900 | { |
---|
[f907fc] | 901 | Warn("wrong order ("); |
---|
[ca7a56] | 902 | wrp(old); |
---|
[87bef42] | 903 | Print(") in %s:%d (pComp=%d)\n",f,l,pComp(old,p)); |
---|
[51c163] | 904 | return FALSE; |
---|
| 905 | } |
---|
[c1489f2] | 906 | if (p != NULL) |
---|
| 907 | { |
---|
| 908 | if (pGetComp(old) == pGetComp(p)) |
---|
| 909 | { |
---|
| 910 | i = pVariables; |
---|
| 911 | for (i=pVariables;i; i--) |
---|
| 912 | { |
---|
| 913 | if (pGetExp(old, i) != pGetExp(p, i)) break; |
---|
| 914 | } |
---|
| 915 | if (i == 0) |
---|
| 916 | { |
---|
| 917 | Warn("different Compare, but same exponent vector for"); |
---|
| 918 | wrp(old); |
---|
| 919 | Print(" in %s%d\n", f, l); |
---|
| 920 | return FALSE; |
---|
| 921 | } |
---|
| 922 | } |
---|
| 923 | } |
---|
[51c163] | 924 | } |
---|
| 925 | return TRUE; |
---|
| 926 | } |
---|
[e95eaa7] | 927 | |
---|
| 928 | #ifdef HAVE_SHIFTED_EXPONENTS |
---|
[8b4b7b] | 929 | int rComp_a(poly p1, poly p2, int i) |
---|
[e95eaa7] | 930 | { |
---|
| 931 | int j; |
---|
| 932 | int o1=0; |
---|
| 933 | int o2=0; |
---|
| 934 | for(j=currRing->block0[i];j<=currRing->block1[i];j++) |
---|
| 935 | { |
---|
| 936 | o1+=pGetExp(p1,j)*currRing->wvhdl[i][j-currRing->block0[i]]; |
---|
| 937 | o2+=pGetExp(p2,j)*currRing->wvhdl[i][j-currRing->block0[i]]; |
---|
| 938 | } |
---|
| 939 | if (o1>o2) |
---|
| 940 | { |
---|
[4301a7] | 941 | return 1; |
---|
[e95eaa7] | 942 | } |
---|
| 943 | if (o1<o2) |
---|
| 944 | { |
---|
[4301a7] | 945 | return -1; |
---|
[e95eaa7] | 946 | } |
---|
| 947 | return 0; |
---|
| 948 | } |
---|
[8b4b7b] | 949 | int rComp_deg(poly p1, poly p2, int i) |
---|
[e95eaa7] | 950 | { |
---|
| 951 | int j; |
---|
| 952 | int o1=0; |
---|
| 953 | int o2=0; |
---|
| 954 | for(j=currRing->block0[i];j<=currRing->block1[i];j++) |
---|
| 955 | { |
---|
| 956 | o1+=pGetExp(p1,j); |
---|
| 957 | o2+=pGetExp(p2,j); |
---|
| 958 | } |
---|
[8b4b7b] | 959 | if (o1>o2) return 1; |
---|
| 960 | if (o1<o2) return -1; |
---|
[e95eaa7] | 961 | return 0; |
---|
| 962 | } |
---|
[8b4b7b] | 963 | int rComp_lex(poly p1, poly p2, int i) |
---|
[e95eaa7] | 964 | { |
---|
| 965 | int j; |
---|
| 966 | for(j=currRing->block0[i];j<=currRing->block1[i];j++) |
---|
| 967 | { |
---|
| 968 | if (pGetExp(p1,j) > pGetExp(p2,j)) |
---|
| 969 | { |
---|
| 970 | return 1; |
---|
| 971 | } |
---|
| 972 | else if (pGetExp(p1,j) < pGetExp(p2,j)) |
---|
| 973 | { |
---|
| 974 | return -1; |
---|
| 975 | } |
---|
| 976 | } |
---|
| 977 | return 0; |
---|
| 978 | } |
---|
[8b4b7b] | 979 | int rComp_revlex(poly p1, poly p2, int i) |
---|
[e95eaa7] | 980 | { |
---|
| 981 | int j; |
---|
| 982 | int e1,e2; |
---|
| 983 | for(j=currRing->block1[i];j>=currRing->block0[i];j--) |
---|
| 984 | { |
---|
| 985 | e1=pGetExp(p1,j); |
---|
| 986 | e2=pGetExp(p2,j); |
---|
| 987 | if (e1 < e2) |
---|
| 988 | { |
---|
| 989 | return 1; |
---|
| 990 | } |
---|
| 991 | else if (e1 > e2) |
---|
| 992 | { |
---|
| 993 | return -1; |
---|
| 994 | } |
---|
| 995 | } |
---|
| 996 | return 0; |
---|
| 997 | } |
---|
| 998 | int rComp0(poly p1, poly p2) |
---|
| 999 | { |
---|
| 1000 | int rr,r; |
---|
| 1001 | _prMonCmp(p1, p2, currRing, {rr=0;goto next_after_comp;}, {rr=1;goto next_after_comp;}, {rr=-1;goto next_after_comp;}); |
---|
| 1002 | next_after_comp: |
---|
| 1003 | int n=rBlocks(currRing)-1; |
---|
| 1004 | int i,j; |
---|
[8b4b7b] | 1005 | // check fur syzcomp - special case, not visible in currRing->order |
---|
| 1006 | if ((currRing->typ[0].ord_typ==ro_syz) |
---|
| 1007 | && (pGetComp(p1)!=pGetComp(p2))) |
---|
| 1008 | { |
---|
| 1009 | i=pGetComp(p1); |
---|
| 1010 | j=pGetComp(p2); |
---|
| 1011 | if (i>0) |
---|
| 1012 | { |
---|
| 1013 | if (i > currRing->typ[0].data.syz.limit) |
---|
| 1014 | i=currRing->typ[0].data.syz.curr_index; |
---|
| 1015 | else |
---|
| 1016 | i=currRing->typ[0].data.syz.syz_index[i]; |
---|
| 1017 | } |
---|
| 1018 | if (j>0) |
---|
| 1019 | { |
---|
| 1020 | if (j > currRing->typ[0].data.syz.limit) |
---|
| 1021 | j=currRing->typ[0].data.syz.curr_index; |
---|
| 1022 | else |
---|
| 1023 | j=currRing->typ[0].data.syz.syz_index[j]; |
---|
| 1024 | } |
---|
| 1025 | if (i >j ) { assume(rr== -1); return -1;} |
---|
| 1026 | if (i <j ) { assume(rr== 1); return 1;} |
---|
| 1027 | } |
---|
[e95eaa7] | 1028 | for(i=0;i<n;i++) |
---|
| 1029 | { |
---|
| 1030 | switch (currRing->order[i]) |
---|
| 1031 | { |
---|
| 1032 | case ringorder_a: |
---|
[8b4b7b] | 1033 | r=rComp_a(p1,p2,i); |
---|
[4301a7] | 1034 | if (r!=0) { assume(r==rr);return rr;} |
---|
[e95eaa7] | 1035 | break; |
---|
| 1036 | |
---|
| 1037 | case ringorder_c: |
---|
| 1038 | if (pGetComp(p1) < pGetComp(p2)) |
---|
| 1039 | { |
---|
| 1040 | assume(rr==1); |
---|
[4301a7] | 1041 | return rr; |
---|
[e95eaa7] | 1042 | } |
---|
| 1043 | if (pGetComp(p1) > pGetComp(p2)) |
---|
| 1044 | { |
---|
| 1045 | assume(rr==-1); |
---|
[4301a7] | 1046 | return rr; |
---|
[e95eaa7] | 1047 | } |
---|
| 1048 | break; |
---|
| 1049 | |
---|
| 1050 | case ringorder_C: |
---|
| 1051 | if (pGetComp(p1) > pGetComp(p2)) |
---|
| 1052 | { |
---|
| 1053 | assume(rr==1); |
---|
[4301a7] | 1054 | return rr; |
---|
[e95eaa7] | 1055 | } |
---|
| 1056 | if (pGetComp(p1) < pGetComp(p2)) |
---|
| 1057 | { |
---|
| 1058 | assume(rr==-1); |
---|
[4301a7] | 1059 | return rr; |
---|
[e95eaa7] | 1060 | } |
---|
| 1061 | break; |
---|
| 1062 | |
---|
| 1063 | case ringorder_M: |
---|
| 1064 | { |
---|
| 1065 | assume(0); // not yet implemented |
---|
[4301a7] | 1066 | return rr; |
---|
[e95eaa7] | 1067 | } |
---|
| 1068 | |
---|
| 1069 | case ringorder_lp: |
---|
[8b4b7b] | 1070 | r=rComp_lex(p1,p2,i); |
---|
[4301a7] | 1071 | if (r!=0) { assume(r==rr);return rr;} |
---|
[e95eaa7] | 1072 | break; |
---|
| 1073 | |
---|
| 1074 | case ringorder_ls: |
---|
[4301a7] | 1075 | r= -rComp_lex(p1,p2,i); |
---|
| 1076 | if (r!=0) { assume(r==rr);return rr;} |
---|
[e95eaa7] | 1077 | break; |
---|
| 1078 | |
---|
| 1079 | case ringorder_dp: |
---|
[8b4b7b] | 1080 | r=rComp_deg(p1,p2,i); |
---|
[4301a7] | 1081 | if (r!=0) { assume(r==rr);return rr;} |
---|
[8b4b7b] | 1082 | r=rComp_revlex(p1,p2,i); |
---|
[4301a7] | 1083 | if (r!=0) { assume(r==rr);return rr;} |
---|
[e95eaa7] | 1084 | break; |
---|
| 1085 | |
---|
| 1086 | case ringorder_Dp: |
---|
[8b4b7b] | 1087 | r=rComp_deg(p1,p2,i); |
---|
[4301a7] | 1088 | if (r!=0) { assume(r==rr);return rr;} |
---|
[8b4b7b] | 1089 | r=rComp_lex(p1,p2,i); |
---|
[4301a7] | 1090 | if (r!=0) { assume(r==rr);return rr;} |
---|
[e95eaa7] | 1091 | break; |
---|
| 1092 | |
---|
| 1093 | case ringorder_ds: |
---|
[8b4b7b] | 1094 | r= -rComp_deg(p1,p2,i); |
---|
[4301a7] | 1095 | if (r!=0) { assume(r==rr);return rr;} |
---|
[8b4b7b] | 1096 | r=rComp_revlex(p1,p2,i); |
---|
[4301a7] | 1097 | if (r!=0) { assume(r==rr);return rr;} |
---|
[e95eaa7] | 1098 | break; |
---|
| 1099 | |
---|
| 1100 | case ringorder_Ds: |
---|
[8b4b7b] | 1101 | r= -rComp_deg(p1,p2,i); |
---|
[4301a7] | 1102 | if (r!=0) { assume(r==rr);return rr;} |
---|
[8b4b7b] | 1103 | r=rComp_lex(p1,p2,i); |
---|
[4301a7] | 1104 | if (r!=0) { assume(r==rr);return rr;} |
---|
[e95eaa7] | 1105 | break; |
---|
| 1106 | |
---|
| 1107 | case ringorder_wp: |
---|
[8b4b7b] | 1108 | r=rComp_a(p1,p2,i); |
---|
[4301a7] | 1109 | if (r!=0) { assume(r==rr);return rr;} |
---|
[8b4b7b] | 1110 | r=rComp_revlex(p1,p2,i); |
---|
[4301a7] | 1111 | if (r!=0) { assume(r==rr);return rr;} |
---|
[e95eaa7] | 1112 | break; |
---|
| 1113 | |
---|
| 1114 | case ringorder_Wp: |
---|
[8b4b7b] | 1115 | r=rComp_a(p1,p2,i); |
---|
[4301a7] | 1116 | if (r!=0) { assume(r==rr);return rr;} |
---|
[8b4b7b] | 1117 | r=rComp_lex(p1,p2,i); |
---|
[4301a7] | 1118 | if (r!=0) { assume(r==rr);return rr;} |
---|
[e95eaa7] | 1119 | break; |
---|
| 1120 | |
---|
| 1121 | case ringorder_ws: |
---|
[8b4b7b] | 1122 | r= -rComp_a(p1,p2,i); |
---|
[4301a7] | 1123 | if (r!=0) { assume(r==rr);return rr;} |
---|
[8b4b7b] | 1124 | r=rComp_revlex(p1,p2,i); |
---|
[4301a7] | 1125 | if (r!=0) { assume(r==rr);return rr;} |
---|
[e95eaa7] | 1126 | break; |
---|
| 1127 | |
---|
| 1128 | case ringorder_Ws: |
---|
[8b4b7b] | 1129 | r= -rComp_a(p1,p2,i); |
---|
[4301a7] | 1130 | if (r!=0) { assume(r==rr);return rr;} |
---|
[8b4b7b] | 1131 | r=rComp_lex(p1,p2,i); |
---|
[4301a7] | 1132 | if (r!=0) { assume(r==rr);return rr;} |
---|
[e95eaa7] | 1133 | break; |
---|
| 1134 | |
---|
| 1135 | case ringorder_S: |
---|
| 1136 | assume(0); |
---|
[4301a7] | 1137 | return rr; |
---|
[e95eaa7] | 1138 | |
---|
| 1139 | case ringorder_s: |
---|
| 1140 | /* ro_syz */ |
---|
[4301a7] | 1141 | // shall not appear: |
---|
| 1142 | assume(0); |
---|
| 1143 | #if 0 |
---|
[e95eaa7] | 1144 | if ((pGetComp(p1) > pDBsyzComp) && (pGetComp(p2) > pDBsyzComp)) break; |
---|
| 1145 | if ((pGetComp(p1) <= pDBsyzComp) && (pGetComp(p2) <= pDBsyzComp)) break; |
---|
[4301a7] | 1146 | if (pGetComp(p1) <= pDBsyzComp) { assume(rr==1); return rr;} |
---|
| 1147 | /* if (pGetComp(p2) <= pDBsyzComp) */ { assume (rr== -1); return rr; } |
---|
| 1148 | #endif |
---|
| 1149 | return rr; |
---|
[e95eaa7] | 1150 | |
---|
| 1151 | case ringorder_unspec: |
---|
| 1152 | case ringorder_no: |
---|
| 1153 | default: |
---|
| 1154 | Print("undef. ringorder used\n"); |
---|
| 1155 | break; |
---|
| 1156 | } |
---|
| 1157 | } |
---|
| 1158 | return rr; |
---|
| 1159 | } |
---|
| 1160 | #endif |
---|
[51c163] | 1161 | #endif // PDEBUG |
---|
| 1162 | |
---|
[a9a7be] | 1163 | static unsigned long GetBitFields(Exponent_t e, |
---|
| 1164 | unsigned int s, unsigned int n) |
---|
[52a5e8] | 1165 | { |
---|
| 1166 | unsigned int i = 0, ev = 0; |
---|
| 1167 | assume(n > 0 && s < BIT_SIZEOF_LONG); |
---|
| 1168 | do |
---|
| 1169 | { |
---|
| 1170 | assume(s+i < BIT_SIZEOF_LONG); |
---|
| 1171 | if (e > (Exponent_t) i) ev |= Sy_bit(s+i); |
---|
| 1172 | else break; |
---|
| 1173 | i++; |
---|
| 1174 | } |
---|
| 1175 | while (i < n); |
---|
| 1176 | return ev; |
---|
| 1177 | } |
---|
| 1178 | |
---|
[a9a7be] | 1179 | // Short Exponent Vectors are used for fast divisibility tests |
---|
| 1180 | // ShortExpVectors "squeeze" an exponent vector into one word as follows: |
---|
[b7b08c] | 1181 | // Let n = BIT_SIZEOF_LONG / pVariables. |
---|
| 1182 | // If n == 0 (i.e. pVariables > BIT_SIZE_OF_LONG), let m == the number |
---|
| 1183 | // of non-zero exponents. If (m>BIT_SIZEOF_LONG), then sev = ~0, else |
---|
| 1184 | // first m bits of sev are set to 1. |
---|
| 1185 | // Otherwise (i.e. pVariables <= BIT_SIZE_OF_LONG) |
---|
[a9a7be] | 1186 | // represented by a bit-field of length n (resp. n+1 for some |
---|
| 1187 | // exponents). If the value of an exponent is greater or equal to n, then |
---|
| 1188 | // all of its respective n bits are set to 1. If the value of an exponent |
---|
| 1189 | // is smaller than n, say m, then only the first m bits of the respective |
---|
[87bef42] | 1190 | // n bits are set to 1, the others are set to 0. |
---|
[b7b08c] | 1191 | // This way, we have: |
---|
[a9a7be] | 1192 | // exp1 / exp2 ==> (ev1 & ~ev2) == 0, i.e., |
---|
| 1193 | // if (ev1 & ~ev2) then exp1 does not divide exp2 |
---|
[52a5e8] | 1194 | unsigned long pGetShortExpVector(poly p) |
---|
| 1195 | { |
---|
[b7b08c] | 1196 | assume(p != NULL); |
---|
| 1197 | if (p == NULL) return 0; |
---|
[52a5e8] | 1198 | unsigned long ev = 0; // short exponent vector |
---|
| 1199 | unsigned int n = BIT_SIZEOF_LONG / pVariables; // number of bits per exp |
---|
| 1200 | unsigned int m1; // highest bit which is filled with (n+1) |
---|
[a9a7be] | 1201 | unsigned int i = 0, j=1; |
---|
| 1202 | |
---|
| 1203 | if (n == 0) |
---|
[52a5e8] | 1204 | { |
---|
[b7b08c] | 1205 | for (; j<=(unsigned long) pVariables; j++) |
---|
| 1206 | { |
---|
| 1207 | if (pGetExp(p,j) > 0) i++; |
---|
| 1208 | if (i == BIT_SIZEOF_LONG) break; |
---|
| 1209 | } |
---|
| 1210 | ev = (unsigned long) ~0 >> ((unsigned long) (BIT_SIZEOF_LONG - i)); |
---|
| 1211 | return ev; |
---|
[52a5e8] | 1212 | } |
---|
| 1213 | else |
---|
| 1214 | { |
---|
| 1215 | m1 = (n+1)*(BIT_SIZEOF_LONG - n*pVariables); |
---|
| 1216 | } |
---|
[a9a7be] | 1217 | |
---|
[52a5e8] | 1218 | n++; |
---|
| 1219 | while (i<m1) |
---|
| 1220 | { |
---|
[a9a7be] | 1221 | ev |= GetBitFields(pGetExp(p, j), i, n); |
---|
[52a5e8] | 1222 | i += n; |
---|
| 1223 | j++; |
---|
| 1224 | } |
---|
| 1225 | |
---|
[a9a7be] | 1226 | n--; |
---|
| 1227 | while (i<BIT_SIZEOF_LONG) |
---|
| 1228 | { |
---|
| 1229 | ev |= GetBitFields(pGetExp(p, j), i, n); |
---|
| 1230 | i += n; |
---|
| 1231 | j++; |
---|
| 1232 | } |
---|
[52a5e8] | 1233 | return ev; |
---|
| 1234 | } |
---|
[51c163] | 1235 | |
---|
[b7b08c] | 1236 | #ifdef PDIV_DEBUG |
---|
| 1237 | static int pDivisibleBy_number = 1; |
---|
| 1238 | static int pDivisibleBy_FALSE = 1; |
---|
| 1239 | static int pDivisibleBy_ShortFalse = 1; |
---|
| 1240 | static int pDivisibleBy_Null = 1; |
---|
[a9a7be] | 1241 | BOOLEAN pDBShortDivisibleBy(poly p1, unsigned long sev_1, |
---|
[87bef42] | 1242 | poly p2, unsigned long not_sev_2, |
---|
[a9a7be] | 1243 | char* f, int l) |
---|
| 1244 | { |
---|
[b7b08c] | 1245 | if (sev_1 != 0 && pGetShortExpVector(p1) != sev_1) |
---|
[a9a7be] | 1246 | { |
---|
[87bef42] | 1247 | Warn("sev1 is %o but should be %o in %s:%d\n", sev_1, |
---|
[a9a7be] | 1248 | pGetShortExpVector(p1), f, l); |
---|
| 1249 | assume(0); |
---|
| 1250 | } |
---|
| 1251 | if (~ pGetShortExpVector(p2) != not_sev_2) |
---|
| 1252 | { |
---|
[87bef42] | 1253 | Warn("not_sev2 is %o but should be %o in %s:%d\n", not_sev_2, |
---|
[a9a7be] | 1254 | ~ pGetShortExpVector(p2), f, l); |
---|
| 1255 | assume(0); |
---|
| 1256 | } |
---|
[b7b08c] | 1257 | if (sev_1 == 0) pDivisibleBy_Null++; |
---|
| 1258 | pDivisibleBy_number++; |
---|
| 1259 | BOOLEAN ret = pDivisibleBy(p1, p2); |
---|
| 1260 | if (! ret) pDivisibleBy_FALSE++; |
---|
[a9a7be] | 1261 | if (sev_1 & not_sev_2) |
---|
| 1262 | { |
---|
[b7b08c] | 1263 | pDivisibleBy_ShortFalse++; |
---|
| 1264 | if (ret) |
---|
[a9a7be] | 1265 | { |
---|
[4f011a] | 1266 | Warn("p1 divides p2, but sev's are wrong in %s:%d\n", f, l); |
---|
[a9a7be] | 1267 | assume(0); |
---|
| 1268 | } |
---|
| 1269 | } |
---|
[b7b08c] | 1270 | return ret; |
---|
| 1271 | } |
---|
| 1272 | |
---|
| 1273 | void pPrintDivisbleByStat() |
---|
| 1274 | { |
---|
| 1275 | Print("#Tests: %d; #FALSE %d(%d); #SHORT %d(%d) #NULL:%d(%d)\n", |
---|
[87bef42] | 1276 | pDivisibleBy_number, |
---|
[b7b08c] | 1277 | pDivisibleBy_FALSE, pDivisibleBy_FALSE*100/pDivisibleBy_number, |
---|
| 1278 | pDivisibleBy_ShortFalse, pDivisibleBy_ShortFalse*100/pDivisibleBy_FALSE, |
---|
| 1279 | pDivisibleBy_Null, pDivisibleBy_Null*100/pDivisibleBy_number); |
---|
| 1280 | |
---|
[a9a7be] | 1281 | } |
---|
| 1282 | #endif |
---|
| 1283 | |
---|
[51c163] | 1284 | #endif // POLYS_IMPL_CC |
---|