Changeset fa167e in git


Ignore:
Timestamp:
Dec 1, 2013, 6:20:33 PM (10 years ago)
Author:
Yue Ren <ren@…>
Branches:
(u'spielwiese', '2a584933abf2a2d3082034c7586d38bb6de1a30a')
Children:
7a4e77045299daf5ac74ab3774c5e1b784ff4e94
Parents:
9bee534198ac4b761325be1435e67598941aac19
git-author:
Yue Ren <ren@mathematik.uni-kl.de>2013-12-01 18:20:33+01:00
git-committer:
Yue Ren <ren@mathematik.uni-kl.de>2015-02-06 13:47:02+01:00
Message:
new: pReduce(poly g, const number p)
Files:
3 added
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/dyn_modules/gfanlib/tropical.cc

    r9bee53 rfa167e  
    55#include <bbcone.h>
    66
     7#include <initialReduction.h>
    78
    89poly initial(poly p)
    … …  
    236237
    237238
    238 /***
    239  * replaces coefficients in g of lowest degree in t
    240  * divisible by p with a suitable power of t
    241  **/
    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 yet
    259     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 new
    266     {
    267       // second, check whether coefficient is divisible by p
    268       if (n_DivBy(p_GetCoeff(pNext(g),currRing->cf),p,currRing->cf))
    269       {
    270         // reduce the term with respect to p-t:
    271         // divide coefficient by p, remove old term,
    272         // add t multiple of old term
    273         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 t-exponent");
    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         else
    289         {
    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 done
    299       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 NDEBUG
    313 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 //NDEBUG
    336 
    337 
    338 /***
    339  * Returns the highest term in g with the matching x-monomial to m
    340  * or, if it does not exist, the NULL pointer
    341  **/
    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 term
    348   // with the same monomial in x as xalpha
    349   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 term
    359   // or is NULL
    360   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_\alpha
    368  **/
    369 static poly powerSeriesCoeff(poly g)
    370 {
    371   // the first term obviously is part of our output
    372   // 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 coeffCache
    381   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 g
    405  * is divisible by t^b, i.e. no divisibility checks
    406  * needed
    407  **/
    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_\alpha
    422   // as alpha runs over all exponent vectors
    423   // of monomials in x occuring in g
    424 
    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 positive
    429   // so that the x won't make the monomial smaller
    430   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 up
    447     // since everything is homogeneous in x and the ordering is local in t,
    448     // this comes down to a divisibility test in two stages
    449     current_not_sev = ~p_GetShortExpVector(current,currRing);
    450     for(i=0; i<=count; i++)
    451     {
    452       // first stage using short exponent vectors
    453       // second stage a proper test
    454       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 entry
    458         // without the x part
    459         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 term
    470       // including the monomial in x
    471       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 terms
    484     // so that gAlpha only contais terms in t
    485     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_\alpha
    495     // and set the input to null as they are deleted in the process
    496     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 gcd
    505     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   else
    514   {
    515     // if g contains only one monomial in x,
    516     // then we can get rid of all the higher t
    517     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^\beta
    531  * 2. Computes the gcd over all (c_\beta + c_{\beta+1}*p + ...) and divides g by it
    532  **/
    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 yet
    551     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 replace
    557     // as many t with p as theoretically feasable
    558     if (doneCurrent)
    559     {
    560       // suppose pNext(gCurrent)=3*t5x and doneCurrent=t3x
    561       // then we want to replace pNext(gCurrent) with 3p2*t3x
    562       // subst = ?*t3x
    563       subst = p_LmInit(doneCurrent,currRing);
    564       // substcoeff = p2
    565       n_Power(p,p_GetExp(subst,1,currRing)-p_GetExp(doneCurrent,1,currRing),&substCoeff,currRing->cf);
    566       // substcoeff = 3p2
    567       n_InpMult(substCoeff,p_GetCoeff(pNext(gCurrent),currRing),currRing->cf);
    568       // subst = 3p2*t3x
    569       p_SetCoeff(subst,substCoeff,currRing);
    570       p_Setm(subst,currRing); pTest(subst);
    571 
    572       // g = g - pNext(gCurrent) + subst
    573       pNext(gCurrent)=p_LmDeleteAndNext(pNext(gCurrent),currRing);
    574       g=p_Add_q(g,subst,currRing);
    575       pTest(pNext(gbeginning));
    576     }
    577     else
    578     {
    579       // if the monomial in x is brand new,
    580       // then we check whether the coefficient is divisible by p
    581       if (n_DivBy(p_GetCoeff(pNext(gCurrent),currRing->cf),p,currRing->cf))
    582       {
    583         // reduce the term with respect to p-t:
    584         // suppose pNext(gCurrent)=4p3*tx
    585         // then we want to replace it with 4*t4x
    586         // divide 4p3 repeatedly by p until it is not divisible anymore,
    587         // keeping track on the final value 4
    588         // and the number of times it has been divided 3
    589         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 overflow
    598         }
    599 
    600         // subst = ?*tx
    601         subst=p_LmInit(pNext(gCurrent),currRing);
    602         // subst = ?*t4x
    603         p_AddExp(subst,1,d,currRing);
    604         // subst = 4*t4x
    605         p_SetCoeff(subst,substCoeffCache,currRing);
    606         p_Setm(subst,currRing); pTest(subst);
    607 
    608         // g = g - pNext(gCurrent) + subst
    609         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 done
    615       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 NDEBUG
    627 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 //NDEBUG
    647 
    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=p-t.
    661        * Step 1: sort elements g_1,...,g_{n-1} such that lm(g_1)>...>lm(g_{n-1})
    662        * (the following algorithm is a bubble sort)
    663        **/
    664       sortingLaterGenerators(I);
    665 
    666       /***
    667        * Step 2: replace coefficient of terms of lowest t-degree divisible by p with t
    668        **/
    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<j
    679        **/
    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<j
    709        **/
    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 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// }
    743744
    744745
    … …  
    969970  p->iiAddCproc("","initial",FALSE,initial);
    970971#ifndef NDEBUG
    971   p->iiAddCproc("","divideByGcd",FALSE,divideByGcd);
    972   p->iiAddCproc("","preduce",FALSE,preduce);
     972  // p->iiAddCproc("","divideByGcd",FALSE,divideByGcd);
     973  p->iiAddCproc("","pReduce",FALSE,pReduce);
    973974#endif //NDEBUG
    974   p->iiAddCproc("","initialReduction",FALSE,initialReduction);
     975  // p->iiAddCproc("","initialReduction",FALSE,initialReduction);
    975976  p->iiAddCproc("","homogeneitySpace",FALSE,homogeneitySpace);
    976977}
Note: See TracChangeset for help on using the changeset viewer.