Changeset 4664b3 in git


Ignore:
Timestamp:
Oct 6, 2014, 11:52:44 AM (10 years ago)
Author:
Yue Ren <ren@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
c9641474cc5811d6eed949bf4509b30db0cbb721
Parents:
dffd1541cc74564f509bd5d284948852484430b0
git-author:
Yue Ren <ren@mathematik.uni-kl.de>2014-10-06 12:52:44+03:00
git-committer:
Yue Ren <ren@mathematik.uni-kl.de>2015-02-06 13:47:05+01:00
Message:
chg: status update 06.10.
Location:
Singular/dyn_modules/gfanlib
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • Singular/dyn_modules/gfanlib/groebnerCone.cc

    rdffd154 r4664b3  
    121121  polyhedralCone.canonicalize();
    122122  interiorPoint = polyhedralCone.getRelativeInteriorPoint();
    123   initialPolynomialIdeal = sloppyInitial(reducedPolynomialIdeal,polynomialRing,interiorPoint);
     123  initialPolynomialIdeal = initial(reducedPolynomialIdeal,polynomialRing,interiorPoint);
    124124  assume(checkOrderingAndCone(polynomialRing,polyhedralCone));
    125125}
    … …  
    186186  polyhedralCone.canonicalize();
    187187  interiorPoint = polyhedralCone.getRelativeInteriorPoint();
    188   initialPolynomialIdeal = sloppyInitial(reducedPolynomialIdeal,polynomialRing,interiorPoint);
     188  initialPolynomialIdeal = initial(reducedPolynomialIdeal,polynomialRing,interiorPoint);
    189189  assume(checkOrderingAndCone(polynomialRing,polyhedralCone));
    190190}
    … …  
    249249  polyhedralCone.canonicalize();
    250250  interiorPoint = polyhedralCone.getRelativeInteriorPoint();
    251   initialPolynomialIdeal = sloppyInitial(reducedPolynomialIdeal,polynomialRing,interiorPoint);
     251  initialPolynomialIdeal = initial(reducedPolynomialIdeal,polynomialRing,interiorPoint);
    252252  assume(checkOrderingAndCone(polynomialRing,polyhedralCone));
    253253}

  • Singular/dyn_modules/gfanlib/initial.cc

    rdffd154 r4664b3  
     1#include <kernel/ideals.h>
     2#include <libpolys/polys/monomials/p_polys.h>
     3
    14#include <gfanlib/gfanlib.h>
    25
    3 #include <kernel/ideals.h>
    4 #include <Singular/subexpr.h>
    5 #include <libpolys/polys/monomials/p_polys.h>
    6 #include <libpolys/polys/simpleideals.h>
    7 
    8 #include <callgfanlib_conversion.h>
    9 #include <gfanlib_exceptions.h>
    10 
    116#include <exception>
    127
    13 /***
    14  * Computes the weighted degree of the leading term of p with respect to w
    15  **/
    168long wDeg(const poly p, const ring r, const gfan::ZVector w)
    179{
    … …  
    2012  {
    2113    if (!w[i].fitsInInt())
    22       throw 0; //weightOverflow;
     14      throw 0; // weightOverflow;
    2315    d += p_GetExp(p,i+1,r)*w[i].toInt();
    2416  }
    25   return d;
    26 }
    27 
    28 /***
    29  * Computes the weighted multidegree of the leading term of p with respect to W.
    30  * The weighted multidegree is a vector whose i-th entry is the weighted degree
    31  * with respect to the i-th row vector of W.
    32  **/
    33 gfan::ZVector WDeg(const poly p, const ring r, const gfan::ZMatrix W)
    34 {
    35   gfan::ZVector d = gfan::ZVector(W.getHeight());
    36   for (int i=0; i<W.getHeight(); i++)
    37     d[i] = wDeg(p,r,W[i]);
    3817  return d;
    3918}
    … …  
    4827}
    4928
    50 /**
    51  * Checks if p is sorted with respect to w.
    52  */
    53 static bool checkSloppyInput(const poly p, const ring r, const gfan::ZVector w)
    54 {
    55   long d = wDeg(p,r,w);
    56   for (poly currentTerm = p->next; currentTerm; pIter(currentTerm))
    57   {
    58     long e = wDeg(currentTerm,r,w);
    59     if (e>d)
    60       return false;
    61   }
    62   return true;
    63 }
    64 
    65 /***
    66  * Returns the terms of p of same weighted degree under w as the leading term.
    67  * Coincides with the initial form of p with respect to w if and only if p was already
    68  * sorted with respect to w in the sense that the leading term is of highest w-degree.
    69  **/
    70 poly sloppyInitial(const poly p, const ring r, const gfan::ZVector w)
    71 {
    72   assume(checkSloppyInput(p,r,w));
    73   int n = rVar(r);
    74   int* expv = (int*) omAlloc(n*sizeof(int));
     29poly initial(const poly p, const ring r, const gfan::ZVector w)
     30{
     31  if (p==NULL)
     32    return NULL;
     33
    7534  poly q0 = p_Head(p,r);
    7635  poly q1 = q0;
    … …  
    7837  for (poly currentTerm = p->next; currentTerm; pIter(currentTerm))
    7938  {
    80     if (wDeg(currentTerm,r,w) == d)
    81     {
    82       pNext(q1) = p_Head(currentTerm,r);
    83       pIter(q1);
    84     }
    85   }
    86   omFreeSize(expv,n*sizeof(int));
    87   return q0;
    88 }
    89 
    90 /***
    91  * Runs the above procedure over all generators of an ideal.
    92  * Coincides with the initial ideal of I with respect to w if and only if
    93  * the elements of I were already sorted with respect to w and
    94  * I is a standard basis form with respect to w.
    95  **/
    96 ideal sloppyInitial(const ideal I, const ring r, const gfan::ZVector w)
    97 {
    98   int k = idSize(I); ideal inI = idInit(k);
    99   for (int i=0; i<k; i++)
    100     inI->m[i] = sloppyInitial(I->m[i],r,w);
    101   return inI;
    102 }
    103 
    104 /***
    105  * Returns the initial form of p with respect to w
    106  **/
    107 poly initial(const poly p, const ring r, const gfan::ZVector w)
    108 {
    109   poly q0 = p_Head(p,r);
    110   poly q1 = q0;
    111   long d = wDeg(p,r,w);
    112   for (poly currentTerm = p->next; currentTerm; pIter(currentTerm))
    113   {
    11439    long e = wDeg(currentTerm,r,w);
    115     if (e>d)
     40    if (d<e)
    11641    {
    11742      p_Delete(&q0,r);
    … …  
    13055}
    13156
    132 /***
    133  * Runs the above procedure over all generators of an ideal.
    134  * Returns the initial ideal if and only if the weight is in the maximal Groebner cone
    135  * of the current ordering.
    136  **/
    13757ideal initial(const ideal I, const ring r, const gfan::ZVector w)
    13858{
    139   idSkipZeroes(I);
    14059  int k = idSize(I); ideal inI = idInit(k);
    14160  for (int i=0; i<k; i++)
    … …  
    14463}
    14564
    146 
    147 /***
    148  * Returns the initial form of p with respect to W,
    149  * i.e. the sum over all terms of p with highest multidegree with respect to W.
    150  **/
    151 poly initial(const poly p, const ring r, const gfan::ZMatrix W)
    152 {
    153   int n = rVar(r);
    154   poly q0 = p_Head(p,r);
    155   poly q1 = q0;
    156   gfan::ZVector d = WDeg(p,r,W);
    157   for (poly currentTerm = p->next; currentTerm; pIter(currentTerm))
    158   {
    159     gfan::ZVector e = WDeg(currentTerm,r,W);
    160     if (d<e)
    161     {
    162       p_Delete(&q0,r);
    163       q0 = p_Head(p,r);
    164       q1 = q0;
    165       d = e;
    166     }
    167     else
    168       if (d==e)
    169       {
    170         pNext(q1) = p_Head(currentTerm,r);
    171         pIter(q1);
    172       }
    173   }
    174   return q0;
    175 }
    176 
    17765poly initial(const poly p, const ring r, const gfan::ZVector w, const gfan::ZMatrix W)
    17866{
    179   int n = rVar(r);
     67  if (p==NULL)
     68    return NULL;
     69
    18070  poly q0 = p_Head(p,r);
    18171  poly q1 = q0;
    … …  
    20191}
    20292
    203 
    204 /***
    205  * Runs the above procedure over all generators of an ideal.
    206  * Returns the initial ideal if and only if the weight is in the maximal Groebner cone
    207  * of the current ordering.
    208  **/
    209 ideal initial(const ideal I, const ring r, const gfan::ZMatrix W)
    210 {
    211   int k = idSize(I); ideal inI = idInit(k);
    212   for (int i=0; i<k; i++)
    213     inI->m[i] = initial(I->m[i],r,W);
    214   return inI;
    215 }
    216 
    21793ideal initial(const ideal I, const ring r, const gfan::ZVector w, const gfan::ZMatrix W)
    21894{
    … …  
    22399}
    224100
    225 #ifndef NDEBUG
    226 BOOLEAN initial0(leftv res, leftv args)
    227 {
    228   leftv u = args;
    229   ideal I = (ideal) u->CopyD();
    230   leftv v = u->next;
    231   bigintmat* w0 = (bigintmat*) v->Data();
    232   gfan::ZVector* w = bigintmatToZVector(w0);
    233   omUpdateInfo();
    234   Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
    235   ideal inI = initial(I,currRing,*w);
    236   id_Delete(&I,currRing);
    237   delete w;
    238   res->rtyp = IDEAL_CMD;
    239   res->data = (char*) inI;
    240   return FALSE;
    241 }
    242 #endif
    243 
     101void initial(poly* pStar, const ring r, const gfan::ZVector w)
     102{
     103  poly p = *pStar;
     104  if (p==NULL)
     105    return;
     106
     107  long d = wDeg(p,r,w);
     108  poly q0 = p;
     109  poly q1 = q0;
     110  pNext(q1) = NULL;
     111  pIter(p);
     112
     113  while(p)
     114  {
     115    long e = wDeg(p,r,w);
     116    if (d<e)
     117    {
     118      p_Delete(&q0,r);
     119      q0 = p;
     120      q1 = q0;
     121      pNext(q1) = NULL;
     122      d = e;
     123      pIter(p);
     124    }
     125    else
     126      if (e==d)
     127      {
     128        pNext(q1) = p;
     129        pIter(q1);
     130        pNext(q1) = NULL;
     131        pIter(p);
     132      }
     133      else
     134        p = p_LmDeleteAndNext(p,r);
     135  }
     136  pStar = &q0;
     137  return;
     138}
     139
     140void initial(ideal* IStar, const ring r, const gfan::ZVector w)
     141{
     142  ideal I = *IStar;
     143  int k = idSize(I);
     144  for (int i=0; i<k; i++)
     145    initial(&I->m[i],r,w);
     146  return;
     147}
     148
     149void initial(poly* pStar, const ring r, const gfan::ZVector w, const gfan::ZMatrix W)
     150{
     151  poly p = *pStar;
     152  if (p==NULL)
     153    return;
     154
     155  gfan::ZVector d = WDeg(p,r,w,W);
     156  poly q0 = p;
     157  poly q1 = q0;
     158  pNext(q1) = NULL;
     159  pIter(p);
     160
     161  while(p)
     162  {
     163    gfan::ZVector e = WDeg(p,r,w,W);
     164    if (d<e)
     165    {
     166      p_Delete(&q0,r);
     167      q0 = p;
     168      q1 = q0;
     169      pNext(q1) = NULL;
     170      d = e;
     171      pIter(p);
     172    }
     173    else
     174      if (d==e)
     175      {
     176        pNext(q1) = p;
     177        pIter(q1);
     178        pNext(q1) = NULL;
     179        pIter(p);
     180      }
     181      else
     182        p = p_LmDeleteAndNext(p,r);
     183  }
     184  pStar = &q0;
     185  return;
     186}
     187
     188void initial(ideal* IStar, const ring r, const gfan::ZVector w, const gfan::ZMatrix W)
     189{
     190  ideal I = *IStar;
     191  int k = idSize(I);
     192  for (int i=0; i<k; i++)
     193    initial(&I->m[i],r,w,W);
     194  return;
     195}
    244196
    245197/***
    246  * Computes the initial form of p with respect to the first row in the order matrix
     198 * throwaway code
    247199 **/
    248 poly initial(const poly p, const ring r)
    249 {
    250   long d = p_Deg(p,r);
    251   poly initialForm0 = p_Head(p,r);
    252   poly initialForm1 = initialForm0;
    253   poly currentTerm = p->next;
    254   while (currentTerm && p_Deg(currentTerm,r)==d)
    255   {
    256     pNext(initialForm1) = p_Head(currentTerm,r);
    257     pIter(currentTerm); pIter(initialForm1);
    258   }
    259   return initialForm0;
    260 }
    261 
    262 /***
    263  * Computes the initial form of all generators of I.
    264  * If I is a standard basis, then this is a standard basis
    265  * of the initial ideal.
    266  **/
    267 ideal initial(const ideal I, const ring r)
    268 {
    269   int k = idSize(I); ideal inI = idInit(k);
    270   for (int i=0; i<k; i++)
    271     inI->m[i] = initial(I->m[i],r);
    272   return inI;
    273 }
    274 
    275 BOOLEAN initial(leftv res, leftv args)
    276 {
    277   leftv u = args;
    278   if ((u != NULL) && (u->Typ() == POLY_CMD) && (u->next == NULL))
    279   {
    280     poly p = (poly) u->Data();
    281     res->rtyp = POLY_CMD;
    282     res->data = (void*) initial(p, currRing);
    283     return FALSE;
    284   }
    285   if ((u != NULL) && (u->Typ() == IDEAL_CMD) && (u->next == NULL))
    286   {
    287     ideal I = (ideal) u->Data();
    288     res->rtyp = IDEAL_CMD;
    289     res->data = (void*) initial(I, currRing);
    290     return FALSE;
    291   }
    292   WerrorS("initial: unexpected parameters");
    293   return TRUE;
    294 }
     200// gfan::ZVector WDeg(const poly p, const ring r, const gfan::ZMatrix W)
     201// {
     202//   gfan::ZVector d = gfan::ZVector(W.getHeight());
     203//   for (int i=0; i<W.getHeight(); i++)
     204//     d[i] = wDeg(p,r,W[i]);
     205//   return d;
     206// }
     207
     208// /**
     209//  * Checks if p is sorted with respect to w.
     210//  */
     211// static bool checkSloppyInput(const poly p, const ring r, const gfan::ZVector w)
     212// {
     213//   long d = wDeg(p,r,w);
     214//   for (poly currentTerm = p->next; currentTerm; pIter(currentTerm))
     215//   {
     216//     long e = wDeg(currentTerm,r,w);
     217//     if (e>d)
     218//       return false;
     219//   }
     220//   return true;
     221// }
     222
     223// /***
     224//  * Returns the terms of p of same weighted degree under w as the leading term.
     225//  * Coincides with the initial form of p with respect to w if and only if p was already
     226//  * sorted with respect to w in the sense that the leading term is of highest w-degree.
     227//  **/
     228// poly sloppyInitial(const poly p, const ring r, const gfan::ZVector w)
     229// {
     230//   assume(checkSloppyInput(p,r,w));
     231//   int n = rVar(r);
     232//   int* expv = (int*) omAlloc(n*sizeof(int));
     233//   poly q0 = p_Head(p,r);
     234//   poly q1 = q0;
     235//   long d = wDeg(p,r,w);
     236//   for (poly currentTerm = p->next; currentTerm; pIter(currentTerm))
     237//   {
     238//     if (wDeg(currentTerm,r,w) == d)
     239//     {
     240//       pNext(q1) = p_Head(currentTerm,r);
     241//       pIter(q1);
     242//     }
     243//   }
     244//   omFreeSize(expv,n*sizeof(int));
     245//   return q0;
     246// }
     247
     248// /***
     249//  * Runs the above procedure over all generators of an ideal.
     250//  * Coincides with the initial ideal of I with respect to w if and only if
     251//  * the elements of I were already sorted with respect to w and
     252//  * I is a standard basis form with respect to w.
     253//  **/
     254// ideal sloppyInitial(const ideal I, const ring r, const gfan::ZVector w)
     255// {
     256//   int k = idSize(I); ideal inI = idInit(k);
     257//   for (int i=0; i<k; i++)
     258//     inI->m[i] = sloppyInitial(I->m[i],r,w);
     259//   return inI;
     260// }
     261
     262// #ifndef NDEBUG
     263// BOOLEAN initial0(leftv res, leftv args)
     264// {
     265//   leftv u = args;
     266//   ideal I = (ideal) u->CopyD();
     267//   leftv v = u->next;
     268//   bigintmat* w0 = (bigintmat*) v->Data();
     269//   gfan::ZVector* w = bigintmatToZVector(w0);
     270//   omUpdateInfo();
     271//   Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
     272//   ideal inI = initial(I,currRing,*w);
     273//   id_Delete(&I,currRing);
     274//   delete w;
     275//   res->rtyp = IDEAL_CMD;
     276//   res->data = (char*) inI;
     277//   return FALSE;
     278// }
     279// #endif
     280
     281// poly initial(const poly p, const ring r, const gfan::ZMatrix W)
     282// {
     283//   poly q0 = p_Head(p,r);
     284//   poly q1 = q0;
     285//   gfan::ZVector d = WDeg(p,r,W);
     286//   for (poly currentTerm = p->next; currentTerm; pIter(currentTerm))
     287//   {
     288//     gfan::ZVector e = WDeg(currentTerm,r,W);
     289//     if (d<e)
     290//     {
     291//       p_Delete(&q0,r);
     292//       q0 = p_Head(p,r);
     293//       q1 = q0;
     294//       d = e;
     295//     }
     296//     else
     297//       if (d==e)
     298//       {
     299//         pNext(q1) = p_Head(currentTerm,r);
     300//         pIter(q1);
     301//       }
     302//   }
     303//   return q0;
     304// }
     305
     306// ideal initial(const ideal I, const ring r, const gfan::ZMatrix W)
     307// {
     308//   int k = idSize(I); ideal inI = idInit(k);
     309//   for (int i=0; i<k; i++)
     310//     inI->m[i] = initial(I->m[i],r,W);
     311//   return inI;
     312// }
     313
     314// /***
     315//  * Computes the initial form of p with respect to the first row in the order matrix
     316//  **/
     317// poly initial(const poly p, const ring r)
     318// {
     319//   long d = p_Deg(p,r);
     320//   poly initialForm0 = p_Head(p,r);
     321//   poly initialForm1 = initialForm0;
     322//   poly currentTerm = p->next;
     323//   while (currentTerm && p_Deg(currentTerm,r)==d)
     324//   {
     325//     pNext(initialForm1) = p_Head(currentTerm,r);
     326//     pIter(currentTerm); pIter(initialForm1);
     327//   }
     328//   return initialForm0;
     329// }
     330
     331// /***
     332//  * Computes the initial form of all generators of I.
     333//  * If I is a standard basis, then this is a standard basis
     334//  * of the initial ideal.
     335//  **/
     336// ideal initial(const ideal I, const ring r)
     337// {
     338//   int k = idSize(I); ideal inI = idInit(k);
     339//   for (int i=0; i<k; i++)
     340//     inI->m[i] = initial(I->m[i],r);
     341//   return inI;
     342// }
     343
     344// BOOLEAN initial(leftv res, leftv args)
     345// {
     346//   leftv u = args;
     347//   if ((u != NULL) && (u->Typ() == POLY_CMD) && (u->next == NULL))
     348//   {
     349//     poly p = (poly) u->Data();
     350//     res->rtyp = POLY_CMD;
     351//     res->data = (void*) initial(p, currRing);
     352//     return FALSE;
     353//   }
     354//   if ((u != NULL) && (u->Typ() == IDEAL_CMD) && (u->next == NULL))
     355//   {
     356//     ideal I = (ideal) u->Data();
     357//     res->rtyp = IDEAL_CMD;
     358//     res->data = (void*) initial(I, currRing);
     359//     return FALSE;
     360//   }
     361//   WerrorS("initial: unexpected parameters");
     362//   return TRUE;
     363// }

  • Singular/dyn_modules/gfanlib/initial.h

    rdffd154 r4664b3  
    22#define INITIAL_H
    33
    4 /***
     4/**
    55 * various functions to compute the initial form of polynomials and ideals
    6  **/
     6 */
     7#include <gfanlib/gfanlib_vector.h>
     8#include <gfanlib/gfanlib_matrix.h>
     9#include <libpolys/polys/monomials/p_polys.h>
    710
    8 #include <gfanlib/gfanlib_vector.h>
    9 #include <libpolys/polys/monomials/p_polys.h>
    10 #include <Singular/ipid.h>
    11 
    12 /***
    13  * Computes the weighted degree of the leading monomial of p with respect to w
    14  **/
     11/**
     12 * Returns the weighted degree of the leading term of p with respect to w
     13 */
    1514long wDeg(const poly p, const ring r, const gfan::ZVector w);
    1615
    17 /***
    18  * Computes the weighted multidegree of the leading term of p with respect to W.
    19  * The weighted multidegree is a vector whose i-th entry is the weighted degree
    20  * with respect to the i-th row vector of W.
    21  **/
    22 gfan::ZVector WDeg(const poly p, const ring r, const gfan::ZMatrix W);
     16/**
     17 * Returns the weighted multidegree of the leading term of p with respect to (w,W).
     18 * The weighted multidegree is a vector of length one higher than the column vectors of W.
     19 * The first entry is the weighted degree with respect to w
     20 * and the i+1st entry is the weighted degree with respect to the i-th row vector of W.
     21 */
     22gfan::ZVector WDeg(const poly p, const ring r, const gfan::ZVector w, const gfan::ZMatrix W);
    2323
    24 /***
    25  * Returns the first terms of p of same weighted degree under w.
    26  * Coincides with the initial form of p with respect to w if and only if p was already
    27  * sorted with respect to w.
    28  **/
    29 poly sloppyInitial(const poly p, const ring r, const gfan::ZVector w);
     24/**
     25 * Returns the initial form of p with respect to w
     26 */
     27poly initial(const poly p, const ring r, const gfan::ZVector w);
    3028
    31 /***
    32  * Runs the above procedure over all generators of an ideal.
    33  * Coincides with the initial ideal of I with respect to w if and only if
    34  * the elements of I were already sorted with respect to w and
    35  * I is a standard basis form with respect to w.
    36  **/
    37 ideal sloppyInitial(const ideal I, const ring r, const gfan::ZVector w);
     29/**
     30 * Returns the initial form of I with respect to w
     31 */
     32ideal initial(const ideal I, const ring r, const gfan::ZVector w);
    3833
    39 poly initial(const poly p, const ring r, const gfan::ZVector w);
    40 ideal initial(const ideal I, const ring r, const gfan::ZVector w);
    41 poly initial(const poly p, const ring r, const gfan::ZMatrix W);
    42 ideal initial(const ideal I, const ring r, const gfan::ZMatrix W);
     34/**
     35 * Returns the initial form of p with respect to w and W
     36 */
    4337poly initial(const poly p, const ring r, const gfan::ZVector w, const gfan::ZMatrix W);
     38
     39/**
     40 * Returns the initial form of I with respect to w and W
     41 */
    4442ideal initial(const ideal I, const ring r, const gfan::ZVector w, const gfan::ZMatrix W);
    4543
     44/**
     45 * Stores the initial form of *pStar with respect to w in pStar
     46 */
     47void initial(poly* pStar, const ring r, const gfan::ZVector w);
    4648
    47 poly initial(const poly p, const ring r);
    48 ideal initial(const ideal I, const ring r);
    49 BOOLEAN initial(leftv res, leftv args);
    50 #ifndef NDEBUG
    51 BOOLEAN initial0(leftv res, leftv args);
    52 #endif
     49/**
     50 * Stores the initial form of *IStar with respect to w in IStar
     51 */
     52void initial(ideal* IStar, const ring r, const gfan::ZVector w);
     53
     54/**
     55 * Stores the initial form of *pStar with respect to w and W in pStar
     56 */
     57void initial(poly* pStar, const ring r, const gfan::ZVector w, const gfan::ZMatrix W);
     58
     59/**
     60 * Stores the initial form of *IStar with respect to w and W in IStar
     61 */
     62void initial(ideal* IStar, const ring r, const gfan::ZVector w, const gfan::ZMatrix W);
     63
     64/* throwaway code */
     65/* #ifndef NDEBUG */
     66/* BOOLEAN initial0(leftv res, leftv args); */
     67/* #endif */
     68/* /\*** */
     69/*  * Returns the first terms of p of same weighted degree under w. */
     70/*  * Coincides with the initial form of p with respect to w if and only if p was already */
     71/*  * sorted with respect to w. */
     72/*  **\/ */
     73/* poly sloppyInitial(const poly p, const ring r, const gfan::ZVector w); */
     74
     75/* /\*** */
     76/*  * Runs the above procedure over all generators of an ideal. */
     77/*  * Coincides with the initial ideal of I with respect to w if and only if */
     78/*  * the elements of I were already sorted with respect to w and */
     79/*  * I is a standard basis form with respect to w. */
     80/*  **\/ */
     81/* ideal sloppyInitial(const ideal I, const ring r, const gfan::ZVector w); */
     82/* poly initial(const poly p, const ring r, const gfan::ZMatrix W); */
     83/* ideal initial(const ideal I, const ring r, const gfan::ZMatrix W); */
     84/* poly initial(const poly p, const ring r); */
     85/* ideal initial(const ideal I, const ring r); */
     86/* BOOLEAN initial(leftv res, leftv args); */
    5387
    5488#endif

  • Singular/dyn_modules/gfanlib/ppinitialReduction.cc

    rdffd154 r4664b3  
    88#include <map>
    99#include <set>
     10#include <iostream>
     11#include <exception>
     12
     13#include <ppinitialReduction.h>
    1014
    1115bool isOrderingLocalInT(const ring r)
    … …  
    2933 * with respect to p-t
    3034 **/
    31 bool pReduce(poly &g, const number p, const ring r)
     35void pReduce(poly &g, const number p, const ring r)
    3236{
    3337  if (g==NULL)
    34     return false;
     38    return;
    3539  p_Test(g,r);
    3640
    … …  
    6872          {
    6973            WerrorS("pReduce: overflow in exponent");
    70             return true;
     74            throw 0;
    7175          }
    7276        }
    … …  
    8993  }
    9094  p_Test(g,r);
    91   return false;
     95  return;
    9296}
    9397
    … …  
    168172#endif //NDEBUG
    169173
    170 bool pReduce0(ideal &I, const number p, const ring r)
     174void pReduce0(ideal &I, const number p, const ring r)
    171175{
    172176  int k = idSize(I);
    … …  
    177181      number c = p_GetCoeff(I->m[i],r);
    178182      if (!n_Equal(p,c,r->cf))
    179         if (pReduce(I->m[i],p,r))
    180           return true;
    181     }
    182   }
    183   return false;
     183        pReduce(I->m[i],p,r);
     184    }
     185  }
     186  return;
    184187}
    185188
    … …  
    215218    h = p_Add_q(q1,q2,r);
    216219    p_Test(h,r);
    217     hStar = &h;
     220    p_Test(g,r);
     221    *hStar = h;
    218222    return true;
    219223  }
     224  p_Test(h,r);
     225  p_Test(g,r);
    220226  return false;
    221227}
    … …  
    279285  } while(n);
    280286  for (int i=0; i<m; i++)
    281     if (pReduce(I->m[i],p,r)) return true;
     287    pReduce(I->m[i],p,r);
    282288
    283289  /***
    … …  
    286292  for (int i=0; i<m-1; i++)
    287293    for (int j=i+1; j<m; j++)
    288       if (ppreduceInitially(&I->m[j], I->m[i], r) && pReduce(I->m[j],p,r)) return true;
     294      if (ppreduceInitially(&I->m[j], I->m[i], r))
     295        pReduce(I->m[j],p,r);
    289296
    290297  /***
    … …  
    293300  for (int i=0; i<m-1; i++)
    294301    for (int j=i+1; j<m; j++)
    295       if (ppreduceInitially(&I->m[i], I->m[j],r) && pReduce(I->m[i],p,r)) return true;
     302      if (ppreduceInitially(&I->m[i], I->m[j],r))
     303        pReduce(I->m[i],p,r);
    296304
    297305  /***
    … …  
    338346/***
    339347 * inserts g into I and reduces I with respect to itself and p-t
     348 * returns the position in I in which g was inserted
    340349 * assumes that I was already sorted and initially reduced in the first place
    341350 **/
    342 bool ppreduceInitially(ideal I, const number p, const poly g, const ring r)
    343 {
     351int ppreduceInitially(ideal I, const number p, const poly g, const ring r)
     352{
     353  id_Test(I,r);
     354  p_Test(g,r);
    344355  idInsertPoly(I,g);
    345356  int n=idSize(I);
    … …  
    362373   **/
    363374  for (int i=0; i<j; i++)
    364     if (ppreduceInitially(&I->m[j], I->m[i], r) && pReduce(I->m[j],p,r)) return true;
     375    if (ppreduceInitially(&I->m[j], I->m[i], r))
     376      pReduce(I->m[j],p,r);
    365377  for (int k=j+1; k<n; k++)
    366     if (ppreduceInitially(&I->m[k], I->m[j], r) && pReduce(I->m[k],p,r)) return true;
     378    if (ppreduceInitially(&I->m[k], I->m[j], r))
     379      pReduce(I->m[k],p,r);
    367380
    368381  /***
    … …  
    372385  for (int i=0; i<j; i++)
    373386    for (int k=j; k<n; k++)
    374       if (ppreduceInitially(&I->m[i], I->m[k], r) && pReduce(I->m[i],p,r)) return true;
     387      if (ppreduceInitially(&I->m[i], I->m[k], r))
     388        pReduce(I->m[i],p,r);
    375389  for (int k=j+1; k<n; k++)
    376     if (ppreduceInitially(&I->m[j], I->m[k], r) && pReduce(I->m[j],p,r)) return true;
     390    if (ppreduceInitially(&I->m[j], I->m[k], r))
     391      pReduce(I->m[j],p,r);
    377392
    378393  /***
    … …  
    380395   **/
    381396  idSkipZeroes(I);
    382   return false;
     397  id_Test(I,r);
     398  return j;
    383399}
    384400
    … …  
    423439
    424440
     441static poly ppNext(poly p, int l)
     442{
     443  poly q = p;
     444  for (int i=0; i<l; i++)
     445  {
     446    if (q==NULL)
     447      break;
     448    pIter(q);
     449  }
     450  return q;
     451}
     452
     453
     454static void sortMarks(const ideal H, const ring r, std::vector<mark> &T)
     455{
     456  std::pair<int,int> pointerToTerm;
     457  int k=T.size();
     458  do
     459  {
     460    int j=0;
     461    for (int i=1; i<k-1; i++)
     462    {
     463      int generatorA = T[i-1].first;
     464      int termA = T[i-1].second;
     465      int generatorB = T[i].first;
     466      int termB = T[i].second;
     467      if (p_LmCmp(ppNext(H->m[generatorA],termA),ppNext(H->m[generatorB],termB),r)<0)
     468      {
     469        mark cache=T[i-1];
     470        T[i-1]=T[i];
     471        T[i]=cache;
     472        j = i;
     473      }
     474    }
     475    k=j;
     476  } while(k);
     477  return;
     478}
     479
     480
     481static poly getTerm(const ideal H, const mark ab)
     482{
     483  int a = ab.first;
     484  int b = ab.second;
     485  return ppNext(H->m[a],b);
     486}
     487
     488
     489static void adjustMarks(std::vector<mark> &T, const int newEntry)
     490{
     491  for (unsigned i=0; i<T.size(); i++)
     492  {
     493    if (T[i].first>=newEntry)
     494      T[i].first = T[i].first+1;
     495  }
     496  return;
     497}
     498
     499
     500static void cleanupMarks(const ideal H, std::vector<mark> &T)
     501{
     502  for (unsigned i=0; i<T.size();)
     503  {
     504    if (getTerm(H,T[i])==NULL)
     505      T.erase(T.begin()+i);
     506    else
     507      i++;
     508  }
     509  return;
     510}
     511
     512
    425513/***
    426514 * reduces H initially with respect to itself, with respect to p-t,
    … …  
    437525
    438526  /***
    439    * Step 2: initialize a working list T and an ideal I in which the reductions will take place
     527   * Step 2: initialize an ideal I in which the reductions will take place-
     528   *   along the reduction it will be enlarged with elements that will be discarded at the end
     529   *   initialize a working list T which keeps track.
     530   *   the working list T is a vector of pairs of integer.
     531   *   if T contains a pair (i,j) then that means that in the i-th element of H
     532   *   term j and subsequent terms need to be checked for reduction.
     533   *   T is sorted by the ordering on the temrs the pairs correspond to.
    440534   **/
    441535  int m=idSize(H),n=0;
    442   ideal I = idInit(m), T = idInit(m);
     536  ideal I = idInit(m);
     537  std::vector<mark> T;
    443538  for (int i=0; i<m; i++)
    444539  {
    445540    I->m[i]=H->m[i];
    446     if (pNext(H->m[i])) T->m[n++]=H->m[i];
    447   }
    448   poly g; int k=n;
    449   do
    450   {
    451     int j=0;
    452     for (int i=1; i<k; i++)
    453     {
    454       if (p_LmCmp(pNext(T->m[i-1]),pNext(T->m[i]),r)<0)
    455       {
    456         g=T->m[i-1];
    457         T->m[i-1]=T->m[i];
    458         T->m[i]=g;
    459         j = i;
    460       }
    461     }
    462     k=j;
    463   } while(k);
     541    if (pNext(I->m[i])!=NULL)
     542      T.push_back(std::make_pair<int,int>(i,1));
     543  }
    464544
    465545  /***
    … …  
    467547   *   by adding suitable elements to I and reducing it initially with respect to itself
    468548   **/
    469   k=idSize(G);
    470   while (n)
    471   {
     549  int k=idSize(G);
     550  while (T.size()>0)
     551  {
     552    sortMarks(I,r,T);
     553    std::cout << "T.size()=" << T.size() << std::endl;
     554    std::cout << "T[0] = (" << T[0].first << "," << T[0].second << ")" << std::endl;
    472555    int i=0; for (; i<k; i++)
    473       if (p_LeadmonomDivisibleBy(G->m[i],pNext(T->m[0]),r)) break;
     556      if (p_LeadmonomDivisibleBy(G->m[i],getTerm(I,T[0]),r)) break;
    474557    if (i<k)
    475558    {
    476       g = p_One(r);
     559      if (T[0].first==10 && T[0].second==3)
     560      {
     561        std::cout << "check this" << std::endl;
     562      }
     563      std::cout << "reducing" << std::endl;
     564      poly g = p_One(r); poly h0 = getTerm(I,T[0]);
     565      assume(h0!=NULL);
    477566      for (int j=2; j<=r->N; j++)
    478         p_SetExp(g,j,p_GetExp(pNext(T->m[0]),j,r)-p_GetExp(G->m[i],j,r),r);
     567        p_SetExp(g,j,p_GetExp(h0,j,r)-p_GetExp(G->m[i],j,r),r);
    479568      p_Setm(g,r);
    480569      g = p_Mult_q(g,p_Copy(G->m[i],r),r);
    481       ppreduceInitially(I,p,g,r);
     570      int newEntry = ppreduceInitially(I,p,g,r);
     571      adjustMarks(T,newEntry);
    482572    }
    483573    else
    484       pIter(T->m[0]);
    485     for (int i=0; i<n;)
    486     {
    487       if (!pNext(T->m[i]))
    488       {
    489         for (int j=i; j<n-1; j++)
    490           T->m[j]=T->m[j+1];
    491         T->m[--n]=NULL;
    492       }
    493       else
    494         i++;
    495     }
    496     int l = n;
    497     do
    498     {
    499       int j=0;
    500       for (int i=1; i<l; i++)
    501       {
    502         if (p_LmCmp(pNext(T->m[i-1]),pNext(T->m[i]),r)<0)
    503         {
    504           g=T->m[i-1];
    505           T->m[i-1]=I->m[i];
    506           T->m[i]=g;
    507           j = i;
    508         }
    509       }
    510       l=j;
    511     } while(l);
     574      T[0].second = T[0].second+1;
     575    cleanupMarks(I,T);
    512576  }
    513577
    … …  
    528592  }
    529593  id_Delete(&I,r);
    530   id_Delete(&T,r);
    531594  return false;
    532595}
    … …  
    701764    m += l; IDELEMS(GG) = m; GG->m = &G->m[n-m];
    702765    // std::vector<int> synch = synchronize(I,it->second);
     766    id_Test(it->second,r);
     767    id_Test(GG,r);
     768    if (ppreduceInitially(it->second,p,r)) return true;
    703769    if (ppreduceInitially(it->second,p,GG,r)) return true;
     770    id_Test(it->second,r);
     771    id_Test(GG,r);
    704772    // synchronize(I,it->second,synch);
    705773    idShallowDelete(&Hi); Hi = it->second;

  • Singular/dyn_modules/gfanlib/ppinitialReduction.h

    rdffd154 r4664b3  
    1414#endif
    1515
     16typedef std::pair<int,int> mark;
     17typedef std::vector<std::pair<int,int> > marks;
     18
    1619bool isOrderingLocalInT(const ring r);
    17 bool pReduce0(ideal &I, const number p, const ring r);
     20void pReduce0(ideal &I, const number p, const ring r);
    1821bool ppreduceInitially(ideal I, const ring r, const number p);
    1922BOOLEAN ppreduceInitially(leftv res, leftv args);

  • Singular/dyn_modules/gfanlib/startingCone.cc

    rdffd154 r4664b3  
    367367
    368368      // compute the initial ideal with respect to the weight
    369       inI = sloppyInitial(inI,s,startingPoint);
     369      inI = initial(inI,s,startingPoint);
    370370      zc = linealitySpaceOfGroebnerFan(inI,s);
    371371    }
    … …  
    417417    inJ = ambientMaximalCone.getPolynomialIdeal();
    418418    s = ambientMaximalCone.getPolynomialRing();
    419     inJ = sloppyInitial(inJ,s,startingPoint);
     419    inJ = initial(inJ,s,startingPoint);
    420420    ideal inI = initial(I,r,startingPoint);
    421421    zc = linealitySpaceOfGroebnerFan(inJ,s);
    … …  
    464464    gfan::ZCone zd = startingConeShortcut.getPolyhedralCone();
    465465    gfan::ZVector interiorPoint = startingConeShortcut.getInteriorPoint();
    466     inJShortcut = sloppyInitial(inJShortcut,sShortcut,interiorPoint);
     466    inJShortcut = initial(inJShortcut,sShortcut,interiorPoint);
    467467    inI = initial(inI,r,interiorPoint);
    468468

  • Singular/dyn_modules/gfanlib/tropical.cc

    rdffd154 r4664b3  
    175175  p->iiAddCproc("","groebnerCone",FALSE,groebnerCone);
    176176  p->iiAddCproc("","maximalGroebnerCone",FALSE,maximalGroebnerCone);
    177   p->iiAddCproc("","initial",FALSE,initial);
     177  // p->iiAddCproc("","initial",FALSE,initial);
    178178  // p->iiAddCproc("","tropicalNeighbours",FALSE,tropicalNeighbours);
    179179#ifndef NDEBUG

  • Singular/dyn_modules/gfanlib/tropicalStrategy.cc

    rdffd154 r4664b3  
    325325}
    326326
    327 bool tropicalStrategy::pReduce(ideal I, const ring r) const
     327void tropicalStrategy::pReduce(ideal I, const ring r) const
    328328{
    329329  rTest(r);
    … …  
    331331
    332332  if (isValuationTrivial())
    333     return false;
     333    return;
    334334
    335335  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
    336336  number p = identity(uniformizingParameter,startingRing->cf,r->cf);
    337   bool b = pReduce0(I,p,r);
     337  pReduce0(I,p,r);
    338338  n_Delete(&p,r->cf);
    339339
    340   return b;
     340  return;
    341341}
    342342

  • Singular/dyn_modules/gfanlib/tropicalStrategy.h

    rdffd154 r4664b3  
    283283  bool reduce(ideal I, const ring r) const;
    284284
    285   bool pReduce(ideal I, const ring r) const;
     285  void pReduce(ideal I, const ring r) const;
    286286
    287287  /**
Note: See TracChangeset for help on using the changeset viewer.