Ignore:
Timestamp:
Apr 28, 2014, 8:50:10 PM (9 years ago)
Author:
Oleksandr Motsak <motsak@…>
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
Children:
1a4c34381788487464977b182596acb57f904d2f
Parents:
c814238a3c69aaf59662652396c4843d8c8c52fa
git-author:
Oleksandr Motsak <motsak@mathematik.uni-kl.de>2014-04-28 20:50:10+02:00
git-committer:
Oleksandr Motsak <motsak@mathematik.uni-kl.de>2014-05-07 04:41:47+02:00
Message:
Traverse API redesign + avoid multiplication as far as possible

add/chg: TraverseNF similar to SchreyerSyzygyNF
add/chg: use TraverseTail(mult, int tail_index) everywhere

add: usage of CReducerFinder::FindReducer controlled by NOPRODUCT macro define
add: _p_LmDivisibleByNoComp, CReducerFinder::FindReducer for products

NOTE: needs p_GetShortExpVector for 'products'
File:
1 edited

Legend:

Unmodified
Added
Removed

  • dyn_modules/syzextra/syzextra.cc

    rc81423 r1cf13b  
    222222}
    223223
    224 
    225224ideal SchreyerSyzygyComputation::Compute1LeadingSyzygyTerms()
    226225{
    … …  
    487486}
    488487
     488poly SchreyerSyzygyComputation::TraverseNF(const poly a, const poly a2) const
     489{
     490  const ideal& L = m_idLeads;
     491  const ideal& T = m_idTails;
     492
     493  const ring& R = m_rBaseRing;
     494
     495  const int r = p_GetComp(a, R) - 1;
     496
     497  assume( r >= 0 && r < IDELEMS(T) );
     498  assume( r >= 0 && r < IDELEMS(L) );
     499
     500  poly aa = leadmonom(a, R); assume( aa != NULL); // :(
     501  poly t = TraverseTail(aa, r);
     502
     503  if( a2 != NULL )
     504  {
     505    assume( __LEAD2SYZ__ );
     506
     507    const int r2 = p_GetComp(a2, R) - 1; poly aa2 = leadmonom(a2, R); // :(
     508
     509    assume( r2 >= 0 && r2 < IDELEMS(T) );
     510
     511    t = p_Add_q(a2, p_Add_q(t, TraverseTail(aa2, r2), R), R);
     512
     513    p_Delete(&aa2, R);
     514  } else
     515    t = p_Add_q(t, ReduceTerm(aa, L->m[r], a), R);
     516
     517  p_Delete(&aa, R);
     518
     519  return t;
     520}
     521
     522
     523
    489524void SchreyerSyzygyComputation::ComputeSyzygy()
    490525{
    491 //  FROM_NAMESPACE(INTERNAL, _ComputeSyzygy(m_idLeads, m_idTails, m_syzLeads, m_syzTails, m_rBaseRing, m_atttributes)); // TODO: just a wrapper for now :/
    492 
    493526  assume( m_idLeads != NULL );
    494527  assume( m_idTails != NULL );
    … …  
    499532  ideal& TT = m_syzTails;
    500533  const ring& R = m_rBaseRing;
    501 //  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
    502 
    503 //  const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
    504 //  const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
    505 //  const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
    506 //  const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
    507 //  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
    508 
    509   assume( R == currRing ); // For attributes :-/
    510534
    511535  assume( IDELEMS(L) == IDELEMS(T) );
    512536
    513537  if( m_syzLeads == NULL )
    514     ComputeLeadingSyzygyTerms( __LEAD2SYZ__ ); // 2 terms OR 1 term!
     538    ComputeLeadingSyzygyTerms( __LEAD2SYZ__ && !__IGNORETAILS__ ); // 2 terms OR 1 term!
    515539
    516540  assume( m_syzLeads != NULL );
    … …  
    541565      TT->m[k] = NULL;
    542566
     567      assume( a2 != NULL );
     568
    543569      if( a2 != NULL )
    544570        p_Delete(&a2, R);
    … …  
    547573    }
    548574   
    549     TT->m[k] = a2;
    550 
    551     const int r = p_GetComp(a, R) - 1;
    552 
    553     assume( r >= 0 && r < IDELEMS(T) );
    554     assume( r >= 0 && r < IDELEMS(L) );         
    555 
    556     poly aa = leadmonom(a, R); assume( aa != NULL); // :(
     575//    TT->m[k] = a2;
    557576
    558577    if( ! __HYBRIDNF__ )
    559578    {
    560       poly t = TraverseTail(aa, T->m[r]);
    561 
    562       if( a2 != NULL )
    563       {
    564         assume( __LEAD2SYZ__ );
    565 
    566         const int r2 = p_GetComp(a2, R) - 1; poly aa2 = leadmonom(a2, R); // :(
    567 
    568         assume( r2 >= 0 && r2 < IDELEMS(T) );
    569 
    570         TT->m[k] = p_Add_q(a2, p_Add_q(t, TraverseTail(aa2, T->m[r2]), R), R);
    571 
    572         p_Delete(&aa2, R);
    573       } else
    574         TT->m[k] = p_Add_q(t, ReduceTerm(aa, L->m[r], a), R);
    575 
     579      TT->m[k] = TraverseNF(a, a2);
    576580    } else
    577581    {
    578       if( a2 == NULL )
    579       {
    580         aa = p_Mult_mm(aa, L->m[r], R); a2 = m_div.FindReducer(aa, a, m_checker);
    581       }
    582       assume( a2 != NULL );
    583 
    584       TT->m[k] = SchreyerSyzygyNF(a, a2); // will copy a2 :(
    585 
    586       p_Delete(&a2, R);
     582      TT->m[k] = SchreyerSyzygyNF(a, a2);
    587583    }
    588584
    589     p_Delete(&aa, R);   
    590585  }
    591586
    … …  
    638633}
    639634
    640 poly SchreyerSyzygyComputation::SchreyerSyzygyNF(poly syz_lead, poly syz_2) const
    641 {
    642 // return FROM_NAMESPACE(INTERNAL, _SchreyerSyzygyNF(syz_lead, syz_2, m_idLeads, m_idTails, m_LS, m_rBaseRing, m_atttributes));
    643 // poly _SchreyerSyzygyNF(poly syz_lead, poly syz_2,
    644 //                       ideal L, ideal T, ideal LS,
    645 //                       const ring r,
    646 //                       const SchreyerSyzygyComputationFlags attributes)
    647 // {
    648 
     635#define NOPRODUCT 1
     636
     637poly SchreyerSyzygyComputation::SchreyerSyzygyNF(const poly syz_lead, poly syz_2) const
     638{
     639 
    649640  assume( !__IGNORETAILS__ );
    650641
    … …  
    653644  const ring& r = m_rBaseRing;
    654645
    655 //   const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
    656 //   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
    657 //   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
    658 //   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
    659 //   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
    660 //   const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
    661 
    662646  assume( syz_lead != NULL );
     647
     648  if( syz_2 == NULL )
     649  {
     650    const int rr = p_GetComp(syz_lead, r) - 1;
     651
     652    assume( rr >= 0 && rr < IDELEMS(T) );
     653    assume( rr >= 0 && rr < IDELEMS(L) );
     654
     655
     656#if NOPRODUCT
     657    syz_2 = m_div.FindReducer(syz_lead, L->m[rr], syz_lead, m_checker);
     658#else   
     659    poly aa = leadmonom(syz_lead, r); assume( aa != NULL); // :(
     660    aa = p_Mult_mm(aa, L->m[rr], r);
     661
     662    syz_2 = m_div.FindReducer(aa, syz_lead, m_checker);
     663
     664    p_Delete(&aa, r);
     665#endif
     666   
     667    assume( syz_2 != NULL ); // by construction of S-Polynomial   
     668  }
     669
     670
     671 
    663672  assume( syz_2 != NULL );
    664673
    … …  
    684693  p_Delete(&p, r);
    685694
    686   poly tail = p_Copy(syz_2, r); // TODO: use bucket!?
     695  poly tail = syz_2; // TODO: use bucket!?
    687696
    688697  while (spoly != NULL)
    … …  
    710719}
    711720
     721poly SchreyerSyzygyComputation::TraverseTail(poly multiplier, const int tail) const
     722{
     723  // TODO: store (multiplier, tail) -.-^-.-^-.--> !
     724  assume(m_idTails !=  NULL && m_idTails->m != NULL);
     725  assume( tail >= 0 && tail < IDELEMS(m_idTails) );
     726
     727  const poly t = m_idTails->m[tail];
     728
     729  if(t != NULL)
     730    return TraverseTail(multiplier, t);
     731
     732  return NULL;
     733}
     734
    712735
    713736poly SchreyerSyzygyComputation::TraverseTail(poly multiplier, poly tail) const
    … …  
    717740  const ideal& L = m_idLeads;
    718741  const ideal& T = m_idTails;
    719 //  const ideal& LS = m_LS;
    720742  const ring& r = m_rBaseRing;
    721 //  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
    722 
    723 //  return FROM_NAMESPACE(INTERNAL, _TraverseTail(multiplier, tail, m_idLeads, m_idTails, m_LS, m_rBaseRing, m_atttributes));
    724 // poly _TraverseTail(poly multiplier, poly tail,
    725 //                    ideal L, ideal T, ideal LS,
    726 //                    const ring r,
    727 //                    const SchreyerSyzygyComputationFlags attributes)
    728 // {
    729 
    730 
    731 //   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
    732 //   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
    733 //   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
    734 //   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
    735 //   const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
    736743
    737744  assume( multiplier != NULL );
    … …  
    758765  const ideal& L = m_idLeads;
    759766  const ideal& T = m_idTails;
    760 //  const ideal& LS = m_LS;
    761767  const ring& r = m_rBaseRing;
    762 //  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
    763 
    764 //  return FROM_NAMESPACE(INTERNAL, _ReduceTerm(multiplier, term4reduction, syztermCheck, m_idLeads, m_idTails, m_LS, m_rBaseRing, m_atttributes));
    765 // poly _ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck,
    766 //                 ideal L, ideal T, ideal LS,
    767 //                 const ring r,
    768 //                 const SchreyerSyzygyComputationFlags attributes)
    769 
    770 
    771 
    772 //   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
    773 //   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
    774 //   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
    775 //   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
    776 //  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
    777768
    778769  assume( multiplier != NULL );
    … …  
    783774  assume( T != NULL );
    784775
    785   // assume(r == currRing); // ?
    786 
    787776  // simple implementation with FindReducer:
    788777  poly s = NULL;
    … …  
    790779  if( (!__TAILREDSYZ__) || m_lcm.Check(multiplier) )
    791780  {
     781#if NOPRODUCT
     782    s = m_div.FindReducer(multiplier, term4reduction, syztermCheck, m_checker);
     783#else   
    792784    // NOTE: only LT(term4reduction) should be used in the following:
    793785    poly product = pp_Mult_mm(multiplier, term4reduction, r);
    794786    s = m_div.FindReducer(product, syztermCheck, m_checker);
    795787    p_Delete(&product, r);
     788#endif
    796789  }
    797790
    … …  
    804797  assume( c >= 0 && c < IDELEMS(T) );
    805798
    806   const poly tail = T->m[c];
    807 
    808   if( tail != NULL )
    809     s = p_Add_q(s, TraverseTail(b, tail), r); 
     799  const poly t = TraverseTail(b, c); // T->m[c];
     800
     801  if( t != NULL )
     802    s = p_Add_q(s, t, r); 
    810803
    811804  return s;
    … …  
    988981
    989982
    990 poly CReducerFinder::FindReducer(const poly product, const poly syzterm, const CReducerFinder& syz_checker) const
    991 {
     983/// _p_LmDivisibleByNoComp for a | b*c
     984static inline BOOLEAN _p_LmDivisibleByNoComp(const poly a, const poly b, const poly c, const ring r)
     985{
     986  int i=r->VarL_Size - 1;
     987  unsigned long divmask = r->divmask;
     988  unsigned long la, lb;
     989
     990  if (r->VarL_LowIndex >= 0)
     991  {
     992    i += r->VarL_LowIndex;
     993    do
     994    {
     995      la = a->exp[i];
     996      lb = b->exp[i] + c->exp[i];
     997      if ((la > lb) ||
     998          (((la & divmask) ^ (lb & divmask)) != ((lb - la) & divmask)))
     999      {
     1000        pDivAssume(p_DebugLmDivisibleByNoComp(a, b, r) == FALSE);
     1001        return FALSE;
     1002      }
     1003      i--;
     1004    }
     1005    while (i>=r->VarL_LowIndex);
     1006  }
     1007  else
     1008  {
     1009    do
     1010    {
     1011      la = a->exp[r->VarL_Offset[i]];
     1012      lb = b->exp[r->VarL_Offset[i]] + c->exp[r->VarL_Offset[i]];
     1013      if ((la > lb) ||
     1014          (((la & divmask) ^ (lb & divmask)) != ((lb - la) & divmask)))
     1015      {
     1016        pDivAssume(p_DebugLmDivisibleByNoComp(a, b, r) == FALSE);
     1017        return FALSE;
     1018      }
     1019      i--;
     1020    }
     1021    while (i>=0);
     1022  }
     1023#ifdef HAVE_RINGS
     1024  assume( !rField_is_Ring(r) ); // not implemented for rings...!
     1025#endif
     1026  return TRUE;
     1027}
     1028
     1029
     1030poly CReducerFinder::FindReducer(const poly multiplier, const poly t,
     1031                                 const poly syzterm, const CReducerFinder& syz_checker) const
     1032{
     1033  // don't case about the module component of multiplier (as it may be
     1034  // the syzygy term)
     1035  // product = multiplier * t?
    9921036  const ring& r = m_rBaseRing;
    9931037
    994   assume( product != NULL );
     1038  assume( multiplier != NULL ); assume( t != NULL );
    9951039
    9961040  const ideal& L = m_L; assume( L != NULL ); // for debug/testing only!
    … …  
    10101054
    10111055    poly lm = p_Mult_mm(leadmonom(syzterm, r), m, r);
    1012     assume( p_EqualPolys(lm, product, r) );
     1056
     1057    poly pr = p_Mult_q( p_LmInit(multiplier, r), p_LmInit(t, r), r);
     1058   
     1059    assume( p_EqualPolys(lm, pr, r) );
    10131060
    10141061    //  def @@c = leadcomp(syzterm); int @@r = int(@@c);
    … …  
    10161063
    10171064    p_Delete(&lm, r);   
    1018   }
    1019 
    1020   const long comp = p_GetComp(product, r);
    1021   const unsigned long not_sev = ~p_GetShortExpVector(product, r);
     1065    p_Delete(&pr, r);   
     1066  }
     1067
     1068  const long comp = p_GetComp(t, r);
     1069 
     1070  const unsigned long not_sev = ~p_GetShortExpVector(multiplier, t, r); // !
    10221071
    10231072  assume( comp >= 0 );
    … …  
    10311080//
    10321081//     const unsigned long p_sev = p_GetShortExpVector(p, r); // to be stored in m_hash!!!
    1033  
     1082
    10341083   // looking for an appropriate diviser p = L[k]...
    10351084  CReducersHash::const_iterator it = m_hash.find(comp); // same module component
    … …  
    10481097  if( __DEBUG__ )
    10491098    p_SetCoeff0(q, 0, r); // for printing q
    1050  
     1099
    10511100  for(TReducers::const_iterator vit = reducers.begin(); vit != reducers.end(); vit++ )
    10521101  {
    … …  
    10631112    assume( p_sev == p_GetShortExpVector(p, r) );     
    10641113
    1065     if( !p_LmShortDivisibleByNoComp(p, p_sev, product, not_sev, r) )
     1114//    if( !p_LmShortDivisibleByNoComp(p, p_sev, product, not_sev, r) )
     1115//      continue;     
     1116
     1117    if (p_sev & not_sev)
     1118      continue;
     1119
     1120    if( !_p_LmDivisibleByNoComp(p, multiplier, t, r) )
    10661121      continue;     
    10671122
    1068 //     // ... which divides the product, looking for the _1st_ appropriate one!
    1069 //     if( !p_LmDivisibleByNoComp(p, product, r) ) // included inside  p_LmShortDivisibleBy!
    1070 //       continue;
    1071 
    1072     p_ExpVectorDiff(q, product, p, r); // (LM(product) / LM(L[k]))
     1123
     1124    p_ExpVectorSum(q, multiplier, t, r); // q == product == multiplier * t       
     1125    p_ExpVectorDiff(q, q, p, r); // (LM(product) / LM(L[k]))
     1126   
    10731127    p_SetComp(q, k + 1, r);
    10741128    p_Setm(q, r);
    … …  
    10941148        PrintS("_FindReducer::Test LS: q is divisible by LS[?] !:((: ");
    10951149      }
     1150
     1151      continue;
     1152    }
     1153
     1154    number n = n_Mult( p_GetCoeff(multiplier, r), p_GetCoeff(t, r), r);
     1155    p_SetCoeff0(q, n_Neg( n_Div(n, p_GetCoeff(p, r), r), r), r);
     1156    n_Delete(&n, r);
     1157   
     1158    return q;
     1159  }
     1160
     1161  p_LmFree(q, r);
     1162
     1163  return NULL;
     1164}
     1165
     1166
     1167poly CReducerFinder::FindReducer(const poly product, const poly syzterm, const CReducerFinder& syz_checker) const
     1168{
     1169  const ring& r = m_rBaseRing;
     1170
     1171  assume( product != NULL );
     1172
     1173  const ideal& L = m_L; assume( L != NULL ); // for debug/testing only!
     1174
     1175  long c = 0;
     1176
     1177  if (syzterm != NULL)
     1178    c = p_GetComp(syzterm, r) - 1;
     1179
     1180  assume( c >= 0 && c < IDELEMS(L) );
     1181
     1182  if (__DEBUG__ && (syzterm != NULL))
     1183  {
     1184    const poly m = L->m[c];
     1185
     1186    assume( m != NULL ); assume( pNext(m) == NULL );
     1187
     1188    poly lm = p_Mult_mm(leadmonom(syzterm, r), m, r);
     1189    assume( p_EqualPolys(lm, product, r) );
     1190
     1191    //  def @@c = leadcomp(syzterm); int @@r = int(@@c);
     1192    //  def @@product = leadmonomial(syzterm) * L[@@r];
     1193
     1194    p_Delete(&lm, r);   
     1195  }
     1196
     1197  const long comp = p_GetComp(product, r);
     1198  const unsigned long not_sev = ~p_GetShortExpVector(product, r);
     1199
     1200  assume( comp >= 0 );
     1201
     1202//   for( int k = IDELEMS(L)-1; k>= 0; k-- )
     1203//   {
     1204//     const poly p = L->m[k];
     1205//
     1206//     if ( p_GetComp(p, r) != comp )
     1207//       continue;
     1208//
     1209//     const unsigned long p_sev = p_GetShortExpVector(p, r); // to be stored in m_hash!!!
     1210 
     1211   // looking for an appropriate diviser p = L[k]...
     1212  CReducersHash::const_iterator it = m_hash.find(comp); // same module component
     1213
     1214  if( it == m_hash.end() )
     1215    return NULL;
     1216
     1217  assume( m_L != NULL );
     1218
     1219  const TReducers& reducers = it->second;
     1220
     1221  const BOOLEAN to_check = (syz_checker.IsNonempty()); // __TAILREDSYZ__ &&
     1222
     1223  const poly q = p_New(r); pNext(q) = NULL;
     1224
     1225  if( __DEBUG__ )
     1226    p_SetCoeff0(q, 0, r); // for printing q
     1227 
     1228  for(TReducers::const_iterator vit = reducers.begin(); vit != reducers.end(); vit++ )
     1229  {
     1230    const poly p = (*vit)->m_lt;
     1231
     1232    assume( p_GetComp(p, r) == comp );
     1233
     1234    const int k = (*vit)->m_label;
     1235
     1236    assume( L->m[k] == p );
     1237
     1238    const unsigned long p_sev = (*vit)->m_sev;
     1239
     1240    assume( p_sev == p_GetShortExpVector(p, r) );     
     1241
     1242    if( !p_LmShortDivisibleByNoComp(p, p_sev, product, not_sev, r) )
     1243      continue;     
     1244
     1245//     // ... which divides the product, looking for the _1st_ appropriate one!
     1246//     if( !p_LmDivisibleByNoComp(p, product, r) ) // included inside  p_LmShortDivisibleBy!
     1247//       continue;
     1248
     1249    p_ExpVectorDiff(q, product, p, r); // (LM(product) / LM(L[k]))
     1250    p_SetComp(q, k + 1, r);
     1251    p_Setm(q, r);
     1252
     1253    // cannot allow something like: a*gen(i) - a*gen(i)
     1254    if (syzterm != NULL && (k == c))
     1255      if (p_ExpVectorEqual(syzterm, q, r))
     1256      {
     1257        if( __DEBUG__ )
     1258        {
     1259          Print("_FindReducer::Test SYZTERM: q == syzterm !:((, syzterm is: ");
     1260          dPrint(syzterm, r, r, 1);
     1261        }   
     1262
     1263        continue;
     1264      }
     1265
     1266    // while the complement (the fraction) is not reducible by leading syzygies
     1267    if( to_check && syz_checker.IsDivisible(q) )
     1268    {
     1269      if( __DEBUG__ )
     1270      {
     1271        PrintS("_FindReducer::Test LS: q is divisible by LS[?] !:((: ");
     1272      }
    10961273     
    10971274      continue;
Note: See TracChangeset for help on using the changeset viewer.