source: git/Singular/singmathic.cc @ 4676d5

spielwiese
Last change on this file since 4676d5 was 4676d5, checked in by Oleksandr Motsak <motsak@…>, 10 years ago
Moved timer.*, feread.cc and febase.* from /kernel/ to /Singular/ (also corrected the option(prot); handling) NOTE: also removes line breakes and timings in the output due to option (prot), which were unaccounted and undocumented
  • Property mode set to 100644
File size: 15.5 KB
Line 
1// include header files
2
3#include <kernel/mod2.h>
4#include <misc/auxiliary.h>
5
6#include <misc/options.h>
7
8#include <Singular/febase.h>
9#include <kernel/ideals.h>
10#include <kernel/polys.h>
11
12#include <Singular/ipid.h>
13#include <Singular/mod_lib.h>
14
15#ifdef HAVE_MATHICGB
16
17#include <mathicgb.h>
18
19typedef mgb::GroebnerConfiguration::Coefficient Coefficient;
20typedef mgb::GroebnerConfiguration::VarIndex VarIndex;
21typedef mgb::GroebnerConfiguration::Exponent Exponent;
22typedef mgb::GroebnerConfiguration::BaseOrder BaseOrder;
23
24// Constructs a Singular ideal.
25class MathicToSingStream {
26public:
27  MathicToSingStream(Coefficient modulus, VarIndex varCount):
28    mModulus(modulus),
29    mVarCount(varCount),
30    mPolyCount(0),
31    mTerm(0),
32    mIdeal(0)
33  {}
34
35  ~MathicToSingStream() {deleteIdeal();}
36
37  // Mathic stream interface
38
39  Coefficient modulus() const {return mModulus;}
40  VarIndex varCount() const {return mModulus;}
41
42  void idealBegin(size_t polyCount) {
43    deleteIdeal();
44    mIdeal = idInit(polyCount);
45    mPolyCount = 0;
46  }
47
48  void appendPolynomialBegin(size_t termCount) {}
49
50  void appendTermBegin() {
51    if (mTerm == 0)
52      mTerm = mIdeal->m[mPolyCount] = pInit();
53    else
54      mTerm = mTerm->next = pInit();
55  }
56
57  void appendExponent(VarIndex index, Exponent exponent) {
58    pSetExp(mTerm, index + 1, exponent);
59  }
60
61  void appendTermDone(Coefficient coefficient) {
62    mTerm->coef = reinterpret_cast<number>(coefficient);
63    pSetm(mTerm);
64  }
65
66  void appendPolynomialDone() {
67    ++mPolyCount;
68    mTerm = 0;
69  }
70
71  void idealDone() {}
72
73
74  // Singular interface
75
76  ::ideal takeIdeal() {
77    ::ideal id = mIdeal;
78    mIdeal = 0;
79    return id;
80  }
81
82private:
83  void deleteIdeal() {
84    if (mIdeal != 0) {
85      idDelete(&mIdeal);
86      mIdeal = 0;
87    }
88  }
89
90  const Coefficient mModulus;
91  const VarIndex mVarCount;
92  size_t mPolyCount;
93  poly mTerm;
94  ::ideal mIdeal;
95};
96
97#include <iostream>
98
99bool setOrder(ring r, mgb::GroebnerConfiguration& conf) {
100  const VarIndex varCount = conf.varCount();
101
102  bool didSetComponentBefore = false;
103  mgb::GroebnerConfiguration::BaseOrder baseOrder =
104    mgb::GroebnerConfiguration::RevLexDescendingBaseOrder;
105
106  std::vector<Exponent> gradings;
107  for (int block = 0; r->order[block] != ringorder_no; ++block) {
108    // *** ringorder_no
109
110    const rRingOrder_t type = static_cast<rRingOrder_t>(r->order[block]);
111    if (r->block0[block] < 0 || r->block1[block] < 0) {
112      WerrorS("Unexpected negative block0/block1 in ring.");
113      return false;
114    }
115    const VarIndex block0 = static_cast<VarIndex>(r->block0[block]);
116    const VarIndex block1 = static_cast<VarIndex>(r->block1[block]);
117    const int* const weights = r->wvhdl[block];
118    if (block0 > block1) {
119      WerrorS("Unexpected block0 > block1 in ring.");
120      return false;
121    }
122
123    // *** ringorder_c and ringorder_C
124    if (type == ringorder_c || type == ringorder_C) {
125      if (block0 != 0 || block1 != 0 || weights != 0) {
126        WerrorS("Unexpected non-zero fields on c/C block in ring.");
127        return false;
128      }
129      if (didSetComponentBefore) {
130        WerrorS("Unexpected two c/C blocks in ring.");
131        return false;
132      }
133      didSetComponentBefore = true;
134      if (r->order[block + 1] == ringorder_no) {
135        conf.setComponentBefore
136          (mgb::GroebnerConfiguration::ComponentAfterBaseOrder);
137      } else
138        conf.setComponentBefore(gradings.size() / varCount);
139      conf.setComponentsAscending(type == ringorder_C);
140      continue;
141    }
142    if (block0 == 0 || block1 == 0) {
143      WerrorS("Expected block0 != 0 and block1 != 0 in ring.");
144      return false;
145    }
146    if (block1 > varCount) {
147      // todo: first handle any block types where this is not true
148      WerrorS("Expected block1 <= #vars in ring.");
149      return false;
150    }
151
152    // dim is how many variables this block concerns.
153    const size_t dim = static_cast<size_t>(block1 - block0 + 1);
154
155    // *** single-graded/ungraded lex/revlex orders
156    // a(w): w-graded and that's it
157    // a64(w): w-graded with 64-bit weights (not supported here)
158    //    lp:               lex from  left (descending)
159    //    Dp:  1-graded,    lex from  left (descending)
160    //    Ds: -1-graded,    lex from  left (descending)
161    // Wp(w):  w-graded,    lex from  left (descending)
162    // Ws(w): -w-graded,    lex from  left (descending)
163    //    rp:               lex from right (ascending)
164    //    rs:            revlex from right (descending)
165    //    dp:  1-graded, revlex from right (descending)
166    //    ds: -1-graded, revlex from right (descending)
167    // wp(w):  w-graded, revlex from right (descending)
168    // ws(w): -w-graded, revlex from right (descending)
169    //    ls:            revlex from  left (ascending)
170
171    if (type == ringorder_a64) {
172      WerrorS("Block type a64 not supported for MathicGB interface.");
173      return false;
174    }
175
176    // * handle the single-grading part
177    const bool oneGrading = (type == ringorder_Dp || type == ringorder_dp);
178    const bool minusOneGrading = (type == ringorder_Ds || type == ringorder_ds);
179    const bool wGrading =
180      (type == ringorder_a || type == ringorder_Wp || type == ringorder_wp);
181    const bool minusWGrading = (type == ringorder_ws || type == ringorder_Ws);
182    if (oneGrading || minusOneGrading || wGrading || minusWGrading) {
183      const VarIndex begin = gradings.size();
184      gradings.resize(begin + varCount);
185      if (oneGrading || minusOneGrading) {
186        if (weights != 0) {
187          WerrorS("Expect wvhdl == 0 in Dp/dp/Ds/ds-block in ring.");
188          return false;
189        }
190        const Exponent value = oneGrading ? 1 : -1;
191        for (int var = block0 - 1; var < block1; ++var)
192          gradings[begin + var] = value;
193      } else {
194        if (weights == 0) {
195          WerrorS("Expect wvhdl != 0 in a/Wp/wp/ws/Ws-block in ring.");
196          return false;
197        }
198        if (wGrading) {
199          for (int var = 0; var < dim; ++var)
200            gradings[begin + (block0 - 1) + var] = weights[var];
201        } else {
202          for (int var = 0; var < dim; ++var)
203            gradings[begin + (block0 - 1) + var] = -weights[var];
204        }
205      }
206    }
207    if (type == ringorder_a)
208      continue; // a has only the grading, so we are done already
209
210    // * handle the lex/revlex part
211    const bool lexFromLeft =
212      type == ringorder_lp ||
213      type == ringorder_Dp ||
214      type == ringorder_Ds ||
215      type == ringorder_Wp ||
216      type == ringorder_Ws;
217    const bool lexFromRight = type == ringorder_rp;
218    const bool revlexFromLeft = type == ringorder_ls;
219    const bool revlexFromRight =
220      type == ringorder_rs ||
221      type == ringorder_dp ||
222      type == ringorder_ds ||
223      type == ringorder_wp ||
224      type == ringorder_ws;
225    if (lexFromLeft || lexFromRight || revlexFromLeft || revlexFromRight) {
226      const int next = r->order[block + 1];
227      bool final = next == ringorder_no;
228      if (!final && r->order[block + 2] == ringorder_no)
229        final = next == ringorder_c || next == ringorder_C;
230      if (final) {
231        if (lexFromRight)
232          baseOrder = mgb::GroebnerConfiguration::LexAscendingBaseOrder;
233        else if (revlexFromRight)
234          baseOrder = mgb::GroebnerConfiguration::RevLexDescendingBaseOrder;
235        else if (lexFromLeft)
236          baseOrder = mgb::GroebnerConfiguration::LexDescendingBaseOrder;
237        else
238          baseOrder = mgb::GroebnerConfiguration::RevLexAscendingBaseOrder;
239        continue;
240      }
241
242      const size_t begin = gradings.size();
243      gradings.resize(begin + dim * varCount);
244      const Exponent value = (lexFromLeft || lexFromRight) ? 1 : -1;
245      if (lexFromLeft || revlexFromLeft) {
246        for (size_t row = 0; row < dim; ++row)
247          gradings[begin + row * varCount + (block0 - 1) + row] = value;
248      } else {
249        for (size_t row = 0; row < dim; ++row)
250          gradings[begin + row * varCount + (block1 - 1) - row] = value;
251      }
252      continue;
253    }
254
255    // *** ringorder_M: a square invertible matrix
256    if (type == ringorder_M) {
257      if (weights == 0) {
258        WerrorS("Expected wvhdl != 0 in M-block in ring.");
259        return false;
260      }
261      const size_t begin = gradings.size();
262      gradings.resize(begin + dim * varCount);
263      for (size_t row = 0; row < dim; ++row)
264        for (size_t col = block0 - 1; col < block1; ++col)
265          gradings[begin + row * varCount + col] = weights[row * dim + col];
266      continue;
267    }
268
269    // *** Miscellaneous unsupported or invalid block types
270    if (
271      type == ringorder_s ||
272      type == ringorder_S ||
273      type == ringorder_IS
274    ) {
275      // todo: Consider supporting this later.
276      WerrorS("Schreyer order s/S/IS not supported in MathicGB interface.");
277      return false;
278    }
279    if (type == ringorder_am) {
280      // This block is a Schreyer-like ordering only used in Spielwiese.
281      // todo: Consider supporting it later.
282      WerrorS("Block type am not supported in MathicGB interface");
283      return false;
284    }
285    if (type == ringorder_L) {
286      WerrorS("Invalid L-block found in order of ring.");
287      return false;
288    }
289    if (type == ringorder_aa) {
290      // I don't know what an aa block is supposed to do.
291      WerrorS("aa ordering not supported by the MathicGB interface.");
292      return false;
293    }
294    if (type == ringorder_unspec) {
295      WerrorS("Invalid unspec-block found in order of ring.");
296      return false;
297    }
298    WerrorS("Unknown block type found in order of ring.");
299    return false;
300  }
301
302  if (!didSetComponentBefore) {
303    WerrorS("Expected to find a c/C block in ring.");
304    return false;
305  }
306
307  if (!conf.setMonomialOrder(baseOrder, gradings)) {
308    WerrorS("MathicGB does not support non-global orders.");
309    return false;
310  }
311  return true;
312}
313
314bool prOrderMatrix(ring r) {
315  const int varCount = r->N;
316  mgb::GroebnerConfiguration conf(101, varCount);
317  if (!setOrder(r, conf))
318    return false;
319  const std::vector<Exponent>& gradings = conf.monomialOrder().second;
320  if (gradings.size() % varCount != 0) {
321    WerrorS("Expected matrix to be a multiple of varCount.");
322    return false;
323  }
324  const size_t rowCount = gradings.size() / varCount;
325  std::cout << "Order matrix:\n";
326  for (size_t row = 0; row < rowCount; ++row) {
327    for (size_t col = 0; col < varCount; ++col)
328      std::cerr << ' ' << gradings[row * varCount + col];
329    std::cerr << '\n';
330  }
331  std::cerr
332    << "Base order: "
333    << mgb::GroebnerConfiguration::baseOrderName(conf.monomialOrder().first)
334    << '\n';
335  std::cerr << "Component before: " << conf.componentBefore() << '\n';
336  std::cerr << "Components ascending: " << conf.componentsAscending() << '\n';
337  std::cerr << "Schreyering: " << conf.schreyering() << '\n';
338}
339
340void prOrder(ring r) {
341  std::cout << "Printing order of ring.\n";
342  for (int block = 0; ; ++block) {
343    switch (r->order[block]) {
344    case ringorder_no: // end of blocks
345      return;
346
347    case ringorder_a:
348      std::cout << "a";
349      break;
350
351    case ringorder_a64: ///< for int64 weights
352      std::cout << "a64";
353      break;
354
355    case ringorder_c:
356      std::cout << "c";
357      break;
358
359    case ringorder_C:
360      std::cout << "C";
361      break;
362
363    case ringorder_M:
364      std::cout << "M";
365      break;
366
367    case ringorder_S: ///< S?
368      std::cout << "S";
369      break;
370
371    case ringorder_s: ///< s?
372      std::cout << "s";
373      break;
374
375    case ringorder_lp:
376      std::cout << "lp";
377      break;
378
379    case ringorder_dp:
380      std::cout << "dp";
381      break;
382
383    case ringorder_rp:
384      std::cout << "rp";
385      break;
386
387    case ringorder_Dp:
388      std::cout << "Dp";
389      break;
390
391    case ringorder_wp:
392      std::cout << "wp";
393      break;
394
395    case ringorder_Wp:
396      std::cout << "Wp";
397      break;
398
399    case ringorder_ls:
400      std::cout << "ls"; // not global
401      break;
402
403    case ringorder_ds:
404      std::cout << "ds"; // not global
405      break;
406
407    case ringorder_Ds:
408      std::cout << "Ds"; // not global
409      break;
410
411    case ringorder_ws:
412      std::cout << "ws"; // not global
413      break;
414
415    case ringorder_Ws:
416      std::cout << "Ws"; // not global
417      break;
418
419    case ringorder_am:
420      std::cout << "am";
421      break;
422
423    case ringorder_L:
424      std::cout << "L";
425      break;
426
427    // the following are only used internally
428    case ringorder_aa: ///< for idElimination, like a, except pFDeg, pWeigths ignore it
429      std::cout << "aa";
430      break;
431
432    case ringorder_rs: ///< opposite of ls
433      std::cout << "rs";
434      break;
435
436    case ringorder_IS: ///< Induced (Schreyer) ordering
437      std::cout << "IS";
438      break;
439
440    case ringorder_unspec:
441      std::cout << "unspec";
442      break;
443    }
444    const int b0 = r->block0[block];
445    const int b1 = r->block1[block];
446    std::cout << ' ' << b0 << ':' << b1 << " (" << r->wvhdl[block] << ")" << std::flush;
447    if (r->wvhdl[block] != 0 && b0 != 0) {
448      for (int v = 0; v <= b1 - b0; ++v)
449        std::cout << ' ' << r->wvhdl[block][v];
450    } else
451      std::cout << " null";
452    std::cout << '\n';
453  }
454}
455
456BOOLEAN prOrderX(leftv result, leftv arg) {
457  if (currRing == 0) {
458    WerrorS("There is no current ring.");
459    return TRUE;
460  }
461  prOrder(currRing);
462  prOrderMatrix(currRing);
463  result->rtyp=NONE;
464  return FALSE;
465}
466
467BOOLEAN setRingGlobal(leftv result, leftv arg) {
468  currRing->OrdSgn = 1;
469  result->rtyp=NONE;
470  return FALSE;
471}
472
473BOOLEAN mathicgb(leftv result, leftv arg)
474{
475  result->rtyp=NONE;
476
477  if (arg == NULL || arg->next != NULL || arg->Typ() != IDEAL_CMD) {
478    WerrorS("Syntax: mathicgb(<ideal>)");
479    return TRUE;
480  }
481  if (!rField_is_Zp(currRing)) {
482    WerrorS("Polynomial ring must be over Zp.");
483    return TRUE;
484  }
485
486  const int characteristic = n_GetChar(currRing);
487  const int varCount = currRing->N;
488  mgb::GroebnerConfiguration conf(characteristic, varCount);
489  if (!setOrder(currRing, conf))
490    return TRUE;
491  if (TEST_OPT_PROT)
492    conf.setLogging("all");
493
494  mgb::GroebnerInputIdealStream toMathic(conf);
495
496  const ideal id = static_cast<const ideal>(arg->Data());
497  const int size = IDELEMS(id);
498  toMathic.idealBegin(size);
499  for (int i = 0; i < size; ++i) {
500    const poly origP = id->m[i];
501    int termCount = 0;
502    for (poly p = origP; p != 0; p = pNext(p))
503      ++termCount;
504    toMathic.appendPolynomialBegin(termCount);
505
506    for (poly p = origP; p != 0; p = pNext(p)) {
507      toMathic.appendTermBegin();
508      for (int i = 1; i <= currRing->N; ++i)
509        toMathic.appendExponent(i - 1, pGetExp(p, i));
510      const long coefLong = reinterpret_cast<long>(pGetCoeff(p));
511      toMathic.appendTermDone(static_cast<int>(coefLong));
512    }
513    toMathic.appendPolynomialDone();
514  }
515  toMathic.idealDone();
516
517  MathicToSingStream fromMathic(characteristic, varCount);
518  mgb::computeGroebnerBasis(toMathic, fromMathic);
519
520  result->rtyp=IDEAL_CMD;
521  result->data = fromMathic.takeIdeal();
522  return FALSE;
523}
524
525template class std::vector<Exponent>;
526template void mgb::computeGroebnerBasis<MathicToSingStream>
527  (mgb::GroebnerInputIdealStream&, MathicToSingStream&);
528
529int SI_MOD_INIT(singmathic)(SModulFunctions* psModulFunctions)
530{
531  PrintS("Initializing Singular-Mathic interface Singmathic.\n");
532  psModulFunctions->iiAddCproc(
533    (currPack->libname ? currPack->libname : ""),
534    "mathicgb",
535    FALSE,
536    mathicgb
537  );
538  psModulFunctions->iiAddCproc(
539    (currPack->libname ? currPack->libname : ""),
540    "mathicgb_prOrder",
541    FALSE,
542    prOrderX
543  );
544  psModulFunctions->iiAddCproc(
545    (currPack->libname ? currPack->libname : ""),
546    "mathicgb_setRingGlobal",
547    FALSE,
548    setRingGlobal
549  );
550  return 0;
551}
552
553#else /* HAVE_MATHICGB */
554
555int SI_MOD_INIT(singmathic)(SModulFunctions* psModulFunctions)
556{
557  WerrorS(
558    "Cannot initialize the Singular interface to MathicGB "
559    "as this Singular executable was built without support "
560    "for MathicGB."
561  );
562  return 1;
563}
564
565#endif /* HAVE_MATHICGB */
Note: See TracBrowser for help on using the repository browser.