Changeset 88c726 in git


Ignore:
Timestamp:
Dec 10, 2013, 6:58:00 PM (10 years ago)
Author:
Yue Ren <ren@…>
Branches:
(u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
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
Message:
new: initialReduction(ideal I)
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • Singular/dyn_modules/gfanlib/tropical.cc

    redb60f r88c726  
    237237
    238238
    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 p-t:
    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 t-exponent");
    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 x-monomial 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 p-t:
    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=p-t.
    662 //        * Step 1: sort elements g_1,...,g_{n-1} such that lm(g_1)>...>lm(g_{n-1})
    663 //        * (the following algorithm is a bubble sort)
    664 //        **/
    665 //       sortingLaterGenerators(I);
    666 
    667 //       /***
    668 //        * Step 2: replace coefficient of terms of lowest t-degree 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 // }
    744 
    745 
    746239#if 0
    747240/***
    … …  
    970463  p->iiAddCproc("","initial",FALSE,initial);
    971464#ifndef NDEBUG
    972   // p->iiAddCproc("","divideByGcd",FALSE,divideByGcd);
    973465  p->iiAddCproc("","pReduce",FALSE,pReduce);
    974466  p->iiAddCproc("","reduceInitially0",FALSE,reduceInitially0);
    975467  p->iiAddCproc("","reduceInitially1",FALSE,reduceInitially1);
    976468  p->iiAddCproc("","reduceInitially2",FALSE,reduceInitially2);
     469  p->iiAddCproc("","reduceInitially3",FALSE,reduceInitially3);
     470  p->iiAddCproc("","reduceInitially4",FALSE,reduceInitially4);
    977471#endif //NDEBUG
    978   // p->iiAddCproc("","initialReduction",FALSE,initialReduction);
     472  p->iiAddCproc("","reduceInitially",FALSE,reduceInitially);
    979473  p->iiAddCproc("","homogeneitySpace",FALSE,homogeneitySpace);
    980474}

  • dyn_modules/callgfanlib/initialReduction.cc

    redb60f r88c726  
    11#include <kernel/polys.h>
     2#include <Singular/ipid.h>
     3
    24#include <libpolys/polys/monomials/p_polys.h>
    35#include <singularWishlist.h>
    46
    5 #include <Singular/ipid.h>
     7#include <map>
     8#include <set>
    69
    710/***
    … …  
    9598
    9699/***
    97  * returns, if it exists, a pointer to the first term in g
    98  * 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 /***
    113100 * reduces h initially with respect to g,
    114101 * returns false if h was initially reduced in the first place,
    … …  
    118105bool reduceInitially(poly &h, const poly g)
    119106{
    120   poly hCache;
     107  pTest(h); pTest(g); poly hCache;
    121108  for (hCache=h; hCache; pIter(hCache))
    122109    if (p_LeadmonomDivisibleBy(g,hCache,currRing)) break;
    … …  
    130117      p_SetExp(hAlphaT,i,0,currRing);
    131118    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);
    135123    pTest(h);
    136124    return true;
    … …  
    173161
    174162
    175 static inline void sortByLeadmonom(ideal I)
    176 {
    177   poly cache; int i, n=IDELEMS(I), newn;
    178   do
    179   {
    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 
    196163/***
    197164 * reduces I initially with respect to itself and with respect to p-t.
    … …  
    201168bool reduceInitially(ideal I, const number p)
    202169{
    203   poly cache; int i,j,n=IDELEMS(I);
     170  int m=idSize(I),n=m; poly cache;
    204171  do
    205172  {
    206     j=0;
    207     for (i=1; i<n; i++)
    208     {
    209       if (pLmCmp(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)
    210177      {
    211178        cache=I->m[i-1];
    … …  
    217184    n=j;
    218185  } while(n);
    219   for (i=1; i<IDELEMS(I); i++)
     186  for (int i=1; i<m; i++)
    220187    if (pReduce(I->m[i],p)) return true;
    221188
    … …  
    223190   * the first pass. removing terms with the same monomials in x as lt(g_i) out of g_j for i<j
    224191   **/
    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++)
    227194      if (reduceInitially(I->m[j], I->m[i]) && pReduce(I->m[j],p)) return true;
    228195
    … …  
    230197   * the second pass. removing terms divisible by lt(g_j) out of g_i for i<j
    231198   **/
    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++)
    234201      if (reduceInitially(I->m[i], I->m[j]) && pReduce(I->m[i],p)) return true;
    235202  return false;
    … …  
    274241 * assumes that I was already sorted and initially reduced in the first place
    275242 **/
    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)
     243bool 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];
    283253      I->m[j] = I->m[j-1];
     254      I->m[j-1] = cache;
     255    }
    284256    else
    285     {
    286       I->m[j] = g;
    287257      break;
    288     }
    289   }
    290   if (j<1) I->m[0] = g;
     258  }
    291259
    292260  /***
    … …  
    294262   * removing terms with the same monomials in x as lt(g_j) out of g_k for j<k
    295263   **/
    296   for (i=0; i<j; i++)
     264  for (int i=0; i<j; i++)
    297265    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++)
    299267    if (reduceInitially(I->m[k], I->m[j]) && pReduce(I->m[k],p)) return true;
    300268
    … …  
    303271   * removing terms divisible by lt(g_k) out of g_j for j<k
    304272   **/
    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++)
    307275      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++)
    309277    if (reduceInitially(I->m[j],I->m[k]) && pReduce(I->m[j],p)) return true;
    310278
    … …  
    350318}
    351319#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 **/
     328bool 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
     433BOOLEAN 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 **/
     477bool 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
     563BOOLEAN 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
     587BOOLEAN 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}

  • dyn_modules/callgfanlib/initialReduction.h

    redb60f r88c726  
    99BOOLEAN reduceInitially1(leftv res, leftv args);
    1010BOOLEAN reduceInitially2(leftv res, leftv args);
     11BOOLEAN reduceInitially3(leftv res, leftv args);
     12BOOLEAN reduceInitially4(leftv res, leftv args);
    1113#endif
    1214
     15BOOLEAN reduceInitially(leftv res, leftv args);
     16
    1317#endif

  • dyn_modules/callgfanlib/singularWishlist.h

    redb60f r88c726  
    5757}
    5858
     59inline void idShallowDelete (ideal *h)
     60{
     61  int j,elems;
     62  if (*h == NULL)
     63    return;
     64  elems=j=(*h)->nrows*(*h)->ncols;
     65  if (j>0)
     66    omFreeSize((ADDRESS)((*h)->m),sizeof(poly)*elems);
     67  omFreeBin((ADDRESS)*h, sip_sideal_bin);
     68  *h=NULL;
     69}
     70
    5971#endif
Note: See TracChangeset for help on using the changeset viewer.