source: git/Singular/singmathic.cc @ 7222f9

spielwiese
Last change on this file since 7222f9 was a9c298, checked in by Hans Schoenemann <hannes@…>, 10 years ago
format stuff
  • Property mode set to 100644
File size: 16.2 KB
Line 
1// include header files
2#ifdef HAVE_CONFIG_H
3#include "singularconfig.h"
4#endif /* HAVE_CONFIG_H */
5
6#include <misc/auxiliary.h>
7#include <kernel/mod2.h>
8
9#include <misc/options.h>
10
11#include <kernel/febase.h>
12#include <kernel/ideals.h>
13#include <kernel/polys.h>
14
15#include <Singular/ipid.h>
16#include <Singular/mod_lib.h>
17
18#ifdef HAVE_MATHICGB
19
20#include <mathicgb.h>
21
22typedef mgb::GroebnerConfiguration::Coefficient Coefficient;
23typedef mgb::GroebnerConfiguration::VarIndex VarIndex;
24typedef mgb::GroebnerConfiguration::Exponent Exponent;
25typedef mgb::GroebnerConfiguration::BaseOrder BaseOrder;
26
27// Constructs a Singular ideal.
28class MathicToSingStream {
29public:
30  MathicToSingStream(Coefficient modulus, VarIndex varCount):
31    mModulus(modulus),
32    mVarCount(varCount),
33    mPolyCount(0),
34    mTerm(0),
35    mIdeal(0)
36  {}
37
38  ~MathicToSingStream() {deleteIdeal();}
39
40  // Mathic stream interface
41
42  Coefficient modulus() const {return mModulus;}
43  VarIndex varCount() const {return mModulus;}
44
45  void idealBegin(size_t polyCount) {
46    deleteIdeal();
47    mIdeal = idInit(polyCount);
48    mPolyCount = 0;
49  }
50
51  void appendPolynomialBegin(size_t termCount) {}
52
53  void appendTermBegin() {
54    if (mTerm == 0)
55      mTerm = mIdeal->m[mPolyCount] = pInit();
56    else
57      mTerm = mTerm->next = pInit();
58  }
59
60  void appendExponent(VarIndex index, Exponent exponent) {
61    pSetExp(mTerm, index + 1, exponent);
62  }
63
64  void appendTermDone(Coefficient coefficient) {
65    mTerm->coef = reinterpret_cast<number>(coefficient);
66    pSetm(mTerm);
67  }
68
69  void appendPolynomialDone() {
70    ++mPolyCount;
71    mTerm = 0;
72  }
73
74  void idealDone() {}
75
76
77  // Singular interface
78
79  ::ideal takeIdeal() {
80    ::ideal id = mIdeal;
81    mIdeal = 0;
82    return id;
83  }
84
85private:
86  void deleteIdeal() {
87    if (mIdeal != 0) {
88      idDelete(&mIdeal);
89      mIdeal = 0;
90    }
91  }
92
93  const Coefficient mModulus;
94  const VarIndex mVarCount;
95  size_t mPolyCount;
96  poly mTerm;
97  ::ideal mIdeal;
98};
99
100#include <iostream>
101
102bool 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
317bool 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
343void 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
459BOOLEAN 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
470BOOLEAN setRingGlobal(leftv result, leftv arg) {
471  currRing->OrdSgn = 1;
472  result->rtyp=NONE;
473  return FALSE;
474}
475
476BOOLEAN mathicgb(leftv result, leftv arg)
477{
478  result->rtyp=NONE;
479
480  if (arg == NULL || arg->next != NULL || arg->Typ() != IDEAL_CMD) {
481    WerrorS("Syntax: mathicgb(<ideal>)");
482    return TRUE;
483  }
484  if (!rField_is_Zp(currRing)) {
485    WerrorS("Polynomial ring must be over Zp.");
486    return TRUE;
487  }
488
489  const int characteristic = n_GetChar(currRing);
490  const int varCount = currRing->N;
491  mgb::GroebnerConfiguration conf(characteristic, varCount);
492  if (!setOrder(currRing, conf))
493    return TRUE;
494  if (TEST_OPT_PROT)
495    conf.setLogging("all");
496
497  mgb::GroebnerInputIdealStream toMathic(conf);
498
499  const ideal id = static_cast<const ideal>(arg->Data());
500  const int size = IDELEMS(id);
501  toMathic.idealBegin(size);
502  for (int i = 0; i < size; ++i) {
503    const poly origP = id->m[i];
504    int termCount = 0;
505    for (poly p = origP; p != 0; p = pNext(p))
506      ++termCount;
507    toMathic.appendPolynomialBegin(termCount);
508
509    for (poly p = origP; p != 0; p = pNext(p)) {
510      toMathic.appendTermBegin();
511      for (int i = 1; i <= currRing->N; ++i)
512        toMathic.appendExponent(i - 1, pGetExp(p, i));
513      const long coefLong = reinterpret_cast<long>(pGetCoeff(p));
514      toMathic.appendTermDone(static_cast<int>(coefLong));
515    }
516    toMathic.appendPolynomialDone();
517  }
518  toMathic.idealDone();
519
520  MathicToSingStream fromMathic(characteristic, varCount);
521  mgb::computeGroebnerBasis(toMathic, fromMathic);
522
523  result->rtyp=IDEAL_CMD;
524  result->data = fromMathic.takeIdeal();
525  return FALSE;
526}
527
528template class std::vector<Exponent>;
529template void mgb::computeGroebnerBasis<MathicToSingStream>
530  (mgb::GroebnerInputIdealStream&, MathicToSingStream&);
531
532int SI_MOD_INIT(singmathic)(SModulFunctions* psModulFunctions)
533{
534  PrintS("Initializing Singular-Mathic interface Singmathic.\n");
535  psModulFunctions->iiAddCproc(
536    (currPack->libname ? currPack->libname : ""),
537    "mathicgb",
538    FALSE,
539    mathicgb
540  );
541  psModulFunctions->iiAddCproc(
542    (currPack->libname ? currPack->libname : ""),
543    "mathicgb_prOrder",
544    FALSE,
545    prOrderX
546  );
547  psModulFunctions->iiAddCproc(
548    (currPack->libname ? currPack->libname : ""),
549    "mathicgb_setRingGlobal",
550    FALSE,
551    setRingGlobal
552  );
553  return 0;
554}
555
556#else /* HAVE_MATHICGB */
557
558int SI_MOD_INIT(singmathic)(SModulFunctions* psModulFunctions)
559{
560  WerrorS(
561    "Cannot initialize the Singular interface to MathicGB "
562    "as this Singular executable was built without support "
563    "for MathicGB."
564  );
565  return 1;
566}
567
568#endif /* HAVE_MATHICGB */
Note: See TracBrowser for help on using the repository browser.