Changeset 3cedf8 in git


Ignore:
Timestamp:
Jun 8, 2022, 12:15:43 PM (23 months ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
1ead6f7a96c0b74c88cce4a314fd1d652c7fb0af
Parents:
d4fe54954b5a73a294f4cada4b4d571d8a1d9c4f
Message:
new Hilbert Series algorithm
Location:
kernel/combinatorics
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • kernel/combinatorics/hdegree.cc

    rd4fe54 r3cedf8  
    939939}
    940940
    941 static void hDegree0(ideal S, ideal Q, const ring tailRing)
    942 {
    943   id_TestTail(S, currRing, tailRing);
    944   if (Q!=NULL) id_TestTail(Q, currRing, tailRing);
     941static void hDegree0(ideal S, ideal Q)
     942{
     943  id_LmTest(S, currRing);
     944  if (Q!=NULL) id_LmTest(Q, currRing);
    945945
    946946  int  mc;
    947   hexist = hInit(S, Q, &hNexist, tailRing);
     947  hexist = hInit(S, Q, &hNexist, currRing);
    948948  if (!hNexist)
    949949  {
    … …  
    10151015int  scMult0Int(ideal S, ideal Q, const ring tailRing)
    10161016{
    1017   id_TestTail(S, currRing, tailRing);
    1018   if (Q!=NULL) id_TestTail(Q, currRing, tailRing);
    1019 
    1020   hDegree0(S, Q, tailRing);
     1017  id_LmTest(S, currRing);
     1018  if (Q!=NULL) id_LmTest(Q, currRing);
     1019
     1020  hDegree0(S, Q);
    10211021  return hMu;
    10221022}
    … …  
    11011101void scComputeHC(ideal S, ideal Q, int ak, poly &hEdge, ring tailRing)
    11021102{
    1103   id_TestTail(S, currRing, tailRing);
    1104   if (Q!=NULL) id_TestTail(Q, currRing, tailRing);
     1103  id_LmTest(S, currRing);
     1104  if (Q!=NULL) id_LmTest(Q, currRing);
    11051105
    11061106  int  i;
    … …  
    11101110  {
    11111111    //consider just monic generators (over rings with zero-divisors)
    1112     ideal SS=id_Copy(S,tailRing);
     1112    ideal SS=id_Head(S,currRing);
    11131113    for(i=0;i<=idElem(S);i++)
    11141114    {
    11151115      if((SS->m[i]!=NULL)
    1116       && ((p_IsPurePower(SS->m[i],tailRing)==0)
    1117         ||(!n_IsUnit(pGetCoeff(SS->m[i]), tailRing->cf))))
    1118       {
    1119         p_Delete(&SS->m[i],tailRing);
    1120       }
    1121     }
    1122     S=id_Copy(SS,tailRing);
     1116      && ((p_IsPurePower(SS->m[i],currRing)==0)
     1117        ||(!n_IsUnit(pGetCoeff(SS->m[i]), currRing->cf))))
     1118      {
     1119        p_Delete(&SS->m[i],currRing);
     1120      }
     1121    }
     1122    S=id_Copy(SS,currRing);
    11231123    idSkipZeroes(S);
    11241124  }
    … …  
    15171517void scComputeHCw(ideal ss, ideal Q, int ak, poly &hEdge, ring tailRing)
    15181518{
    1519   id_TestTail(ss, currRing, tailRing);
    1520   if (Q!=NULL) id_TestTail(Q, currRing, tailRing);
     1519  id_LmTest(ss, currRing);
     1520  if (Q!=NULL) id_LmTest(Q, currRing);
    15211521
    15221522  int  i, di;

  • kernel/combinatorics/hilb.cc

    rd4fe54 r3cedf8  
    1818#include "polys/monomials/p_polys.h"
    1919#include "polys/simpleideals.h"
     20#include "polys/weight.h"
    2021
    2122#if SIZEOF_LONG == 8
    … …  
    4950#endif
    5051
    51 STATIC_VAR int64  **Qpol;
    52 STATIC_VAR int64  *Q0, *Ql;
    53 STATIC_VAR int  hLength;
    54 
    55 
    56 static int hMinModulweight(intvec *modulweight)
    57 {
    58   int i,j,k;
    59 
    60   if(modulweight==NULL) return 0;
    61   j=(*modulweight)[0];
    62   for(i=modulweight->rows()-1;i!=0;i--)
    63   {
    64     k=(*modulweight)[i];
    65     if(k<j) j=k;
    66   }
    67   return j;
    68 }
    69 
    70 static void hHilbEst(scfmon stc, int Nstc, varset var, int Nvar)
    71 {
    72   int  i, j;
    73   int  x, y, z = 1;
    74   int64  *p;
    75   for (i = Nvar; i>0; i--)
    76   {
    77     x = 0;
    78     for (j = 0; j < Nstc; j++)
    79     {
    80       y = stc[j][var[i]];
    81       if (y > x)
    82         x = y;
    83     }
    84     z += x;
    85     j = i - 1;
    86     if (z > Ql[j])
    87     {
    88       if (z>(MAX_INT_VAL)/2)
    89       {
    90        WerrorS("internal arrays too big");
    91        return;
    92       }
    93       p = (int64 *)omAlloc((unsigned long)z * sizeof(int64));
    94       if (Ql[j]!=0)
    95       {
    96         if (j==0)
    97           memcpy(p, Qpol[j], Ql[j] * sizeof(int64));
    98         omFreeSize((ADDRESS)Qpol[j], Ql[j] * sizeof(int64));
    99       }
    100       if (j==0)
    101       {
    102         for (x = Ql[j]; x < z; x++)
    103           p[x] = 0;
    104       }
    105       Ql[j] = z;
    106       Qpol[j] = p;
    107     }
    108   }
    109 }
    110 
    111 static int64 *hAddHilb(int Nv, int x, int64 *pol, int *lp)
    112 {
    113   int  l = *lp, ln, i;
    114   int64  *pon;
    115   *lp = ln = l + x;
    116   pon = Qpol[Nv];
    117   memcpy(pon, pol, l * sizeof(int64));
    118   if (l > x)
    119   {/*pon[i] -= pol[i - x];*/
    120     for (i = x; i < l; i++)
    121     {
    122       #ifndef __SIZEOF_INT128__
    123       int64 t=pon[i];
    124       int64 t2=pol[i - x];
    125       t-=t2;
    126       if ((t>=OVERFLOW_MIN)&&(t<=OVERFLOW_MAX)) pon[i]=t;
    127       else if (!errorreported) WerrorS("int overflow in hilb 1");
    128       #else
    129       __int128 t=pon[i];
    130       __int128 t2=pol[i - x];
    131       t-=t2;
    132       if ((t>=LONG_MIN)&&(t<=LONG_MAX)) pon[i]=t;
    133       else if (!errorreported) WerrorS("long int overflow in hilb 1");
    134       #endif
    135     }
    136     for (i = l; i < ln; i++)
    137     { /*pon[i] = -pol[i - x];*/
    138       #ifndef __SIZEOF_INT128__
    139       int64 t= -pol[i - x];
    140       if ((t>=OVERFLOW_MIN)&&(t<=OVERFLOW_MAX)) pon[i]=t;
    141       else if (!errorreported) WerrorS("int overflow in hilb 2");
    142       #else
    143       __int128 t= -pol[i - x];
    144       if ((t>=LONG_MIN)&&(t<=LONG_MAX)) pon[i]=t;
    145       else if (!errorreported) WerrorS("long int overflow in hilb 2");
    146       #endif
    147     }
    148   }
    149   else
    150   {
    151     for (i = l; i < x; i++)
    152       pon[i] = 0;
    153     for (i = x; i < ln; i++)
    154       pon[i] = -pol[i - x];
    155   }
    156   return pon;
    157 }
    158 
    159 static void hLastHilb(scmon pure, int Nv, varset var, int64 *pol, int lp)
    160 {
    161   int  l = lp, x, i, j;
    162   int64  *pl;
    163   int64  *p;
    164   p = pol;
    165   for (i = Nv; i>0; i--)
    166   {
    167     x = pure[var[i + 1]];
    168     if (x!=0)
    169       p = hAddHilb(i, x, p, &l);
    170   }
    171   pl = *Qpol;
    172   j = Q0[Nv + 1];
    173   for (i = 0; i < l; i++)
    174   { /* pl[i + j] += p[i];*/
    175     #ifndef __SIZEOF_INT128__
    176     int64 t=pl[i+j];
    177     int64 t2=p[i];
    178     t+=t2;
    179     if ((t>=OVERFLOW_MIN)&&(t<=OVERFLOW_MAX)) pl[i+j]=t;
    180     else if (!errorreported) WerrorS("int overflow in hilb 3");
    181     #else
    182     __int128 t=pl[i+j];
    183     __int128 t2=p[i];
    184     t+=t2;
    185     if ((t>=LONG_MIN)&&(t<=LONG_MAX)) pl[i+j]=t;
    186     else if (!errorreported) WerrorS("long int overflow in hilb 3");
    187     #endif
    188   }
    189   x = pure[var[1]];
    190   if (x!=0)
    191   {
    192     j += x;
    193     for (i = 0; i < l; i++)
    194     { /* pl[i + j] -= p[i];*/
    195       #ifndef __SIZEOF_INT128__
    196       int64 t=pl[i+j];
    197       int64 t2=p[i];
    198       t-=t2;
    199       if ((t>=OVERFLOW_MIN)&&(t<=OVERFLOW_MAX)) pl[i+j]=t;
    200       else if (!errorreported) WerrorS("int overflow in hilb 4");
    201       #else
    202       __int128 t=pl[i+j];
    203       __int128 t2=p[i];
    204       t-=t2;
    205       if ((t>=LONG_MIN)&&(t<=LONG_MAX)) pl[i+j]=t;
    206       else if (!errorreported) WerrorS("long int overflow in hilb 4");
    207       #endif
    208     }
    209   }
    210   j += l;
    211   if (j > hLength)
    212     hLength = j;
    213 }
    214 
    215 static void hHilbStep(scmon pure, scfmon stc, int Nstc, varset var,
    216  int Nvar, int64 *pol, int Lpol)
    217 {
    218   int  iv = Nvar -1, ln, a, a0, a1, b, i;
    219   int  x, x0;
    220   scmon pn;
    221   scfmon sn;
    222   int64  *pon;
    223   if (Nstc==0)
    224   {
    225     hLastHilb(pure, iv, var, pol, Lpol);
    226     return;
    227   }
    228   x = a = 0;
    229   pn = hGetpure(pure);
    230   sn = hGetmem(Nstc, stc, stcmem[iv]);
    231   hStepS(sn, Nstc, var, Nvar, &a, &x);
    232   Q0[iv] = Q0[Nvar];
    233   ln = Lpol;
    234   pon = pol;
    235   if (a == Nstc)
    236   {
    237     x = pure[var[Nvar]];
    238     if (x!=0)
    239       pon = hAddHilb(iv, x, pon, &ln);
    240     hHilbStep(pn, sn, a, var, iv, pon, ln);
    241     return;
    242   }
    243   else
    244   {
    245     pon = hAddHilb(iv, x, pon, &ln);
    246     hHilbStep(pn, sn, a, var, iv, pon, ln);
    247   }
    248   b = a;
    249   x0 = 0;
    250   loop
    251   {
    252     Q0[iv] += (x - x0);
    253     a0 = a;
    254     x0 = x;
    255     hStepS(sn, Nstc, var, Nvar, &a, &x);
    256     hElimS(sn, &b, a0, a, var, iv);
    257     a1 = a;
    258     hPure(sn, a0, &a1, var, iv, pn, &i);
    259     hLex2S(sn, b, a0, a1, var, iv, hwork);
    260     b += (a1 - a0);
    261     ln = Lpol;
    262     if (a < Nstc)
    263     {
    264       pon = hAddHilb(iv, x - x0, pol, &ln);
    265       hHilbStep(pn, sn, b, var, iv, pon, ln);
    266     }
    267     else
    268     {
    269       x = pure[var[Nvar]];
    270       if (x!=0)
    271         pon = hAddHilb(iv, x - x0, pol, &ln);
    272       else
    273         pon = pol;
    274       hHilbStep(pn, sn, b, var, iv, pon, ln);
    275       return;
    276     }
    277   }
    278 }
    279 
    28052/*
    28153*basic routines
    28254*/
    283 static void hWDegree(intvec *wdegree)
    284 {
    285   int i, k;
    286   int x;
    287 
    288   for (i=(currRing->N); i; i--)
    289   {
    290     x = (*wdegree)[i-1];
    291     if (x != 1)
    292     {
    293       for (k=hNexist-1; k>=0; k--)
    294       {
    295         hexist[k][i] *= x;
    296       }
    297     }
    298   }
    299 }
    300 // ---------------------------------- ADICHANGES ---------------------------------------------
    301 //!!!!!!!!!!!!!!!!!!!!! Just for Monomial Ideals !!!!!!!!!!!!!!!!!!!!!!!!!!!!
    302 
    303 #if 0 // unused
    304 //Tests if the ideal is sorted by degree
    305 static bool idDegSortTest(ideal I)
    306 {
    307     if((I == NULL)||(idIs0(I)))
    308     {
    309         return(TRUE);
    310     }
    311     for(int i = 0; i<IDELEMS(I)-1; i++)
    312     {
    313         if(p_Totaldegree(I->m[i],currRing)>p_Totaldegree(I->m[i+1],currRing))
    314         {
    315             idPrint(I);
    316             WerrorS("Ideal is not deg sorted!!");
    317             return(FALSE);
    318         }
    319     }
    320     return(TRUE);
    321 }
    322 #endif
    32355
    32456//adds the new polynomial at the coresponding position
    … …  
    546278}
    547279
    548 #if 0 // unused
    549 //choice XL: last entry divided by x (xy10z15 -> y9z14)
    550 static poly ChoosePXL(ideal I)
    551 {
    552     int i,j,dummy=0;
    553     poly m;
    554     for(i = IDELEMS(I)-1; (i>=0) && (dummy == 0); i--)
    555     {
    556         for(j = 1; (j<=currRing->N) && (dummy == 0); j++)
    557         {
    558             if(p_GetExp(I->m[i],j, currRing)>1)
    559             {
    560                 dummy = 1;
    561             }
    562         }
    563     }
    564     m = p_Copy(I->m[i+1],currRing);
    565     for(j = 1; j<=currRing->N; j++)
    566     {
    567         dummy = p_GetExp(m,j,currRing);
    568         if(dummy >= 1)
    569         {
    570             p_SetExp(m, j, dummy-1, currRing);
    571         }
    572     }
    573     if(!p_IsOne(m, currRing))
    574     {
    575         p_Setm(m, currRing);
    576         return(m);
    577     }
    578     m = ChoosePVar(I);
    579     return(m);
    580 }
    581 #endif
    582 
    583 #if 0 // unused
    584 //choice XF: first entry divided by x (xy10z15 -> y9z14)
    585 static poly ChoosePXF(ideal I)
    586 {
    587     int i,j,dummy=0;
    588     poly m;
    589     for(i =0 ; (i<=IDELEMS(I)-1) && (dummy == 0); i++)
    590     {
    591         for(j = 1; (j<=currRing->N) && (dummy == 0); j++)
    592         {
    593             if(p_GetExp(I->m[i],j, currRing)>1)
    594             {
    595                 dummy = 1;
    596             }
    597         }
    598     }
    599     m = p_Copy(I->m[i-1],currRing);
    600     for(j = 1; j<=currRing->N; j++)
    601     {
    602         dummy = p_GetExp(m,j,currRing);
    603         if(dummy >= 1)
    604         {
    605             p_SetExp(m, j, dummy-1, currRing);
    606         }
    607     }
    608     if(!p_IsOne(m, currRing))
    609     {
    610         p_Setm(m, currRing);
    611         return(m);
    612     }
    613     m = ChoosePVar(I);
    614     return(m);
    615 }
    616 #endif
    617 
    618 #if 0 // unused
    619 //choice OL: last entry the first power (xy10z15 -> xy9z15)
    620 static poly ChoosePOL(ideal I)
    621 {
    622     int i,j,dummy;
    623     poly m;
    624     for(i = IDELEMS(I)-1;i>=0;i--)
    625     {
    626         m = p_Copy(I->m[i],currRing);
    627         for(j=1;j<=currRing->N;j++)
    628         {
    629             dummy = p_GetExp(m,j,currRing);
    630             if(dummy > 0)
    631             {
    632                 p_SetExp(m,j,dummy-1,currRing);
    633                 p_Setm(m,currRing);
    634             }
    635         }
    636         if(!p_IsOne(m, currRing))
    637         {
    638             return(m);
    639         }
    640         else
    641         {
    642             p_Delete(&m,currRing);
    643         }
    644     }
    645     m = ChoosePVar(I);
    646     return(m);
    647 }
    648 #endif
    649 
    650 #if 0 // unused
    651 //choice OF: first entry the first power (xy10z15 -> xy9z15)
    652 static poly ChoosePOF(ideal I)
    653 {
    654     int i,j,dummy;
    655     poly m;
    656     for(i = 0 ;i<=IDELEMS(I)-1;i++)
    657     {
    658         m = p_Copy(I->m[i],currRing);
    659         for(j=1;j<=currRing->N;j++)
    660         {
    661             dummy = p_GetExp(m,j,currRing);
    662             if(dummy > 0)
    663             {
    664                 p_SetExp(m,j,dummy-1,currRing);
    665                 p_Setm(m,currRing);
    666             }
    667         }
    668         if(!p_IsOne(m, currRing))
    669         {
    670             return(m);
    671         }
    672         else
    673         {
    674             p_Delete(&m,currRing);
    675         }
    676     }
    677     m = ChoosePVar(I);
    678     return(m);
    679 }
    680 #endif
    681 
    682 #if 0 // unused
    683 //choice VL: last entry the first variable with power (xy10z15 -> y)
    684 static poly ChoosePVL(ideal I)
    685 {
    686     int i,j,dummy;
    687     bool flag = TRUE;
    688     poly m = p_ISet(1,currRing);
    689     for(i = IDELEMS(I)-1;(i>=0) && (flag);i--)
    690     {
    691         flag = TRUE;
    692         for(j=1;(j<=currRing->N) && (flag);j++)
    693         {
    694             dummy = p_GetExp(I->m[i],j,currRing);
    695             if(dummy >= 2)
    696             {
    697                 p_SetExp(m,j,1,currRing);
    698                 p_Setm(m,currRing);
    699                 flag = FALSE;
    700             }
    701         }
    702         if(!p_IsOne(m, currRing))
    703         {
    704             return(m);
    705         }
    706     }
    707     m = ChoosePVar(I);
    708     return(m);
    709 }
    710 #endif
    711 
    712 #if 0 // unused
    713 //choice VF: first entry the first variable with power (xy10z15 -> y)
    714 static poly ChoosePVF(ideal I)
    715 {
    716     int i,j,dummy;
    717     bool flag = TRUE;
    718     poly m = p_ISet(1,currRing);
    719     for(i = 0;(i<=IDELEMS(I)-1) && (flag);i++)
    720     {
    721         flag = TRUE;
    722         for(j=1;(j<=currRing->N) && (flag);j++)
    723         {
    724             dummy = p_GetExp(I->m[i],j,currRing);
    725             if(dummy >= 2)
    726             {
    727                 p_SetExp(m,j,1,currRing);
    728                 p_Setm(m,currRing);
    729                 flag = FALSE;
    730             }
    731         }
    732         if(!p_IsOne(m, currRing))
    733         {
    734             return(m);
    735         }
    736     }
    737     m = ChoosePVar(I);
    738     return(m);
    739 }
    740 #endif
    741 
    742280//choice JL: last entry just variable with power (xy10z15 -> y10)
    743281static poly ChoosePJL(ideal I)
    … …  
    769307}
    770308
    771 #if 0
    772 //choice JF: last entry just variable with power -1 (xy10z15 -> y9)
    773 static poly ChoosePJF(ideal I)
    774 {
    775     int i,j,dummy;
    776     bool flag = TRUE;
    777     poly m = p_ISet(1,currRing);
    778     for(i = 0;(i<=IDELEMS(I)-1) && (flag);i++)
    779     {
    780         flag = TRUE;
    781         for(j=1;(j<=currRing->N) && (flag);j++)
    782         {
    783             dummy = p_GetExp(I->m[i],j,currRing);
    784             if(dummy >= 2)
    785             {
    786                 p_SetExp(m,j,dummy-1,currRing);
    787                 p_Setm(m,currRing);
    788                 flag = FALSE;
    789             }
    790         }
    791         if(!p_IsOne(m, currRing))
    792         {
    793             return(m);
    794         }
    795     }
    796     m = ChoosePVar(I);
    797     return(m);
    798 }
    799 #endif
    800 
    801309//chooses 1 \neq p \not\in S. This choice should be made optimal
    802310static poly ChooseP(ideal I)
    803311{
    804312  poly m;
    805   //  TEST TO SEE WHICH ONE IS BETTER
    806   //m = ChoosePXL(I);
    807   //m = ChoosePXF(I);
    808   //m = ChoosePOL(I);
    809   //m = ChoosePOF(I);
    810   //m = ChoosePVL(I);
    811   //m = ChoosePVF(I);
    812313  m = ChoosePJL(I);
    813   //m = ChoosePJF(I);
    814314  return(m);
    815315}
    … …  
    844344static bool JustVar(ideal I)
    845345{
    846     #if 0
    847     int i,j;
    848     bool foundone;
    849     for(i=0;i<=IDELEMS(I)-1;i++)
    850     {
    851         foundone = FALSE;
    852         for(j = 1;j<=currRing->N;j++)
    853         {
    854             if(p_GetExp(I->m[i], j, currRing)>0)
    855             {
    856                 if(foundone == TRUE)
    857                 {
    858                     return(FALSE);
    859                 }
    860                 foundone = TRUE;
    861             }
    862         }
     346    if(p_Totaldegree(I->m[IDELEMS(I)-1],currRing)>1)
     347    {
     348        return(FALSE);
    863349    }
    864350    return(TRUE);
    865     #else
    866     if(p_Totaldegree(I->m[IDELEMS(I)-1],currRing)>1)
    867     {
    868         return(FALSE);
    869     }
    870     return(TRUE);
    871     #endif
    872351}
    873352
    … …  
    1010489
    1011490//the Roune Slice Algorithm
    1012 void rouneslice(ideal I, ideal S, poly q, poly x, int &prune, int &moreprune, int &steps, int &NNN, mpz_ptr &hilbertcoef, int* &hilbpower)
     491static void rouneslice(ideal I, ideal S, poly q, poly x, int &prune, int &moreprune, int &steps, int &NNN, mpz_ptr &hilbertcoef, int* &hilbpower)
    1013492{
    1014493  loop
    … …  
    1207686}
    1208687
    1209 // -------------------------------- END OF CHANGES -------------------------------------------
    1210 static intvec * hSeries(ideal S, intvec *modulweight,
    1211                 int /*notstc*/, intvec *wdegree, ideal Q, ring tailRing)
    1212 {
    1213 //  id_TestTail(S, currRing, tailRing);
    1214 
    1215   intvec *work, *hseries1=NULL;
    1216   int  mc;
    1217   int64  p0;
    1218   int  i, j, k, l, ii, mw;
    1219   hexist = hInit(S, Q, &hNexist, tailRing);
    1220   if (hNexist==0)
    1221   {
    1222     hseries1=new intvec(2);
    1223     (*hseries1)[0]=1;
    1224     (*hseries1)[1]=0;
    1225     return hseries1;
    1226   }
    1227 
    1228   #if 0
    1229   if (wdegree == NULL)
    1230     hWeight();
    1231   else
    1232     hWDegree(wdegree);
    1233   #else
    1234   if (wdegree != NULL) hWDegree(wdegree);
    1235   #endif
    1236 
    1237   p0 = 1;
    1238   hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
    1239   hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
    1240   hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
    1241   stcmem = hCreate((currRing->N) - 1);
    1242   Qpol = (int64 **)omAlloc(((currRing->N) + 1) * sizeof(int64 *));
    1243   Ql = (int64 *)omAlloc0(((currRing->N) + 1) * sizeof(int64));
    1244   Q0 = (int64 *)omAlloc(((currRing->N) + 1) * sizeof(int64));
    1245   *Qpol = NULL;
    1246   hLength = k = j = 0;
    1247   mc = hisModule;
    1248   if (mc!=0)
    1249   {
    1250     mw = hMinModulweight(modulweight);
    1251     hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
    1252   }
    1253   else
    1254   {
    1255     mw = 0;
    1256     hstc = hexist;
    1257     hNstc = hNexist;
    1258   }
    1259   loop
    1260   {
    1261     if (mc!=0)
    1262     {
    1263       hComp(hexist, hNexist, mc, hstc, &hNstc);
    1264       if (modulweight != NULL)
    1265         j = (*modulweight)[mc-1]-mw;
    1266     }
    1267     if (hNstc!=0)
    1268     {
    1269       hNvar = (currRing->N);
    1270       for (i = hNvar; i>=0; i--)
    1271         hvar[i] = i;
    1272       //if (notstc) // TODO: no mon divides another
    1273         hStaircase(hstc, &hNstc, hvar, hNvar);
    1274       hSupp(hstc, hNstc, hvar, &hNvar);
    1275       if (hNvar!=0)
    1276       {
    1277         if ((hNvar > 2) && (hNstc > 10))
    1278           hOrdSupp(hstc, hNstc, hvar, hNvar);
    1279         hHilbEst(hstc, hNstc, hvar, hNvar);
    1280         memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
    1281         hPure(hstc, 0, &hNstc, hvar, hNvar, hpure, &hNpure);
    1282         hLexS(hstc, hNstc, hvar, hNvar);
    1283         Q0[hNvar] = 0;
    1284         hHilbStep(hpure, hstc, hNstc, hvar, hNvar, &p0, 1);
    1285       }
    1286     }
    1287     else
    1288     {
    1289       if(*Qpol!=NULL)
    1290         (**Qpol)++;
    1291       else
    1292       {
    1293         *Qpol = (int64 *)omAlloc(sizeof(int64));
    1294         hLength = *Ql = **Qpol = 1;
    1295       }
    1296     }
    1297     if (*Qpol!=NULL)
    1298     {
    1299       i = hLength;
    1300       while ((i > 0) && ((*Qpol)[i - 1] == 0))
    1301         i--;
    1302       if (i > 0)
    1303       {
    1304         l = i + j;
    1305         if (l > k)
    1306         {
    1307           work = new intvec(l);
    1308           for (ii=0; ii<k; ii++)
    1309             (*work)[ii] = (*hseries1)[ii];
    1310           if (hseries1 != NULL)
    1311             delete hseries1;
    1312           hseries1 = work;
    1313           k = l;
    1314         }
    1315         while (i > 0)
    1316         {
    1317           (*hseries1)[i + j - 1] += (*Qpol)[i - 1];
    1318           (*Qpol)[i - 1] = 0;
    1319           i--;
    1320         }
    1321       }
    1322     }
    1323     mc--;
    1324     if (mc <= 0)
    1325       break;
    1326   }
    1327   if (k==0)
    1328   {
    1329     hseries1=new intvec(2);
    1330     (*hseries1)[0]=0;
    1331     (*hseries1)[1]=0;
    1332   }
    1333   else
    1334   {
    1335     l = k+1;
    1336     while ((*hseries1)[l-2]==0) l--;
    1337     if (l!=k)
    1338     {
    1339       work = new intvec(l);
    1340       for (ii=l-2; ii>=0; ii--)
    1341         (*work)[ii] = (*hseries1)[ii];
    1342       delete hseries1;
    1343       hseries1 = work;
    1344     }
    1345     (*hseries1)[l-1] = mw;
    1346   }
    1347   for (i = 0; i <= (currRing->N); i++)
    1348   {
    1349     if (Ql[i]!=0)
    1350       omFreeSize((ADDRESS)Qpol[i], Ql[i] * sizeof(int64));
    1351   }
    1352   omFreeSize((ADDRESS)Q0, ((currRing->N) + 1) * sizeof(int64));
    1353   omFreeSize((ADDRESS)Ql, ((currRing->N) + 1) * sizeof(int64));
    1354   omFreeSize((ADDRESS)Qpol, ((currRing->N) + 1) * sizeof(int64 *));
    1355   hKill(stcmem, (currRing->N) - 1);
    1356   omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
    1357   omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
    1358   omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
    1359   hDelete(hexist, hNexist);
    1360   if (hisModule!=0)
    1361     omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
    1362   return hseries1;
    1363 }
    1364 
    1365 
    1366 intvec * hHstdSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
    1367 {
    1368   id_TestTail(S, currRing, tailRing);
    1369   if (Q!=NULL) id_TestTail(Q, currRing, tailRing);
    1370   return hSeries(S, modulweight, 0, wdegree, Q, tailRing);
    1371 }
    1372 
    1373 intvec * hFirstSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
    1374 {
    1375   id_TestTail(S, currRing, tailRing);
    1376   if (Q!= NULL) id_TestTail(Q, currRing, tailRing);
    1377 
    1378   intvec *hseries1= hSeries(S, modulweight, 1, wdegree, Q, tailRing);
    1379   if (errorreported) { delete hseries1; hseries1=NULL; }
    1380   return hseries1;
    1381 }
    1382 
    1383688intvec * hSecondSeries(intvec *hseries1)
    1384689{
    … …  
    1460765*caller
    1461766*/
    1462 void hLookSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
    1463 {
    1464   id_TestTail(S, currRing, tailRing);
    1465 
    1466   intvec *hseries1 = hFirstSeries(S, modulweight, Q, wdegree, tailRing);
     767void hLookSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring src)
     768{
     769  id_LmTest(S, currRing);
     770
     771  intvec *hseries1 = hFirstSeries(S, modulweight, Q, wdegree, src);
    1467772  if (errorreported) return;
    1468773
    … …  
    18021107}
    18031108
    1804 void sortMonoIdeal_pCompare(ideal I)
     1109static void sortMonoIdeal_pCompare(ideal I)
    18051110{
    18061111  /*
    … …  
    20121317void HilbertSeries_OrbitData(ideal S, int lV, bool IG_CASE, bool mgrad, bool odp, int trunDegHs)
    20131318{
    2014        
    2015         /* new story:
    2016         no lV is needed, i.e.  it is to be determined
    2017         the rest is extracted from the interface input list in extra.cc and makes the input of this proc
    2018         called from extra.cc
    2019         */
    2020        
     1319
     1320        /* new story:
     1321        no lV is needed, i.e.  it is to be determined
     1322        the rest is extracted from the interface input list in extra.cc and makes the input of this proc
     1323        called from extra.cc
     1324        */
     1325
    20211326  /*
    20221327   * This is based on iterative right colon operations on a
    … …  
    23701675}
    23711676
     1677/* ------------------------------------------------------------------------ */
     1678
     1679/****************************************
     1680*  Computer Algebra System SINGULAR     *
     1681****************************************/
     1682/*
     1683*  ABSTRACT -  Hilbert series
     1684*/
     1685
     1686#include "kernel/mod2.h"
     1687
     1688#include "misc/mylimits.h"
     1689#include "misc/intvec.h"
     1690
     1691#include "polys/monomials/ring.h"
     1692#include "polys/monomials/p_polys.h"
     1693#include "polys/simpleideals.h"
     1694#include "coeffs/coeffs.h"
     1695
     1696#include "kernel/ideals.h"
     1697
     1698static void p_Div_hi(poly p, const int* exp_q, const ring src)
     1699{
     1700  // e=max(0,p-q) for all exps
     1701  for(int i=1;i<=src->N;i++)
     1702  {
     1703    p_SetExp(p,i,si_max(0L,p_GetExp(p,i,src)-exp_q[i]),src);
     1704  }
     1705  p_Setm(p,src);
     1706}
     1707
     1708poly hilbert_series(ideal A, const ring src, const intvec* wdegree, const ring Qt)
     1709// accoding to:
     1710// Algorithm 2.6 of
     1711// Dave Bayer, Mike Stillman - Computation of Hilbert Function
     1712// J.Symbolic Computaion (1992) 14, 31-50
     1713{
     1714  int r=id_Elem(A,src);
     1715  poly h=NULL;
     1716  if (r==0)
     1717    return p_One(Qt);
     1718  else
     1719  {
     1720    if (wdegree!=NULL)
     1721    {
     1722      int* exp=(int*)omAlloc((src->N+1)*sizeof(int));
     1723      for(int i=IDELEMS(A)-1; i>=0;i--)
     1724      {
     1725        if (A->m[i]!=NULL)
     1726        {
     1727          p_GetExpV(A->m[i],exp,src);
     1728          for(int j=src->N;j>0;j--)
     1729            exp[j]*=ABS((*wdegree)[j-1]);
     1730          p_SetExpV(A->m[i],exp,src);
     1731          p_Setm(A->m[i],src);
     1732        }
     1733      }
     1734      omFreeSize(exp,(src->N+1)*sizeof(int));
     1735    }
     1736    //Print("start hilbert_series, r=%d\n",r);
     1737    h=p_One(Qt);
     1738    p_SetExp(h,1,p_Totaldegree(A->m[0],src),Qt);
     1739    p_Setm(h,Qt);
     1740    h=p_Neg(h,Qt);
     1741    h=p_Add_q(h,p_One(Qt),Qt);
     1742  }
     1743  int i;
     1744  int *exp_q=(int*)omAlloc((src->N+1)*sizeof(int));
     1745  for (i=1;i<r;i++)
     1746  {
     1747    ideal J=id_Copy(A,src);
     1748    for (int ii=i;ii<r;ii++) p_Delete(&J->m[ii],src);
     1749    idSkipZeroes(J);
     1750    for(int ii=1;ii<=src->N;ii++)
     1751      exp_q[ii]=p_GetExp(A->m[i],ii,src);
     1752    for(int ii=0;ii<i;ii++) p_Div_hi(J->m[ii],exp_q,src);
     1753    id_DelDiv(J,src);
     1754    // search linear elems:
     1755    int k=0;
     1756    for (int ii=0;ii<IDELEMS(J);ii++)
     1757    {
     1758      if((J->m[ii]!=NULL) && (p_Totaldegree(J->m[ii],src)==1))
     1759      {
     1760        k++;
     1761        p_Delete(&J->m[ii],src);
     1762      }
     1763    }
     1764    idSkipZeroes(J);
     1765    poly h_J=hilbert_series(J,src,NULL,Qt);// J_1
     1766    poly tmp;
     1767    if (k>0)
     1768    {
     1769      // hilbert_series of unmodified J:
     1770      tmp=p_One(Qt);
     1771      p_SetExp(tmp,1,1,Qt);
     1772      tmp=p_Neg(tmp,Qt);
     1773      tmp=p_Add_q(tmp,p_One(Qt),Qt);
     1774      if (k>1)
     1775      {
     1776        tmp=p_Power(tmp,k,Qt); // (1-t)^k
     1777      }
     1778      h_J=p_Mult_q(h_J,tmp,Qt);
     1779    }
     1780    // forget about J:
     1781    id_Delete(&J,src);
     1782    // t^|A_i|
     1783    tmp=p_One(Qt);
     1784    p_SetExp(tmp,1,p_Totaldegree(A->m[i],src),Qt);
     1785    tmp=p_Neg(tmp,Qt);
     1786    tmp=p_Mult_q(tmp,h_J,Qt);
     1787    h=p_Add_q(h,tmp,Qt);
     1788  }
     1789  omFreeSize(exp_q,(src->N+1)*sizeof(int));
     1790  //Print("end hilbert_series, r=%d\n",r);
     1791  return h;
     1792}
     1793
     1794static ring makeQt()
     1795{
     1796  ring Qt=(ring) omAlloc0Bin(sip_sring_bin);
     1797  Qt->cf = nInitChar(n_Q, NULL);
     1798  Qt->N=1;
     1799  Qt->names=(char**)omAlloc(sizeof(char_ptr));
     1800  Qt->names[0]=omStrDup("t");
     1801  Qt->wvhdl=(int **)omAlloc0(3 * sizeof(int_ptr));
     1802  Qt->order = (rRingOrder_t *) omAlloc(3 * sizeof(rRingOrder_t *));
     1803  Qt->block0 = (int *)omAlloc0(3 * sizeof(int *));
     1804  Qt->block1 = (int *)omAlloc0(3 * sizeof(int *));
     1805  /* ringorder lp for the first block: var 1 */
     1806  Qt->order[0]  = ringorder_lp;
     1807  Qt->block0[0] = 1;
     1808  Qt->block1[0] = 1;
     1809  /* ringorder C for the second block: no vars */
     1810  Qt->order[1]  = ringorder_C;
     1811  /* the last block: everything is 0 */
     1812  Qt->order[2]  = (rRingOrder_t)0;
     1813  rComplete(Qt);
     1814  return Qt;
     1815}
     1816
     1817intvec* hFirstSeries0(ideal A,ideal Q, intvec *wdegree, const ring src, const ring Qt)
     1818{
     1819  A=id_Head(A,src);
     1820  id_Test(A,src);
     1821  ideal AA;
     1822  if (Q!=NULL)
     1823  {
     1824    ideal QQ=id_Head(Q,src);
     1825    AA=id_SimpleAdd(A,QQ,src);
     1826    id_Delete(&QQ,src);
     1827    id_Delete(&A,src);
     1828  }
     1829  else AA=A;
     1830  id_DelDiv(AA,src);
     1831  idSkipZeroes(AA);
     1832  poly s=hilbert_series(AA,src,wdegree,Qt);
     1833  id_Delete(&AA,src);
     1834  intvec *ss;
     1835  if (s==NULL)
     1836    ss=new intvec(2);
     1837  else
     1838  {
     1839    ss=new intvec(p_Totaldegree(s,Qt)+2);
     1840    while(s!=NULL)
     1841    {
     1842      int i=p_Totaldegree(s,Qt);
     1843      (*ss)[i]=n_Int(pGetCoeff(s),Qt->cf);
     1844      p_LmDelete(&s,Qt);
     1845    }
     1846  }
     1847  return ss;
     1848}
     1849
     1850static ideal getModuleComp(ideal A, int c, const ring src)
     1851{
     1852  ideal res=idInit(IDELEMS(A),A->rank);
     1853  for (int i=0;i<IDELEMS(A);i++)
     1854  {
     1855    if ((A->m[i]!=NULL) && (p_GetComp(A->m[i],src)==c))
     1856      res->m[i]=p_Head(A->m[i],src);
     1857  }
     1858  return res;
     1859}
     1860
     1861static BOOLEAN isModule(ideal A, const ring src)
     1862{
     1863  for (int i=0;i<IDELEMS(A);i++)
     1864  {
     1865    if (A->m[i]!=NULL)
     1866    {
     1867      if (p_GetComp(A->m[i],src)>0)
     1868        return TRUE;
     1869      else
     1870        return FALSE;
     1871    }
     1872  }
     1873  return FALSE;
     1874}
     1875
     1876
     1877intvec* hFirstSeries(ideal A,intvec *module_w,ideal Q, intvec *wdegree, ring src)
     1878{
     1879  static ring hilb_Qt=NULL;
     1880  if (hilb_Qt==NULL) hilb_Qt=makeQt();
     1881  if (!isModule(A,src))
     1882    return hFirstSeries0(A,Q,wdegree,src,hilb_Qt);
     1883  intvec *res=NULL;
     1884  int w_max=0,w_min=0;
     1885  if (module_w!=NULL)
     1886  {
     1887    w_max=module_w->max_in();
     1888    w_min=module_w->min_in();
     1889  }
     1890  for(int c=1;c<=A->rank;c++)
     1891  {
     1892    ideal Ac=getModuleComp(A,c,src);
     1893    intvec *res_c=hFirstSeries0(Ac,Q,wdegree,src,hilb_Qt);
     1894    intvec *tmp=NULL;
     1895    if (res==NULL)
     1896      res=new intvec(res_c->length()+(w_max-w_min));
     1897    if ((module_w==NULL) || ((*module_w)[c-1]==0)) tmp=ivAdd(res,res_c);
     1898    else tmp=ivAddShift(res, res_c,(*module_w)[c-1]-w_min);
     1899    delete res_c;
     1900    if (tmp!=NULL)
     1901    {
     1902      delete res;
     1903      res=tmp;
     1904    }
     1905  }
     1906  (*res)[res->length()-1]=w_min;
     1907  return res;
     1908}

  • kernel/combinatorics/hilb.h

    rd4fe54 r3cedf8  
    1212#include "misc/intvec.h"
    1313
    14 intvec * hHstdSeries(ideal S, intvec *modulweight, ideal Q=NULL, intvec *wdegree=NULL, ring tailRing = currRing);
    15 intvec * hFirstSeries(ideal S, intvec *modulweight, ideal Q=NULL, intvec *wdegree=NULL, ring tailRing = currRing);
     14intvec * hFirstSeries(ideal S, intvec *modulweight, ideal Q=NULL, intvec *wdegree=NULL, ring src = currRing);
     15intvec * hFirstSeries0(ideal S, ideal Q, intvec *wdegree, const ring src, const ring Qt);
     16poly hilbert_series(ideal A, const ring src,  const intvec *wdegree, const ring Qt);
    1617
    1718intvec * hSecondSeries(intvec *hseries1);
    1819
    19 void hLookSeries(ideal S, intvec *modulweight, ideal Q=NULL, intvec *wdegree=NULL, ring tailRing = currRing);
    20 
    21 void sortMonoIdeal_pCompare(ideal I);
     20void hLookSeries(ideal S, intvec *modulweight, ideal Q=NULL, intvec *wdegree=NULL, ring src = currRing);
    2221#endif
    23 
    24 

  • kernel/combinatorics/hutil.cc

    rd4fe54 r3cedf8  
    3131scfmon hInit(ideal S, ideal Q, int *Nexist, ring tailRing)
    3232{
    33   id_TestTail(S, currRing, tailRing);
    34   if (Q!=NULL) id_TestTail(Q, currRing, tailRing);
     33  id_LmTest(S, currRing);
     34  if (Q!=NULL) id_LmTest(Q, currRing);
    3535
    3636//   if (tailRing != currRing)
Note: See TracChangeset for help on using the changeset viewer.