Changeset 204092 in git


Ignore:
Timestamp:
Jul 3, 2012, 7:51:13 PM (12 years ago)
Author:
Oleksandr Motsak <motsak@…>
Branches:
(u'spielwiese', '4a9821a93ffdc22a6696668bd4f6b8c9de3e6c5f')
Children:
7088f18ba4c47a3e87c604a59adbacd8ef35e979
Parents:
cd5fefcfcb3762114c087ec176782c753f2a14d8
git-author:
Oleksandr Motsak <motsak@mathematik.uni-kl.de>2012-07-03 19:51:13+02:00
git-committer:
Oleksandr Motsak <motsak@mathematik.uni-kl.de>2014-05-07 04:41:45+02:00
Message:
moved/separated new functions related to Schreyer Syzygy computation

chg: prefixed corresponding wrappers with underscore due to interpreter registratrar' general expectation
Location:
dyn_modules/syzextra
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • dyn_modules/syzextra/mod_main.cc

    rcd5fefc r204092  
    3434
    3535#include "singularxx_defs.h"
     36
    3637#include "DebugPrint.h"
    3738#include "myNF.h"
    … …  
    5152#include <string.h>
    5253
    53 
    54 #ifdef _GNU_SOURCE
    55 #define qsort_my(m, s, ss, r, cmp) qsort_r(m, s, ss, cmp, r)
    56 #else
    57 #define qsort_my(m, s, ss, r, cmp) qsort_r(m, s, ss, cmp)
    58 #endif
    5954
    6055
    … …  
    6459extern ring rAssure_InducedSchreyerOrdering(const ring r, BOOLEAN complete, int sign);
    6560extern int rGetISPos(const int p, const ring r);
     61
     62// USING_NAMESPACE_SINGULARXX;
    6663
    6764USING_NAMESPACE( SINGULARXXNAME :: DEBUG )
    … …  
    509506
    510507
    511 #ifdef _GNU_SOURCE
    512 static int cmp_c_ds(const void *p1, const void *p2, void *R)
    513 {
    514 #else
    515 static int cmp_c_ds(const void *p1, const void *p2)
    516 {
    517   void *R = currRing;
    518 #endif
    519 
    520   const int YES = 1;
    521   const int NO = -1;
    522 
    523   const ring r =  (const ring) R; // TODO/NOTE: the structure is known: C, lp!!!
    524 
    525   assume( r == currRing );
    526 
    527   const poly a = *(const poly*)p1;
    528   const poly b = *(const poly*)p2;
    529 
    530   assume( a != NULL );
    531   assume( b != NULL );
    532  
    533   assume( p_LmTest(a, r) );
    534   assume( p_LmTest(b, r) );
    535 
    536 
    537   const signed long iCompDiff = p_GetComp(a, r) - p_GetComp(b, r);
    538 
    539   // TODO: test this!!!!!!!!!!!!!!!!
    540 
    541   //return -( compare (c, qsorts) )
    542 
    543   const int __DEBUG__ = 0;
    544 
    545 #ifndef NDEBUG
    546   if( __DEBUG__ )
    547   {
    548     PrintS("cmp_c_ds: a, b: \np1: "); dPrint(a, r, r, 2);
    549     PrintS("b: "); dPrint(b, r, r, 2);
    550     PrintLn();
    551   }
    552 #endif
    553 
    554 
    555   if( iCompDiff > 0 )
    556     return YES;
    557 
    558   if( iCompDiff < 0 )
    559     return  NO;
    560 
    561   assume( iCompDiff == 0 );
    562 
    563   const signed long iDegDiff = p_Totaldegree(a, r) - p_Totaldegree(b, r);
    564 
    565   if( iDegDiff > 0 )
    566     return YES;
    567 
    568   if( iDegDiff < 0 )
    569     return  NO;
    570 
    571   assume( iDegDiff == 0 );
    572 
    573 #ifndef NDEBUG
    574   if( __DEBUG__ )
    575   {
    576     PrintS("cmp_c_ds: a & b have the same comp & deg! "); PrintLn();
    577   }
    578 #endif
    579 
    580   for (int v = rVar(r); v > 0; v--)
    581   {
    582     assume( v > 0 );
    583     assume( v <= rVar(r) );
    584 
    585     const signed int d = p_GetExp(a, v, r) - p_GetExp(b, v, r);
    586 
    587     if( d > 0 )
    588       return YES;
    589 
    590     if( d < 0 )
    591       return NO;
    592 
    593     assume( d == 0 );
    594   }
    595 
    596   return 0; 
    597 }
    598 
    599 
    600 
    601 static ideal ComputeLeadingSyzygyTerms(const ideal& id, const ring r)
    602 {
    603   // 1. set of components S?
    604   // 2. for each component c from S: set of indices of leading terms
    605   // with this component?
    606   // 3. short exp. vectors for each leading term?
    607 
    608   const int size = IDELEMS(id);
    609 
    610   if( size < 2 )
    611   {
    612     const ideal newid = idInit(1, 0); newid->m[0] = NULL; // zero ideal...       
    613     return newid;
    614   }
    615 
    616 
    617   // TODO/NOTE: input is supposed to be (reverse-) sorted wrt "(c,ds)"!??
    618 
    619   // components should come in groups: count elements in each group
    620   // && estimate the real size!!!
    621 
    622 
    623   // use just a vector instead???
    624   const ideal newid = idInit( (size * (size-1))/2, size); // maximal size: ideal case!
    625 
    626   int k = 0;
    627 
    628   for (int j = 0; j < size; j++)
    629   {
    630     const poly p = id->m[j];
    631     assume( p != NULL );
    632     const int  c = p_GetComp(p, r);
    633 
    634     for (int i = j - 1; i >= 0; i--)
    635     {
    636       const poly pp = id->m[i];
    637       assume( pp != NULL );
    638       const int  cc = p_GetComp(pp, r);
    639 
    640       if( c != cc )
    641         continue;
    642 
    643       const poly m = p_Init(r); // p_New???
    644 
    645       // m = LCM(p, pp) / p! // TODO: optimize: knowing the ring structure: (C/lp)!
    646       for (int v = rVar(r); v > 0; v--)
    647       {
    648         assume( v > 0 );
    649         assume( v <= rVar(r) );
    650 
    651         const short e1 = p_GetExp(p , v, r);
    652         const short e2 = p_GetExp(pp, v, r);
    653 
    654         if( e1 >= e2 )
    655           p_SetExp(m, v, 0, r);
    656         else
    657           p_SetExp(m, v, e2 - e1, r);
    658 
    659       }
    660 
    661       assume( (j > i) && (i >= 0) );
    662 
    663       p_SetComp(m, j + 1, r);
    664       pNext(m) = NULL;
    665       p_SetCoeff0(m, n_Init(1, r->cf), r); // for later...
    666 
    667       p_Setm(m, r); // should not do anything!!!
    668 
    669       newid->m[k++] = m;
    670     }
    671   }
    672 
    673 //   if( __DEBUG__ && FALSE )
    674 //   {
    675 //     PrintS("ComputeLeadingSyzygyTerms::Temp0: \n");
    676 //     dPrint(newid, r, r, 1);
    677 //   }
    678 
    679   // the rest of newid is assumed to be zeroes...
    680 
    681   // simplify(newid, 2 + 32)??
    682   // sort(newid, "C,ds")[1]???
    683   id_DelDiv(newid, r); // #define SIMPL_LMDIV 32
    684 
    685 //   if( __DEBUG__ && FALSE )
    686 //   {
    687 //     PrintS("ComputeLeadingSyzygyTerms::Temp1: \n");
    688 //     dPrint(newid, r, r, 1);
    689 //   }
    690 
    691   idSkipZeroes(newid); // #define SIMPL_NULL 2
    692 
    693 //   if( __DEBUG__ )
    694 //   {
    695 //     PrintS("ComputeLeadingSyzygyTerms::Output: \n");
    696 //     dPrint(newid, r, r, 1);
    697 //   }
    698 
    699   const int sizeNew = IDELEMS(newid);
    700 
    701   if( sizeNew >= 2 )
    702     qsort_my(newid->m, sizeNew, sizeof(poly), r, cmp_c_ds);
    703 
    704   newid->rank = id_RankFreeModule(newid, r);
    705 
    706   return newid;
    707 }
    708 
    709 
    710 static BOOLEAN ComputeLeadingSyzygyTerms(leftv res, leftv h)
     508static BOOLEAN _ComputeLeadingSyzygyTerms(leftv res, leftv h)
    711509{
    712510  const ring r = currRing;
    … …  
    767565/// change the input inplace!!!
    768566// TODO: use a ring with >_{c, ds}!???
    769 static BOOLEAN Sort_c_ds(leftv res, leftv h)
     567static BOOLEAN _Sort_c_ds(leftv res, leftv h)
    770568{
    771569  NoReturn(res);
    … …  
    817615    id_Test(id, r);
    818616
    819     const int size = IDELEMS(id);
    820 
    821     const ideal newid = id; // id_Copy(id, r); // copy???
    822 
    823     if( size >= 2 )
    824       qsort_my(newid->m, size, sizeof(poly), r, cmp_c_ds);
    825    
    826 //    res->data = newid;
     617    Sort_c_ds(id, r); // NOT A COPY! inplace sorting!!!
     618
     619//    res->data = id;
    827620//    res->rtyp = h->Typ();
    828621   
    … …  
    830623    {
    831624      PrintS("Sort_c_ds::Output: \n");
    832       dPrint(newid, r, r, 1);
     625      dPrint(id, r, r, 1);
    833626    }
    834627
     628    // NOTE: nothing is to be returned!!!
    835629    return FALSE;
    836630  }
    … …  
    841635
    842636
    843 
    844 static ideal Compute2LeadingSyzygyTerms(const ideal& id, const ring r)
    845 {
     637static BOOLEAN _Compute2LeadingSyzygyTerms(leftv res, leftv h)
     638{
     639  const ring r = currRing;
     640  NoReturn(res);
     641
     642  if( h == NULL )
     643  {
     644    WarnS("Compute2LeadingSyzygyTerms needs an argument...");
     645    return TRUE;
     646  }
     647
     648  assume( h != NULL ); 
     649
    846650#ifndef NDEBUG
    847651  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE)));
    … …  
    850654#endif
    851655
    852   const BOOLEAN __TAILREDSYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"TAILREDSYZ",INT_CMD, (void*)0)));
    853   const BOOLEAN __SYZCHECK__ = (BOOLEAN)((long)(atGet(currRingHdl,"SYZCHECK",INT_CMD, (void*)0)));   
    854 
    855 
    856 
    857   // 1. set of components S?
    858   // 2. for each component c from S: set of indices of leading terms
    859   // with this component?
    860   // 3. short exp. vectors for each leading term?
    861 
    862   const int size = IDELEMS(id);
    863 
    864   if( size < 2 )
    865   {
    866     const ideal newid = idInit(1, 1); newid->m[0] = NULL; // zero module...   
    867     return newid;
    868   }
    869 
    870 
    871     // TODO/NOTE: input is supposed to be sorted wrt "C,ds"!??
    872 
    873     // components should come in groups: count elements in each group
    874     // && estimate the real size!!!
    875 
    876 
    877     // use just a vector instead???
    878   ideal newid = idInit( (size * (size-1))/2, size); // maximal size: ideal case!
    879 
    880   int k = 0;
    881 
    882   for (int j = 0; j < size; j++)
    883   {
    884     const poly p = id->m[j];
    885     assume( p != NULL );
    886     const int  c = p_GetComp(p, r);
    887 
    888     for (int i = j - 1; i >= 0; i--)
    889     {
    890       const poly pp = id->m[i];
    891       assume( pp != NULL );
    892       const int  cc = p_GetComp(pp, r);
    893 
    894       if( c != cc )
    895         continue;
    896 
    897         // allocate memory & zero it out!
    898       const poly m = p_Init(r); const poly mm = p_Init(r);
    899 
    900 
    901         // m = LCM(p, pp) / p! mm = LCM(p, pp) / pp!
    902         // TODO: optimize: knowing the ring structure: (C/lp)!
    903 
    904       for (int v = rVar(r); v > 0; v--)
    905       {
    906         assume( v > 0 );
    907         assume( v <= rVar(r) );
    908 
    909         const short e1 = p_GetExp(p , v, r);
    910         const short e2 = p_GetExp(pp, v, r);
    911 
    912         if( e1 >= e2 )
    913           p_SetExp(mm, v, e1 - e2, r); //            p_SetExp(m, v, 0, r);
    914         else
    915           p_SetExp(m, v, e2 - e1, r); //            p_SetExp(mm, v, 0, r);
    916 
    917       }
    918 
    919       assume( (j > i) && (i >= 0) );
    920 
    921       p_SetComp(m, j + 1, r);
    922       p_SetComp(mm, i + 1, r);
    923 
    924       const number& lc1 = p_GetCoeff(p , r);
    925       const number& lc2 = p_GetCoeff(pp, r);
    926 
    927       number g = n_Lcm( lc1, lc2, r );
    928 
    929       p_SetCoeff0(m ,       n_Div(g, lc1, r), r);
    930       p_SetCoeff0(mm, n_Neg(n_Div(g, lc2, r), r), r);
    931 
    932       n_Delete(&g, r);
    933 
    934       p_Setm(m, r); // should not do anything!!!
    935       p_Setm(mm, r); // should not do anything!!!
    936 
    937       pNext(m) = mm; //        pNext(mm) = NULL;
    938 
    939       newid->m[k++] = m;
    940     }
    941   }
    942 
    943 //   if( __DEBUG__ && FALSE )
    944 //   {
    945 //     PrintS("Compute2LeadingSyzygyTerms::Temp0: \n");
    946 //     dPrint(newid, r, r, 1);
    947 //   }
    948 
    949   if( !__TAILREDSYZ__ )
    950   {
    951       // simplify(newid, 2 + 32)??
    952       // sort(newid, "C,ds")[1]???
    953     id_DelDiv(newid, r); // #define SIMPL_LMDIV 32
    954 
    955 //     if( __DEBUG__ && FALSE )
    956 //     {
    957 //       PrintS("Compute2LeadingSyzygyTerms::Temp1 (deldiv): \n");
    958 //       dPrint(newid, r, r, 1);
    959 //     }
    960   }
    961   else
    962   {
    963       //      option(redSB); option(redTail);
    964       //      TEST_OPT_REDSB
    965       //      TEST_OPT_REDTAIL
    966     assume( r == currRing );
    967     BITSET _save_test = test;
    968     test |= (Sy_bit(OPT_REDTAIL) | Sy_bit(OPT_REDSB));
    969 
    970     intvec* w=new intvec(IDELEMS(newid));
    971     ideal tmp = kStd(newid, currQuotient, isHomog, &w);
    972     delete w;
    973 
    974     test = _save_test;
    975 
    976     id_Delete(&newid, r);
    977     newid = tmp;
    978 
    979 //     if( __DEBUG__ && FALSE )
    980 //     {
    981 //       PrintS("Compute2LeadingSyzygyTerms::Temp1 (std): \n");
    982 //       dPrint(newid, r, r, 1);
    983 //     }
    984 
    985   }
    986 
    987   idSkipZeroes(newid);
    988 
    989   const int sizeNew = IDELEMS(newid);
    990 
    991   if( sizeNew >= 2 )
    992     qsort_my(newid->m, sizeNew, sizeof(poly), r, cmp_c_ds);
    993 
    994   newid->rank = id_RankFreeModule(newid, r);
    995  
    996   return newid;
    997 }
    998 
    999 
    1000 
    1001 static BOOLEAN Compute2LeadingSyzygyTerms(leftv res, leftv h)
    1002 {
    1003   const ring r = currRing;
    1004   NoReturn(res);
    1005 
    1006   if( h == NULL )
    1007   {
    1008     WarnS("Compute2LeadingSyzygyTerms needs an argument...");
    1009     return TRUE;
    1010   }
    1011 
    1012   assume( h != NULL ); 
    1013 
    1014 #ifndef NDEBUG
    1015   const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE)));
    1016 #else
    1017   const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE)));
    1018 #endif
    1019 
    1020656  const BOOLEAN __LEAD2SYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"LEAD2SYZ",INT_CMD, (void*)0)));
    1021657 
    … …  
    1047683
    1048684
    1049 /// return a new term: leading coeff * leading monomial of p
    1050 /// with 0 leading component!
    1051 static inline poly leadmonom(const poly p, const ring r)
    1052 {
    1053   poly m = NULL;
    1054 
    1055   if( p != NULL )
    1056   {
    1057     assume( p != NULL );
    1058     assume( p_LmTest(p, r) );
    1059 
    1060     m = p_LmInit(p, r);
    1061     p_SetCoeff0(m, n_Copy(p_GetCoeff(p, r), r), r);
    1062 
    1063     p_SetComp(m, 0, r);
    1064     p_Setm(m, r);
    1065 
    1066     assume( p_GetComp(m, r) == 0 );
    1067     assume( m != NULL );
    1068     assume( pNext(m) == NULL );
    1069     assume( p_LmTest(m, r) );
    1070   }
    1071 
    1072   return m;
    1073 }
    1074 
    1075 
    1076 
    1077 static poly FindReducer(poly product, poly syzterm,
    1078                         ideal L, ideal LS,
    1079                         const ring r)
    1080 {
    1081   assume( product != NULL );
    1082   assume( L != NULL );
    1083 
    1084   int c = 0;
    1085 
    1086   if (syzterm != NULL)
    1087     c = p_GetComp(syzterm, r) - 1;
    1088 
    1089   assume( c >= 0 && c < IDELEMS(L) );
    1090  
    1091 #ifndef NDEBUG
    1092   const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE)));
    1093 #else
    1094   const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE)));
    1095 #endif
    1096 
    1097   long debug = __DEBUG__;
    1098   const BOOLEAN __SYZCHECK__ = (BOOLEAN)((long)(atGet(currRingHdl,"SYZCHECK",INT_CMD, (void*)debug)));   
    1099  
    1100   if (__SYZCHECK__ && syzterm != NULL)
    1101   {
    1102     const poly m = L->m[c];
    1103 
    1104     assume( m != NULL ); assume( pNext(m) == NULL );
    1105 
    1106     poly lm = p_Mult_mm(leadmonom(syzterm, r), m, r);
    1107     assume( p_EqualPolys(lm, product, r) );
    1108 
    1109     //  def @@c = leadcomp(syzterm); int @@r = int(@@c);
    1110     //  def @@product = leadmonomial(syzterm) * L[@@r];
    1111 
    1112     p_Delete(&lm, r);   
    1113   }
    1114  
    1115   // looking for an appropriate diviser q = L[k]...
    1116   for( int k = IDELEMS(L)-1; k>= 0; k-- )
    1117   {
    1118     const poly p = L->m[k];   
    1119 
    1120     // ... which divides the product, looking for the _1st_ appropriate one!
    1121     if( !p_LmDivisibleBy(p, product, r) )
    1122       continue;
    1123 
    1124 
    1125     const poly q = p_New(r);
    1126     pNext(q) = NULL;
    1127     p_ExpVectorDiff(q, product, p, r); // (LM(product) / LM(L[k]))
    1128     p_SetComp(q, k + 1, r);
    1129     p_Setm(q, r);
    1130 
    1131     // cannot allow something like: a*gen(i) - a*gen(i)
    1132     if (syzterm != NULL && (k == c))
    1133       if (p_ExpVectorEqual(syzterm, q, r))
    1134       {
    1135         if( __DEBUG__ )
    1136         {
    1137           Print("FindReducer::Test SYZTERM: q == syzterm !:((, syzterm is: ");
    1138           dPrint(syzterm, r, r, 1);
    1139         }   
    1140 
    1141         p_LmFree(q, r);
    1142         continue;
    1143       }
    1144 
    1145     // while the complement (the fraction) is not reducible by leading syzygies
    1146     if( LS != NULL )
    1147     {
    1148       BOOLEAN ok = TRUE;
    1149 
    1150       for(int kk = IDELEMS(LS)-1; kk>= 0; kk-- )
    1151       {
    1152         const poly pp = LS->m[kk];
    1153 
    1154         if( p_LmDivisibleBy(pp, q, r) )
    1155         {
    1156          
    1157           if( __DEBUG__ )
    1158           {
    1159             Print("FindReducer::Test LS: q is divisible by LS[%d] !:((, diviser is: ", kk+1);
    1160             dPrint(pp, r, r, 1);
    1161           }   
    1162 
    1163           ok = FALSE; // q in <LS> :((
    1164           break;                 
    1165         }
    1166       }
    1167 
    1168       if(!ok)
    1169       {
    1170         p_LmFree(q, r);
    1171         continue;
    1172       }
    1173     }
    1174 
    1175     p_SetCoeff0(q, n_Neg( n_Div( p_GetCoeff(product, r), p_GetCoeff(p, r), r), r), r);
    1176     return q;
    1177 
    1178   }
    1179 
    1180  
    1181   return NULL;
    1182 }
    1183 
    1184685
    1185686/// TODO: save shortcut (syz: |-.->) LM(LM(m) * "t") -> syz?
    1186687/// proc SSFindReducer(def product, def syzterm, def L, def T, list #)
    1187 static BOOLEAN FindReducer(leftv res, leftv h)
     688static BOOLEAN _FindReducer(leftv res, leftv h)
    1188689{
    1189690  const char* usage = "`FindReducer(<poly/vector>, <vector/0>, <ideal/module>[,<module>])` expected";
    … …  
    1278779}
    1279780
    1280 static poly ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck,
    1281                        ideal L, ideal T, ideal LS, const ring r);
    1282 
    1283 static poly TraverseTail(poly multiplier, poly tail,
    1284                          ideal L, ideal T, ideal LS,
    1285                          const ring r)
    1286 {
    1287   assume( multiplier != NULL );
    1288 
    1289   assume( L != NULL );
    1290   assume( T != NULL );
    1291 
    1292   poly s = NULL;
    1293 
    1294   // iterate over the tail
    1295   for(poly p = tail; p != NULL; p = pNext(p))
    1296     s = p_Add_q(s, ReduceTerm(multiplier, p, NULL, L, T, LS, r), r); 
    1297    
    1298   return s;
    1299 }
    1300 
    1301 
    1302 static poly ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck,
    1303                        ideal L, ideal T, ideal LS, const ring r)
    1304 {
    1305   assume( multiplier != NULL );
    1306   assume( term4reduction != NULL );
    1307 
    1308 
    1309   assume( L != NULL );
    1310   assume( T != NULL );
    1311  
    1312   // assume(r == currRing); // ?
    1313 
    1314   // simple implementation with FindReducer:
    1315   poly s = NULL;
    1316  
    1317   if(1)
    1318   {
    1319     // NOTE: only LT(term4reduction) should be used in the following:
    1320     poly product = pp_Mult_mm(multiplier, term4reduction, r);
    1321     s = FindReducer(product, syztermCheck, L, LS, r);
    1322     p_Delete(&product, r);
    1323   }
    1324 
    1325   if( s == NULL ) // No Reducer?
    1326     return s;
    1327 
    1328   poly b = leadmonom(s, r);
    1329 
    1330   const int c = p_GetComp(s, r) - 1;
    1331   assume( c >= 0 && c < IDELEMS(T) );
    1332 
    1333   const poly tail = T->m[c];
    1334 
    1335   if( tail != NULL )
    1336     s = p_Add_q(s, TraverseTail(b, tail, L, T, LS, r), r); 
    1337  
    1338   return s;
    1339 }
    1340 
    1341 
    1342 static poly SchreyerSyzygyNF(poly syz_lead, poly syz_2, ideal L, ideal T, ideal LS, const ring r)
    1343 {
    1344   assume( syz_lead != NULL );
    1345   assume( syz_2 != NULL );
    1346 
    1347   assume( L != NULL );
    1348   assume( T != NULL );
    1349 
    1350   assume( IDELEMS(L) == IDELEMS(T) );
    1351 
    1352   int  c = p_GetComp(syz_lead, r) - 1;
    1353 
    1354   assume( c >= 0 && c < IDELEMS(T) );
    1355 
    1356   poly p = leadmonom(syz_lead, r); // :( 
    1357   poly spoly = pp_Mult_qq(p, T->m[c], r);
    1358   p_Delete(&p, r);
    1359 
    1360  
    1361   c = p_GetComp(syz_2, r) - 1;
    1362   assume( c >= 0 && c < IDELEMS(T) );
    1363 
    1364   p = leadmonom(syz_2, r); // :(
    1365   spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);
    1366   p_Delete(&p, r);
    1367 
    1368   poly tail = p_Copy(syz_2, r); // TODO: use bucket!?
    1369 
    1370   while (spoly != NULL)
    1371   {
    1372     poly t = FindReducer(spoly, NULL, L, LS, r);
    1373 
    1374     p_LmDelete(&spoly, r);
    1375    
    1376     if( t != NULL )
    1377     {
    1378       p = leadmonom(t, r); // :(
    1379       c = p_GetComp(t, r) - 1;
    1380 
    1381       assume( c >= 0 && c < IDELEMS(T) );
    1382      
    1383       spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);
    1384 
    1385       p_Delete(&p, r);
    1386      
    1387       tail = p_Add_q(tail, t, r);
    1388     }   
    1389   }
    1390  
    1391   return tail;
    1392 }
    1393 
    1394781// proc SchreyerSyzygyNF(vector syz_lead, vector syz_2, def L, def T, list #)
    1395 static BOOLEAN SchreyerSyzygyNF(leftv res, leftv h)
     782static BOOLEAN _SchreyerSyzygyNF(leftv res, leftv h)
    1396783{
    1397784  const char* usage = "`SchreyerSyzygyNF(<vector>, <vector>, <ideal/module>, <ideal/module>[,<module>])` expected";
    … …  
    1503890/// TODO: save shortcut (syz: |-.->) LM(m) * "t" -> ?
    1504891/// proc SSReduceTerm(poly m, def t, def syzterm, def L, def T, list #)
    1505 static BOOLEAN ReduceTerm(leftv res, leftv h)
     892static BOOLEAN _ReduceTerm(leftv res, leftv h)
    1506893{
    1507894  const char* usage = "`ReduceTerm(<poly>, <poly/vector>, <vector/0>, <ideal/module>, <ideal/module>[,<module>])` expected";
    … …  
    16481035// TODO: store m * @tail -.-^-.-^-.--> ?
    16491036// proc SSTraverseTail(poly m, def @tail, def L, def T, list #)
    1650 static BOOLEAN TraverseTail(leftv res, leftv h)
     1037static BOOLEAN _TraverseTail(leftv res, leftv h)
    16511038{
    16521039  const char* usage = "`TraverseTail(<poly>, <poly/vector>, <ideal/module>, <ideal/module>[,<module>])` expected";
    … …  
    17501137
    17511138
    1752 
    1753 static void ComputeSyzygy(const ideal L, const ideal T, ideal& LL, ideal& TT, const ring R)
    1754 {
    1755   assume( R == currRing ); // For attributes :-/
    1756 
    1757   assume( IDELEMS(L) == IDELEMS(T) );
    1758  
     1139// module (N, LL, TT) = SSComputeSyzygy(L, T);
     1140// Compute Syz(L ++ T) = N = LL ++ TT
     1141// proc SSComputeSyzygy(def L, def T)
     1142static BOOLEAN _ComputeSyzygy(leftv res, leftv h)
     1143{
     1144  const char* usage = "`ComputeSyzygy(<ideal/module>, <ideal/module>])` expected";
     1145  const ring r = currRing;
     1146
     1147  NoReturn(res);
     1148
    17591149#ifndef NDEBUG
    17601150  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE)));
    … …  
    17631153#endif
    17641154
    1765   const BOOLEAN __LEAD2SYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"LEAD2SYZ",INT_CMD, (void*)1)));
    1766   const BOOLEAN __TAILREDSYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"TAILREDSYZ",INT_CMD, (void*)1)));
    1767   const BOOLEAN __SYZCHECK__ = (BOOLEAN)((long)(atGet(currRingHdl,"SYZCHECK",INT_CMD, (void*)__DEBUG__)));
    1768 
    1769   const BOOLEAN __HYBRIDNF__ = (BOOLEAN)((long)(atGet(currRingHdl,"HYBRIDNF",INT_CMD, (void*)0)));
    1770 
    1771 
    1772   if( __LEAD2SYZ__ )
    1773     LL = Compute2LeadingSyzygyTerms(L, R); // 2 terms!
    1774   else
    1775     LL = ComputeLeadingSyzygyTerms(L, R); // 1 term!
    1776 
    1777   const int size = IDELEMS(LL);
    1778 
    1779   TT = idInit(size, 0);
    1780 
    1781   if( size == 1 && LL->m[0] == NULL )
    1782     return;
    1783  
    1784 
    1785   ideal LS = NULL;
    1786  
    1787   if( __TAILREDSYZ__ )
    1788     LS = LL;
    1789 
    1790   for( int k = size - 1; k >= 0; k-- )
    1791   {
    1792     const poly a = LL->m[k]; assume( a != NULL );
    1793    
    1794     const int r = p_GetComp(a, R) - 1;
    1795    
    1796     assume( r >= 0 && r < IDELEMS(T) );
    1797     assume( r >= 0 && r < IDELEMS(L) );
    1798 
    1799     poly aa = leadmonom(a, R); assume( aa != NULL); // :(   
    1800    
    1801     poly a2 = pNext(a);   
    1802 
    1803     if( a2 != NULL )
    1804     {
    1805       TT->m[k] = a2; pNext(a) = NULL;
    1806     }
    1807 
    1808     if( ! __HYBRIDNF__ )
    1809     {
    1810       poly t = TraverseTail(aa, T->m[r], L, T, LS, R);
    1811 
    1812       if( a2 != NULL )
    1813       {
    1814         assume( __LEAD2SYZ__ );
    1815 
    1816         const int r2 = p_GetComp(a2, R) - 1; poly aa2 = leadmonom(a2, R); // :(
    1817 
    1818         assume( r2 >= 0 && r2 < IDELEMS(T) );
    1819        
    1820         TT->m[k] = p_Add_q(a2, p_Add_q(t, TraverseTail(aa2, T->m[r2], L, T, LS, R), R), R);
    1821 
    1822         p_Delete(&aa2, R);
    1823       } else
    1824         TT->m[k] = p_Add_q(t, ReduceTerm(aa, L->m[r], a, L, T, LS, R), R);
    1825 
    1826     } else
    1827     {
    1828       if( a2 == NULL )
    1829       {
    1830         aa = p_Mult_mm(aa, L->m[r], R);
    1831         a2 = FindReducer(aa, a, L, LS, R);
    1832       }
    1833       assume( a2 != NULL );
    1834        
    1835       TT->m[k] = SchreyerSyzygyNF(a, a2, L, T, LS, R); // will copy a2 :(
    1836      
    1837       p_Delete(&a2, R);
    1838     }
    1839     p_Delete(&aa, R);   
    1840   }
    1841 
    1842   TT->rank = id_RankFreeModule(TT, R);
    1843 }
    1844 
    1845 
    1846 
    1847 // module (N, LL, TT) = SSComputeSyzygy(L, T);
    1848 // Compute Syz(L ++ T) = N = LL ++ TT
    1849 // proc SSComputeSyzygy(def L, def T)
    1850 static BOOLEAN ComputeSyzygy(leftv res, leftv h)
    1851 {
    1852   const char* usage = "`ComputeSyzygy(<ideal/module>, <ideal/module>])` expected";
    1853   const ring r = currRing;
    1854 
    1855   NoReturn(res);
    1856 
    1857 #ifndef NDEBUG
    1858   const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE)));
    1859 #else
    1860   const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE)));
    1861 #endif
    1862 
    18631155  if ((h==NULL) || (h->Typ()!=IDEAL_CMD && h->Typ() !=MODUL_CMD) || (h->Data() == NULL))
    18641156  {
    … …  
    19151207}
    19161208
    1917 
    1918 
    1919 
    1920 
    1921 
    1922 
    19231209/// Get leading term without a module component
    1924 static BOOLEAN leadmonom(leftv res, leftv h)
     1210static BOOLEAN _leadmonom(leftv res, leftv h)
    19251211{
    19261212  NoReturn(res);
    … …  
    19311217    const poly p = (poly)(h->Data());
    19321218
    1933     const poly m = leadmonom(p, r);
    1934 
     1219    res->data = reinterpret_cast<void *>(  leadmonom(p, r) );
    19351220    res->rtyp = POLY_CMD;
    1936     res->data = reinterpret_cast<void *>(m);
    19371221
    19381222    return FALSE;
    … …  
    25301814
    25311815  ADD(psModulFunctions, currPack->libname, "DetailedPrint", FALSE, DetailedPrint);
    2532   ADD(psModulFunctions, currPack->libname, "leadmonomial", FALSE, leadmonom);
     1816  ADD(psModulFunctions, currPack->libname, "leadmonomial", FALSE, _leadmonom);
    25331817  ADD(psModulFunctions, currPack->libname, "leadcomp", FALSE, leadcomp);
    25341818  ADD(psModulFunctions, currPack->libname, "leadrawexp", FALSE, leadrawexp);
    … …  
    25511835  ADD(psModulFunctions, currPack->libname, "Tail", FALSE, Tail);
    25521836 
    2553   ADD(psModulFunctions, currPack->libname, "ComputeLeadingSyzygyTerms", FALSE, ComputeLeadingSyzygyTerms);
    2554   ADD(psModulFunctions, currPack->libname, "Compute2LeadingSyzygyTerms", FALSE, Compute2LeadingSyzygyTerms);
    2555  
    2556   ADD(psModulFunctions, currPack->libname, "Sort_c_ds", FALSE, Sort_c_ds);
    2557   ADD(psModulFunctions, currPack->libname, "FindReducer", FALSE, FindReducer);
    2558 
    2559 
    2560   ADD(psModulFunctions, currPack->libname, "ReduceTerm", FALSE, ReduceTerm);
    2561   ADD(psModulFunctions, currPack->libname, "TraverseTail", FALSE, TraverseTail);
     1837  ADD(psModulFunctions, currPack->libname, "ComputeLeadingSyzygyTerms", FALSE, _ComputeLeadingSyzygyTerms);
     1838  ADD(psModulFunctions, currPack->libname, "Compute2LeadingSyzygyTerms", FALSE, _Compute2LeadingSyzygyTerms);
     1839 
     1840  ADD(psModulFunctions, currPack->libname, "Sort_c_ds", FALSE, _Sort_c_ds);
     1841  ADD(psModulFunctions, currPack->libname, "FindReducer", FALSE, _FindReducer);
     1842
     1843
     1844  ADD(psModulFunctions, currPack->libname, "ReduceTerm", FALSE, _ReduceTerm);
     1845  ADD(psModulFunctions, currPack->libname, "TraverseTail", FALSE, _TraverseTail);
    25621846
    25631847   
    2564   ADD(psModulFunctions, currPack->libname, "SchreyerSyzygyNF", FALSE, SchreyerSyzygyNF);
    2565   ADD(psModulFunctions, currPack->libname, "ComputeSyzygy", FALSE, ComputeSyzygy);
     1848  ADD(psModulFunctions, currPack->libname, "SchreyerSyzygyNF", FALSE, _SchreyerSyzygyNF);
     1849  ADD(psModulFunctions, currPack->libname, "ComputeSyzygy", FALSE, _ComputeSyzygy);
    25661850 
    25671851  //  ADD(psModulFunctions, currPack->libname, "GetAMData", FALSE, GetAMData);

  • dyn_modules/syzextra/syzextra.cc

    rcd5fefc r204092  
    1919#include "syzextra.h"
    2020
     21#include "DebugPrint.h"
     22
    2123#include <omalloc/omalloc.h>
     24
     25#include <misc/intvec.h>
     26#include <misc/options.h>
     27
     28#include <coeffs/coeffs.h>
     29
    2230#include <polys/monomials/p_polys.h>
    23 
     31#include <polys/monomials/ring.h>
     32
     33#include <kernel/kstd1.h>
     34#include <kernel/polys.h>
     35#include <kernel/syz.h>
    2436#include <kernel/ideals.h>
    2537
    2638
     39
     40#include <Singular/tok.h>
     41#include <Singular/ipid.h>
     42#include <Singular/lists.h>
     43#include <Singular/attrib.h>
     44
     45#include <Singular/ipid.h>
     46#include <Singular/ipshell.h> // For iiAddCproc
     47
     48#include <stdio.h>
     49#include <stdlib.h>
     50#include <string.h>
     51
     52// USING_NAMESPACE_SINGULARXX;
     53USING_NAMESPACE( SINGULARXXNAME :: DEBUG )
     54
     55
    2756BEGIN_NAMESPACE_SINGULARXX     BEGIN_NAMESPACE(SYZEXTRA)
    2857
    29 
     58BEGIN_NAMESPACE(SORT_c_ds)
     59
     60
     61#ifdef _GNU_SOURCE
     62static int cmp_c_ds(const void *p1, const void *p2, void *R)
     63{
     64#else
     65static int cmp_c_ds(const void *p1, const void *p2)
     66{
     67  void *R = currRing;
     68#endif
     69
     70  const int YES = 1;
     71  const int NO = -1;
     72
     73  const ring r =  (const ring) R; // TODO/NOTE: the structure is known: C, lp!!!
     74
     75  assume( r == currRing );
     76
     77  const poly a = *(const poly*)p1;
     78  const poly b = *(const poly*)p2;
     79
     80  assume( a != NULL );
     81  assume( b != NULL );
     82
     83  assume( p_LmTest(a, r) );
     84  assume( p_LmTest(b, r) );
     85
     86
     87  const signed long iCompDiff = p_GetComp(a, r) - p_GetComp(b, r);
     88
     89  // TODO: test this!!!!!!!!!!!!!!!!
     90
     91  //return -( compare (c, qsorts) )
     92
     93  const int __DEBUG__ = 0;
     94
     95#ifndef NDEBUG
     96  if( __DEBUG__ )
     97  {
     98    PrintS("cmp_c_ds: a, b: \np1: "); dPrint(a, r, r, 2);
     99    PrintS("b: "); dPrint(b, r, r, 2);
     100    PrintLn();
     101  }
     102#endif
     103
     104
     105  if( iCompDiff > 0 )
     106    return YES;
     107
     108  if( iCompDiff < 0 )
     109    return  NO;
     110
     111  assume( iCompDiff == 0 );
     112
     113  const signed long iDegDiff = p_Totaldegree(a, r) - p_Totaldegree(b, r);
     114
     115  if( iDegDiff > 0 )
     116    return YES;
     117
     118  if( iDegDiff < 0 )
     119    return  NO;
     120
     121  assume( iDegDiff == 0 );
     122
     123#ifndef NDEBUG
     124  if( __DEBUG__ )
     125  {
     126    PrintS("cmp_c_ds: a & b have the same comp & deg! "); PrintLn();
     127  }
     128#endif
     129
     130  for (int v = rVar(r); v > 0; v--)
     131  {
     132    assume( v > 0 );
     133    assume( v <= rVar(r) );
     134
     135    const signed int d = p_GetExp(a, v, r) - p_GetExp(b, v, r);
     136
     137    if( d > 0 )
     138      return YES;
     139
     140    if( d < 0 )
     141      return NO;
     142
     143    assume( d == 0 );
     144  }
     145
     146  return 0; 
     147}
     148
     149END_NAMESPACE
     150/* namespace SORT_c_ds */
     151
     152/// return a new term: leading coeff * leading monomial of p
     153/// with 0 leading component!
     154poly leadmonom(const poly p, const ring r)
     155{
     156  poly m = NULL;
     157
     158  if( p != NULL )
     159  {
     160    assume( p != NULL );
     161    assume( p_LmTest(p, r) );
     162
     163    m = p_LmInit(p, r);
     164    p_SetCoeff0(m, n_Copy(p_GetCoeff(p, r), r), r);
     165
     166    p_SetComp(m, 0, r);
     167    p_Setm(m, r);
     168
     169    assume( p_GetComp(m, r) == 0 );
     170    assume( m != NULL );
     171    assume( pNext(m) == NULL );
     172    assume( p_LmTest(m, r) );
     173  }
     174
     175  return m;
     176}
     177
     178
     179   
    30180poly p_Tail(const poly p, const ring r)
    31181{
    … …  
    54204
    55205
     206void Sort_c_ds(const ideal id, const ring r)
     207{
     208  const int sizeNew = IDELEMS(id);
     209
     210#ifdef _GNU_SOURCE
     211#define qsort_my(m, s, ss, r, cmp) qsort_r(m, s, ss, cmp, r)
     212#else
     213#define qsort_my(m, s, ss, r, cmp) qsort_r(m, s, ss, cmp)
     214#endif
     215 
     216  if( sizeNew >= 2 )
     217    qsort_my(id->m, sizeNew, sizeof(poly), r, SORT_c_ds::cmp_c_ds);
     218
     219#undef qsort_my
     220
     221  id->rank = id_RankFreeModule(id, r);
     222}
     223
     224
     225
     226ideal ComputeLeadingSyzygyTerms(const ideal& id, const ring r)
     227{
     228  // 1. set of components S?
     229  // 2. for each component c from S: set of indices of leading terms
     230  // with this component?
     231  // 3. short exp. vectors for each leading term?
     232
     233  const int size = IDELEMS(id);
     234
     235  if( size < 2 )
     236  {
     237    const ideal newid = idInit(1, 0); newid->m[0] = NULL; // zero ideal...       
     238    return newid;
     239  }
     240
     241
     242  // TODO/NOTE: input is supposed to be (reverse-) sorted wrt "(c,ds)"!??
     243
     244  // components should come in groups: count elements in each group
     245  // && estimate the real size!!!
     246
     247
     248  // use just a vector instead???
     249  const ideal newid = idInit( (size * (size-1))/2, size); // maximal size: ideal case!
     250
     251  int k = 0;
     252
     253  for (int j = 0; j < size; j++)
     254  {
     255    const poly p = id->m[j];
     256    assume( p != NULL );
     257    const int  c = p_GetComp(p, r);
     258
     259    for (int i = j - 1; i >= 0; i--)
     260    {
     261      const poly pp = id->m[i];
     262      assume( pp != NULL );
     263      const int  cc = p_GetComp(pp, r);
     264
     265      if( c != cc )
     266        continue;
     267
     268      const poly m = p_Init(r); // p_New???
     269
     270      // m = LCM(p, pp) / p! // TODO: optimize: knowing the ring structure: (C/lp)!
     271      for (int v = rVar(r); v > 0; v--)
     272      {
     273        assume( v > 0 );
     274        assume( v <= rVar(r) );
     275
     276        const short e1 = p_GetExp(p , v, r);
     277        const short e2 = p_GetExp(pp, v, r);
     278
     279        if( e1 >= e2 )
     280          p_SetExp(m, v, 0, r);
     281        else
     282          p_SetExp(m, v, e2 - e1, r);
     283
     284      }
     285
     286      assume( (j > i) && (i >= 0) );
     287
     288      p_SetComp(m, j + 1, r);
     289      pNext(m) = NULL;
     290      p_SetCoeff0(m, n_Init(1, r->cf), r); // for later...
     291
     292      p_Setm(m, r); // should not do anything!!!
     293
     294      newid->m[k++] = m;
     295    }
     296  }
     297
     298//   if( __DEBUG__ && FALSE )
     299//   {
     300//     PrintS("ComputeLeadingSyzygyTerms::Temp0: \n");
     301//     dPrint(newid, r, r, 1);
     302//   }
     303
     304  // the rest of newid is assumed to be zeroes...
     305
     306  // simplify(newid, 2 + 32)??
     307  // sort(newid, "C,ds")[1]???
     308  id_DelDiv(newid, r); // #define SIMPL_LMDIV 32
     309
     310//   if( __DEBUG__ && FALSE )
     311//   {
     312//     PrintS("ComputeLeadingSyzygyTerms::Temp1: \n");
     313//     dPrint(newid, r, r, 1);
     314//   }
     315
     316  idSkipZeroes(newid); // #define SIMPL_NULL 2
     317
     318//   if( __DEBUG__ )
     319//   {
     320//     PrintS("ComputeLeadingSyzygyTerms::Output: \n");
     321//     dPrint(newid, r, r, 1);
     322//   }
     323
     324  Sort_c_ds(newid, r);
     325
     326  return newid;
     327}
     328
     329
     330
     331
     332ideal Compute2LeadingSyzygyTerms(const ideal& id, const ring r)
     333{
     334#ifndef NDEBUG
     335  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE)));
     336#else
     337  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE)));
     338#endif
     339
     340  const BOOLEAN __TAILREDSYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"TAILREDSYZ",INT_CMD, (void*)0)));
     341  const BOOLEAN __SYZCHECK__ = (BOOLEAN)((long)(atGet(currRingHdl,"SYZCHECK",INT_CMD, (void*)0)));   
     342
     343
     344
     345  // 1. set of components S?
     346  // 2. for each component c from S: set of indices of leading terms
     347  // with this component?
     348  // 3. short exp. vectors for each leading term?
     349
     350  const int size = IDELEMS(id);
     351
     352  if( size < 2 )
     353  {
     354    const ideal newid = idInit(1, 1); newid->m[0] = NULL; // zero module...   
     355    return newid;
     356  }
     357
     358
     359    // TODO/NOTE: input is supposed to be sorted wrt "C,ds"!??
     360
     361    // components should come in groups: count elements in each group
     362    // && estimate the real size!!!
     363
     364
     365    // use just a vector instead???
     366  ideal newid = idInit( (size * (size-1))/2, size); // maximal size: ideal case!
     367
     368  int k = 0;
     369
     370  for (int j = 0; j < size; j++)
     371  {
     372    const poly p = id->m[j];
     373    assume( p != NULL );
     374    const int  c = p_GetComp(p, r);
     375
     376    for (int i = j - 1; i >= 0; i--)
     377    {
     378      const poly pp = id->m[i];
     379      assume( pp != NULL );
     380      const int  cc = p_GetComp(pp, r);
     381
     382      if( c != cc )
     383        continue;
     384
     385        // allocate memory & zero it out!
     386      const poly m = p_Init(r); const poly mm = p_Init(r);
     387
     388
     389        // m = LCM(p, pp) / p! mm = LCM(p, pp) / pp!
     390        // TODO: optimize: knowing the ring structure: (C/lp)!
     391
     392      for (int v = rVar(r); v > 0; v--)
     393      {
     394        assume( v > 0 );
     395        assume( v <= rVar(r) );
     396
     397        const short e1 = p_GetExp(p , v, r);
     398        const short e2 = p_GetExp(pp, v, r);
     399
     400        if( e1 >= e2 )
     401          p_SetExp(mm, v, e1 - e2, r); //            p_SetExp(m, v, 0, r);
     402        else
     403          p_SetExp(m, v, e2 - e1, r); //            p_SetExp(mm, v, 0, r);
     404
     405      }
     406
     407      assume( (j > i) && (i >= 0) );
     408
     409      p_SetComp(m, j + 1, r);
     410      p_SetComp(mm, i + 1, r);
     411
     412      const number& lc1 = p_GetCoeff(p , r);
     413      const number& lc2 = p_GetCoeff(pp, r);
     414
     415      number g = n_Lcm( lc1, lc2, r );
     416
     417      p_SetCoeff0(m ,       n_Div(g, lc1, r), r);
     418      p_SetCoeff0(mm, n_Neg(n_Div(g, lc2, r), r), r);
     419
     420      n_Delete(&g, r);
     421
     422      p_Setm(m, r); // should not do anything!!!
     423      p_Setm(mm, r); // should not do anything!!!
     424
     425      pNext(m) = mm; //        pNext(mm) = NULL;
     426
     427      newid->m[k++] = m;
     428    }
     429  }
     430
     431//   if( __DEBUG__ && FALSE )
     432//   {
     433//     PrintS("Compute2LeadingSyzygyTerms::Temp0: \n");
     434//     dPrint(newid, r, r, 1);
     435//   }
     436
     437  if( !__TAILREDSYZ__ )
     438  {
     439      // simplify(newid, 2 + 32)??
     440      // sort(newid, "C,ds")[1]???
     441    id_DelDiv(newid, r); // #define SIMPL_LMDIV 32
     442
     443//     if( __DEBUG__ && FALSE )
     444//     {
     445//       PrintS("Compute2LeadingSyzygyTerms::Temp1 (deldiv): \n");
     446//       dPrint(newid, r, r, 1);
     447//     }
     448  }
     449  else
     450  {
     451      //      option(redSB); option(redTail);
     452      //      TEST_OPT_REDSB
     453      //      TEST_OPT_REDTAIL
     454    assume( r == currRing );
     455    BITSET _save_test = test;
     456    test |= (Sy_bit(OPT_REDTAIL) | Sy_bit(OPT_REDSB));
     457
     458    intvec* w=new intvec(IDELEMS(newid));
     459    ideal tmp = kStd(newid, currQuotient, isHomog, &w);
     460    delete w;
     461
     462    test = _save_test;
     463
     464    id_Delete(&newid, r);
     465    newid = tmp;
     466
     467//     if( __DEBUG__ && FALSE )
     468//     {
     469//       PrintS("Compute2LeadingSyzygyTerms::Temp1 (std): \n");
     470//       dPrint(newid, r, r, 1);
     471//     }
     472
     473  }
     474
     475  idSkipZeroes(newid);
     476
     477  Sort_c_ds(newid, r);
     478 
     479  return newid;
     480}
     481
     482poly FindReducer(poly product, poly syzterm,
     483                        ideal L, ideal LS,
     484                        const ring r)
     485{
     486  assume( product != NULL );
     487  assume( L != NULL );
     488
     489  int c = 0;
     490
     491  if (syzterm != NULL)
     492    c = p_GetComp(syzterm, r) - 1;
     493
     494  assume( c >= 0 && c < IDELEMS(L) );
     495
     496#ifndef NDEBUG
     497  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE)));
     498#else
     499  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE)));
     500#endif
     501
     502  long debug = __DEBUG__;
     503  const BOOLEAN __SYZCHECK__ = (BOOLEAN)((long)(atGet(currRingHdl,"SYZCHECK",INT_CMD, (void*)debug)));   
     504
     505  if (__SYZCHECK__ && syzterm != NULL)
     506  {
     507    const poly m = L->m[c];
     508
     509    assume( m != NULL ); assume( pNext(m) == NULL );
     510
     511    poly lm = p_Mult_mm(leadmonom(syzterm, r), m, r);
     512    assume( p_EqualPolys(lm, product, r) );
     513
     514    //  def @@c = leadcomp(syzterm); int @@r = int(@@c);
     515    //  def @@product = leadmonomial(syzterm) * L[@@r];
     516
     517    p_Delete(&lm, r);   
     518  }
     519
     520  // looking for an appropriate diviser q = L[k]...
     521  for( int k = IDELEMS(L)-1; k>= 0; k-- )
     522  {
     523    const poly p = L->m[k];   
     524
     525    // ... which divides the product, looking for the _1st_ appropriate one!
     526    if( !p_LmDivisibleBy(p, product, r) )
     527      continue;
     528
     529
     530    const poly q = p_New(r);
     531    pNext(q) = NULL;
     532    p_ExpVectorDiff(q, product, p, r); // (LM(product) / LM(L[k]))
     533    p_SetComp(q, k + 1, r);
     534    p_Setm(q, r);
     535
     536    // cannot allow something like: a*gen(i) - a*gen(i)
     537    if (syzterm != NULL && (k == c))
     538      if (p_ExpVectorEqual(syzterm, q, r))
     539      {
     540        if( __DEBUG__ )
     541        {
     542          Print("FindReducer::Test SYZTERM: q == syzterm !:((, syzterm is: ");
     543          dPrint(syzterm, r, r, 1);
     544        }   
     545
     546        p_LmFree(q, r);
     547        continue;
     548      }
     549
     550    // while the complement (the fraction) is not reducible by leading syzygies
     551    if( LS != NULL )
     552    {
     553      BOOLEAN ok = TRUE;
     554
     555      for(int kk = IDELEMS(LS)-1; kk>= 0; kk-- )
     556      {
     557        const poly pp = LS->m[kk];
     558
     559        if( p_LmDivisibleBy(pp, q, r) )
     560        {
     561
     562          if( __DEBUG__ )
     563          {
     564            Print("FindReducer::Test LS: q is divisible by LS[%d] !:((, diviser is: ", kk+1);
     565            dPrint(pp, r, r, 1);
     566          }   
     567
     568          ok = FALSE; // q in <LS> :((
     569          break;                 
     570        }
     571      }
     572
     573      if(!ok)
     574      {
     575        p_LmFree(q, r);
     576        continue;
     577      }
     578    }
     579
     580    p_SetCoeff0(q, n_Neg( n_Div( p_GetCoeff(product, r), p_GetCoeff(p, r), r), r), r);
     581    return q;
     582
     583  }
     584
     585
     586  return NULL;
     587}
     588
     589poly TraverseTail(poly multiplier, poly tail,
     590                         ideal L, ideal T, ideal LS,
     591                         const ring r)
     592{
     593  assume( multiplier != NULL );
     594
     595  assume( L != NULL );
     596  assume( T != NULL );
     597
     598  poly s = NULL;
     599
     600  // iterate over the tail
     601  for(poly p = tail; p != NULL; p = pNext(p))
     602    s = p_Add_q(s, ReduceTerm(multiplier, p, NULL, L, T, LS, r), r); 
     603
     604  return s;
     605}
     606
     607
     608poly ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck,
     609                       ideal L, ideal T, ideal LS, const ring r)
     610{
     611  assume( multiplier != NULL );
     612  assume( term4reduction != NULL );
     613
     614
     615  assume( L != NULL );
     616  assume( T != NULL );
     617
     618  // assume(r == currRing); // ?
     619
     620  // simple implementation with FindReducer:
     621  poly s = NULL;
     622
     623  if(1)
     624  {
     625    // NOTE: only LT(term4reduction) should be used in the following:
     626    poly product = pp_Mult_mm(multiplier, term4reduction, r);
     627    s = FindReducer(product, syztermCheck, L, LS, r);
     628    p_Delete(&product, r);
     629  }
     630
     631  if( s == NULL ) // No Reducer?
     632    return s;
     633
     634  poly b = leadmonom(s, r);
     635
     636  const int c = p_GetComp(s, r) - 1;
     637  assume( c >= 0 && c < IDELEMS(T) );
     638
     639  const poly tail = T->m[c];
     640
     641  if( tail != NULL )
     642    s = p_Add_q(s, TraverseTail(b, tail, L, T, LS, r), r); 
     643
     644  return s;
     645}
     646
     647
     648poly SchreyerSyzygyNF(poly syz_lead, poly syz_2, ideal L, ideal T, ideal LS, const ring r)
     649{
     650  assume( syz_lead != NULL );
     651  assume( syz_2 != NULL );
     652
     653  assume( L != NULL );
     654  assume( T != NULL );
     655
     656  assume( IDELEMS(L) == IDELEMS(T) );
     657
     658  int  c = p_GetComp(syz_lead, r) - 1;
     659
     660  assume( c >= 0 && c < IDELEMS(T) );
     661
     662  poly p = leadmonom(syz_lead, r); // :( 
     663  poly spoly = pp_Mult_qq(p, T->m[c], r);
     664  p_Delete(&p, r);
     665
     666
     667  c = p_GetComp(syz_2, r) - 1;
     668  assume( c >= 0 && c < IDELEMS(T) );
     669
     670  p = leadmonom(syz_2, r); // :(
     671  spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);
     672  p_Delete(&p, r);
     673
     674  poly tail = p_Copy(syz_2, r); // TODO: use bucket!?
     675
     676  while (spoly != NULL)
     677  {
     678    poly t = FindReducer(spoly, NULL, L, LS, r);
     679
     680    p_LmDelete(&spoly, r);
     681
     682    if( t != NULL )
     683    {
     684      p = leadmonom(t, r); // :(
     685      c = p_GetComp(t, r) - 1;
     686
     687      assume( c >= 0 && c < IDELEMS(T) );
     688
     689      spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);
     690
     691      p_Delete(&p, r);
     692
     693      tail = p_Add_q(tail, t, r);
     694    }   
     695  }
     696
     697  return tail;
     698}
     699
     700
     701void ComputeSyzygy(const ideal L, const ideal T, ideal& LL, ideal& TT, const ring R)
     702{
     703  assume( R == currRing ); // For attributes :-/
     704
     705  assume( IDELEMS(L) == IDELEMS(T) );
     706
     707#ifndef NDEBUG
     708  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)TRUE)));
     709#else
     710  const BOOLEAN __DEBUG__ = (BOOLEAN)((long)(atGet(currRingHdl,"DEBUG",INT_CMD, (void*)FALSE)));
     711#endif
     712
     713  const BOOLEAN __LEAD2SYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"LEAD2SYZ",INT_CMD, (void*)1)));
     714  const BOOLEAN __TAILREDSYZ__ = (BOOLEAN)((long)(atGet(currRingHdl,"TAILREDSYZ",INT_CMD, (void*)1)));
     715  const BOOLEAN __SYZCHECK__ = (BOOLEAN)((long)(atGet(currRingHdl,"SYZCHECK",INT_CMD, (void*)__DEBUG__)));
     716
     717  const BOOLEAN __HYBRIDNF__ = (BOOLEAN)((long)(atGet(currRingHdl,"HYBRIDNF",INT_CMD, (void*)0)));
     718
     719
     720  if( __LEAD2SYZ__ )
     721    LL = Compute2LeadingSyzygyTerms(L, R); // 2 terms!
     722  else
     723    LL = ComputeLeadingSyzygyTerms(L, R); // 1 term!
     724
     725  const int size = IDELEMS(LL);
     726
     727  TT = idInit(size, 0);
     728
     729  if( size == 1 && LL->m[0] == NULL )
     730    return;
     731
     732
     733  ideal LS = NULL;
     734
     735  if( __TAILREDSYZ__ )
     736    LS = LL;
     737
     738  for( int k = size - 1; k >= 0; k-- )
     739  {
     740    const poly a = LL->m[k]; assume( a != NULL );
     741
     742    const int r = p_GetComp(a, R) - 1;
     743
     744    assume( r >= 0 && r < IDELEMS(T) );
     745    assume( r >= 0 && r < IDELEMS(L) );
     746
     747    poly aa = leadmonom(a, R); assume( aa != NULL); // :(   
     748
     749    poly a2 = pNext(a);   
     750
     751    if( a2 != NULL )
     752    {
     753      TT->m[k] = a2; pNext(a) = NULL;
     754    }
     755
     756    if( ! __HYBRIDNF__ )
     757    {
     758      poly t = TraverseTail(aa, T->m[r], L, T, LS, R);
     759
     760      if( a2 != NULL )
     761      {
     762        assume( __LEAD2SYZ__ );
     763
     764        const int r2 = p_GetComp(a2, R) - 1; poly aa2 = leadmonom(a2, R); // :(
     765
     766        assume( r2 >= 0 && r2 < IDELEMS(T) );
     767
     768        TT->m[k] = p_Add_q(a2, p_Add_q(t, TraverseTail(aa2, T->m[r2], L, T, LS, R), R), R);
     769
     770        p_Delete(&aa2, R);
     771      } else
     772        TT->m[k] = p_Add_q(t, ReduceTerm(aa, L->m[r], a, L, T, LS, R), R);
     773
     774    } else
     775    {
     776      if( a2 == NULL )
     777      {
     778        aa = p_Mult_mm(aa, L->m[r], R);
     779        a2 = FindReducer(aa, a, L, LS, R);
     780      }
     781      assume( a2 != NULL );
     782
     783      TT->m[k] = SchreyerSyzygyNF(a, a2, L, T, LS, R); // will copy a2 :(
     784
     785      p_Delete(&a2, R);
     786    }
     787    p_Delete(&aa, R);   
     788  }
     789
     790  TT->rank = id_RankFreeModule(TT, R);
     791}
     792
     793
     794
    56795
    57796

  • dyn_modules/syzextra/syzextra.h

    rcd5fefc r204092  
    3333BEGIN_NAMESPACE_SINGULARXX    BEGIN_NAMESPACE(SYZEXTRA)
    3434
     35poly leadmonom(const poly p, const ring r);
    3536
    3637/// return the tail of a given polynomial or vector
    … …  
    4748
    4849
     50/// inplace sorting of the module (ideal) id wrt >_(c,ds)
     51void Sort_c_ds(const ideal id, const ring r);
     52
     53ideal ComputeLeadingSyzygyTerms(const ideal& id, const ring r);
     54
     55ideal Compute2LeadingSyzygyTerms(const ideal& id, const ring r);
     56
     57
     58poly FindReducer(poly product, poly syzterm,
     59                 ideal L, ideal LS,
     60                 const ring r);
     61
     62
     63poly TraverseTail(poly multiplier, poly tail,
     64                  ideal L, ideal T, ideal LS,
     65                  const ring r);
     66
     67poly ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck,
     68                ideal L, ideal T, ideal LS, const ring r);
     69
     70
     71
     72poly SchreyerSyzygyNF(poly syz_lead, poly syz_2, ideal L, ideal T, ideal LS, const ring r);
     73
     74
     75void ComputeSyzygy(const ideal L, const ideal T, ideal& LL, ideal& TT, const ring R);
     76
    4977
    5078END_NAMESPACE               END_NAMESPACE_SINGULARXX
Note: See TracChangeset for help on using the changeset viewer.