My Project  debian-1:4.1.1-p2+ds-4build4
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/feOpt.h"
14 #include "Singular/mod_lib.h"
15 
16 #include <mathicgb.h>
17 
22 
23 // Constructs a Singular ideal.
25 public:
29  mPolyCount(0),
30  mTerm(0),
31  mIdeal(0)
32  {}
33 
35 
36  // Mathic stream interface
37 
38  Coefficient modulus() const {return mModulus;}
39  VarIndex varCount() const {return mModulus;}
40 
41  void idealBegin(size_t polyCount) {
42  deleteIdeal();
43  mIdeal = idInit(polyCount);
44  mPolyCount = 0;
45  }
46 
47  void appendPolynomialBegin(size_t termCount) {}
48 
49  void appendTermBegin(const mgb::GroebnerConfiguration::Component c) {
50  if (mTerm == 0)
51  mTerm = mIdeal->m[mPolyCount] = pInit();
52  else
53  mTerm = mTerm->next = pInit();
54  pSetComp(mTerm,c);
55  }
56 
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 
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 
82 private:
83  void deleteIdeal() {
84  if (mIdeal != 0) {
85  idDelete(&mIdeal);
86  mIdeal = 0;
87  }
88  }
89 
92  size_t mPolyCount;
93  poly mTerm;
94  ::ideal mIdeal;
95 };
96 
97 #include <iostream>
98 
99 bool setOrder(ring r, mgb::GroebnerConfiguration& conf) {
100  const VarIndex varCount = conf.varCount();
101 
102  bool didSetComponentBefore = false;
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 
314 bool prOrderMatrix(ring r) {
315  const int varCount = r->N;
316  mgb::GroebnerConfiguration conf(101, varCount,0);
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 
340 void 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 
457  if (currRing == 0) {
458  WerrorS("There is no current ring.");
459  return TRUE;
460  }
461  prOrder(currRing);
463  result->rtyp=NONE;
464  return FALSE;
465 }
466 
468 {
469  result->rtyp=NONE;
470 
471  if (arg == NULL || arg->next != NULL ||
472  ((arg->Typ() != IDEAL_CMD) &&(arg->Typ() != MODUL_CMD))){
473  WerrorS("Syntax: mathicgb(<ideal>/<module>)");
474  return TRUE;
475  }
476  if (!rField_is_Zp(currRing)) {
477  WerrorS("Polynomial ring must be over Zp.");
478  return TRUE;
479  }
480 
481  const int characteristic = n_GetChar(currRing);
482  const int varCount = currRing->N;
483  const ideal I=(ideal) arg->Data();
484  mgb::GroebnerConfiguration conf(characteristic, varCount,I->rank);
485  feOptIndex fno=feGetOptIndex(FE_OPT_THREADS);
486  //conf.setMaxThreadCount(0); // default number of cores
487  conf.setMaxThreadCount((int)(long)feOptSpec[fno].value);
488  if (!setOrder(currRing, conf))
489  return TRUE;
490  if (TEST_OPT_PROT)
491  conf.setLogging("all");
492 
493  mgb::GroebnerInputIdealStream toMathic(conf);
494 
495  const ideal id = static_cast<const ideal>(arg->Data());
496  const int size = IDELEMS(id);
497  toMathic.idealBegin(size);
498  for (int i = 0; i < size; ++i) {
499  const poly origP = id->m[i];
500  int termCount = 0;
501  for (poly p = origP; p != 0; p = pNext(p))
502  ++termCount;
503  toMathic.appendPolynomialBegin(termCount);
504 
505  for (poly p = origP; p != 0; p = pNext(p)) {
506  toMathic.appendTermBegin(pGetComp(p));
507  for (int i = 1; i <= currRing->N; ++i)
508  toMathic.appendExponent(i - 1, pGetExp(p, i));
509  const long coefLong = reinterpret_cast<long>(pGetCoeff(p));
510  toMathic.appendTermDone(static_cast<int>(coefLong));
511  }
512  toMathic.appendPolynomialDone();
513  }
514  toMathic.idealDone();
515 
516  MathicToSingStream fromMathic(characteristic, varCount);
517  mgb::computeGroebnerBasis(toMathic, fromMathic);
518 
519  result->rtyp = arg->Typ();
520  result->data = fromMathic.takeIdeal();
521  return FALSE;
522 }
523 
524 template class std::vector<Exponent>;
525 template void mgb::computeGroebnerBasis<MathicToSingStream>
526  (mgb::GroebnerInputIdealStream&, MathicToSingStream&);
527 
528 extern "C" int SI_MOD_INIT(singmathic)(SModulFunctions* psModulFunctions)
529 {
530  psModulFunctions->iiAddCproc(
531  (currPack->libname ? currPack->libname : ""),
532  "mathicgb",
533  FALSE,
534  mathicgb
535  );
536  psModulFunctions->iiAddCproc(
537  (currPack->libname ? currPack->libname : ""),
538  "mathicgb_prOrder",
539  FALSE,
540  prOrderX
541  );
542  return MAX_TOK;
543 }
544 
545 /* #else
546 
547 int SI_MOD_INIT(singmathic)(SModulFunctions* psModulFunctions)
548 {
549  WerrorS(
550  "Cannot initialize the Singular interface to MathicGB "
551  "as this Singular executable was built without support "
552  "for MathicGB."
553  );
554  return 1;
555 }
556 */
557 
558 /* ressources: ------------------------------------------------------------
559 
560 http://stackoverflow.com/questions/3786408/number-of-threads-used-by-intel-tbb
561 When you create the scheduler, you can specify the number of threads as
562 tbb::task_scheduler_init init(nthread);
563 
564  How do I know how many threads are available?
565 
566  Do not ask!
567 
568  Not even the scheduler knows how many threads really are available
569  There may be other processes running on the machine
570  Routine may be nested inside other parallel routines
571 
572  conf.setMaxThreadCount(0); // default number of cores
573 */
574 #endif /* HAVE_MATHICGB */
int BOOLEAN
Definition: auxiliary.h:85
#define TRUE
Definition: auxiliary.h:98
#define FALSE
Definition: auxiliary.h:94
int size(const CanonicalForm &f, const Variable &v)
int size ( const CanonicalForm & f, const Variable & v )
Definition: cf_ops.cc:600
int i
Definition: cfEzgcd.cc:125
int p
Definition: cfModGcd.cc:4019
::ideal takeIdeal()
Definition: singmathic.cc:76
void appendPolynomialDone()
Definition: singmathic.cc:66
void appendTermBegin(const mgb::GroebnerConfiguration::Component c)
Definition: singmathic.cc:49
VarIndex varCount() const
Definition: singmathic.cc:39
Coefficient modulus() const
Definition: singmathic.cc:38
void idealBegin(size_t polyCount)
Definition: singmathic.cc:41
MathicToSingStream(Coefficient modulus, VarIndex varCount)
Definition: singmathic.cc:26
void appendTermDone(Coefficient coefficient)
Definition: singmathic.cc:61
void appendExponent(VarIndex index, Exponent exponent)
Definition: singmathic.cc:57
const VarIndex mVarCount
Definition: singmathic.cc:91
void appendPolynomialBegin(size_t termCount)
Definition: singmathic.cc:47
const Coefficient mModulus
Definition: singmathic.cc:90
Class used for (list of) interpreter objects.
Definition: subexpr.h:83
int Typ()
Definition: subexpr.cc:992
void * Data()
Definition: subexpr.cc:1134
leftv next
Definition: subexpr.h:86
static FORCE_INLINE int n_GetChar(const coeffs r)
Return the characteristic of the coeff. domain.
Definition: coeffs.h:445
return result
Definition: facAbsBiFact.cc:76
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
void WerrorS(const char *s)
Definition: feFopen.cc:24
feOptIndex
Definition: feOptGen.h:15
feOptIndex feGetOptIndex(const char *name)
Definition: feOpt.cc:101
struct fe_option feOptSpec[]
int exponent(const CanonicalForm &f, int q)
int exponent ( const CanonicalForm & f, int q )
@ IDEAL_CMD
Definition: grammar.cc:283
@ MODUL_CMD
Definition: grammar.cc:286
#define idDelete(H)
delete an ideal
Definition: ideals.h:29
package currPack
Definition: ipid.cc:59
ListNode * next
Definition: janet.h:31
#define pNext(p)
Definition: monomials.h:43
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
#define NULL
Definition: omList.c:10
#define TEST_OPT_PROT
Definition: options.h:102
static int index(p_Length length, p_Ord ord)
Definition: p_Procs_Impl.h:592
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
Compatiblity layer for legacy polynomial operations (over currRing)
#define pSetm(p)
Definition: polys.h:257
#define pGetComp(p)
Component.
Definition: polys.h:37
#define pSetComp(p, v)
Definition: polys.h:38
#define pGetExp(p, i)
Exponent.
Definition: polys.h:41
#define pInit()
allocates a new monomial and initializes everything to 0
Definition: polys.h:61
#define pSetExp(p, i, v)
Definition: polys.h:42
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:495
rRingOrder_t
order stuff
Definition: ring.h:75
@ ringorder_lp
Definition: ring.h:84
@ ringorder_a
Definition: ring.h:77
@ ringorder_am
Definition: ring.h:95
@ ringorder_a64
for int64 weights
Definition: ring.h:78
@ ringorder_rs
opposite of ls
Definition: ring.h:99
@ ringorder_C
Definition: ring.h:80
@ ringorder_S
S?
Definition: ring.h:82
@ ringorder_ds
Definition: ring.h:91
@ ringorder_Dp
Definition: ring.h:87
@ ringorder_unspec
Definition: ring.h:101
@ ringorder_L
Definition: ring.h:96
@ ringorder_Ds
Definition: ring.h:92
@ ringorder_dp
Definition: ring.h:85
@ ringorder_c
Definition: ring.h:79
@ ringorder_rp
Definition: ring.h:86
@ ringorder_aa
for idElimination, like a, except pFDeg, pWeigths ignore it
Definition: ring.h:98
@ ringorder_no
Definition: ring.h:76
@ ringorder_Wp
Definition: ring.h:89
@ ringorder_ws
Definition: ring.h:93
@ ringorder_Ws
Definition: ring.h:94
@ ringorder_IS
Induced (Schreyer) ordering.
Definition: ring.h:100
@ ringorder_ls
Definition: ring.h:90
@ ringorder_s
s?
Definition: ring.h:83
@ ringorder_wp
Definition: ring.h:88
@ ringorder_M
Definition: ring.h:81
#define block
Definition: scanner.cc:665
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:37
#define IDELEMS(i)
Definition: simpleideals.h:24
mgb::GroebnerConfiguration::Coefficient Coefficient
Definition: singmathic.cc:18
int SI_MOD_INIT() singmathic(SModulFunctions *psModulFunctions)
Definition: singmathic.cc:528
BOOLEAN mathicgb(leftv result, leftv arg)
Definition: singmathic.cc:467
BOOLEAN prOrderX(leftv result, leftv arg)
Definition: singmathic.cc:456
mgb::GroebnerConfiguration::Exponent Exponent
Definition: singmathic.cc:20
bool prOrderMatrix(ring r)
Definition: singmathic.cc:314
bool setOrder(ring r, mgb::GroebnerConfiguration &conf)
Definition: singmathic.cc:99
mgb::GroebnerConfiguration::VarIndex VarIndex
Definition: singmathic.cc:19
mgb::GroebnerConfiguration::BaseOrder BaseOrder
Definition: singmathic.cc:21
void prOrder(ring r)
Definition: singmathic.cc:340
@ MAX_TOK
Definition: tok.h:215
#define NONE
Definition: tok.h:218
int dim(ideal I, ring r)