singmathic.cc
Go to the documentation of this file.
1 #include <kernel/mod2.h>
2 
3 #ifdef HAVE_MATHICGB
4 
5 #include <kernel/mod2.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 
21 
22 // Constructs a Singular ideal.
24 public:
26  mModulus(modulus),
27  mVarCount(varCount),
28  mPolyCount(0),
29  mTerm(0),
30  mIdeal(0)
31  {}
32 
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 
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 
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 
81 private:
82  void deleteIdeal() {
83  if (mIdeal != 0) {
84  idDelete(&mIdeal);
85  mIdeal = 0;
86  }
87  }
88 
91  size_t mPolyCount;
93  ::ideal mIdeal;
94 };
95 
96 #include <iostream>
97 
98 bool setOrder(ring r, mgb::GroebnerConfiguration& conf) {
99  const VarIndex varCount = conf.varCount();
100 
101  bool didSetComponentBefore = false;
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 
313 bool 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 
339 void 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 
456  if (currRing == 0) {
457  WerrorS("There is no current ring.");
458  return TRUE;
459  }
460  prOrder(currRing);
462  result->rtyp=NONE;
463  return FALSE;
464 }
465 
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 
521 template class std::vector<Exponent>;
522 template void mgb::computeGroebnerBasis<MathicToSingStream>
523  (mgb::GroebnerInputIdealStream&, MathicToSingStream&);
524 
525 extern "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 
544 int 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 
557 http://stackoverflow.com/questions/3786408/number-of-threads-used-by-intel-tbb
558 When you create the scheduler, you can specify the number of threads as
559 tbb::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 */
for idElimination, like a, except pFDeg, pWeigths ignore it
Definition: ring.h:99
for int64 weights
Definition: ring.h:79
Class used for (list of) interpreter objects.
Definition: subexpr.h:83
#define pSetm(p)
Definition: polys.h:253
BOOLEAN prOrderX(leftv result, leftv arg)
Definition: singmathic.cc:455
#define block
Definition: scanner.cc:662
#define idDelete(H)
delete an ideal
Definition: ideals.h:29
bool setOrder(ring r, mgb::GroebnerConfiguration &conf)
Definition: singmathic.cc:98
int SI_MOD_INIT() singmathic(SModulFunctions *psModulFunctions)
Definition: singmathic.cc:525
#define TEST_OPT_PROT
Definition: options.h:98
#define pSetExp(p, i, v)
Definition: polys.h:42
mgb::GroebnerConfiguration::VarIndex VarIndex
Definition: singmathic.cc:18
#define FALSE
Definition: auxiliary.h:94
Compatiblity layer for legacy polynomial operations (over currRing)
return P p
Definition: myNF.cc:203
opposite of ls
Definition: ring.h:100
MathicToSingStream(Coefficient modulus, VarIndex varCount)
Definition: singmathic.cc:25
Definition: tok.h:213
static FORCE_INLINE int n_GetChar(const coeffs r)
Return the characteristic of the coeff. domain.
Definition: coeffs.h:448
#define TRUE
Definition: auxiliary.h:98
void appendTermBegin(const mgb::GroebnerConfiguration::Component c)
Definition: singmathic.cc:48
void WerrorS(const char *s)
Definition: feFopen.cc:24
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy ...
Definition: monomials.h:51
int Typ()
Definition: subexpr.cc:996
void idealBegin(size_t polyCount)
Definition: singmathic.cc:40
#define pGetComp(p)
Component.
Definition: polys.h:37
void * data
Definition: subexpr.h:89
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:10
#define pGetExp(p, i)
Exponent.
Definition: polys.h:41
const ring r
Definition: syzextra.cc:208
VarIndex varCount() const
Definition: singmathic.cc:38
void appendExponent(VarIndex index, Exponent exponent)
Definition: singmathic.cc:56
mgb::GroebnerConfiguration::Exponent Exponent
Definition: singmathic.cc:19
rRingOrder_t
order stuff
Definition: ring.h:75
void prOrder(ring r)
Definition: singmathic.cc:339
#define pSetComp(p, v)
Definition: polys.h:38
int dim(ideal I, ring r)
const VarIndex mVarCount
Definition: singmathic.cc:90
int i
Definition: cfEzgcd.cc:123
Induced (Schreyer) ordering.
Definition: ring.h:101
void appendPolynomialBegin(size_t termCount)
Definition: singmathic.cc:46
S?
Definition: ring.h:83
BOOLEAN mathicgb(leftv result, leftv arg)
Definition: singmathic.cc:466
#define IDELEMS(i)
Definition: simpleideals.h:24
leftv next
Definition: subexpr.h:87
static int index(p_Length length, p_Ord ord)
Definition: p_Procs_Impl.h:592
int size(const CanonicalForm &f, const Variable &v)
int size ( const CanonicalForm & f, const Variable & v )
Definition: cf_ops.cc:600
bool prOrderMatrix(ring r)
Definition: singmathic.cc:313
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:495
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:38
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
void appendTermDone(Coefficient coefficient)
Definition: singmathic.cc:60
#define NULL
Definition: omList.c:10
const Coefficient mModulus
Definition: singmathic.cc:89
Coefficient modulus() const
Definition: singmathic.cc:37
#define pInit()
allocates a new monomial and initializes everything to 0
Definition: polys.h:61
package currPack
Definition: ipid.cc:63
int rtyp
Definition: subexpr.h:92
#define pNext(p)
Definition: monomials.h:43
void * Data()
Definition: subexpr.cc:1138
mgb::GroebnerConfiguration::BaseOrder BaseOrder
Definition: singmathic.cc:20
int exponent(const CanonicalForm &f, int q)
int exponent ( const CanonicalForm & f, int q )
polyrec * poly
Definition: hilb.h:10
mgb::GroebnerConfiguration::Coefficient Coefficient
Definition: singmathic.cc:17
void appendPolynomialDone()
Definition: singmathic.cc:65
int BOOLEAN
Definition: auxiliary.h:85
s?
Definition: ring.h:84
#define NONE
Definition: tok.h:216
::ideal takeIdeal()
Definition: singmathic.cc:75
return result
Definition: facAbsBiFact.cc:76
ListNode * next
Definition: janet.h:31