source: git/Singular/dyn_modules/singmathic/singmathic.cc @ 00d5f4

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