# Changeset fa167e in git for Singular/dyn_modules/gfanlib/tropical.cc

Ignore:
Timestamp:
Dec 1, 2013, 6:20:33 PM (9 years ago)
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
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)
File:
1 edited

Unmodified
Added
Removed
• ## Singular/dyn_modules/gfanlib/tropical.cc

 r9bee53 #include #include poly initial(poly p) /*** * replaces coefficients in g of lowest degree in t * divisible by p with a suitable power of t **/ static bool preduce(poly g, const number p) { // go along g and look for relevant terms. // monomials in x which have been checked are tracked in done. // because we assume the leading coefficient of g is 1, // the leading term does not need to be considered. poly done = p_LmInit(g,currRing); p_SetExp(done,1,0,currRing); p_SetCoeff(done,n_Init(1,currRing->cf),currRing); p_Setm(done,currRing); poly doneEnd = done; poly doneCache; number dnumber; long d; poly subst; number ppower; while(pNext(g)) { // check whether next term needs to be reduced: // first, check whether monomial in x has come up yet for (doneCache=done; doneCache; pIter(doneCache)) { if (p_LmDivisibleBy(doneCache,pNext(g),currRing)) break; } if (!doneCache) // if for loop did not terminate prematurely, // then the monomial in x is new { // second, check whether coefficient is divisible by p if (n_DivBy(p_GetCoeff(pNext(g),currRing->cf),p,currRing->cf)) { // reduce the term with respect to p-t: // divide coefficient by p, remove old term, // add t multiple of old term dnumber = n_Div(p_GetCoeff(pNext(g),currRing->cf),p,currRing->cf); d = n_Int(dnumber,currRing->cf); n_Delete(&dnumber,currRing->cf); if (!d) { p_Delete(&done,currRing); WerrorS("initialReduction: overflow in a t-exponent"); return true; } subst=p_LmInit(pNext(g),currRing); if (d>0) { p_AddExp(subst,1,d,currRing); p_SetCoeff(subst,n_Init(1,currRing->cf),currRing); } else { p_AddExp(subst,1,-d,currRing); p_SetCoeff(subst,n_Init(-1,currRing->cf),currRing); } p_Setm(subst,currRing); pTest(subst); pNext(g)=p_LmDeleteAndNext(pNext(g),currRing); pNext(g)=p_Add_q(pNext(g),subst,currRing); pTest(pNext(g)); } // either way, add monomial in x to done pNext(doneEnd)=p_LmInit(pNext(g),currRing); pIter(doneEnd); p_SetExp(doneEnd,1,0,currRing); p_SetCoeff(doneEnd,n_Init(1,currRing->cf),currRing); p_Setm(doneEnd,currRing); } pIter(g); } p_Delete(&done,currRing); return false; } #ifndef NDEBUG BOOLEAN preduce(leftv res, leftv args) { leftv u = args; if ((u != NULL) && (u->Typ() == POLY_CMD)) { poly g; bool b; number p = n_Init(3,currRing->cf); omUpdateInfo(); Print("usedBytes=%ld\n",om_Info.UsedBytes); g = (poly) u->CopyD(); b = preduce(g,p); p_Delete(&g,currRing); if (b) return TRUE; omUpdateInfo(); Print("usedBytes=%ld\n",om_Info.UsedBytes); n_Delete(&p,currRing->cf); res->rtyp = NONE; res->data = NULL; return FALSE; } return TRUE; } #endif //NDEBUG /*** * Returns the highest term in g with the matching x-monomial to m * or, if it does not exist, the NULL pointer **/ static poly highestMatchingX(poly g, const poly m) { pTest(g); pTest(m); poly xalpha=p_LmInit(m,currRing); // go along g and find the first term // with the same monomial in x as xalpha while (g) { p_SetExp(xalpha,1,p_GetExp(g,1,currRing),currRing); p_Setm(xalpha,currRing); if (p_LmEqual(g,xalpha,currRing)) break; pIter(g); } // gCache now either points at the wanted term // or is NULL p_Delete(&xalpha,currRing); pTest(g); return g; } /*** * Given g with lm(g)=t^\beta x^\alpha, returns g_\alpha **/ static poly powerSeriesCoeff(poly g) { // the first term obviously is part of our output // so we copy it... poly xalpha=p_LmInit(g,currRing); poly coeff=p_One(currRing); p_SetCoeff(coeff,n_Copy(p_GetCoeff(g,currRing),currRing->cf),currRing); p_SetExp(coeff,1,p_GetExp(g,1,currRing),currRing); p_Setm(coeff,currRing); // ..and proceed with the remaining terms, // appending the relevant terms to coeff via coeffCache poly coeffCache=coeff; pIter(g); while (g) { p_SetExp(xalpha,1,p_GetExp(g,1,currRing),currRing); p_Setm(xalpha,currRing); if (p_LmEqual(g,xalpha,currRing)) { pNext(coeffCache)=p_Init(currRing); pIter(coeffCache); p_SetCoeff(coeffCache,n_Copy(p_GetCoeff(g,currRing),currRing->cf),currRing); p_SetExp(coeffCache,1,p_GetExp(g,1,currRing),currRing); p_Setm(coeffCache,currRing); } pIter(g); } p_Delete(&xalpha,currRing); return coeff; } /*** * divides g by t^b knowing that each term of g * is divisible by t^b, i.e. no divisibility checks * needed **/ static void divideByT(poly g, const long b) { while (g) { p_SetExp(g,1,p_GetExp(g,1,currRing)-b,currRing); p_Setm(g,currRing); pIter(g); } } static void divideByGcd(poly &g) { // first determine all g_\alpha // as alpha runs over all exponent vectors // of monomials in x occuring in g // gAlpha will store all g_\alpha, // the first term will, for comparison purposes for now, // also keep their monomial in x. // we assume that the weight on the x are positive // so that the x won't make the monomial smaller ideal gAlphaFront = idInit(pLength(g)); gAlphaFront->m[0] = p_Head(g,currRing); p_SetExp(gAlphaFront->m[0],1,0,currRing); p_Setm(gAlphaFront->m[0],currRing); ideal gAlphaEnd = idInit(pLength(g)); gAlphaEnd->m[0] = gAlphaFront->m[0]; unsigned long gAlpha_sev[pLength(g)]; gAlpha_sev[0] = p_GetShortExpVector(g,currRing); long beta[pLength(g)]; beta[0] = p_GetExp(g,1,currRing); int count=0; poly current = pNext(g); unsigned long current_not_sev; int i; while (current) { // see whether the monomial in x of current already came up // since everything is homogeneous in x and the ordering is local in t, // this comes down to a divisibility test in two stages current_not_sev = ~p_GetShortExpVector(current,currRing); for(i=0; i<=count; i++) { // first stage using short exponent vectors // second stage a proper test if (p_LmShortDivisibleBy(gAlphaFront->m[i],gAlpha_sev[i],current,current_not_sev, currRing) && p_LmDivisibleBy(gAlphaFront->m[i],current,currRing)) { // if it already occured, add the term to the respective entry // without the x part pNext(gAlphaEnd->m[i])=p_Init(currRing); pIter(gAlphaEnd->m[i]); p_SetExp(gAlphaEnd->m[i],1,p_GetExp(current,1,currRing)-beta[i],currRing); p_SetCoeff(gAlphaEnd->m[i],n_Copy(p_GetCoeff(current,currRing),currRing->cf),currRing); p_Setm(gAlphaEnd->m[i],currRing); break; } } if (i>count) { // if it is new, create a new entry for the term // including the monomial in x count++; gAlphaFront->m[count] = p_Head(current,currRing); beta[count] = p_GetExp(current,1,currRing); gAlphaEnd->m[count] = gAlphaFront->m[count]; gAlpha_sev[count] = p_GetShortExpVector(current,currRing); } pIter(current); } if (count) { // now remove all the x in the leading terms // so that gAlpha only contais terms in t int j; long d; for (i=0; i<=count; i++) { for (j=2; j<=currRing->N; j++) p_SetExp(gAlphaFront->m[i],j,0,currRing); p_Setm(gAlphaFront->m[i],currRing); gAlphaEnd->m[i]=NULL; } // now compute the gcd over all g_\alpha // and set the input to null as they are deleted in the process poly gcd = singclap_gcd(gAlphaFront->m[0],gAlphaFront->m[1],currRing); gAlphaFront->m[0] = NULL; gAlphaFront->m[1] = NULL; for(i=2; i<=count; i++) { gcd = singclap_gcd(gcd,gAlphaFront->m[i],currRing); gAlphaFront->m[i] = NULL; } // divide g by the gcd poly h = singclap_pdivide(g,gcd,currRing); p_Delete(&gcd,currRing); p_Delete(&g,currRing); g = h; id_Delete(&gAlphaFront,currRing); id_Delete(&gAlphaEnd,currRing); } else { // if g contains only one monomial in x, // then we can get rid of all the higher t p_Delete(&g,currRing); g = gAlphaFront->m[0]; pIter(gAlphaFront->m[0]); pNext(g)=NULL; gAlphaEnd->m[0] = NULL; id_Delete(&gAlphaFront,currRing); id_Delete(&gAlphaEnd,currRing); } } /*** * 1. For each \alpha\in\NN^n, changes (c_\beta t^\beta + c_{\beta+1} t^{\beta+1} + ...) *    to (c_\beta + c_{\beta+1}*p + ...) t^\beta * 2. Computes the gcd over all (c_\beta + c_{\beta+1}*p + ...) and divides g by it **/ static void simplify(poly g, const number p) { // go along g and look for relevant terms. // monomials in x which have been checked are tracked in done. poly done = p_LmInit(g,currRing); p_SetCoeff(done,n_Init(1,currRing->cf),currRing); p_Setm(done,currRing); poly doneEnd = done; poly doneCurrent; poly subst; number substCoeff, substCoeffCache; unsigned long d; poly gCurrent = g; while(pNext(gCurrent)) { // check whether next term needs to be reduced: // first, check whether monomial in x has come up yet for (doneCurrent=done; doneCurrent; pIter(doneCurrent)) { if (p_LmDivisibleBy(doneCurrent,pNext(gCurrent),currRing)) break; } // if the monomial in x already occured, then we want to replace // as many t with p as theoretically feasable if (doneCurrent) { // suppose pNext(gCurrent)=3*t5x and doneCurrent=t3x // then we want to replace pNext(gCurrent) with 3p2*t3x // subst = ?*t3x subst = p_LmInit(doneCurrent,currRing); // substcoeff = p2 n_Power(p,p_GetExp(subst,1,currRing)-p_GetExp(doneCurrent,1,currRing),&substCoeff,currRing->cf); // substcoeff = 3p2 n_InpMult(substCoeff,p_GetCoeff(pNext(gCurrent),currRing),currRing->cf); // subst = 3p2*t3x p_SetCoeff(subst,substCoeff,currRing); p_Setm(subst,currRing); pTest(subst); // g = g - pNext(gCurrent) + subst pNext(gCurrent)=p_LmDeleteAndNext(pNext(gCurrent),currRing); g=p_Add_q(g,subst,currRing); pTest(pNext(gbeginning)); } else { // if the monomial in x is brand new, // then we check whether the coefficient is divisible by p if (n_DivBy(p_GetCoeff(pNext(gCurrent),currRing->cf),p,currRing->cf)) { // reduce the term with respect to p-t: // suppose pNext(gCurrent)=4p3*tx // then we want to replace it with 4*t4x // divide 4p3 repeatedly by p until it is not divisible anymore, // keeping track on the final value 4 // and the number of times it has been divided 3 substCoeff = n_Div(p_GetCoeff(pNext(gCurrent),currRing->cf),p,currRing->cf); d = 1; while (n_DivBy(substCoeff,p,currRing->cf)) { substCoeffCache = n_Div(p_GetCoeff(pNext(gCurrent),currRing->cf),p,currRing->cf); n_Delete(&substCoeff,currRing->cf); substCoeff = substCoeffCache; d++; assume(d>0); // check for overflow } // subst = ?*tx subst=p_LmInit(pNext(gCurrent),currRing); // subst = ?*t4x p_AddExp(subst,1,d,currRing); // subst = 4*t4x p_SetCoeff(subst,substCoeffCache,currRing); p_Setm(subst,currRing); pTest(subst); // g = g - pNext(gCurrent) + subst pNext(gCurrent)=p_LmDeleteAndNext(pNext(gCurrent),currRing); pNext(gCurrent)=p_Add_q(pNext(gCurrent),subst,currRing); pTest(pNext(gCurrent)); } // either way, add monomial in x with minimal t to done pNext(doneEnd)=p_LmInit(pNext(gCurrent),currRing); pIter(doneEnd); p_SetCoeff(doneEnd,n_Init(1,currRing->cf),currRing); p_Setm(doneEnd,currRing); } pIter(gCurrent); } p_Delete(&done,currRing); } #ifndef NDEBUG BOOLEAN divideByGcd(leftv res, leftv args) { leftv u = args; if ((u != NULL) && (u->Typ() == POLY_CMD)) { poly g; omUpdateInfo(); Print("usedBytes1=%ld\n",om_Info.UsedBytes); g = (poly) u->CopyD(); divideByGcd(g); p_Delete(&g,currRing); omUpdateInfo(); Print("usedBytes1=%ld\n",om_Info.UsedBytes); res->rtyp = NONE; res->data = NULL; return FALSE; } return TRUE; } #endif //NDEBUG BOOLEAN initialReduction(leftv res, leftv args) { leftv u = args; if ((u != NULL) && (u->Typ() == IDEAL_CMD)) { leftv v = u->next; if (v == NULL) { ideal I = (ideal) u->Data(); /*** * Suppose I=g_0,...,g_n and g_0=p-t. * Step 1: sort elements g_1,...,g_{n-1} such that lm(g_1)>...>lm(g_{n-1}) * (the following algorithm is a bubble sort) **/ sortingLaterGenerators(I); /*** * Step 2: replace coefficient of terms of lowest t-degree divisible by p with t **/ number p = p_GetCoeff(I->m[0],currRing); for (int i=1; im[i],p)) return TRUE; } /*** * Step 3: the first pass. removing terms with the same monomials in x as lt(g_i) * out of g_j for im[j],I->m[i]); if (unitj) { unitj = powerSeriesCoeff(unitj); divideByT(unitj,p_GetExp(I->m[i],1,currRing)); uniti = powerSeriesCoeff(I->m[i]); divideByT(uniti,p_GetExp(I->m[i],1,currRing)); pTest(unitj); pTest(uniti); pTest(I->m[j]); pTest(I->m[i]); I->m[j]=p_Add_q(p_Mult_q(uniti,I->m[j],currRing), p_Neg(p_Mult_q(unitj,p_Copy(I->m[i],currRing),currRing),currRing), currRing); divideByGcd(I->m[j]); } } } for (int i=1; im[i],p)) return TRUE; } /*** * Step 4: the second pass. removing terms divisible by lt(g_j) out of g_i for im[i],I->m[j]); if (uniti && p_GetExp(uniti,1,currRing)>=p_GetExp(I->m[j],1,currRing)) { uniti = powerSeriesCoeff(uniti); divideByT(uniti,p_GetExp(I->m[j],1,currRing)); unitj = powerSeriesCoeff(I->m[j]); divideByT(unitj,p_GetExp(I->m[j],1,currRing)); I->m[i] = p_Add_q(p_Mult_q(unitj,I->m[i],currRing), p_Neg(p_Mult_q(uniti,p_Copy(I->m[j],currRing),currRing),currRing), currRing); divideByGcd(I->m[j]); } } } for (int i=1; im[i],p)) return TRUE; } res->rtyp = NONE; res->data = NULL; IDDATA((idhdl)u->data) = (char*) I; return FALSE; } } WerrorS("initialReduction: unexpected parameters"); return TRUE; } // /*** // * replaces coefficients in g of lowest degree in t // * divisible by p with a suitable power of t // **/ // static bool preduce(poly g, const number p) // { //   // go along g and look for relevant terms. //   // monomials in x which have been checked are tracked in done. //   // because we assume the leading coefficient of g is 1, //   // the leading term does not need to be considered. //   poly done = p_LmInit(g,currRing); //   p_SetExp(done,1,0,currRing); p_SetCoeff(done,n_Init(1,currRing->cf),currRing); //   p_Setm(done,currRing); //   poly doneEnd = done; //   poly doneCache; //   number dnumber; long d; //   poly subst; number ppower; //   while(pNext(g)) //   { //     // check whether next term needs to be reduced: //     // first, check whether monomial in x has come up yet //     for (doneCache=done; doneCache; pIter(doneCache)) //     { //       if (p_LmDivisibleBy(doneCache,pNext(g),currRing)) //         break; //     } //     if (!doneCache) // if for loop did not terminate prematurely, //                     // then the monomial in x is new //     { //       // second, check whether coefficient is divisible by p //       if (n_DivBy(p_GetCoeff(pNext(g),currRing->cf),p,currRing->cf)) //       { //         // reduce the term with respect to p-t: //         // divide coefficient by p, remove old term, //         // add t multiple of old term //         dnumber = n_Div(p_GetCoeff(pNext(g),currRing->cf),p,currRing->cf); //         d = n_Int(dnumber,currRing->cf); //         n_Delete(&dnumber,currRing->cf); //         if (!d) //         { //           p_Delete(&done,currRing); //           WerrorS("initialReduction: overflow in a t-exponent"); //           return true; //         } //         subst=p_LmInit(pNext(g),currRing); //         if (d>0) //         { //           p_AddExp(subst,1,d,currRing); //           p_SetCoeff(subst,n_Init(1,currRing->cf),currRing); //         } //         else //         { //           p_AddExp(subst,1,-d,currRing); //           p_SetCoeff(subst,n_Init(-1,currRing->cf),currRing); //         } //         p_Setm(subst,currRing); pTest(subst); //         pNext(g)=p_LmDeleteAndNext(pNext(g),currRing); //         pNext(g)=p_Add_q(pNext(g),subst,currRing); //         pTest(pNext(g)); //       } //       // either way, add monomial in x to done //       pNext(doneEnd)=p_LmInit(pNext(g),currRing); //       pIter(doneEnd); //       p_SetExp(doneEnd,1,0,currRing); //       p_SetCoeff(doneEnd,n_Init(1,currRing->cf),currRing); //       p_Setm(doneEnd,currRing); //     } //     pIter(g); //   } //   p_Delete(&done,currRing); //   return false; // } // #ifndef NDEBUG // BOOLEAN preduce(leftv res, leftv args) // { //   leftv u = args; //   if ((u != NULL) && (u->Typ() == POLY_CMD)) //   { //     poly g; bool b; //     number p = n_Init(3,currRing->cf); //     omUpdateInfo(); //     Print("usedBytes=%ld\n",om_Info.UsedBytes); //     g = (poly) u->CopyD(); //     b = preduce(g,p); //     p_Delete(&g,currRing); //     if (b) return TRUE; //     omUpdateInfo(); //     Print("usedBytes=%ld\n",om_Info.UsedBytes); //     n_Delete(&p,currRing->cf); //     res->rtyp = NONE; //     res->data = NULL; //     return FALSE; //   } //   return TRUE; // } // #endif //NDEBUG // /*** // * Returns the highest term in g with the matching x-monomial to m // * or, if it does not exist, the NULL pointer // **/ // static poly highestMatchingX(poly g, const poly m) // { //   pTest(g); pTest(m); //   poly xalpha=p_LmInit(m,currRing); //   // go along g and find the first term //   // with the same monomial in x as xalpha //   while (g) //   { //     p_SetExp(xalpha,1,p_GetExp(g,1,currRing),currRing); //     p_Setm(xalpha,currRing); //     if (p_LmEqual(g,xalpha,currRing)) //       break; //     pIter(g); //   } //   // gCache now either points at the wanted term //   // or is NULL //   p_Delete(&xalpha,currRing); //   pTest(g); //   return g; // } // /*** // * Given g with lm(g)=t^\beta x^\alpha, returns g_\alpha // **/ // static poly powerSeriesCoeff(poly g) // { //   // the first term obviously is part of our output //   // so we copy it... //   poly xalpha=p_LmInit(g,currRing); //   poly coeff=p_One(currRing); //   p_SetCoeff(coeff,n_Copy(p_GetCoeff(g,currRing),currRing->cf),currRing); //   p_SetExp(coeff,1,p_GetExp(g,1,currRing),currRing); //   p_Setm(coeff,currRing); //   // ..and proceed with the remaining terms, //   // appending the relevant terms to coeff via coeffCache //   poly coeffCache=coeff; //   pIter(g); //   while (g) //   { //     p_SetExp(xalpha,1,p_GetExp(g,1,currRing),currRing); //     p_Setm(xalpha,currRing); //     if (p_LmEqual(g,xalpha,currRing)) //     { //       pNext(coeffCache)=p_Init(currRing); //       pIter(coeffCache); //       p_SetCoeff(coeffCache,n_Copy(p_GetCoeff(g,currRing),currRing->cf),currRing); //       p_SetExp(coeffCache,1,p_GetExp(g,1,currRing),currRing); //       p_Setm(coeffCache,currRing); //     } //     pIter(g); //   } //   p_Delete(&xalpha,currRing); //   return coeff; // } // /*** // * divides g by t^b knowing that each term of g // * is divisible by t^b, i.e. no divisibility checks // * needed // **/ // static void divideByT(poly g, const long b) // { //   while (g) //   { //     p_SetExp(g,1,p_GetExp(g,1,currRing)-b,currRing); //     p_Setm(g,currRing); //     pIter(g); //   } // } // static void divideByGcd(poly &g) // { //   // first determine all g_\alpha //   // as alpha runs over all exponent vectors //   // of monomials in x occuring in g //   // gAlpha will store all g_\alpha, //   // the first term will, for comparison purposes for now, //   // also keep their monomial in x. //   // we assume that the weight on the x are positive //   // so that the x won't make the monomial smaller //   ideal gAlphaFront = idInit(pLength(g)); //   gAlphaFront->m[0] = p_Head(g,currRing); //   p_SetExp(gAlphaFront->m[0],1,0,currRing); //   p_Setm(gAlphaFront->m[0],currRing); //   ideal gAlphaEnd = idInit(pLength(g)); //   gAlphaEnd->m[0] = gAlphaFront->m[0]; //   unsigned long gAlpha_sev[pLength(g)]; //   gAlpha_sev[0] = p_GetShortExpVector(g,currRing); //   long beta[pLength(g)]; //   beta[0] = p_GetExp(g,1,currRing); //   int count=0; //   poly current = pNext(g); unsigned long current_not_sev; //   int i; //   while (current) //   { //     // see whether the monomial in x of current already came up //     // since everything is homogeneous in x and the ordering is local in t, //     // this comes down to a divisibility test in two stages //     current_not_sev = ~p_GetShortExpVector(current,currRing); //     for(i=0; i<=count; i++) //     { //       // first stage using short exponent vectors //       // second stage a proper test //       if (p_LmShortDivisibleBy(gAlphaFront->m[i],gAlpha_sev[i],current,current_not_sev, currRing) //             && p_LmDivisibleBy(gAlphaFront->m[i],current,currRing)) //       { //         // if it already occured, add the term to the respective entry //         // without the x part //         pNext(gAlphaEnd->m[i])=p_Init(currRing); //         pIter(gAlphaEnd->m[i]); //         p_SetExp(gAlphaEnd->m[i],1,p_GetExp(current,1,currRing)-beta[i],currRing); //         p_SetCoeff(gAlphaEnd->m[i],n_Copy(p_GetCoeff(current,currRing),currRing->cf),currRing); //         p_Setm(gAlphaEnd->m[i],currRing); //         break; //       } //     } //     if (i>count) //     { //       // if it is new, create a new entry for the term //       // including the monomial in x //       count++; //       gAlphaFront->m[count] = p_Head(current,currRing); //       beta[count] = p_GetExp(current,1,currRing); //       gAlphaEnd->m[count] = gAlphaFront->m[count]; //       gAlpha_sev[count] = p_GetShortExpVector(current,currRing); //     } //     pIter(current); //   } //   if (count) //   { //     // now remove all the x in the leading terms //     // so that gAlpha only contais terms in t //     int j; long d; //     for (i=0; i<=count; i++) //     { //       for (j=2; j<=currRing->N; j++) //         p_SetExp(gAlphaFront->m[i],j,0,currRing); //       p_Setm(gAlphaFront->m[i],currRing); //       gAlphaEnd->m[i]=NULL; //     } //     // now compute the gcd over all g_\alpha //     // and set the input to null as they are deleted in the process //     poly gcd = singclap_gcd(gAlphaFront->m[0],gAlphaFront->m[1],currRing); //     gAlphaFront->m[0] = NULL; //     gAlphaFront->m[1] = NULL; //     for(i=2; i<=count; i++) //     { //       gcd = singclap_gcd(gcd,gAlphaFront->m[i],currRing); //       gAlphaFront->m[i] = NULL; //     } //     // divide g by the gcd //     poly h = singclap_pdivide(g,gcd,currRing); //     p_Delete(&gcd,currRing); //     p_Delete(&g,currRing); //     g = h; //     id_Delete(&gAlphaFront,currRing); //     id_Delete(&gAlphaEnd,currRing); //   } //   else //   { //     // if g contains only one monomial in x, //     // then we can get rid of all the higher t //     p_Delete(&g,currRing); //     g = gAlphaFront->m[0]; //     pIter(gAlphaFront->m[0]); //     pNext(g)=NULL; //     gAlphaEnd->m[0] = NULL; //     id_Delete(&gAlphaFront,currRing); //     id_Delete(&gAlphaEnd,currRing); //   } // } // /*** // * 1. For each \alpha\in\NN^n, changes (c_\beta t^\beta + c_{\beta+1} t^{\beta+1} + ...) // *    to (c_\beta + c_{\beta+1}*p + ...) t^\beta // * 2. Computes the gcd over all (c_\beta + c_{\beta+1}*p + ...) and divides g by it // **/ // static void simplify(poly g, const number p) // { //   // go along g and look for relevant terms. //   // monomials in x which have been checked are tracked in done. //   poly done = p_LmInit(g,currRing); //   p_SetCoeff(done,n_Init(1,currRing->cf),currRing); //   p_Setm(done,currRing); //   poly doneEnd = done; //   poly doneCurrent; //   poly subst; number substCoeff, substCoeffCache; //   unsigned long d; //   poly gCurrent = g; //   while(pNext(gCurrent)) //   { //     // check whether next term needs to be reduced: //     // first, check whether monomial in x has come up yet //     for (doneCurrent=done; doneCurrent; pIter(doneCurrent)) //     { //       if (p_LmDivisibleBy(doneCurrent,pNext(gCurrent),currRing)) //         break; //     } //     // if the monomial in x already occured, then we want to replace //     // as many t with p as theoretically feasable //     if (doneCurrent) //     { //       // suppose pNext(gCurrent)=3*t5x and doneCurrent=t3x //       // then we want to replace pNext(gCurrent) with 3p2*t3x //       // subst = ?*t3x //       subst = p_LmInit(doneCurrent,currRing); //       // substcoeff = p2 //       n_Power(p,p_GetExp(subst,1,currRing)-p_GetExp(doneCurrent,1,currRing),&substCoeff,currRing->cf); //       // substcoeff = 3p2 //       n_InpMult(substCoeff,p_GetCoeff(pNext(gCurrent),currRing),currRing->cf); //       // subst = 3p2*t3x //       p_SetCoeff(subst,substCoeff,currRing); //       p_Setm(subst,currRing); pTest(subst); //       // g = g - pNext(gCurrent) + subst //       pNext(gCurrent)=p_LmDeleteAndNext(pNext(gCurrent),currRing); //       g=p_Add_q(g,subst,currRing); //       pTest(pNext(gbeginning)); //     } //     else //     { //       // if the monomial in x is brand new, //       // then we check whether the coefficient is divisible by p //       if (n_DivBy(p_GetCoeff(pNext(gCurrent),currRing->cf),p,currRing->cf)) //       { //         // reduce the term with respect to p-t: //         // suppose pNext(gCurrent)=4p3*tx //         // then we want to replace it with 4*t4x //         // divide 4p3 repeatedly by p until it is not divisible anymore, //         // keeping track on the final value 4 //         // and the number of times it has been divided 3 //         substCoeff = n_Div(p_GetCoeff(pNext(gCurrent),currRing->cf),p,currRing->cf); //         d = 1; //         while (n_DivBy(substCoeff,p,currRing->cf)) //         { //           substCoeffCache = n_Div(p_GetCoeff(pNext(gCurrent),currRing->cf),p,currRing->cf); //           n_Delete(&substCoeff,currRing->cf); //           substCoeff = substCoeffCache; //           d++; //           assume(d>0); // check for overflow //         } //         // subst = ?*tx //         subst=p_LmInit(pNext(gCurrent),currRing); //         // subst = ?*t4x //         p_AddExp(subst,1,d,currRing); //         // subst = 4*t4x //         p_SetCoeff(subst,substCoeffCache,currRing); //         p_Setm(subst,currRing); pTest(subst); //         // g = g - pNext(gCurrent) + subst //         pNext(gCurrent)=p_LmDeleteAndNext(pNext(gCurrent),currRing); //         pNext(gCurrent)=p_Add_q(pNext(gCurrent),subst,currRing); //         pTest(pNext(gCurrent)); //       } //       // either way, add monomial in x with minimal t to done //       pNext(doneEnd)=p_LmInit(pNext(gCurrent),currRing); //       pIter(doneEnd); //       p_SetCoeff(doneEnd,n_Init(1,currRing->cf),currRing); //       p_Setm(doneEnd,currRing); //     } //     pIter(gCurrent); //   } //   p_Delete(&done,currRing); // } // #ifndef NDEBUG // BOOLEAN divideByGcd(leftv res, leftv args) // { //   leftv u = args; //   if ((u != NULL) && (u->Typ() == POLY_CMD)) //   { //     poly g; //     omUpdateInfo(); //     Print("usedBytes1=%ld\n",om_Info.UsedBytes); //     g = (poly) u->CopyD(); //     divideByGcd(g); //     p_Delete(&g,currRing); //     omUpdateInfo(); //     Print("usedBytes1=%ld\n",om_Info.UsedBytes); //     res->rtyp = NONE; //     res->data = NULL; //     return FALSE; //   } //   return TRUE; // } // #endif //NDEBUG // BOOLEAN initialReduction(leftv res, leftv args) // { //   leftv u = args; //   if ((u != NULL) && (u->Typ() == IDEAL_CMD)) //   { //     leftv v = u->next; //     if (v == NULL) //     { //       ideal I = (ideal) u->Data(); //       /*** //        * Suppose I=g_0,...,g_n and g_0=p-t. //        * Step 1: sort elements g_1,...,g_{n-1} such that lm(g_1)>...>lm(g_{n-1}) //        * (the following algorithm is a bubble sort) //        **/ //       sortingLaterGenerators(I); //       /*** //        * Step 2: replace coefficient of terms of lowest t-degree divisible by p with t //        **/ //       number p = p_GetCoeff(I->m[0],currRing); //       for (int i=1; im[i],p)) //           return TRUE; //       } //       /*** //        * Step 3: the first pass. removing terms with the same monomials in x as lt(g_i) //        * out of g_j for im[j],I->m[i]); //           if (unitj) //           { //             unitj = powerSeriesCoeff(unitj); //             divideByT(unitj,p_GetExp(I->m[i],1,currRing)); //             uniti = powerSeriesCoeff(I->m[i]); //             divideByT(uniti,p_GetExp(I->m[i],1,currRing)); //             pTest(unitj); pTest(uniti); pTest(I->m[j]); pTest(I->m[i]); //             I->m[j]=p_Add_q(p_Mult_q(uniti,I->m[j],currRing), //                             p_Neg(p_Mult_q(unitj,p_Copy(I->m[i],currRing),currRing),currRing), //                             currRing); //             divideByGcd(I->m[j]); //           } //         } //       } //       for (int i=1; im[i],p)) //           return TRUE; //       } //       /*** //        * Step 4: the second pass. removing terms divisible by lt(g_j) out of g_i for im[i],I->m[j]); //           if (uniti && p_GetExp(uniti,1,currRing)>=p_GetExp(I->m[j],1,currRing)) //           { //             uniti = powerSeriesCoeff(uniti); //             divideByT(uniti,p_GetExp(I->m[j],1,currRing)); //             unitj = powerSeriesCoeff(I->m[j]); //             divideByT(unitj,p_GetExp(I->m[j],1,currRing)); //             I->m[i] = p_Add_q(p_Mult_q(unitj,I->m[i],currRing), //                               p_Neg(p_Mult_q(uniti,p_Copy(I->m[j],currRing),currRing),currRing), //                               currRing); //             divideByGcd(I->m[j]); //           } //         } //       } //       for (int i=1; im[i],p)) //           return TRUE; //       } //       res->rtyp = NONE; //       res->data = NULL; //       IDDATA((idhdl)u->data) = (char*) I; //       return FALSE; //     } //   } //   WerrorS("initialReduction: unexpected parameters"); //   return TRUE; // } p->iiAddCproc("","initial",FALSE,initial); #ifndef NDEBUG p->iiAddCproc("","divideByGcd",FALSE,divideByGcd); p->iiAddCproc("","preduce",FALSE,preduce); // p->iiAddCproc("","divideByGcd",FALSE,divideByGcd); p->iiAddCproc("","pReduce",FALSE,pReduce); #endif //NDEBUG p->iiAddCproc("","initialReduction",FALSE,initialReduction); // p->iiAddCproc("","initialReduction",FALSE,initialReduction); p->iiAddCproc("","homogeneitySpace",FALSE,homogeneitySpace); }
Note: See TracChangeset for help on using the changeset viewer.