source: git/Singular/dyn_modules/gfanlib/tropical.cc @ d770e6

#include <kernel/polys.h>
#include <kernel/kstd1.h>
#include <libpolys/coeffs/longrat.h>
#include <libpolys/polys/clapsing.h>
#include <bbcone.h>


poly initial(poly p)
{
  poly g = p;
  poly h = p_Head(g,currRing);
  poly f = h;
  long d = p_Deg(g,currRing);
  pIter(g);
  while ((g != NULL) && (p_Deg(g,currRing) == d))
  {
    pNext(h) = p_Head(g,currRing);
    pIter(h);
    pIter(g);
  }
  return(f);
}


BOOLEAN initial(leftv res, leftv args)
{
  leftv u = args;
  if ((u != NULL) && (u->Typ() == POLY_CMD))
  {
    leftv v = u->next;
    if (v == NULL)
    {
      poly p = (poly) u->Data();
      res->rtyp = POLY_CMD;
      res->data = (void*) initial(p);
      return FALSE;
    }
  }
  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
  {
    leftv v = u->next;
    if (v == NULL)
    {
      ideal I = (ideal) u->Data();
      ideal inI = idInit(IDELEMS(I));
      for (int i=0; i<IDELEMS(I); i++)
        inI->m[i]=initial(I->m[i]);
      res->rtyp = IDEAL_CMD;
      res->data = (void*) inI;
      return FALSE;
    }
  }
  WerrorS("initial: unexpected parameters");
  return TRUE;
}


BOOLEAN homogeneitySpace(leftv res, leftv args)
{
  leftv u = args;
  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
  {
    leftv v = u->next;
    if (v == NULL)
    {
      int n = currRing->N;
      ideal I = (ideal) u->Data();
      poly g;
      int* leadexpv = (int*) omAlloc((n+1)*sizeof(int));
      int* tailexpv = (int*) omAlloc((n+1)*sizeof(int));
      gfan::ZVector leadexpw = gfan::ZVector(n);
      gfan::ZVector tailexpw = gfan::ZVector(n);
      gfan::ZMatrix equations = gfan::ZMatrix(0,n);
      for (int i=0; i<IDELEMS(I); i++)
      {
        g = (poly) I->m[i]; pGetExpV(g,leadexpv);
        leadexpw = intStar2ZVector(n, leadexpv);
        pIter(g);
        while (g != NULL)
        {
          pGetExpV(g,tailexpv);
          tailexpw = intStar2ZVector(n, tailexpv);
          equations.appendRow(leadexpw-tailexpw);
          pIter(g);
        }
      }
      gfan::ZCone* gCone = new gfan::ZCone(gfan::ZMatrix(0, equations.getWidth()),equations);
      omFreeSize(leadexpv,(n+1)*sizeof(int));
      omFreeSize(tailexpv,(n+1)*sizeof(int));

      res->rtyp = coneID;
      res->data = (void*) gCone;
      return FALSE;
    }
  }
  WerrorS("homogeneitySpace: unexpected parameters");
  return TRUE;
}


BOOLEAN groebnerCone(leftv res, leftv args)
{
  leftv u = args;
  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
  {
    leftv v = u->next;
    if (v == NULL)
    {
      int n = currRing->N;
      ideal I = (ideal) u->Data();
      poly g = NULL;
      int* leadexpv = (int*) omAlloc((n+1)*sizeof(int));
      int* tailexpv = (int*) omAlloc((n+1)*sizeof(int));
      gfan::ZVector leadexpw = gfan::ZVector(n);
      gfan::ZVector tailexpw = gfan::ZVector(n);
      gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
      gfan::ZMatrix equations = gfan::ZMatrix(0,n);
      long d;
      for (int i=0; i<IDELEMS(I); i++)
      {
        g = (poly) I->m[i]; pGetExpV(g,leadexpv);
        leadexpw = intStar2ZVector(n, leadexpv);
        pIter(g);
        d = p_Deg(g,currRing);
        while ((g != NULL) && (p_Deg(g,currRing) == d))
        {
          pGetExpV(g,tailexpv);
          tailexpw = intStar2ZVector(n, tailexpv);
          equations.appendRow(leadexpw-tailexpw);
          pIter(g);
        }

        if (g != NULL)
        {
          while (g != NULL)
          {
            pGetExpV(g,tailexpv);
            tailexpw = intStar2ZVector(n, tailexpv);
            inequalities.appendRow(leadexpw-tailexpw);
            pIter(g);
          }
        }
      }
      gfan::ZCone* gCone = new gfan::ZCone(inequalities,equations);
      omFreeSize(leadexpv,(n+1)*sizeof(int));
      omFreeSize(tailexpv,(n+1)*sizeof(int));

      res->rtyp = coneID;
      res->data = (void*) gCone;
      return FALSE;
    }
  }
  WerrorS("groebnerCone: unexpected parameters");
  return TRUE;
}


gfan::ZCone* maximalGroebnerCone(const ring &r, const ideal &I)
{
  int n = rVar(r);
  poly g = NULL;
  int* leadexpv = (int*) omAlloc((n+1)*sizeof(int));
  int* tailexpv = (int*) omAlloc((n+1)*sizeof(int));
  gfan::ZVector leadexpw = gfan::ZVector(n);
  gfan::ZVector tailexpw = gfan::ZVector(n);
  gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
  for (int i=0; i<IDELEMS(I); i++)
  {
    g = (poly) I->m[i]; pGetExpV(g,leadexpv);
    leadexpw = intStar2ZVector(n, leadexpv);
    pIter(g);
    while (g != NULL)
    {
      pGetExpV(g,tailexpv);
      tailexpw = intStar2ZVector(n, tailexpv);
      inequalities.appendRow(leadexpw-tailexpw);
      pIter(g);
    }
  }
  omFreeSize(leadexpv,(n+1)*sizeof(int));
  omFreeSize(tailexpv,(n+1)*sizeof(int));
  return new gfan::ZCone(inequalities,gfan::ZMatrix(0, inequalities.getWidth()));
}


BOOLEAN maximalGroebnerCone(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();
      res->rtyp = coneID;
      res->data = (void*) maximalGroebnerCone(currRing, I);
      return FALSE;
    }
  }
  WerrorS("maximalGroebnerCone: unexpected parameters");
  return TRUE;
}

/***
 * If 1, replaces all occuring t with prime p,
 * where theoretically feasable.
 * Also computes GCD over integers than
 * over univariate polynomials
 **/
#define T_TO_P 0

/***
 * Suppose I=g_0,...,g_{n-1}.
 * This function uses bubble sort to sorts g_1,...,g_{n-1}
 * such that lm(g_1)>...>lm(g_{n-1})
 **/
static inline void sortingLaterGenerators(ideal I)
{
  poly cache; int i; int n=IDELEMS(I); int newn;
  do
  {
    newn=0;
    for (i=2; i<n; i++)
    {
      if (pLmCmp(I->m[i-1],I->m[i])<0)
      {
        cache=I->m[i-1];
        I->m[i-1]=I->m[i];
        I->m[i]=cache;
        newn = i;
      }
    }
    n=newn;
  } while(n);
}


/***
 * 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; i<IDELEMS(I); i++)
      {
        if (preduce(I->m[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 i<j
       **/
      int i,j;
      poly uniti, unitj;
      for (i=1; i<IDELEMS(I)-1; i++)
      {
        for (j=i+1; j<IDELEMS(I); j++)
        {
          unitj = highestMatchingX(I->m[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; i<IDELEMS(I); i++)
      {
        if (preduce(I->m[i],p))
          return TRUE;
      }

      /***
       * Step 4: the second pass. removing terms divisible by lt(g_j) out of g_i for i<j
       **/
      for (i=1; i<IDELEMS(I)-1; i++)
      {
        for (j=i+1; j<IDELEMS(I); j++)
        {
          uniti = highestMatchingX(I->m[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; i<IDELEMS(I); i++)
      {
        if (preduce(I->m[i],p))
          return TRUE;
      }

      res->rtyp = NONE;
      res->data = NULL;
      IDDATA((idhdl)u->data) = (char*) I;
      return FALSE;
    }
  }
  WerrorS("initialReduction: unexpected parameters");
  return TRUE;
}


#if 0
/***
 * Given a general ring r with any ordering,
 * changes the ordering to a(v),ws(-w)
 **/
bool changetoAWSRing(ring r, gfan::ZVector v, gfan::ZVector w)
{
  omFree(r->order);
  r->order  = (int*) omAlloc0(4*sizeof(int));
  omFree(r->block0);
  r->block0 = (int*) omAlloc0(4*sizeof(int));
  omFree(r->block1);
  r->block1 = (int*) omAlloc0(4*sizeof(int));
  for (int i=0; r->wvhdl[i]; i++)
  { omFree(r->wvhdl[i]); }
  omFree(r->wvhdl);
  r->wvhdl  = (int**) omAlloc0(4*sizeof(int*));

  bool ok = false;
  r->order[0]  = ringorder_a;
  r->block0[0] = 1;
  r->block1[0] = r->N;
  r->wvhdl[0]  = ZVectorToIntStar(v,ok);
  r->order[1]  = ringorder_ws;
  r->block0[1] = 1;
  r->block1[1] = r->N;
  r->wvhdl[1]  = ZVectorToIntStar(w,ok);
  r->order[2]=ringorder_C;
  return ok;
}


/***
 * Given a ring with ordering a(v'),ws(w'),
 * changes the weights to v,w
 **/
bool changeAWSWeights(ring r, gfan::ZVector v, gfan::ZVector w)
{
  omFree(r->wvhdl[0]);
  omFree(r->wvhdl[1]);
  bool ok = false;
  r->wvhdl[0]  = ZVectorToIntStar(v,ok);
  r->wvhdl[1]  = ZVectorToIntStar(w,ok);
  return ok;
}


// /***
//  * Creates an int* representing the transposition of the last two variables
//  **/
// static inline int* createPermutationVectorForSaturation(static const ring &r)
// {
//   int* w = (int*) omAlloc0((rVar(r)+1)*sizeof(int));
//   for (int i=1; i<=rVar(r)-2; i++)
//     w[i] = i;
//   w[rVar(r)-1] = rVar(r);
//   w[rVar(r)] = rVar(r)-1;
// }


/***
 * Creates an int* representing the permutation
 * 1 -> 1, ..., i-1 -> i-1, i -> n, i+1 -> n-1, ... , n -> i
 **/
static inline int* createPermutationVectorForSaturation(const ring &r, const int i)
{
  int* sigma = (int*) omAlloc0((rVar(r)+1)*sizeof(int));
  int j;
  for (j=1; j<i; j++)
    sigma[j] = j;
  for (; j<=rVar(r); j++)
    sigma[j] = rVar(r)-j+i;
  return(sigma);
}


/***
 * Changes the int* representing the permutation
 * 1 -> 1, ..., i -> i, i+1 -> n, i+2 -> n-1, ... , n -> i+1
 * to an int* representing the permutation
 * 1 -> 1, ..., i-1 -> i-1, i -> n, i+1 -> n-1, ... , n -> i
 **/
static void changePermutationVectorForSaturation(int* sigma, const ring &r, const int i)
{
  for (int j=i; j<rVar(r); j++)
    sigma[j] = rVar(r)-j+i;
  sigma[rVar(r)] = i;
}


/***
 * returns a ring in which the weights of the ring variables are permuted
 * if handed over a poly in which the variables are permuted, this is basically
 * as good as permuting the variables of the ring itself.
 **/
static ring permuteWeighstOfRingVariables(const ring &r, const int* const sigma)
{
  ring s = rCopy0(r);
  for (int j=0; j<rVar(r); j++)
  {
    s->wvhdl[0][j] = r->wvhdl[0][sigma[j+1]];
    s->wvhdl[1][j] = r->wvhdl[1][sigma[j+1]];
  }
  rComplete(s,1);
  return s;
}


/***
 * creates a ring s that is a copy of r except with ordering wp(w)
 **/
static inline ring createInitialRingForSaturation(const ring &r, const gfan::ZVector &w, bool &ok)
{
  assume(rVar(r) == (int) w.size());

  ring s = rCopy0(r); int i;
  for (i=0; s->order[i]; i++)
    omFreeSize(s->wvhdl[i],rVar(r)*sizeof(int));
  i++;
  omFreeSize(s->order,i*sizeof(int));
  s->order  = (int*) omAlloc0(3*sizeof(int));
  omFreeSize(s->block0,i*sizeof(int));
  s->block0 = (int*) omAlloc0(3*sizeof(int));
  omFreeSize(s->block1,i*sizeof(int));
  s->block1 = (int*) omAlloc0(3*sizeof(int));
  omFreeSize(s->wvhdl,i*sizeof(int*));
  s->wvhdl  = (int**) omAlloc0(3*sizeof(int*));

  s->order[0]  = ringorder_wp;
  s->block0[0] = 1;
  s->block1[0] = rVar(r);
  s->wvhdl[0]  = ZVectorToIntStar(w,ok);
  s->order[1]=ringorder_C;

  rComplete(s,1);
  return s;
}


/***
 * Given an weighted homogeneous ideal I with respect to weight w
 * that in standard basis form with respect to the ordering ws(-w),
 * derives the standard basis of I:<x_n>^\infty
 * and returns a long k such that I:<x_n>^\infty=I:<x_n>^k
 **/
static long deriveStandardBasisOfSaturation(ideal &I, ring &r)
{
  long k=0, l; poly current;
  for (int i=0; i<IDELEMS(I); i++)
  {
    current = I->m[i];
    l = p_GetExp(current,rVar(r),r);
    if (k<l) k=l;
    while (current)
    {
      p_SubExp(current,rVar(r),l,r); p_Setm(current,r);
      pIter(current);
    }
  }
  return k;
}


/***
 * Given a weighted homogeneous ideal I with respect to weight w
 * with constant first element,
 * returns NULL if I does not contain a monomial
 * otherwise returns the monomial contained in I
 **/
poly containsMonomial(const ideal &I, const gfan::ZVector &w)
{
  assume(rField_is_Ring_Z(currRing));

  // first we switch to the ground field currRing->cf / I->m[0]
  ring r = rCopy0(currRing);
  nKillChar(r->cf);
  r->cf = nInitChar(n_Zp,(void*)(long)n_Int(p_GetCoeff(I->m[0],currRing),currRing->cf));
  rComplete(r);

  ideal J = id_Copy(I, currRing); poly cache; number temp;
  for (int i=0; i<IDELEMS(I); i++)
  {
    cache = J->m[i];
    while (cache)
    {
      // TODO: temp = npMapGMP(p_GetCoeff(cache,currRing),currRing->cf,r->cf);
      p_SetCoeff(cache,temp,r); pIter(cache);
    }
  }


  J = kStd(J,NULL,isHomog,NULL);

  bool b = false;
  ring s = createInitialRingForSaturation(currRing, w, b);
  if (b)
  {
    WerrorS("containsMonomial: overflow in weight vector");
    return NULL;
  }

  return NULL;
}


gfan::ZCone* startingCone(ideal I)
{
  I = kStd(I,NULL,isNotHomog,NULL);
  gfan::ZCone* zc = maximalGroebnerCone(currRing,I);
  gfan::ZMatrix rays = zc->extremeRays();
  gfan::ZVector v;
  for (int i=0; i<rays.getHeight(); i++)
  {
    v = rays[i];
  }
  return zc;
}
#endif


void tropical_setup(SModulFunctions* p)
{
  p->iiAddCproc("","groebnerCone",FALSE,groebnerCone);
  p->iiAddCproc("","maximalGroebnerCone",FALSE,maximalGroebnerCone);
  p->iiAddCproc("","initial",FALSE,initial);
#ifndef NDEBUG
  p->iiAddCproc("","divideByGcd",FALSE,divideByGcd);
  p->iiAddCproc("","preduce",FALSE,preduce);
#endif //NDEBUG
  p->iiAddCproc("","initialReduction",FALSE,initialReduction);
  p->iiAddCproc("","homogeneitySpace",FALSE,homogeneitySpace);
}