Changeset 8f761ae in git


Ignore:
Timestamp:
Aug 28, 2013, 4:56:50 PM (11 years ago)
Author:
Oleksandr Motsak <motsak@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
79502e381504b0ec8fe16939b38060226ecebf6f
Parents:
be561246e7064a6d59d9b79c1eff7305b7d7aa00
git-author:
Oleksandr Motsak <motsak@mathematik.uni-kl.de>2013-08-28 16:56:50+02:00
git-committer:
Oleksandr Motsak <motsak@mathematik.uni-kl.de>2013-08-28 17:05:28+02:00
Message:
Update of MATHICGB according to the PullRequest #351
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • Singular/singmathic.cc

    • Property mode changed from 100755 to 100644
    rbe56124 r8f761ae  
    2020#include <mathicgb.h>
    2121
    22 typedef mgb::GroebnerConfiguration::Coefficient Coefficient;
    23 typedef mgb::GroebnerConfiguration::VarIndex VarIndex;
    24 typedef mgb::GroebnerConfiguration::Exponent Exponent;
    25 typedef mgb::GroebnerConfiguration::BaseOrder BaseOrder;
    26 
    2722// Constructs a Singular ideal.
    2823class MathicToSingStream {
    2924public:
     25  typedef mgb::GroebnerConfiguration::Coefficient Coefficient;
     26  typedef mgb::GroebnerConfiguration::VarIndex VarIndex;
     27  typedef mgb::GroebnerConfiguration::Exponent Exponent;
     28
    3029  MathicToSingStream(Coefficient modulus, VarIndex varCount):
    3130    mModulus(modulus),
    3231    mVarCount(varCount),
    33     mPolyCount(0),
     32    polyCount(0),
    3433    mTerm(0),
    3534    mIdeal(0)
     
    4645    deleteIdeal();
    4746    mIdeal = idInit(polyCount);
    48     mPolyCount = 0;
     47    polyCount = 0;
    4948  }
    5049
     
    5352  void appendTermBegin() {
    5453    if (mTerm == 0)
    55       mTerm = mIdeal->m[mPolyCount] = pInit();
     54      mTerm = mIdeal->m[polyCount] = pInit();
    5655    else
    5756      mTerm = mTerm->next = pInit();
     
    6867
    6968  void appendPolynomialDone() {
    70     ++mPolyCount;
     69    ++polyCount;
    7170    mTerm = 0;
    7271  }
     
    9392  const Coefficient mModulus;
    9493  const VarIndex mVarCount;
    95   size_t mPolyCount;
     94  size_t polyCount;
    9695  poly mTerm;
    9796  ::ideal mIdeal;
    9897};
    9998
    100 #include <iostream>
    101 
    102 bool setOrder(ring r, mgb::GroebnerConfiguration& conf) {
    103   const VarIndex varCount = conf.varCount();
    104 
    105   bool didSetComponentBefore = false;
    106   mgb::GroebnerConfiguration::BaseOrder baseOrder =
    107     mgb::GroebnerConfiguration::RevLexDescendingBaseOrder;
    108 
    109   std::vector<Exponent> gradings;
    110   for (int block = 0; r->order[block] != ringorder_no; ++block) {
    111     // *** ringorder_no
    112 
    113     const rRingOrder_t type = static_cast<rRingOrder_t>(r->order[block]);
    114     if (r->block0[block] < 0 || r->block1[block] < 0) {
    115       WerrorS("Unexpected negative block0/block1 in ring.");
    116       return false;
    117     }
    118     const VarIndex block0 = static_cast<VarIndex>(r->block0[block]);
    119     const VarIndex block1 = static_cast<VarIndex>(r->block1[block]);
    120     const int* const weights = r->wvhdl[block];
    121     if (block0 > block1) {
    122       WerrorS("Unexpected block0 > block1 in ring.");
    123       return false;
    124     }
    125 
    126     // *** ringorder_c and ringorder_C
    127     if (type == ringorder_c || type == ringorder_C) {
    128       if (block0 != 0 || block1 != 0 || weights != 0) {
    129         WerrorS("Unexpected non-zero fields on c/C block in ring.");
    130         return false;
    131       }
    132       if (didSetComponentBefore) {
    133         WerrorS("Unexpected two c/C blocks in ring.");
    134         return false;
    135       }
    136       didSetComponentBefore = true;
    137       if (r->order[block + 1] == ringorder_no) {
    138         conf.setComponentBefore
    139           (mgb::GroebnerConfiguration::ComponentAfterBaseOrder);
    140       } else
    141         conf.setComponentBefore(gradings.size() / varCount);
    142       conf.setComponentsAscending(type == ringorder_C);
    143       continue;
    144     }
    145     if (block0 == 0 || block1 == 0) {
    146       WerrorS("Expected block0 != 0 and block1 != 0 in ring.");
    147       return false;
    148     }
    149     if (block1 > varCount) {
    150       // todo: first handle any block types where this is not true
    151       WerrorS("Expected block1 <= #vars in ring.");
    152       return false;
    153     }
    154 
    155     // dim is how many variables this block concerns.
    156     const size_t dim = static_cast<size_t>(block1 - block0 + 1);
    157 
    158     // *** single-graded/ungraded lex/revlex orders
    159     // a(w): w-graded and that's it
    160     // a64(w): w-graded with 64-bit weights (not supported here)
    161     //    lp:               lex from  left (descending)
    162     //    Dp:  1-graded,    lex from  left (descending)
    163     //    Ds: -1-graded,    lex from  left (descending)
    164     // Wp(w):  w-graded,    lex from  left (descending)
    165     // Ws(w): -w-graded,    lex from  left (descending)
    166     //    rp:               lex from right (ascending)
    167     //    rs:            revlex from right (descending)
    168     //    dp:  1-graded, revlex from right (descending)
    169     //    ds: -1-graded, revlex from right (descending)
    170     // wp(w):  w-graded, revlex from right (descending)
    171     // ws(w): -w-graded, revlex from right (descending)
    172     //    ls:            revlex from  left (ascending)
    173 
    174     if (type == ringorder_a64) {
    175       WerrorS("Block type a64 not supported for MathicGB interface.");
    176       return false;
    177     }
    178 
    179     // * handle the single-grading part
    180     const bool oneGrading = (type == ringorder_Dp || type == ringorder_dp);
    181     const bool minusOneGrading = (type == ringorder_Ds || type == ringorder_ds);
    182     const bool wGrading =
    183       (type == ringorder_a || type == ringorder_Wp || type == ringorder_wp);
    184     const bool minusWGrading = (type == ringorder_ws || type == ringorder_Ws);
    185     if (oneGrading || minusOneGrading || wGrading || minusWGrading) {
    186       const VarIndex begin = gradings.size();
    187       gradings.resize(begin + varCount);
    188       if (oneGrading || minusOneGrading) {
    189         if (weights != 0) {
    190           WerrorS("Expect wvhdl == 0 in Dp/dp/Ds/ds-block in ring.");
    191           return false;
    192         }
    193         const Exponent value = oneGrading ? 1 : -1;
    194         for (int var = block0 - 1; var < block1; ++var)
    195           gradings[begin + var] = value;
    196       } else {
    197         if (weights == 0) {
    198           WerrorS("Expect wvhdl != 0 in a/Wp/wp/ws/Ws-block in ring.");
    199           return false;
    200         }
    201         if (wGrading) {
    202           for (int var = 0; var < dim; ++var)
    203             gradings[begin + (block0 - 1) + var] = weights[var];
    204         } else {
    205           for (int var = 0; var < dim; ++var)
    206             gradings[begin + (block0 - 1) + var] = -weights[var];
    207         }
    208       }
    209     }
    210     if (type == ringorder_a)
    211       continue; // a has only the grading, so we are done already
    212 
    213     // * handle the lex/revlex part
    214     const bool lexFromLeft =
    215       type == ringorder_lp ||
    216       type == ringorder_Dp ||
    217       type == ringorder_Ds ||
    218       type == ringorder_Wp ||
    219       type == ringorder_Ws;
    220     const bool lexFromRight = type == ringorder_rp;
    221     const bool revlexFromLeft = type == ringorder_ls;
    222     const bool revlexFromRight =
    223       type == ringorder_rs ||
    224       type == ringorder_dp ||
    225       type == ringorder_ds ||
    226       type == ringorder_wp ||
    227       type == ringorder_ws;
    228     if (lexFromLeft || lexFromRight || revlexFromLeft || revlexFromRight) {
    229       const int next = r->order[block + 1];
    230       bool final = next == ringorder_no;
    231       if (!final && r->order[block + 2] == ringorder_no)
    232         final = next == ringorder_c || next == ringorder_C;
    233       if (final) {
    234         if (lexFromRight)
    235           baseOrder = mgb::GroebnerConfiguration::LexAscendingBaseOrder;
    236         else if (revlexFromRight)
    237           baseOrder = mgb::GroebnerConfiguration::RevLexDescendingBaseOrder;
    238         else if (lexFromLeft)
    239           baseOrder = mgb::GroebnerConfiguration::LexDescendingBaseOrder;
    240         else
    241           baseOrder = mgb::GroebnerConfiguration::RevLexAscendingBaseOrder;
    242         continue;
    243       }
    244 
    245       const size_t begin = gradings.size();
    246       gradings.resize(begin + dim * varCount);
    247       const Exponent value = (lexFromLeft || lexFromRight) ? 1 : -1;
    248       if (lexFromLeft || revlexFromLeft) {
    249         for (size_t row = 0; row < dim; ++row)
    250           gradings[begin + row * varCount + (block0 - 1) + row] = value;
    251       } else {
    252         for (size_t row = 0; row < dim; ++row)
    253           gradings[begin + row * varCount + (block1 - 1) - row] = value;
    254       }
    255       continue;
    256     }
    257 
    258     // *** ringorder_M: a square invertible matrix
    259     if (type == ringorder_M) {
    260       if (weights == 0) {
    261         WerrorS("Expected wvhdl != 0 in M-block in ring.");
    262         return false;
    263       }
    264       const size_t begin = gradings.size();
    265       gradings.resize(begin + dim * varCount);
    266       for (size_t row = 0; row < dim; ++row)
    267         for (size_t col = block0 - 1; col < block1; ++col)
    268           gradings[begin + row * varCount + col] = weights[row * dim + col];
    269       continue;
    270     }
    271 
    272     // *** Miscellaneous unsupported or invalid block types
    273     if (
    274       type == ringorder_s ||
    275       type == ringorder_S ||
    276       type == ringorder_IS
    277     ) {
    278       // todo: Consider supporting this later.
    279       WerrorS("Schreyer order s/S/IS not supported in MathicGB interface.");
    280       return false;
    281     }
    282     if (type == ringorder_am) {
    283       // This block is a Schreyer-like ordering only used in Spielwiese.
    284       // todo: Consider supporting it later.
    285       WerrorS("Block type am not supported in MathicGB interface");
    286       return false;
    287     }
    288     if (type == ringorder_L) {
    289       WerrorS("Invalid L-block found in order of ring.");
    290       return false;
    291     }
    292     if (type == ringorder_aa) {
    293       // I don't know what an aa block is supposed to do.
    294       WerrorS("aa ordering not supported by the MathicGB interface.");
    295       return false;
    296     }
    297     if (type == ringorder_unspec) {
    298       WerrorS("Invalid unspec-block found in order of ring.");
    299       return false;
    300     }
    301     WerrorS("Unknown block type found in order of ring.");
    302     return false;
    303   }
    304 
    305   if (!didSetComponentBefore) {
    306     WerrorS("Expected to find a c/C block in ring.");
    307     return false;
    308   }
    309 
    310   if (!conf.setMonomialOrder(baseOrder, gradings)) {
    311     WerrorS("MathicGB does not support non-global orders.");
    312     return false;
    313   }
    314   return true;
    315 }
    316 
    317 bool prOrderMatrix(ring r) {
    318   const int varCount = r->N;
    319   mgb::GroebnerConfiguration conf(101, varCount);
    320   if (!setOrder(r, conf))
    321     return false;
    322   const std::vector<Exponent>& gradings = conf.monomialOrder().second;
    323   if (gradings.size() % varCount != 0) {
    324     WerrorS("Expected matrix to be a multiple of varCount.");
    325     return false;
    326   }
    327   const size_t rowCount = gradings.size() / varCount;
    328   std::cout << "Order matrix:\n";
    329   for (size_t row = 0; row < rowCount; ++row) {
    330     for (size_t col = 0; col < varCount; ++col)
    331       std::cerr << ' ' << gradings[row * varCount + col];
    332     std::cerr << '\n';
    333   }
    334   std::cerr
    335     << "Base order: "
    336     << mgb::GroebnerConfiguration::baseOrderName(conf.monomialOrder().first)
    337     << '\n';
    338   std::cerr << "Component before: " << conf.componentBefore() << '\n';
    339   std::cerr << "Components ascending: " << conf.componentsAscending() << '\n';
    340   std::cerr << "Schreyering: " << conf.schreyering() << '\n';
    341 }
    342 
    343 void prOrder(ring r) {
    344   std::cout << "Printing order of ring.\n";
    345   for (int block = 0; ; ++block) {
    346     switch (r->order[block]) {
    347     case ringorder_no: // end of blocks
    348       return;
    349 
    350     case ringorder_a:
    351       std::cout << "a";
    352       break;
    353 
    354     case ringorder_a64: ///< for int64 weights
    355       std::cout << "a64";
    356       break;
    357 
    358     case ringorder_c:
    359       std::cout << "c";
    360       break;
    361 
    362     case ringorder_C:
    363       std::cout << "C";
    364       break;
    365 
    366     case ringorder_M:
    367       std::cout << "M";
    368       break;
    369      
    370     case ringorder_S: ///< S?
    371       std::cout << "S";
    372       break;
    373      
    374     case ringorder_s: ///< s?
    375       std::cout << "s";
    376       break;
    377 
    378     case ringorder_lp:
    379       std::cout << "lp";
    380       break;
    381      
    382     case ringorder_dp:
    383       std::cout << "dp";
    384       break;
    385 
    386     case ringorder_rp:
    387       std::cout << "rp";
    388       break;
    389 
    390     case ringorder_Dp:
    391       std::cout << "Dp";
    392       break;
    393 
    394     case ringorder_wp:
    395       std::cout << "wp";
    396       break;
    397 
    398     case ringorder_Wp:
    399       std::cout << "Wp";
    400       break;
    401 
    402     case ringorder_ls:
    403       std::cout << "ls"; // not global
    404       break;
    405 
    406     case ringorder_ds:
    407       std::cout << "ds"; // not global
    408       break;
    409 
    410     case ringorder_Ds:
    411       std::cout << "Ds"; // not global
    412       break;
    413 
    414     case ringorder_ws:
    415       std::cout << "ws"; // not global
    416       break;
    417 
    418     case ringorder_Ws:
    419       std::cout << "Ws"; // not global
    420       break;
    421 
    422     case ringorder_am:
    423       std::cout << "am";
    424       break;
    425 
    426     case ringorder_L:
    427       std::cout << "L";
    428       break;
    429 
    430     // the following are only used internally
    431     case ringorder_aa: ///< for idElimination, like a, except pFDeg, pWeigths ignore it
    432       std::cout << "aa";
    433       break;
    434 
    435     case ringorder_rs: ///< opposite of ls
    436       std::cout << "rs";
    437       break;
    438 
    439     case ringorder_IS: ///< Induced (Schreyer) ordering
    440       std::cout << "IS";
    441       break;
    442 
    443     case ringorder_unspec:
    444       std::cout << "unspec";
    445       break;
    446     }
    447     const int b0 = r->block0[block];
    448     const int b1 = r->block1[block];
    449     std::cout << ' ' << b0 << ':' << b1 << " (" << r->wvhdl[block] << ")" << std::flush;
    450     if (r->wvhdl[block] != 0 && b0 != 0) {
    451       for (int v = 0; v <= b1 - b0; ++v)
    452         std::cout << ' ' << r->wvhdl[block][v];
    453     } else
    454       std::cout << " null";
    455     std::cout << '\n';
    456   }
    457 }
    458 
    459 BOOLEAN prOrderX(leftv result, leftv arg) {
    460   if (currRing == 0) {
    461     WerrorS("There is no current ring.");
    462     return TRUE;
    463   }
    464   prOrder(currRing);
    465   prOrderMatrix(currRing);
    466   result->rtyp=NONE;
    467   return FALSE;
    468 }
    469 
    470 BOOLEAN setRingGlobal(leftv result, leftv arg) {
    471   currRing->OrdSgn = 1;
    472   result->rtyp=NONE;
    473   return FALSE;
    474 }
    475 
    47699BOOLEAN mathicgb(leftv result, leftv arg)
    477100{
    478   result->rtyp=NONE;
    479 
    480101  if (arg == NULL || arg->next != NULL || arg->Typ() != IDEAL_CMD) {
    481102    WerrorS("Syntax: mathicgb(<ideal>)");
     
    490111  const int varCount = currRing->N;
    491112  mgb::GroebnerConfiguration conf(characteristic, varCount);
    492   if (!setOrder(currRing, conf))
    493     return TRUE;
    494   if (TEST_OPT_PROT)
    495     conf.setLogging("all");
    496 
    497113  mgb::GroebnerInputIdealStream toMathic(conf);
    498114
     
    511127      for (int i = 1; i <= currRing->N; ++i)
    512128        toMathic.appendExponent(i - 1, pGetExp(p, i));
    513       const long coefLong = reinterpret_cast<long>(pGetCoeff(p));
    514       toMathic.appendTermDone(static_cast<int>(coefLong));
     129      toMathic.appendTermDone(reinterpret_cast<int>(pGetCoeff(p)));
    515130    }
    516131    toMathic.appendPolynomialDone();
     
    538153    mathicgb
    539154  );
    540   psModulFunctions->iiAddCproc(
    541     (currPack->libname ? currPack->libname : ""),
    542     "mathicgb_prOrder",
    543     FALSE,
    544     prOrderX
    545   );
    546   psModulFunctions->iiAddCproc(
    547     (currPack->libname ? currPack->libname : ""),
    548     "mathicgb_setRingGlobal",
    549     FALSE,
    550     setRingGlobal
    551   );
    552155  return 1;
    553156}
  • m4/mathic-check.m4

    rbe56124 r8f761ae  
    33#
    44 AC_ARG_WITH(mathicgb,
    5    [
    6     AS_HELP_STRING(
     5   [AS_HELP_STRING(
    76     [--with-mathicgb=yes|no],
    87     [Use the MathicGB library. Default is no.]
    9     )
    10    ],
     8   )],
    119   [],
    1210   [with_mathicgb="no"]
     
    1816 AS_IF( [test "x$with_mathicgb" = xyes],
    1917 [
    20   AC_CHECK_LIB([tbb], [TBB_runtime_interface_version], [], [])
     18  AC_LANG_PUSH([C++])
    2119  AC_CHECK_LIB(memtailor, libmemtailorIsPresent, [],
    2220    [AC_MSG_ERROR([Cannot find libmemtailor, which is required for MathicGB.])])
     
    2422    [AC_MSG_ERROR([Cannot find libmathic, which is required for MathicGB.])])
    2523  AC_CHECK_LIB(mathicgb, libmathicgbIsPresent, [],
    26     [AC_MSG_ERROR([Cannot find the MathicGB library.])], [-lstdc++])
     24    [AC_MSG_ERROR([Cannot find the MathicGB library.])])
    2725  AC_CHECK_HEADER([mathicgb.h])
     26  AC_LANG_POP([C++])
    2827  AC_DEFINE(HAVE_MATHICGB,1,[Define if mathicgb is to be used])
     28#  AC_SUBST(HAVE_MATHICGB_VALUE, 1) 
    2929 ])
    3030 #
Note: See TracChangeset for help on using the changeset viewer.