Changeset 88c726 in git for dyn_modules/callgfanlib/initialReduction.cc
- Timestamp:
- Dec 10, 2013, 6:58:00 PM (10 years ago)
- Branches:
- (u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
- Children:
- 5222d2aaa879330d2c8daa004f176ee41aed4b78
- Parents:
- edb60f7c79d1b2f53a09d2650752f01a8d8da321
- git-author:
- Yue Ren <ren@mathematik.uni-kl.de>2013-12-10 18:58:00+01:00
- git-committer:
- Yue Ren <ren@mathematik.uni-kl.de>2015-02-06 13:47:02+01:00
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
dyn_modules/callgfanlib/initialReduction.cc
redb60f r88c726 1 1 #include <kernel/polys.h> 2 #include <Singular/ipid.h> 3 2 4 #include <libpolys/polys/monomials/p_polys.h> 3 5 #include <singularWishlist.h> 4 6 5 #include <Singular/ipid.h> 7 #include <map> 8 #include <set> 6 9 7 10 /*** … … 95 98 96 99 /*** 97 * returns, if it exists, a pointer to the first term in g98 * that is divisible by (the leading monomial of) m or, if it does not exist, the NULL pointer.99 * if g is homogeneous in x with the same degree as m,100 * then it returns the first term with the same monomial in x as m,101 * if the t-degree of the term is higher than the t-degree of m, or NULL otherwise.102 **/103 static inline poly firstTermDivisibleBy(const poly g, const poly m)104 {105 poly gCache = NULL;106 for (gCache=g; gCache; pIter(gCache))107 if (p_LeadmonomDivisibleBy(m,gCache,currRing)) break;108 return gCache;109 }110 111 112 /***113 100 * reduces h initially with respect to g, 114 101 * returns false if h was initially reduced in the first place, … … 118 105 bool reduceInitially(poly &h, const poly g) 119 106 { 120 p oly hCache;107 pTest(h); pTest(g); poly hCache; 121 108 for (hCache=h; hCache; pIter(hCache)) 122 109 if (p_LeadmonomDivisibleBy(g,hCache,currRing)) break; … … 130 117 p_SetExp(hAlphaT,i,0,currRing); 131 118 p_Setm(hAlphaT,currRing); pTest(hAlphaT); 132 h = p_Add_q(p_Mult_nn(h,gAlpha,currRing), 133 p_Neg(p_Mult_q(p_Copy(g,currRing),hAlphaT,currRing),currRing), 134 currRing); 119 poly q1 = p_Mult_nn(h,gAlpha,currRing); pTest(q1); 120 poly q2 = p_Mult_q(p_Copy(g,currRing),hAlphaT,currRing); pTest(q2); 121 q2 = p_Neg(q2,currRing); pTest(q2); 122 h = p_Add_q(q1,q2,currRing); 135 123 pTest(h); 136 124 return true; … … 173 161 174 162 175 static inline void sortByLeadmonom(ideal I)176 {177 poly cache; int i, n=IDELEMS(I), newn;178 do179 {180 newn=0;181 for (i=1; i<n; i++)182 {183 if (pLmCmp(I->m[i-1],I->m[i])<0)184 {185 cache=I->m[i-1];186 I->m[i-1]=I->m[i];187 I->m[i]=cache;188 newn = i;189 }190 }191 n=newn;192 } while(n);193 }194 195 196 163 /*** 197 164 * reduces I initially with respect to itself and with respect to p-t. … … 201 168 bool reduceInitially(ideal I, const number p) 202 169 { 203 poly cache; int i,j,n=IDELEMS(I);170 int m=idSize(I),n=m; poly cache; 204 171 do 205 172 { 206 j=0;207 for (i =1; i<n; i++)208 { 209 if (p LmCmp(I->m[i-1],I->m[i])<0)173 int j=0; 174 for (int i=1; i<n; i++) 175 { 176 if (p_LmCmp(I->m[i-1],I->m[i],currRing)<0) 210 177 { 211 178 cache=I->m[i-1]; … … 217 184 n=j; 218 185 } while(n); 219 for (i =1; i<IDELEMS(I); i++)186 for (int i=1; i<m; i++) 220 187 if (pReduce(I->m[i],p)) return true; 221 188 … … 223 190 * the first pass. removing terms with the same monomials in x as lt(g_i) out of g_j for i<j 224 191 **/ 225 for (i =0; i<IDELEMS(I)-1; i++)226 for ( j=i+1; j<IDELEMS(I); j++)192 for (int i=0; i<m-1; i++) 193 for (int j=i+1; j<m; j++) 227 194 if (reduceInitially(I->m[j], I->m[i]) && pReduce(I->m[j],p)) return true; 228 195 … … 230 197 * the second pass. removing terms divisible by lt(g_j) out of g_i for i<j 231 198 **/ 232 for (i =0; i<IDELEMS(I)-1; i++)233 for ( j=i+1; j<IDELEMS(I); j++)199 for (int i=0; i<m-1; i++) 200 for (int j=i+1; j<m; j++) 234 201 if (reduceInitially(I->m[i], I->m[j]) && pReduce(I->m[i],p)) return true; 235 202 return false; … … 274 241 * assumes that I was already sorted and initially reduced in the first place 275 242 **/ 276 bool reduceInitially(ideal I, const number p, poly g) 277 { 278 pEnlargeSet(&(I->m),IDELEMS(I),1); 279 IDELEMS(I)++; int i,j,k; 280 for (j=IDELEMS(I)-1; j>0; j--) 281 { 282 if (pLmCmp(I->m[j-1],g)<0) 243 bool reduceInitially(ideal I, const number p, const poly g) 244 { 245 int n=idSize(I); 246 idInsertPoly(I,g); 247 poly cache; n++; int j; 248 for (j=n-1; j>0; j--) 249 { 250 if (p_LmCmp(I->m[j], I->m[j-1],currRing)>0) 251 { 252 cache = I->m[j]; 283 253 I->m[j] = I->m[j-1]; 254 I->m[j-1] = cache; 255 } 284 256 else 285 {286 I->m[j] = g;287 257 break; 288 } 289 } 290 if (j<1) I->m[0] = g; 258 } 291 259 292 260 /*** … … 294 262 * removing terms with the same monomials in x as lt(g_j) out of g_k for j<k 295 263 **/ 296 for (i =0; i<j; i++)264 for (int i=0; i<j; i++) 297 265 if (reduceInitially(I->m[j], I->m[i]) && pReduce(I->m[j],p)) return true; 298 for ( k=j+1; k<IDELEMS(I); k++)266 for (int k=j+1; k<n; k++) 299 267 if (reduceInitially(I->m[k], I->m[j]) && pReduce(I->m[k],p)) return true; 300 268 … … 303 271 * removing terms divisible by lt(g_k) out of g_j for j<k 304 272 **/ 305 for (i =0; i<j; i++)306 for ( k=j; k<IDELEMS(I); k++)273 for (int i=0; i<j; i++) 274 for (int k=j; k<n; k++) 307 275 if (reduceInitially(I->m[i], I->m[k]) && pReduce(I->m[i],p)) return true; 308 for ( k=j+1; k<IDELEMS(I); k++)276 for (int k=j+1; k<n; k++) 309 277 if (reduceInitially(I->m[j],I->m[k]) && pReduce(I->m[j],p)) return true; 310 278 … … 350 318 } 351 319 #endif //NDEBUG 320 321 322 /*** 323 * reduces H initially with respect to itself, with respect to p-t, 324 * and with respect to G. 325 * assumes that the generators of H are homogeneous in x of the same degree, 326 * assumes that the generators of G are homogeneous in x of lesser degree. 327 **/ 328 bool reduceInitially(ideal H, const number p, const ideal G) 329 { 330 /*** 331 * Step 1: reduce H initially with respect to itself and with respect to p-t 332 **/ 333 if (reduceInitially(H,p)) return true; 334 335 /*** 336 * Step 2: initialize a working list T and an ideal I in which the reductions will take place 337 **/ 338 int m=idSize(H),n=0; 339 ideal I = idInit(m), T = idInit(m); 340 for (int i=0; i<m; i++) 341 { 342 I->m[i]=H->m[i]; 343 if (pNext(H->m[i])) T->m[n++]=H->m[i]; 344 } 345 poly g; int k=n; 346 do 347 { 348 int j=0; 349 for (int i=1; i<k; i++) 350 { 351 if (p_LmCmp(pNext(T->m[i-1]),pNext(T->m[i]),currRing)<0) 352 { 353 g=T->m[i-1]; 354 T->m[i-1]=I->m[i]; 355 T->m[i]=g; 356 j = i; 357 } 358 } 359 k=j; 360 } while(k); 361 362 /*** 363 * Step 3: as long as the working list is not empty, successively reduce terms in it 364 * by adding suitable elements to I and reducing it initially with respect to itself 365 **/ 366 k=idSize(G); 367 while (n) 368 { 369 int i=0; for (; i<k; i++) 370 if (p_LeadmonomDivisibleBy(G->m[i],pNext(T->m[0]),currRing)) break; 371 if (i<k) 372 { 373 g = p_Init(currRing); 374 for (int j=2; j<=currRing->N; j++) 375 p_SetExp(g,j,p_GetExp(pNext(T->m[0]),j,currRing)-p_GetExp(G->m[i],j,currRing),currRing); 376 p_SetCoeff(g,n_Init(1,currRing->cf),currRing); p_Setm(g,currRing); 377 g = p_Mult_q(g,p_Copy(G->m[i],currRing),currRing); 378 reduceInitially(I,p,g); 379 } 380 else 381 pIter(T->m[0]); 382 for (int i=0; i<n;) 383 { 384 if (!pNext(T->m[i])) 385 { 386 for (int j=i; j<n-1; j++) 387 T->m[j]=T->m[j+1]; 388 T->m[--n]=NULL; 389 } 390 else 391 i++; 392 } 393 int l = n; 394 do 395 { 396 int j=0; 397 for (int i=1; i<l; i++) 398 { 399 if (p_LmCmp(pNext(T->m[i-1]),pNext(T->m[i]),currRing)<0) 400 { 401 g=T->m[i-1]; 402 T->m[i-1]=I->m[i]; 403 T->m[i]=g; 404 j = i; 405 } 406 } 407 l=j; 408 } while(l); 409 } 410 411 /*** 412 * Step 4: cleanup, delete all polynomials in I which have been added in Step 3 413 **/ 414 k=idSize(I); 415 for (int i=0; i<k; i++) 416 { 417 for (int j=0; j<m; j++) 418 { 419 if (p_LeadmonomDivisibleBy(H->m[j],I->m[i],currRing)) 420 { 421 I->m[i]=NULL; 422 break; 423 } 424 } 425 } 426 id_Delete(&I,currRing); 427 id_Delete(&T,currRing); 428 return false; 429 } 430 431 432 #ifndef NDEBUG 433 BOOLEAN reduceInitially3(leftv res, leftv args) 434 { 435 leftv u = args; 436 if ((u != NULL) && (u->Typ() == IDEAL_CMD)) 437 { 438 leftv v = u->next; 439 if ((v != NULL) && (v->Typ() == NUMBER_CMD)) 440 { 441 leftv w = v->next; 442 if ((w != NULL) && (w->Typ() == IDEAL_CMD)) 443 { 444 ideal H,G; number p; 445 omUpdateInfo(); 446 Print("usedBytesBefore=%ld\n",om_Info.UsedBytes); 447 H = (ideal) u->CopyD(); 448 p = (number) v->CopyD(); 449 G = (ideal) w->CopyD(); 450 (void) reduceInitially(H,p,G); 451 id_Delete(&H,currRing); 452 id_Delete(&G,currRing); 453 n_Delete(&p,currRing->cf); 454 omUpdateInfo(); 455 Print("usedBytesAfter=%ld\n",om_Info.UsedBytes); 456 H = (ideal) u->CopyD(); 457 p = (number) v->CopyD(); 458 G = (ideal) w->CopyD(); 459 (void) reduceInitially(H,p,G); 460 n_Delete(&p,currRing->cf); 461 id_Delete(&G,currRing); 462 res->rtyp = IDEAL_CMD; 463 res->data = (char*) H; 464 return FALSE; 465 } 466 } 467 } 468 return TRUE; 469 } 470 #endif //NDEBUG 471 472 473 /*** 474 * reduces I initially with respect to itself. 475 * assumes that the generators of I are homogeneous in x and that p-t is in I. 476 **/ 477 bool reduceInitially(ideal I) 478 { 479 /*** 480 * Step 1: split up I into components of same degree in x 481 * the lowest component should only contain p-t 482 **/ 483 std::map<long,ideal> H; int n = idSize(I); 484 for (int i=0; i<n; i++) 485 { 486 long d = 0; 487 for (int j=2; j<=currRing->N; j++) 488 d += p_GetExp(I->m[i],j,currRing); 489 std::map<long,ideal>::iterator it = H.find(d); 490 if (it != H.end()) 491 idInsertPoly(it->second,I->m[i]); 492 else 493 { 494 std::pair<std::map<long,ideal>::iterator,bool> ret; 495 ret = H.insert(std::pair<long,ideal>(d,idInit(16))); 496 idInsertPoly(ret.first->second,I->m[i]); 497 } 498 } 499 500 std::map<long,ideal>::iterator it=H.begin(); 501 ideal Hi = it->second; 502 assume(idSize(Hi)==1); 503 assume(pLength(Hi->m[0])==2 && p_GetExp(Hi->m[0],1,currRing)==0 504 && p_GetExp(Hi->m[0]->next,1,currRing)==1); 505 number p = p_GetCoeff(Hi->m[0],currRing); 506 assume(!n_IsUnit(p,currRing->cf)); 507 idShallowDelete(&it->second); 508 509 /*** 510 * Step 2: reduce each component initially with respect to itself 511 * and all lower components 512 **/ 513 it++; Hi = it->second; n--; 514 if (reduceInitially(Hi,p)) return true; 515 516 ideal G = idInit(n); int m=0; 517 ideal GG = (ideal) omAllocBin(sip_sideal_bin); 518 GG->nrows = 1; GG->rank = 1; GG->m=NULL; 519 520 for (it++; it!=H.end(); it++) 521 { 522 int l=idSize(Hi); int k=l; poly cache; 523 do 524 { 525 int j=0; 526 for (int i=1; i<k; i++) 527 { 528 if (p_GetExp(Hi->m[i-1],1,currRing)<p_GetExp(Hi->m[i],1,currRing)) 529 { 530 cache=Hi->m[i-1]; 531 Hi->m[i-1]=Hi->m[i]; 532 Hi->m[i]=cache; 533 j = i; 534 } 535 } 536 k=j; 537 } while(k); 538 int kG=n-m, kH=0; 539 for (int i=n-m-l; i<n; i++) 540 { 541 if (kG==n) 542 { 543 memcpy(&(G->m[i]),&(Hi->m[kH]),(n-i)*sizeof(poly)); 544 break; 545 } 546 if (p_GetExp(G->m[kG],1,currRing)>p_GetExp(Hi->m[kH],1,currRing)) 547 G->m[i] = G->m[kG++]; 548 else 549 G->m[i] = Hi->m[kH++]; 550 } 551 m += l; IDELEMS(GG) = m; GG->m = &G->m[n-m]; 552 if (reduceInitially(it->second,p,GG)) return true; 553 idShallowDelete(&Hi); Hi = it->second; 554 } 555 idShallowDelete(&Hi); 556 557 omFreeBin((ADDRESS)GG, sip_sideal_bin); idShallowDelete(&G); 558 return false; 559 } 560 561 562 #ifndef NDEBUG 563 BOOLEAN reduceInitially4(leftv res, leftv args) 564 { 565 leftv u = args; 566 if ((u != NULL) && (u->Typ() == IDEAL_CMD)) 567 { 568 ideal I; 569 omUpdateInfo(); 570 Print("usedBytesBefore=%ld\n",om_Info.UsedBytes); 571 I = (ideal) u->CopyD(); 572 (void) reduceInitially(I); 573 id_Delete(&I,currRing); 574 omUpdateInfo(); 575 Print("usedBytesAfter=%ld\n",om_Info.UsedBytes); 576 I = (ideal) u->CopyD(); 577 (void) reduceInitially(I); 578 res->rtyp = IDEAL_CMD; 579 res->data = (char*) I; 580 return FALSE; 581 } 582 return TRUE; 583 } 584 #endif 585 586 587 BOOLEAN reduceInitially(leftv res, leftv args) 588 { 589 leftv u = args; 590 if ((u != NULL) && (u->Typ() == IDEAL_CMD)) 591 { 592 ideal I = (ideal) u->CopyD(); 593 (void) reduceInitially(I); 594 res->rtyp = IDEAL_CMD; 595 res->data = (char*) I; 596 return FALSE; 597 } 598 return TRUE; 599 }
Note: See TracChangeset
for help on using the changeset viewer.