Changeset fa167e in git for Singular/dyn_modules/gfanlib/tropical.cc
 Timestamp:
 Dec 1, 2013, 6:20:33 PM (9 years ago)
 Branches:
 (u'jengelhdatetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
 Children:
 7a4e77045299daf5ac74ab3774c5e1b784ff4e94
 Parents:
 9bee534198ac4b761325be1435e67598941aac19
 gitauthor:
 Yue Ren <ren@mathematik.unikl.de>20131201 18:20:33+01:00
 gitcommitter:
 Yue Ren <ren@mathematik.unikl.de>20150206 13:47:02+01:00
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/dyn_modules/gfanlib/tropical.cc
r9bee53 rfa167e 5 5 #include <bbcone.h> 6 6 7 #include <initialReduction.h> 7 8 8 9 poly initial(poly p) … … 236 237 237 238 238 / ***239 * replaces coefficients in g of lowest degree in t240 * divisible by p with a suitable power of t241 **/242 static bool preduce(poly g, const number p)243 {244 // go along g and look for relevant terms.245 // monomials in x which have been checked are tracked in done.246 // because we assume the leading coefficient of g is 1,247 // the leading term does not need to be considered.248 poly done = p_LmInit(g,currRing);249 p_SetExp(done,1,0,currRing); p_SetCoeff(done,n_Init(1,currRing>cf),currRing);250 p_Setm(done,currRing);251 poly doneEnd = done;252 poly doneCache;253 number dnumber; long d;254 poly subst; number ppower;255 while(pNext(g))256 {257 // check whether next term needs to be reduced:258 // first, check whether monomial in x has come up yet259 for (doneCache=done; doneCache; pIter(doneCache))260 {261 if (p_LmDivisibleBy(doneCache,pNext(g),currRing))262 break;263 }264 if (!doneCache) // if for loop did not terminate prematurely,265 // then the monomial in x is new266 {267 // second, check whether coefficient is divisible by p268 if (n_DivBy(p_GetCoeff(pNext(g),currRing>cf),p,currRing>cf))269 {270 // reduce the term with respect to pt:271 // divide coefficient by p, remove old term,272 // add t multiple of old term273 dnumber = n_Div(p_GetCoeff(pNext(g),currRing>cf),p,currRing>cf);274 d = n_Int(dnumber,currRing>cf);275 n_Delete(&dnumber,currRing>cf);276 if (!d)277 {278 p_Delete(&done,currRing);279 WerrorS("initialReduction: overflow in a texponent");280 return true;281 }282 subst=p_LmInit(pNext(g),currRing);283 if (d>0)284 {285 p_AddExp(subst,1,d,currRing);286 p_SetCoeff(subst,n_Init(1,currRing>cf),currRing);287 }288 else289 {290 p_AddExp(subst,1,d,currRing);291 p_SetCoeff(subst,n_Init(1,currRing>cf),currRing);292 }293 p_Setm(subst,currRing); pTest(subst);294 pNext(g)=p_LmDeleteAndNext(pNext(g),currRing);295 pNext(g)=p_Add_q(pNext(g),subst,currRing);296 pTest(pNext(g));297 }298 // either way, add monomial in x to done299 pNext(doneEnd)=p_LmInit(pNext(g),currRing);300 pIter(doneEnd);301 p_SetExp(doneEnd,1,0,currRing);302 p_SetCoeff(doneEnd,n_Init(1,currRing>cf),currRing);303 p_Setm(doneEnd,currRing);304 }305 pIter(g);306 }307 p_Delete(&done,currRing);308 return false;309 }310 311 312 #ifndef NDEBUG313 BOOLEAN preduce(leftv res, leftv args)314 {315 leftv u = args;316 if ((u != NULL) && (u>Typ() == POLY_CMD))317 {318 poly g; bool b;319 number p = n_Init(3,currRing>cf);320 omUpdateInfo();321 Print("usedBytes=%ld\n",om_Info.UsedBytes);322 g = (poly) u>CopyD();323 b = preduce(g,p);324 p_Delete(&g,currRing);325 if (b) return TRUE;326 omUpdateInfo();327 Print("usedBytes=%ld\n",om_Info.UsedBytes);328 n_Delete(&p,currRing>cf);329 res>rtyp = NONE;330 res>data = NULL;331 return FALSE;332 }333 return TRUE;334 }335 #endif //NDEBUG336 337 338 / ***339 * Returns the highest term in g with the matching xmonomial to m340 * or, if it does not exist, the NULL pointer341 **/342 static poly highestMatchingX(poly g, const poly m)343 {344 pTest(g); pTest(m);345 poly xalpha=p_LmInit(m,currRing);346 347 // go along g and find the first term348 // with the same monomial in x as xalpha349 while (g)350 {351 p_SetExp(xalpha,1,p_GetExp(g,1,currRing),currRing);352 p_Setm(xalpha,currRing);353 if (p_LmEqual(g,xalpha,currRing))354 break;355 pIter(g);356 }357 358 // gCache now either points at the wanted term359 // or is NULL360 p_Delete(&xalpha,currRing);361 pTest(g);362 return g;363 }364 365 366 / ***367 * Given g with lm(g)=t^\beta x^\alpha, returns g_\alpha368 **/369 static poly powerSeriesCoeff(poly g)370 {371 // the first term obviously is part of our output372 // so we copy it...373 poly xalpha=p_LmInit(g,currRing);374 poly coeff=p_One(currRing);375 p_SetCoeff(coeff,n_Copy(p_GetCoeff(g,currRing),currRing>cf),currRing);376 p_SetExp(coeff,1,p_GetExp(g,1,currRing),currRing);377 p_Setm(coeff,currRing);378 379 // ..and proceed with the remaining terms,380 // appending the relevant terms to coeff via coeffCache381 poly coeffCache=coeff;382 pIter(g);383 while (g)384 {385 p_SetExp(xalpha,1,p_GetExp(g,1,currRing),currRing);386 p_Setm(xalpha,currRing);387 if (p_LmEqual(g,xalpha,currRing))388 {389 pNext(coeffCache)=p_Init(currRing);390 pIter(coeffCache);391 p_SetCoeff(coeffCache,n_Copy(p_GetCoeff(g,currRing),currRing>cf),currRing);392 p_SetExp(coeffCache,1,p_GetExp(g,1,currRing),currRing);393 p_Setm(coeffCache,currRing);394 }395 pIter(g);396 }397 398 p_Delete(&xalpha,currRing);399 return coeff;400 }401 402 403 / ***404 * divides g by t^b knowing that each term of g405 * is divisible by t^b, i.e. no divisibility checks406 * needed407 **/408 static void divideByT(poly g, const long b)409 {410 while (g)411 {412 p_SetExp(g,1,p_GetExp(g,1,currRing)b,currRing);413 p_Setm(g,currRing);414 pIter(g);415 }416 }417 418 419 static void divideByGcd(poly &g)420 {421 // first determine all g_\alpha422 // as alpha runs over all exponent vectors423 // of monomials in x occuring in g424 425 // gAlpha will store all g_\alpha,426 // the first term will, for comparison purposes for now,427 // also keep their monomial in x.428 // we assume that the weight on the x are positive429 // so that the x won't make the monomial smaller430 ideal gAlphaFront = idInit(pLength(g));431 gAlphaFront>m[0] = p_Head(g,currRing);432 p_SetExp(gAlphaFront>m[0],1,0,currRing);433 p_Setm(gAlphaFront>m[0],currRing);434 ideal gAlphaEnd = idInit(pLength(g));435 gAlphaEnd>m[0] = gAlphaFront>m[0];436 unsigned long gAlpha_sev[pLength(g)];437 gAlpha_sev[0] = p_GetShortExpVector(g,currRing);438 long beta[pLength(g)];439 beta[0] = p_GetExp(g,1,currRing);440 int count=0;441 442 poly current = pNext(g); unsigned long current_not_sev;443 int i;444 while (current)445 {446 // see whether the monomial in x of current already came up447 // since everything is homogeneous in x and the ordering is local in t,448 // this comes down to a divisibility test in two stages449 current_not_sev = ~p_GetShortExpVector(current,currRing);450 for(i=0; i<=count; i++)451 {452 // first stage using short exponent vectors453 // second stage a proper test454 if (p_LmShortDivisibleBy(gAlphaFront>m[i],gAlpha_sev[i],current,current_not_sev, currRing)455 && p_LmDivisibleBy(gAlphaFront>m[i],current,currRing))456 {457 // if it already occured, add the term to the respective entry458 // without the x part459 pNext(gAlphaEnd>m[i])=p_Init(currRing);460 pIter(gAlphaEnd>m[i]);461 p_SetExp(gAlphaEnd>m[i],1,p_GetExp(current,1,currRing)beta[i],currRing);462 p_SetCoeff(gAlphaEnd>m[i],n_Copy(p_GetCoeff(current,currRing),currRing>cf),currRing);463 p_Setm(gAlphaEnd>m[i],currRing);464 break;465 }466 }467 if (i>count)468 {469 // if it is new, create a new entry for the term470 // including the monomial in x471 count++;472 gAlphaFront>m[count] = p_Head(current,currRing);473 beta[count] = p_GetExp(current,1,currRing);474 gAlphaEnd>m[count] = gAlphaFront>m[count];475 gAlpha_sev[count] = p_GetShortExpVector(current,currRing);476 }477 478 pIter(current);479 }480 481 if (count)482 {483 // now remove all the x in the leading terms484 // so that gAlpha only contais terms in t485 int j; long d;486 for (i=0; i<=count; i++)487 {488 for (j=2; j<=currRing>N; j++)489 p_SetExp(gAlphaFront>m[i],j,0,currRing);490 p_Setm(gAlphaFront>m[i],currRing);491 gAlphaEnd>m[i]=NULL;492 }493 494 // now compute the gcd over all g_\alpha495 // and set the input to null as they are deleted in the process496 poly gcd = singclap_gcd(gAlphaFront>m[0],gAlphaFront>m[1],currRing);497 gAlphaFront>m[0] = NULL;498 gAlphaFront>m[1] = NULL;499 for(i=2; i<=count; i++)500 {501 gcd = singclap_gcd(gcd,gAlphaFront>m[i],currRing);502 gAlphaFront>m[i] = NULL;503 }504 // divide g by the gcd505 poly h = singclap_pdivide(g,gcd,currRing);506 p_Delete(&gcd,currRing);507 p_Delete(&g,currRing);508 g = h;509 510 id_Delete(&gAlphaFront,currRing);511 id_Delete(&gAlphaEnd,currRing);512 }513 else514 {515 // if g contains only one monomial in x,516 // then we can get rid of all the higher t517 p_Delete(&g,currRing);518 g = gAlphaFront>m[0];519 pIter(gAlphaFront>m[0]);520 pNext(g)=NULL;521 gAlphaEnd>m[0] = NULL;522 id_Delete(&gAlphaFront,currRing);523 id_Delete(&gAlphaEnd,currRing);524 }525 }526 527 528 / ***529 * 1. For each \alpha\in\NN^n, changes (c_\beta t^\beta + c_{\beta+1} t^{\beta+1} + ...)530 * to (c_\beta + c_{\beta+1}*p + ...) t^\beta531 * 2. Computes the gcd over all (c_\beta + c_{\beta+1}*p + ...) and divides g by it532 **/533 static void simplify(poly g, const number p)534 {535 // go along g and look for relevant terms.536 // monomials in x which have been checked are tracked in done.537 poly done = p_LmInit(g,currRing);538 p_SetCoeff(done,n_Init(1,currRing>cf),currRing);539 p_Setm(done,currRing);540 poly doneEnd = done;541 poly doneCurrent;542 543 poly subst; number substCoeff, substCoeffCache;544 unsigned long d;545 546 poly gCurrent = g;547 while(pNext(gCurrent))548 {549 // check whether next term needs to be reduced:550 // first, check whether monomial in x has come up yet551 for (doneCurrent=done; doneCurrent; pIter(doneCurrent))552 {553 if (p_LmDivisibleBy(doneCurrent,pNext(gCurrent),currRing))554 break;555 }556 // if the monomial in x already occured, then we want to replace557 // as many t with p as theoretically feasable558 if (doneCurrent)559 {560 // suppose pNext(gCurrent)=3*t5x and doneCurrent=t3x561 // then we want to replace pNext(gCurrent) with 3p2*t3x562 // subst = ?*t3x563 subst = p_LmInit(doneCurrent,currRing);564 // substcoeff = p2565 n_Power(p,p_GetExp(subst,1,currRing)p_GetExp(doneCurrent,1,currRing),&substCoeff,currRing>cf);566 // substcoeff = 3p2567 n_InpMult(substCoeff,p_GetCoeff(pNext(gCurrent),currRing),currRing>cf);568 // subst = 3p2*t3x569 p_SetCoeff(subst,substCoeff,currRing);570 p_Setm(subst,currRing); pTest(subst);571 572 // g = g  pNext(gCurrent) + subst573 pNext(gCurrent)=p_LmDeleteAndNext(pNext(gCurrent),currRing);574 g=p_Add_q(g,subst,currRing);575 pTest(pNext(gbeginning));576 }577 else578 {579 // if the monomial in x is brand new,580 // then we check whether the coefficient is divisible by p581 if (n_DivBy(p_GetCoeff(pNext(gCurrent),currRing>cf),p,currRing>cf))582 {583 // reduce the term with respect to pt:584 // suppose pNext(gCurrent)=4p3*tx585 // then we want to replace it with 4*t4x586 // divide 4p3 repeatedly by p until it is not divisible anymore,587 // keeping track on the final value 4588 // and the number of times it has been divided 3589 substCoeff = n_Div(p_GetCoeff(pNext(gCurrent),currRing>cf),p,currRing>cf);590 d = 1;591 while (n_DivBy(substCoeff,p,currRing>cf))592 {593 substCoeffCache = n_Div(p_GetCoeff(pNext(gCurrent),currRing>cf),p,currRing>cf);594 n_Delete(&substCoeff,currRing>cf);595 substCoeff = substCoeffCache;596 d++;597 assume(d>0); // check for overflow598 }599 600 // subst = ?*tx601 subst=p_LmInit(pNext(gCurrent),currRing);602 // subst = ?*t4x603 p_AddExp(subst,1,d,currRing);604 // subst = 4*t4x605 p_SetCoeff(subst,substCoeffCache,currRing);606 p_Setm(subst,currRing); pTest(subst);607 608 // g = g  pNext(gCurrent) + subst609 pNext(gCurrent)=p_LmDeleteAndNext(pNext(gCurrent),currRing);610 pNext(gCurrent)=p_Add_q(pNext(gCurrent),subst,currRing);611 pTest(pNext(gCurrent));612 }613 614 // either way, add monomial in x with minimal t to done615 pNext(doneEnd)=p_LmInit(pNext(gCurrent),currRing);616 pIter(doneEnd);617 p_SetCoeff(doneEnd,n_Init(1,currRing>cf),currRing);618 p_Setm(doneEnd,currRing);619 }620 pIter(gCurrent);621 }622 p_Delete(&done,currRing);623 }624 625 626 #ifndef NDEBUG627 BOOLEAN divideByGcd(leftv res, leftv args)628 {629 leftv u = args;630 if ((u != NULL) && (u>Typ() == POLY_CMD))631 {632 poly g;633 omUpdateInfo();634 Print("usedBytes1=%ld\n",om_Info.UsedBytes);635 g = (poly) u>CopyD();636 divideByGcd(g);637 p_Delete(&g,currRing);638 omUpdateInfo();639 Print("usedBytes1=%ld\n",om_Info.UsedBytes);640 res>rtyp = NONE;641 res>data = NULL;642 return FALSE;643 }644 return TRUE;645 }646 #endif //NDEBUG647 648 649 BOOLEAN initialReduction(leftv res, leftv args)650 {651 leftv u = args;652 if ((u != NULL) && (u>Typ() == IDEAL_CMD))653 {654 leftv v = u>next;655 if (v == NULL)656 {657 ideal I = (ideal) u>Data();658 659 /***660 * Suppose I=g_0,...,g_n and g_0=pt.661 * Step 1: sort elements g_1,...,g_{n1} such that lm(g_1)>...>lm(g_{n1})662 * (the following algorithm is a bubble sort)663 **/664 sortingLaterGenerators(I);665 666 /***667 * Step 2: replace coefficient of terms of lowest tdegree divisible by p with t668 **/669 number p = p_GetCoeff(I>m[0],currRing);670 for (int i=1; i<IDELEMS(I); i++)671 {672 if (preduce(I>m[i],p))673 return TRUE;674 }675 676 /***677 * Step 3: the first pass. removing terms with the same monomials in x as lt(g_i)678 * out of g_j for i<j679 **/680 int i,j;681 poly uniti, unitj;682 for (i=1; i<IDELEMS(I)1; i++)683 {684 for (j=i+1; j<IDELEMS(I); j++)685 {686 unitj = highestMatchingX(I>m[j],I>m[i]);687 if (unitj)688 {689 unitj = powerSeriesCoeff(unitj);690 divideByT(unitj,p_GetExp(I>m[i],1,currRing));691 uniti = powerSeriesCoeff(I>m[i]);692 divideByT(uniti,p_GetExp(I>m[i],1,currRing));693 pTest(unitj); pTest(uniti); pTest(I>m[j]); pTest(I>m[i]);694 I>m[j]=p_Add_q(p_Mult_q(uniti,I>m[j],currRing),695 p_Neg(p_Mult_q(unitj,p_Copy(I>m[i],currRing),currRing),currRing),696 currRing);697 divideByGcd(I>m[j]);698 }699 }700 }701 for (int i=1; i<IDELEMS(I); i++)702 {703 if (preduce(I>m[i],p))704 return TRUE;705 }706 707 /***708 * Step 4: the second pass. removing terms divisible by lt(g_j) out of g_i for i<j709 **/710 for (i=1; i<IDELEMS(I)1; i++)711 {712 for (j=i+1; j<IDELEMS(I); j++)713 {714 uniti = highestMatchingX(I>m[i],I>m[j]);715 if (uniti && p_GetExp(uniti,1,currRing)>=p_GetExp(I>m[j],1,currRing))716 {717 uniti = powerSeriesCoeff(uniti);718 divideByT(uniti,p_GetExp(I>m[j],1,currRing));719 unitj = powerSeriesCoeff(I>m[j]);720 divideByT(unitj,p_GetExp(I>m[j],1,currRing));721 I>m[i] = p_Add_q(p_Mult_q(unitj,I>m[i],currRing),722 p_Neg(p_Mult_q(uniti,p_Copy(I>m[j],currRing),currRing),currRing),723 currRing);724 divideByGcd(I>m[j]);725 }726 }727 }728 for (int i=1; i<IDELEMS(I); i++)729 {730 if (preduce(I>m[i],p))731 return TRUE;732 }733 734 res>rtyp = NONE;735 res>data = NULL;736 IDDATA((idhdl)u>data) = (char*) I;737 return FALSE;738 }739 }740 WerrorS("initialReduction: unexpected parameters");741 return TRUE;742 }239 // /*** 240 // * replaces coefficients in g of lowest degree in t 241 // * divisible by p with a suitable power of t 242 // **/ 243 // static bool preduce(poly g, const number p) 244 // { 245 // // go along g and look for relevant terms. 246 // // monomials in x which have been checked are tracked in done. 247 // // because we assume the leading coefficient of g is 1, 248 // // the leading term does not need to be considered. 249 // poly done = p_LmInit(g,currRing); 250 // p_SetExp(done,1,0,currRing); p_SetCoeff(done,n_Init(1,currRing>cf),currRing); 251 // p_Setm(done,currRing); 252 // poly doneEnd = done; 253 // poly doneCache; 254 // number dnumber; long d; 255 // poly subst; number ppower; 256 // while(pNext(g)) 257 // { 258 // // check whether next term needs to be reduced: 259 // // first, check whether monomial in x has come up yet 260 // for (doneCache=done; doneCache; pIter(doneCache)) 261 // { 262 // if (p_LmDivisibleBy(doneCache,pNext(g),currRing)) 263 // break; 264 // } 265 // if (!doneCache) // if for loop did not terminate prematurely, 266 // // then the monomial in x is new 267 // { 268 // // second, check whether coefficient is divisible by p 269 // if (n_DivBy(p_GetCoeff(pNext(g),currRing>cf),p,currRing>cf)) 270 // { 271 // // reduce the term with respect to pt: 272 // // divide coefficient by p, remove old term, 273 // // add t multiple of old term 274 // dnumber = n_Div(p_GetCoeff(pNext(g),currRing>cf),p,currRing>cf); 275 // d = n_Int(dnumber,currRing>cf); 276 // n_Delete(&dnumber,currRing>cf); 277 // if (!d) 278 // { 279 // p_Delete(&done,currRing); 280 // WerrorS("initialReduction: overflow in a texponent"); 281 // return true; 282 // } 283 // subst=p_LmInit(pNext(g),currRing); 284 // if (d>0) 285 // { 286 // p_AddExp(subst,1,d,currRing); 287 // p_SetCoeff(subst,n_Init(1,currRing>cf),currRing); 288 // } 289 // else 290 // { 291 // p_AddExp(subst,1,d,currRing); 292 // p_SetCoeff(subst,n_Init(1,currRing>cf),currRing); 293 // } 294 // p_Setm(subst,currRing); pTest(subst); 295 // pNext(g)=p_LmDeleteAndNext(pNext(g),currRing); 296 // pNext(g)=p_Add_q(pNext(g),subst,currRing); 297 // pTest(pNext(g)); 298 // } 299 // // either way, add monomial in x to done 300 // pNext(doneEnd)=p_LmInit(pNext(g),currRing); 301 // pIter(doneEnd); 302 // p_SetExp(doneEnd,1,0,currRing); 303 // p_SetCoeff(doneEnd,n_Init(1,currRing>cf),currRing); 304 // p_Setm(doneEnd,currRing); 305 // } 306 // pIter(g); 307 // } 308 // p_Delete(&done,currRing); 309 // return false; 310 // } 311 312 313 // #ifndef NDEBUG 314 // BOOLEAN preduce(leftv res, leftv args) 315 // { 316 // leftv u = args; 317 // if ((u != NULL) && (u>Typ() == POLY_CMD)) 318 // { 319 // poly g; bool b; 320 // number p = n_Init(3,currRing>cf); 321 // omUpdateInfo(); 322 // Print("usedBytes=%ld\n",om_Info.UsedBytes); 323 // g = (poly) u>CopyD(); 324 // b = preduce(g,p); 325 // p_Delete(&g,currRing); 326 // if (b) return TRUE; 327 // omUpdateInfo(); 328 // Print("usedBytes=%ld\n",om_Info.UsedBytes); 329 // n_Delete(&p,currRing>cf); 330 // res>rtyp = NONE; 331 // res>data = NULL; 332 // return FALSE; 333 // } 334 // return TRUE; 335 // } 336 // #endif //NDEBUG 337 338 339 // /*** 340 // * Returns the highest term in g with the matching xmonomial to m 341 // * or, if it does not exist, the NULL pointer 342 // **/ 343 // static poly highestMatchingX(poly g, const poly m) 344 // { 345 // pTest(g); pTest(m); 346 // poly xalpha=p_LmInit(m,currRing); 347 348 // // go along g and find the first term 349 // // with the same monomial in x as xalpha 350 // while (g) 351 // { 352 // p_SetExp(xalpha,1,p_GetExp(g,1,currRing),currRing); 353 // p_Setm(xalpha,currRing); 354 // if (p_LmEqual(g,xalpha,currRing)) 355 // break; 356 // pIter(g); 357 // } 358 359 // // gCache now either points at the wanted term 360 // // or is NULL 361 // p_Delete(&xalpha,currRing); 362 // pTest(g); 363 // return g; 364 // } 365 366 367 // /*** 368 // * Given g with lm(g)=t^\beta x^\alpha, returns g_\alpha 369 // **/ 370 // static poly powerSeriesCoeff(poly g) 371 // { 372 // // the first term obviously is part of our output 373 // // so we copy it... 374 // poly xalpha=p_LmInit(g,currRing); 375 // poly coeff=p_One(currRing); 376 // p_SetCoeff(coeff,n_Copy(p_GetCoeff(g,currRing),currRing>cf),currRing); 377 // p_SetExp(coeff,1,p_GetExp(g,1,currRing),currRing); 378 // p_Setm(coeff,currRing); 379 380 // // ..and proceed with the remaining terms, 381 // // appending the relevant terms to coeff via coeffCache 382 // poly coeffCache=coeff; 383 // pIter(g); 384 // while (g) 385 // { 386 // p_SetExp(xalpha,1,p_GetExp(g,1,currRing),currRing); 387 // p_Setm(xalpha,currRing); 388 // if (p_LmEqual(g,xalpha,currRing)) 389 // { 390 // pNext(coeffCache)=p_Init(currRing); 391 // pIter(coeffCache); 392 // p_SetCoeff(coeffCache,n_Copy(p_GetCoeff(g,currRing),currRing>cf),currRing); 393 // p_SetExp(coeffCache,1,p_GetExp(g,1,currRing),currRing); 394 // p_Setm(coeffCache,currRing); 395 // } 396 // pIter(g); 397 // } 398 399 // p_Delete(&xalpha,currRing); 400 // return coeff; 401 // } 402 403 404 // /*** 405 // * divides g by t^b knowing that each term of g 406 // * is divisible by t^b, i.e. no divisibility checks 407 // * needed 408 // **/ 409 // static void divideByT(poly g, const long b) 410 // { 411 // while (g) 412 // { 413 // p_SetExp(g,1,p_GetExp(g,1,currRing)b,currRing); 414 // p_Setm(g,currRing); 415 // pIter(g); 416 // } 417 // } 418 419 420 // static void divideByGcd(poly &g) 421 // { 422 // // first determine all g_\alpha 423 // // as alpha runs over all exponent vectors 424 // // of monomials in x occuring in g 425 426 // // gAlpha will store all g_\alpha, 427 // // the first term will, for comparison purposes for now, 428 // // also keep their monomial in x. 429 // // we assume that the weight on the x are positive 430 // // so that the x won't make the monomial smaller 431 // ideal gAlphaFront = idInit(pLength(g)); 432 // gAlphaFront>m[0] = p_Head(g,currRing); 433 // p_SetExp(gAlphaFront>m[0],1,0,currRing); 434 // p_Setm(gAlphaFront>m[0],currRing); 435 // ideal gAlphaEnd = idInit(pLength(g)); 436 // gAlphaEnd>m[0] = gAlphaFront>m[0]; 437 // unsigned long gAlpha_sev[pLength(g)]; 438 // gAlpha_sev[0] = p_GetShortExpVector(g,currRing); 439 // long beta[pLength(g)]; 440 // beta[0] = p_GetExp(g,1,currRing); 441 // int count=0; 442 443 // poly current = pNext(g); unsigned long current_not_sev; 444 // int i; 445 // while (current) 446 // { 447 // // see whether the monomial in x of current already came up 448 // // since everything is homogeneous in x and the ordering is local in t, 449 // // this comes down to a divisibility test in two stages 450 // current_not_sev = ~p_GetShortExpVector(current,currRing); 451 // for(i=0; i<=count; i++) 452 // { 453 // // first stage using short exponent vectors 454 // // second stage a proper test 455 // if (p_LmShortDivisibleBy(gAlphaFront>m[i],gAlpha_sev[i],current,current_not_sev, currRing) 456 // && p_LmDivisibleBy(gAlphaFront>m[i],current,currRing)) 457 // { 458 // // if it already occured, add the term to the respective entry 459 // // without the x part 460 // pNext(gAlphaEnd>m[i])=p_Init(currRing); 461 // pIter(gAlphaEnd>m[i]); 462 // p_SetExp(gAlphaEnd>m[i],1,p_GetExp(current,1,currRing)beta[i],currRing); 463 // p_SetCoeff(gAlphaEnd>m[i],n_Copy(p_GetCoeff(current,currRing),currRing>cf),currRing); 464 // p_Setm(gAlphaEnd>m[i],currRing); 465 // break; 466 // } 467 // } 468 // if (i>count) 469 // { 470 // // if it is new, create a new entry for the term 471 // // including the monomial in x 472 // count++; 473 // gAlphaFront>m[count] = p_Head(current,currRing); 474 // beta[count] = p_GetExp(current,1,currRing); 475 // gAlphaEnd>m[count] = gAlphaFront>m[count]; 476 // gAlpha_sev[count] = p_GetShortExpVector(current,currRing); 477 // } 478 479 // pIter(current); 480 // } 481 482 // if (count) 483 // { 484 // // now remove all the x in the leading terms 485 // // so that gAlpha only contais terms in t 486 // int j; long d; 487 // for (i=0; i<=count; i++) 488 // { 489 // for (j=2; j<=currRing>N; j++) 490 // p_SetExp(gAlphaFront>m[i],j,0,currRing); 491 // p_Setm(gAlphaFront>m[i],currRing); 492 // gAlphaEnd>m[i]=NULL; 493 // } 494 495 // // now compute the gcd over all g_\alpha 496 // // and set the input to null as they are deleted in the process 497 // poly gcd = singclap_gcd(gAlphaFront>m[0],gAlphaFront>m[1],currRing); 498 // gAlphaFront>m[0] = NULL; 499 // gAlphaFront>m[1] = NULL; 500 // for(i=2; i<=count; i++) 501 // { 502 // gcd = singclap_gcd(gcd,gAlphaFront>m[i],currRing); 503 // gAlphaFront>m[i] = NULL; 504 // } 505 // // divide g by the gcd 506 // poly h = singclap_pdivide(g,gcd,currRing); 507 // p_Delete(&gcd,currRing); 508 // p_Delete(&g,currRing); 509 // g = h; 510 511 // id_Delete(&gAlphaFront,currRing); 512 // id_Delete(&gAlphaEnd,currRing); 513 // } 514 // else 515 // { 516 // // if g contains only one monomial in x, 517 // // then we can get rid of all the higher t 518 // p_Delete(&g,currRing); 519 // g = gAlphaFront>m[0]; 520 // pIter(gAlphaFront>m[0]); 521 // pNext(g)=NULL; 522 // gAlphaEnd>m[0] = NULL; 523 // id_Delete(&gAlphaFront,currRing); 524 // id_Delete(&gAlphaEnd,currRing); 525 // } 526 // } 527 528 529 // /*** 530 // * 1. For each \alpha\in\NN^n, changes (c_\beta t^\beta + c_{\beta+1} t^{\beta+1} + ...) 531 // * to (c_\beta + c_{\beta+1}*p + ...) t^\beta 532 // * 2. Computes the gcd over all (c_\beta + c_{\beta+1}*p + ...) and divides g by it 533 // **/ 534 // static void simplify(poly g, const number p) 535 // { 536 // // go along g and look for relevant terms. 537 // // monomials in x which have been checked are tracked in done. 538 // poly done = p_LmInit(g,currRing); 539 // p_SetCoeff(done,n_Init(1,currRing>cf),currRing); 540 // p_Setm(done,currRing); 541 // poly doneEnd = done; 542 // poly doneCurrent; 543 544 // poly subst; number substCoeff, substCoeffCache; 545 // unsigned long d; 546 547 // poly gCurrent = g; 548 // while(pNext(gCurrent)) 549 // { 550 // // check whether next term needs to be reduced: 551 // // first, check whether monomial in x has come up yet 552 // for (doneCurrent=done; doneCurrent; pIter(doneCurrent)) 553 // { 554 // if (p_LmDivisibleBy(doneCurrent,pNext(gCurrent),currRing)) 555 // break; 556 // } 557 // // if the monomial in x already occured, then we want to replace 558 // // as many t with p as theoretically feasable 559 // if (doneCurrent) 560 // { 561 // // suppose pNext(gCurrent)=3*t5x and doneCurrent=t3x 562 // // then we want to replace pNext(gCurrent) with 3p2*t3x 563 // // subst = ?*t3x 564 // subst = p_LmInit(doneCurrent,currRing); 565 // // substcoeff = p2 566 // n_Power(p,p_GetExp(subst,1,currRing)p_GetExp(doneCurrent,1,currRing),&substCoeff,currRing>cf); 567 // // substcoeff = 3p2 568 // n_InpMult(substCoeff,p_GetCoeff(pNext(gCurrent),currRing),currRing>cf); 569 // // subst = 3p2*t3x 570 // p_SetCoeff(subst,substCoeff,currRing); 571 // p_Setm(subst,currRing); pTest(subst); 572 573 // // g = g  pNext(gCurrent) + subst 574 // pNext(gCurrent)=p_LmDeleteAndNext(pNext(gCurrent),currRing); 575 // g=p_Add_q(g,subst,currRing); 576 // pTest(pNext(gbeginning)); 577 // } 578 // else 579 // { 580 // // if the monomial in x is brand new, 581 // // then we check whether the coefficient is divisible by p 582 // if (n_DivBy(p_GetCoeff(pNext(gCurrent),currRing>cf),p,currRing>cf)) 583 // { 584 // // reduce the term with respect to pt: 585 // // suppose pNext(gCurrent)=4p3*tx 586 // // then we want to replace it with 4*t4x 587 // // divide 4p3 repeatedly by p until it is not divisible anymore, 588 // // keeping track on the final value 4 589 // // and the number of times it has been divided 3 590 // substCoeff = n_Div(p_GetCoeff(pNext(gCurrent),currRing>cf),p,currRing>cf); 591 // d = 1; 592 // while (n_DivBy(substCoeff,p,currRing>cf)) 593 // { 594 // substCoeffCache = n_Div(p_GetCoeff(pNext(gCurrent),currRing>cf),p,currRing>cf); 595 // n_Delete(&substCoeff,currRing>cf); 596 // substCoeff = substCoeffCache; 597 // d++; 598 // assume(d>0); // check for overflow 599 // } 600 601 // // subst = ?*tx 602 // subst=p_LmInit(pNext(gCurrent),currRing); 603 // // subst = ?*t4x 604 // p_AddExp(subst,1,d,currRing); 605 // // subst = 4*t4x 606 // p_SetCoeff(subst,substCoeffCache,currRing); 607 // p_Setm(subst,currRing); pTest(subst); 608 609 // // g = g  pNext(gCurrent) + subst 610 // pNext(gCurrent)=p_LmDeleteAndNext(pNext(gCurrent),currRing); 611 // pNext(gCurrent)=p_Add_q(pNext(gCurrent),subst,currRing); 612 // pTest(pNext(gCurrent)); 613 // } 614 615 // // either way, add monomial in x with minimal t to done 616 // pNext(doneEnd)=p_LmInit(pNext(gCurrent),currRing); 617 // pIter(doneEnd); 618 // p_SetCoeff(doneEnd,n_Init(1,currRing>cf),currRing); 619 // p_Setm(doneEnd,currRing); 620 // } 621 // pIter(gCurrent); 622 // } 623 // p_Delete(&done,currRing); 624 // } 625 626 627 // #ifndef NDEBUG 628 // BOOLEAN divideByGcd(leftv res, leftv args) 629 // { 630 // leftv u = args; 631 // if ((u != NULL) && (u>Typ() == POLY_CMD)) 632 // { 633 // poly g; 634 // omUpdateInfo(); 635 // Print("usedBytes1=%ld\n",om_Info.UsedBytes); 636 // g = (poly) u>CopyD(); 637 // divideByGcd(g); 638 // p_Delete(&g,currRing); 639 // omUpdateInfo(); 640 // Print("usedBytes1=%ld\n",om_Info.UsedBytes); 641 // res>rtyp = NONE; 642 // res>data = NULL; 643 // return FALSE; 644 // } 645 // return TRUE; 646 // } 647 // #endif //NDEBUG 648 649 650 // BOOLEAN initialReduction(leftv res, leftv args) 651 // { 652 // leftv u = args; 653 // if ((u != NULL) && (u>Typ() == IDEAL_CMD)) 654 // { 655 // leftv v = u>next; 656 // if (v == NULL) 657 // { 658 // ideal I = (ideal) u>Data(); 659 660 // /*** 661 // * Suppose I=g_0,...,g_n and g_0=pt. 662 // * Step 1: sort elements g_1,...,g_{n1} such that lm(g_1)>...>lm(g_{n1}) 663 // * (the following algorithm is a bubble sort) 664 // **/ 665 // sortingLaterGenerators(I); 666 667 // /*** 668 // * Step 2: replace coefficient of terms of lowest tdegree divisible by p with t 669 // **/ 670 // number p = p_GetCoeff(I>m[0],currRing); 671 // for (int i=1; i<IDELEMS(I); i++) 672 // { 673 // if (preduce(I>m[i],p)) 674 // return TRUE; 675 // } 676 677 // /*** 678 // * Step 3: the first pass. removing terms with the same monomials in x as lt(g_i) 679 // * out of g_j for i<j 680 // **/ 681 // int i,j; 682 // poly uniti, unitj; 683 // for (i=1; i<IDELEMS(I)1; i++) 684 // { 685 // for (j=i+1; j<IDELEMS(I); j++) 686 // { 687 // unitj = highestMatchingX(I>m[j],I>m[i]); 688 // if (unitj) 689 // { 690 // unitj = powerSeriesCoeff(unitj); 691 // divideByT(unitj,p_GetExp(I>m[i],1,currRing)); 692 // uniti = powerSeriesCoeff(I>m[i]); 693 // divideByT(uniti,p_GetExp(I>m[i],1,currRing)); 694 // pTest(unitj); pTest(uniti); pTest(I>m[j]); pTest(I>m[i]); 695 // I>m[j]=p_Add_q(p_Mult_q(uniti,I>m[j],currRing), 696 // p_Neg(p_Mult_q(unitj,p_Copy(I>m[i],currRing),currRing),currRing), 697 // currRing); 698 // divideByGcd(I>m[j]); 699 // } 700 // } 701 // } 702 // for (int i=1; i<IDELEMS(I); i++) 703 // { 704 // if (preduce(I>m[i],p)) 705 // return TRUE; 706 // } 707 708 // /*** 709 // * Step 4: the second pass. removing terms divisible by lt(g_j) out of g_i for i<j 710 // **/ 711 // for (i=1; i<IDELEMS(I)1; i++) 712 // { 713 // for (j=i+1; j<IDELEMS(I); j++) 714 // { 715 // uniti = highestMatchingX(I>m[i],I>m[j]); 716 // if (uniti && p_GetExp(uniti,1,currRing)>=p_GetExp(I>m[j],1,currRing)) 717 // { 718 // uniti = powerSeriesCoeff(uniti); 719 // divideByT(uniti,p_GetExp(I>m[j],1,currRing)); 720 // unitj = powerSeriesCoeff(I>m[j]); 721 // divideByT(unitj,p_GetExp(I>m[j],1,currRing)); 722 // I>m[i] = p_Add_q(p_Mult_q(unitj,I>m[i],currRing), 723 // p_Neg(p_Mult_q(uniti,p_Copy(I>m[j],currRing),currRing),currRing), 724 // currRing); 725 // divideByGcd(I>m[j]); 726 // } 727 // } 728 // } 729 // for (int i=1; i<IDELEMS(I); i++) 730 // { 731 // if (preduce(I>m[i],p)) 732 // return TRUE; 733 // } 734 735 // res>rtyp = NONE; 736 // res>data = NULL; 737 // IDDATA((idhdl)u>data) = (char*) I; 738 // return FALSE; 739 // } 740 // } 741 // WerrorS("initialReduction: unexpected parameters"); 742 // return TRUE; 743 // } 743 744 744 745 … … 969 970 p>iiAddCproc("","initial",FALSE,initial); 970 971 #ifndef NDEBUG 971 p>iiAddCproc("","divideByGcd",FALSE,divideByGcd);972 p>iiAddCproc("","p reduce",FALSE,preduce);972 // p>iiAddCproc("","divideByGcd",FALSE,divideByGcd); 973 p>iiAddCproc("","pReduce",FALSE,pReduce); 973 974 #endif //NDEBUG 974 p>iiAddCproc("","initialReduction",FALSE,initialReduction);975 // p>iiAddCproc("","initialReduction",FALSE,initialReduction); 975 976 p>iiAddCproc("","homogeneitySpace",FALSE,homogeneitySpace); 976 977 }
Note: See TracChangeset
for help on using the changeset viewer.