Changeset d770e6 in git


Ignore:
Timestamp:
Oct 29, 2013, 1:54:17 PM (9 years ago)
Author:
Yue Ren <ren@…>
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'ad2543eab51733612ba7d118afc77edca719600e')
Children:
7395dd26103d1edd6fc4505722010d2bba012d96
Parents:
a73fae59ddd878f87c48e12088f6b0046549aaf7
git-author:
Yue Ren <ren@mathematik.uni-kl.de>2013-10-29 13:54:17+01:00
git-committer:
Yue Ren <ren@mathematik.uni-kl.de>2015-02-06 13:47:01+01:00
Message:
new: initialReduction
Location:
Singular/dyn_modules/gfanlib
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • Singular/dyn_modules/gfanlib/bbcone.cc

    ra73fae rd770e6  
    113113    zv[j]=i[j+1];
    114114  return zv;
     115}
     116
     117int* ZVectorToIntStar(const gfan::ZVector &v, bool &overflow)
     118{
     119  int* w = (int*) omAlloc(v.size()*sizeof(int));
     120  for (unsigned i=0; i<v.size(); i++)
     121  {
     122    if (!v[i].fitsInInt())
     123    {
     124      omFree(w);
     125      WerrorS("intoverflow converting gfan:ZVector to int*");
     126      overflow = true;
     127      return NULL;
     128    }
     129    w[i]=v[i].toInt();
     130  }
     131  return w;
    115132}
    116133
  • Singular/dyn_modules/gfanlib/bbcone.h

    ra73fae rd770e6  
    2727
    2828gfan::ZVector intStar2ZVector(const int d, const int* i);
     29int* ZVectorToIntStar(const gfan::ZVector &v, bool &overflow);
    2930char* toString(gfan::ZMatrix const &m);
    3031std::string toString(const gfan::ZCone* const c);
  • Singular/dyn_modules/gfanlib/tropical.cc

    ra73fae rd770e6  
    11#include <kernel/polys.h>
     2#include <kernel/kstd1.h>
    23#include <libpolys/coeffs/longrat.h>
     4#include <libpolys/polys/clapsing.h>
    35#include <bbcone.h>
    46
     
    154156
    155157
     158gfan::ZCone* maximalGroebnerCone(const ring &r, const ideal &I)
     159{
     160  int n = rVar(r);
     161  poly g = NULL;
     162  int* leadexpv = (int*) omAlloc((n+1)*sizeof(int));
     163  int* tailexpv = (int*) omAlloc((n+1)*sizeof(int));
     164  gfan::ZVector leadexpw = gfan::ZVector(n);
     165  gfan::ZVector tailexpw = gfan::ZVector(n);
     166  gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
     167  for (int i=0; i<IDELEMS(I); i++)
     168  {
     169    g = (poly) I->m[i]; pGetExpV(g,leadexpv);
     170    leadexpw = intStar2ZVector(n, leadexpv);
     171    pIter(g);
     172    while (g != NULL)
     173    {
     174      pGetExpV(g,tailexpv);
     175      tailexpw = intStar2ZVector(n, tailexpv);
     176      inequalities.appendRow(leadexpw-tailexpw);
     177      pIter(g);
     178    }
     179  }
     180  omFreeSize(leadexpv,(n+1)*sizeof(int));
     181  omFreeSize(tailexpv,(n+1)*sizeof(int));
     182  return new gfan::ZCone(inequalities,gfan::ZMatrix(0, inequalities.getWidth()));
     183}
     184
     185
    156186BOOLEAN maximalGroebnerCone(leftv res, leftv args)
    157187{
     
    162192    if (v == NULL)
    163193    {
    164       int n = currRing->N;
    165194      ideal I = (ideal) u->Data();
    166       poly g = NULL;
    167       int* leadexpv = (int*) omAlloc((n+1)*sizeof(int));
    168       int* tailexpv = (int*) omAlloc((n+1)*sizeof(int));
    169       gfan::ZVector leadexpw = gfan::ZVector(n);
    170       gfan::ZVector tailexpw = gfan::ZVector(n);
    171       gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
    172       for (int i=0; i<IDELEMS(I); i++)
    173       {
    174         g = (poly) I->m[i]; pGetExpV(g,leadexpv);
    175         leadexpw = intStar2ZVector(n, leadexpv);
    176         pIter(g);
    177         while (g != NULL)
    178         {
    179           pGetExpV(g,tailexpv);
    180           tailexpw = intStar2ZVector(n, tailexpv);
    181           inequalities.appendRow(leadexpw-tailexpw);
    182           pIter(g);
    183         }
    184       }
    185       gfan::ZCone* gCone = new gfan::ZCone(inequalities,gfan::ZMatrix(0, inequalities.getWidth()));
    186       omFreeSize(leadexpv,(n+1)*sizeof(int));
    187       omFreeSize(tailexpv,(n+1)*sizeof(int));
    188 
    189195      res->rtyp = coneID;
    190       res->data = (void*) gCone;
     196      res->data = (void*) maximalGroebnerCone(currRing, I);
    191197      return FALSE;
    192198    }
     
    195201  return TRUE;
    196202}
     203
     204/***
     205 * If 1, replaces all occuring t with prime p,
     206 * where theoretically feasable.
     207 * Also computes GCD over integers than
     208 * over univariate polynomials
     209 **/
     210#define T_TO_P 0
     211
     212/***
     213 * Suppose I=g_0,...,g_{n-1}.
     214 * This function uses bubble sort to sorts g_1,...,g_{n-1}
     215 * such that lm(g_1)>...>lm(g_{n-1})
     216 **/
     217static inline void sortingLaterGenerators(ideal I)
     218{
     219  poly cache; int i; int n=IDELEMS(I); int newn;
     220  do
     221  {
     222    newn=0;
     223    for (i=2; i<n; i++)
     224    {
     225      if (pLmCmp(I->m[i-1],I->m[i])<0)
     226      {
     227        cache=I->m[i-1];
     228        I->m[i-1]=I->m[i];
     229        I->m[i]=cache;
     230        newn = i;
     231      }
     232    }
     233    n=newn;
     234  } while(n);
     235}
     236
     237
     238/***
     239 * replaces coefficients in g of lowest degree in t
     240 * divisible by p with a suitable power of t
     241 **/
     242static 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
     313BOOLEAN 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 **/
     342static 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 **/
     369static 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 **/
     408static 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
     419static 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 **/
     533static 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
     627BOOLEAN 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
     649BOOLEAN 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}
     743
     744
     745#if 0
     746/***
     747 * Given a general ring r with any ordering,
     748 * changes the ordering to a(v),ws(-w)
     749 **/
     750bool changetoAWSRing(ring r, gfan::ZVector v, gfan::ZVector w)
     751{
     752  omFree(r->order);
     753  r->order  = (int*) omAlloc0(4*sizeof(int));
     754  omFree(r->block0);
     755  r->block0 = (int*) omAlloc0(4*sizeof(int));
     756  omFree(r->block1);
     757  r->block1 = (int*) omAlloc0(4*sizeof(int));
     758  for (int i=0; r->wvhdl[i]; i++)
     759  { omFree(r->wvhdl[i]); }
     760  omFree(r->wvhdl);
     761  r->wvhdl  = (int**) omAlloc0(4*sizeof(int*));
     762
     763  bool ok = false;
     764  r->order[0]  = ringorder_a;
     765  r->block0[0] = 1;
     766  r->block1[0] = r->N;
     767  r->wvhdl[0]  = ZVectorToIntStar(v,ok);
     768  r->order[1]  = ringorder_ws;
     769  r->block0[1] = 1;
     770  r->block1[1] = r->N;
     771  r->wvhdl[1]  = ZVectorToIntStar(w,ok);
     772  r->order[2]=ringorder_C;
     773  return ok;
     774}
     775
     776
     777/***
     778 * Given a ring with ordering a(v'),ws(w'),
     779 * changes the weights to v,w
     780 **/
     781bool changeAWSWeights(ring r, gfan::ZVector v, gfan::ZVector w)
     782{
     783  omFree(r->wvhdl[0]);
     784  omFree(r->wvhdl[1]);
     785  bool ok = false;
     786  r->wvhdl[0]  = ZVectorToIntStar(v,ok);
     787  r->wvhdl[1]  = ZVectorToIntStar(w,ok);
     788  return ok;
     789}
     790
     791
     792// /***
     793//  * Creates an int* representing the transposition of the last two variables
     794//  **/
     795// static inline int* createPermutationVectorForSaturation(static const ring &r)
     796// {
     797//   int* w = (int*) omAlloc0((rVar(r)+1)*sizeof(int));
     798//   for (int i=1; i<=rVar(r)-2; i++)
     799//     w[i] = i;
     800//   w[rVar(r)-1] = rVar(r);
     801//   w[rVar(r)] = rVar(r)-1;
     802// }
     803
     804
     805/***
     806 * Creates an int* representing the permutation
     807 * 1 -> 1, ..., i-1 -> i-1, i -> n, i+1 -> n-1, ... , n -> i
     808 **/
     809static inline int* createPermutationVectorForSaturation(const ring &r, const int i)
     810{
     811  int* sigma = (int*) omAlloc0((rVar(r)+1)*sizeof(int));
     812  int j;
     813  for (j=1; j<i; j++)
     814    sigma[j] = j;
     815  for (; j<=rVar(r); j++)
     816    sigma[j] = rVar(r)-j+i;
     817  return(sigma);
     818}
     819
     820
     821/***
     822 * Changes the int* representing the permutation
     823 * 1 -> 1, ..., i -> i, i+1 -> n, i+2 -> n-1, ... , n -> i+1
     824 * to an int* representing the permutation
     825 * 1 -> 1, ..., i-1 -> i-1, i -> n, i+1 -> n-1, ... , n -> i
     826 **/
     827static void changePermutationVectorForSaturation(int* sigma, const ring &r, const int i)
     828{
     829  for (int j=i; j<rVar(r); j++)
     830    sigma[j] = rVar(r)-j+i;
     831  sigma[rVar(r)] = i;
     832}
     833
     834
     835/***
     836 * returns a ring in which the weights of the ring variables are permuted
     837 * if handed over a poly in which the variables are permuted, this is basically
     838 * as good as permuting the variables of the ring itself.
     839 **/
     840static ring permuteWeighstOfRingVariables(const ring &r, const int* const sigma)
     841{
     842  ring s = rCopy0(r);
     843  for (int j=0; j<rVar(r); j++)
     844  {
     845    s->wvhdl[0][j] = r->wvhdl[0][sigma[j+1]];
     846    s->wvhdl[1][j] = r->wvhdl[1][sigma[j+1]];
     847  }
     848  rComplete(s,1);
     849  return s;
     850}
     851
     852
     853/***
     854 * creates a ring s that is a copy of r except with ordering wp(w)
     855 **/
     856static inline ring createInitialRingForSaturation(const ring &r, const gfan::ZVector &w, bool &ok)
     857{
     858  assume(rVar(r) == (int) w.size());
     859
     860  ring s = rCopy0(r); int i;
     861  for (i=0; s->order[i]; i++)
     862    omFreeSize(s->wvhdl[i],rVar(r)*sizeof(int));
     863  i++;
     864  omFreeSize(s->order,i*sizeof(int));
     865  s->order  = (int*) omAlloc0(3*sizeof(int));
     866  omFreeSize(s->block0,i*sizeof(int));
     867  s->block0 = (int*) omAlloc0(3*sizeof(int));
     868  omFreeSize(s->block1,i*sizeof(int));
     869  s->block1 = (int*) omAlloc0(3*sizeof(int));
     870  omFreeSize(s->wvhdl,i*sizeof(int*));
     871  s->wvhdl  = (int**) omAlloc0(3*sizeof(int*));
     872
     873  s->order[0]  = ringorder_wp;
     874  s->block0[0] = 1;
     875  s->block1[0] = rVar(r);
     876  s->wvhdl[0]  = ZVectorToIntStar(w,ok);
     877  s->order[1]=ringorder_C;
     878
     879  rComplete(s,1);
     880  return s;
     881}
     882
     883
     884/***
     885 * Given an weighted homogeneous ideal I with respect to weight w
     886 * that in standard basis form with respect to the ordering ws(-w),
     887 * derives the standard basis of I:<x_n>^\infty
     888 * and returns a long k such that I:<x_n>^\infty=I:<x_n>^k
     889 **/
     890static long deriveStandardBasisOfSaturation(ideal &I, ring &r)
     891{
     892  long k=0, l; poly current;
     893  for (int i=0; i<IDELEMS(I); i++)
     894  {
     895    current = I->m[i];
     896    l = p_GetExp(current,rVar(r),r);
     897    if (k<l) k=l;
     898    while (current)
     899    {
     900      p_SubExp(current,rVar(r),l,r); p_Setm(current,r);
     901      pIter(current);
     902    }
     903  }
     904  return k;
     905}
     906
     907
     908/***
     909 * Given a weighted homogeneous ideal I with respect to weight w
     910 * with constant first element,
     911 * returns NULL if I does not contain a monomial
     912 * otherwise returns the monomial contained in I
     913 **/
     914poly containsMonomial(const ideal &I, const gfan::ZVector &w)
     915{
     916  assume(rField_is_Ring_Z(currRing));
     917
     918  // first we switch to the ground field currRing->cf / I->m[0]
     919  ring r = rCopy0(currRing);
     920  nKillChar(r->cf);
     921  r->cf = nInitChar(n_Zp,(void*)(long)n_Int(p_GetCoeff(I->m[0],currRing),currRing->cf));
     922  rComplete(r);
     923
     924  ideal J = id_Copy(I, currRing); poly cache; number temp;
     925  for (int i=0; i<IDELEMS(I); i++)
     926  {
     927    cache = J->m[i];
     928    while (cache)
     929    {
     930      // TODO: temp = npMapGMP(p_GetCoeff(cache,currRing),currRing->cf,r->cf);
     931      p_SetCoeff(cache,temp,r); pIter(cache);
     932    }
     933  }
     934
     935
     936  J = kStd(J,NULL,isHomog,NULL);
     937
     938  bool b = false;
     939  ring s = createInitialRingForSaturation(currRing, w, b);
     940  if (b)
     941  {
     942    WerrorS("containsMonomial: overflow in weight vector");
     943    return NULL;
     944  }
     945
     946  return NULL;
     947}
     948
     949
     950gfan::ZCone* startingCone(ideal I)
     951{
     952  I = kStd(I,NULL,isNotHomog,NULL);
     953  gfan::ZCone* zc = maximalGroebnerCone(currRing,I);
     954  gfan::ZMatrix rays = zc->extremeRays();
     955  gfan::ZVector v;
     956  for (int i=0; i<rays.getHeight(); i++)
     957  {
     958    v = rays[i];
     959  }
     960  return zc;
     961}
     962#endif
    197963
    198964
     
    202968  p->iiAddCproc("","maximalGroebnerCone",FALSE,maximalGroebnerCone);
    203969  p->iiAddCproc("","initial",FALSE,initial);
     970#ifndef NDEBUG
     971  p->iiAddCproc("","divideByGcd",FALSE,divideByGcd);
     972  p->iiAddCproc("","preduce",FALSE,preduce);
     973#endif //NDEBUG
     974  p->iiAddCproc("","initialReduction",FALSE,initialReduction);
    204975  p->iiAddCproc("","homogeneitySpace",FALSE,homogeneitySpace);
    205976}
Note: See TracChangeset for help on using the changeset viewer.