sparsmat.cc
Go to the documentation of this file.
1 /****************************************
2 * Computer Algebra System SINGULAR *
3 ****************************************/
4 
5 /*
6 * ABSTRACT: operations with sparse matrices (bareiss, ...)
7 */
8 
9 
10 
11 
12 #include <misc/auxiliary.h>
13 
14 #include <omalloc/omalloc.h>
15 #include <misc/mylimits.h>
16 
17 #include <misc/options.h>
18 
19 #include <reporter/reporter.h>
20 #include <misc/intvec.h>
21 
22 #include <coeffs/numbers.h>
23 
24 #include "monomials/ring.h"
25 #include "monomials/p_polys.h"
26 
27 #include "simpleideals.h"
28 
29 
30 #include "sparsmat.h"
31 #include "prCopy.h"
32 
33 #include "templates/p_Procs.h"
34 
35 #include "kbuckets.h"
36 #include "operations/p_Mult_q.h"
37 
38 // define SM_NO_BUCKETS, if sparsemat stuff should not use buckets
39 // #define SM_NO_BUCKETS
40 
41 // this is also influenced by TEST_OPT_NOTBUCKETS
42 #ifndef SM_NO_BUCKETS
43 // buckets do carry a small additional overhead: only use them if
44 // min-length of polys is >= SM_MIN_LENGTH_BUCKET
45 #define SM_MIN_LENGTH_BUCKET MIN_LENGTH_BUCKET - 5
46 #else
47 #define SM_MIN_LENGTH_BUCKET MAX_INT_VAL
48 #endif
49 
50 typedef struct smprec sm_prec;
51 typedef sm_prec * smpoly;
52 struct smprec
53 {
54  smpoly n; // the next element
55  int pos; // position
56  int e; // level
57  poly m; // the element
58  float f; // complexity of the element
59 };
60 
61 
62 /* declare internal 'C' stuff */
63 static void sm_ExactPolyDiv(poly, poly, const ring);
64 static BOOLEAN sm_IsNegQuot(poly, const poly, const poly, const ring);
65 static void sm_ExpMultDiv(poly, const poly, const poly, const ring);
66 static void sm_PolyDivN(poly, const number, const ring);
67 static BOOLEAN smSmaller(poly, poly);
68 static void sm_CombineChain(poly *, poly, const ring);
69 static void sm_FindRef(poly *, poly *, poly, const ring);
70 
71 static void sm_ElemDelete(smpoly *, const ring);
72 static smpoly smElemCopy(smpoly);
73 static float sm_PolyWeight(smpoly, const ring);
74 static smpoly sm_Poly2Smpoly(poly, const ring);
75 static poly sm_Smpoly2Poly(smpoly, const ring);
76 static BOOLEAN sm_HaveDenom(poly, const ring);
77 static number sm_Cleardenom(ideal, const ring);
78 
80 
82  poly a, poly b, const ring currRing)
83 {
84  if (rOrd_is_Comp_dp(currRing) && currRing->ExpL_Size > 2)
85  {
86  // pp_Mult_Coeff_mm_DivSelectMult only works for (c/C,dp) and
87  // ExpL_Size > 2
88  // should be generalized, at least to dp with ExpL_Size == 2
89  // (is the case for 1 variable)
90  int shorter;
91  p = currRing->p_Procs->pp_Mult_Coeff_mm_DivSelectMult(p, m, a, b,
92  shorter, currRing);
93  lp -= shorter;
94  }
95  else
96  {
97  p = pp_Mult_Coeff_mm_DivSelect(p, lp, m, currRing);
98  sm_ExpMultDiv(p, a, b, currRing);
99  }
100  return p;
101 }
102 
104 {
105  int lp = 0;
106  return pp_Mult_Coeff_mm_DivSelect_MultDiv(p, lp, m, a, b, currRing);
107 }
108 
109 
110 /* class for sparse matrix:
111 * 3 parts of matrix during the algorithm
112 * m_act[cols][pos(rows)] => m_row[rows][pos(cols)] => m_res[cols][pos(rows)]
113 * input pivotcols as rows result
114 * pivot like a stack from pivot and pivotcols
115 * elimination rows reordered according
116 * to pivot choise
117 * stored in perm
118 * a step is as follows
119 * - search of pivot (smPivot)
120 * - swap pivot column and last column (smSwapC)
121 * select pivotrow to piv and red (smSelectPR)
122 * consider sign
123 * - elimination (smInitElim, sm0Elim, sm1Elim)
124 * clear zero column as result of elimination (smZeroElim)
125 * - tranfer from
126 * piv and m_row to m_res (smRowToCol)
127 * m_act to m_row (smColToRow)
128 */
130 private:
131  int nrows, ncols; // dimension of the problem
132  int sign; // for determinant (start: 1)
133  int act; // number of unreduced columns (start: ncols)
134  int crd; // number of reduced columns (start: 0)
135  int tored; // border for rows to reduce
136  int inred; // unreducable part
137  int rpiv, cpiv; // position of the pivot
138  int normalize; // Normalization flag
139  int *perm; // permutation of rows
140  float wpoints; // weight of all points
141  float *wrw, *wcl; // weights of rows and columns
142  smpoly * m_act; // unreduced columns
143  smpoly * m_res; // reduced columns (result)
144  smpoly * m_row; // reduced part of rows
145  smpoly red; // row to reduce
146  smpoly piv, oldpiv; // pivot and previous pivot
147  smpoly dumm; // allocated dummy
148  ring _R;
149 
150  void smColToRow();
151  void smRowToCol();
152  void smFinalMult();
153  void smSparseHomog();
154  void smWeights();
155  void smPivot();
156  void smNewWeights();
157  void smNewPivot();
158  void smZeroElim();
159  void smToredElim();
160  void smCopToRes();
161  void smSelectPR();
162  void sm1Elim();
163  void smHElim();
164  void smMultCol();
165  poly smMultPoly(smpoly);
166  void smActDel();
167  void smColDel();
168  void smPivDel();
169  void smSign();
170  void smInitPerm();
171  int smCheckNormalize();
172  void smNormalize();
173 public:
174  sparse_mat(ideal, const ring);
175  ~sparse_mat();
176  int smGetSign() { return sign; }
177  smpoly * smGetAct() { return m_act; }
178  int smGetRed() { return tored; }
179  ideal smRes2Mod();
180  poly smDet();
181  void smNewBareiss(int, int);
182  void smToIntvec(intvec *);
183 };
184 
185 /*
186 * estimate maximal exponent for det or minors,
187 * the module m has di vectors and maximal rank ra,
188 * estimate yields for the tXt minors,
189 * we have di,ra >= t
190 */
191 static void smMinSelect(long *, int, int);
192 
193 long sm_ExpBound( ideal m, int di, int ra, int t, const ring currRing)
194 {
195  poly p;
196  long kr, kc;
197  long *r, *c;
198  int al, bl, i, j, k;
199 
200  if (ra==0) ra=1;
201  al = di*sizeof(long);
202  c = (long *)omAlloc(al);
203  bl = ra*sizeof(long);
204  r = (long *)omAlloc0(bl);
205  for (i=di-1;i>=0;i--)
206  {
207  kc = 0;
208  p = m->m[i];
209  while(p!=NULL)
210  {
211  k = p_GetComp(p, currRing)-1;
212  kr = r[k];
213  for (j=rVar(currRing);j>0;j--)
214  {
215  if(p_GetExp(p,j, currRing)>kc)
216  kc=p_GetExp(p,j, currRing);
217  if(p_GetExp(p,j, currRing)>kr)
218  kr=p_GetExp(p,j, currRing);
219  }
220  r[k] = kr;
221  pIter(p);
222  }
223  c[i] = kc;
224  }
225  if (t<di) smMinSelect(c, t, di);
226  if (t<ra) smMinSelect(r, t, ra);
227  kr = kc = 0;
228  for (j=t-1;j>=0;j--)
229  {
230  kr += r[j];
231  kc += c[j];
232  }
233  omFreeSize((ADDRESS)c, al);
234  omFreeSize((ADDRESS)r, bl);
235  if (kr<kc) kc = kr;
236  if (kr<1) kr = 1;
237  return kr;
238 }
239 
240 static void smMinSelect(long *c, int t, int d)
241 {
242  long m;
243  int pos, i;
244  do
245  {
246  d--;
247  pos = d;
248  m = c[pos];
249  for (i=d-1;i>=0;i--)
250  {
251  if(c[i]<m)
252  {
253  pos = i;
254  m = c[i];
255  }
256  }
257  for (i=pos;i<d;i++) c[i] = c[i+1];
258  } while (d>t);
259 }
260 
261 /* ----------------- ops with rings ------------------ */
262 ring sm_RingChange(const ring origR, long bound)
263 {
264 // *origR =currRing;
265  ring tmpR=rCopy0(origR,FALSE,FALSE);
266  int *ord=(int*)omAlloc0(3*sizeof(int));
267  int *block0=(int*)omAlloc(3*sizeof(int));
268  int *block1=(int*)omAlloc(3*sizeof(int));
269  ord[0]=ringorder_c;
270  ord[1]=ringorder_dp;
271  tmpR->order=ord;
272  tmpR->OrdSgn=1;
273  block0[1]=1;
274  tmpR->block0=block0;
275  block1[1]=tmpR->N;
276  tmpR->block1=block1;
277  tmpR->bitmask = 2*bound;
278  tmpR->wvhdl = (int **)omAlloc0((3) * sizeof(int*));
279 
280 // ???
281 // if (tmpR->bitmask > currRing->bitmask) tmpR->bitmask = currRing->bitmask;
282 
283  rComplete(tmpR,1);
284  if (origR->qideal!=NULL)
285  {
286  tmpR->qideal= idrCopyR_NoSort(origR->qideal, origR, tmpR);
287  }
288  if (TEST_OPT_PROT)
289  Print("[%ld:%d]", (long) tmpR->bitmask, tmpR->ExpL_Size);
290  return tmpR;
291 }
292 
294 {
295  if (r->qideal!=NULL) id_Delete(&(r->qideal),r);
297 }
298 
299 /*2
300 * Bareiss or Chinese remainder ?
301 * I is dXd
302 * sw = TRUE -> I Matrix
303 * FALSE -> I Module
304 * return True -> change Type
305 * FALSE -> same Type
306 */
307 BOOLEAN sm_CheckDet(ideal I, int d, BOOLEAN sw, const ring r)
308 {
309  int s,t,i;
310  poly p;
311 
312  if (d>100)
313  return sw;
314  if (!rField_is_Q(r))
315  return sw;
316  s = t = 0;
317  if (sw)
318  {
319  for(i=IDELEMS(I)-1;i>=0;i--)
320  {
321  p=I->m[i];
322  if (p!=NULL)
323  {
324  if(!p_IsConstant(p,r))
325  return sw;
326  s++;
327  t+=n_Size(pGetCoeff(p),r->cf);
328  }
329  }
330  }
331  else
332  {
333  for(i=IDELEMS(I)-1;i>=0;i--)
334  {
335  p=I->m[i];
336  if (!p_IsConstantPoly(p,r))
337  return sw;
338  while (p!=NULL)
339  {
340  s++;
341  t+=n_Size(pGetCoeff(p),r->cf);
342  pIter(p);
343  }
344  }
345  }
346  s*=15;
347  if (t>s)
348  return !sw;
349  else
350  return sw;
351 }
352 
353 /* ----------------- basics (used from 'C') ------------------ */
354 /*2
355 *returns the determinant of the module I;
356 *uses Bareiss algorithm
357 */
358 poly sm_CallDet(ideal I,const ring R)
359 {
360  if (I->ncols != I->rank)
361  {
362  Werror("det of %ld x %d module (matrix)",I->rank,I->ncols);
363  return NULL;
364  }
365  int r=id_RankFreeModule(I,R);
366  if (I->ncols != r) // some 0-lines at the end
367  {
368  return NULL;
369  }
370  long bound=sm_ExpBound(I,r,r,r,R);
371  number diag,h=n_Init(1,R->cf);
372  poly res;
373  ring tmpR;
374  sparse_mat *det;
375  ideal II;
376 
377  tmpR=sm_RingChange(R,bound);
378  II = idrCopyR(I, R, tmpR);
379  diag = sm_Cleardenom(II,tmpR);
380  det = new sparse_mat(II,tmpR);
381  id_Delete(&II,tmpR);
382  if (det->smGetAct() == NULL)
383  {
384  delete det;
385  sm_KillModifiedRing(tmpR);
386  return NULL;
387  }
388  res=det->smDet();
389  if(det->smGetSign()<0) res=p_Neg(res,tmpR);
390  delete det;
391  res = prMoveR(res, tmpR, R);
392  sm_KillModifiedRing(tmpR);
393  if (!n_Equal(diag,h,R->cf))
394  {
395  p_Mult_nn(res,diag,R);
396  p_Normalize(res,R);
397  }
398  n_Delete(&diag,R->cf);
399  n_Delete(&h,R->cf);
400  return res;
401 }
402 
403 void sm_CallBareiss(ideal I, int x, int y, ideal & M, intvec **iv, const ring R)
404 {
405  int r=id_RankFreeModule(I,R),t=r;
406  int c=IDELEMS(I),s=c;
407  long bound;
408  ring tmpR;
409  sparse_mat *bareiss;
410 
411  if ((x>0) && (x<t))
412  t-=x;
413  if ((y>1) && (y<s))
414  s-=y;
415  if (t>s) t=s;
416  bound=sm_ExpBound(I,c,r,t,R);
417  tmpR=sm_RingChange(R,bound);
418  ideal II = idrCopyR(I, R, tmpR);
419  bareiss = new sparse_mat(II,tmpR);
420  if (bareiss->smGetAct() == NULL)
421  {
422  delete bareiss;
423  *iv=new intvec(1,rVar(tmpR));
424  }
425  else
426  {
427  id_Delete(&II,tmpR);
428  bareiss->smNewBareiss(x, y);
429  II = bareiss->smRes2Mod();
430  *iv = new intvec(bareiss->smGetRed());
431  bareiss->smToIntvec(*iv);
432  delete bareiss;
433  II = idrMoveR(II,tmpR,R);
434  }
435  sm_KillModifiedRing(tmpR);
436  M=II;
437 }
438 
439 /*
440 * constructor
441 */
442 sparse_mat::sparse_mat(ideal smat, const ring RR)
443 {
444  int i;
445  poly* pmat;
446  _R=RR;
447 
448  ncols = smat->ncols;
449  nrows = id_RankFreeModule(smat,RR);
450  if (nrows <= 0)
451  {
452  m_act = NULL;
453  return;
454  }
455  sign = 1;
456  inred = act = ncols;
457  crd = 0;
458  tored = nrows; // without border
459  i = tored+1;
460  perm = (int *)omAlloc(sizeof(int)*(i+1));
461  perm[i] = 0;
462  m_row = (smpoly *)omAlloc0(sizeof(smpoly)*i);
463  wrw = (float *)omAlloc(sizeof(float)*i);
464  i = ncols+1;
465  wcl = (float *)omAlloc(sizeof(float)*i);
466  m_act = (smpoly *)omAlloc(sizeof(smpoly)*i);
467  m_res = (smpoly *)omAlloc0(sizeof(smpoly)*i);
468  dumm = (smpoly)omAllocBin(smprec_bin);
469  m_res[0] = (smpoly)omAllocBin(smprec_bin);
470  m_res[0]->m = NULL;
471  pmat = smat->m;
472  for(i=ncols; i; i--)
473  {
474  m_act[i] = sm_Poly2Smpoly(pmat[i-1], RR);
475  pmat[i-1] = NULL;
476  }
477  this->smZeroElim();
478  oldpiv = NULL;
479 }
480 
481 /*
482 * destructor
483 */
485 {
486  int i;
487  if (m_act == NULL) return;
488  omFreeBin((ADDRESS)m_res[0], smprec_bin);
489  omFreeBin((ADDRESS)dumm, smprec_bin);
490  i = ncols+1;
491  omFreeSize((ADDRESS)m_res, sizeof(smpoly)*i);
492  omFreeSize((ADDRESS)m_act, sizeof(smpoly)*i);
493  omFreeSize((ADDRESS)wcl, sizeof(float)*i);
494  i = nrows+1;
495  omFreeSize((ADDRESS)wrw, sizeof(float)*i);
496  omFreeSize((ADDRESS)m_row, sizeof(smpoly)*i);
497  omFreeSize((ADDRESS)perm, sizeof(int)*(i+1));
498 }
499 
500 /*
501 * transform the result to a module
502 */
503 
505 {
506  ideal res = idInit(crd, crd);
507  int i;
508 
509  for (i=crd; i; i--)
510  {
511  res->m[i-1] = sm_Smpoly2Poly(m_res[i],_R);
512  res->rank=si_max(res->rank, p_MaxComp(res->m[i-1],_R));
513  }
514  return res;
515 }
516 
517 /*
518 * permutation of rows
519 */
521 {
522  int i;
523 
524  for (i=v->rows()-1; i>=0; i--)
525  (*v)[i] = perm[i+1];
526 }
527 /* ---------------- the algorithm's ------------------ */
528 /*
529 * the determinant (up to sign), uses new Bareiss elimination
530 */
532 {
533  poly res = NULL;
534 
535  if (sign == 0)
536  {
537  this->smActDel();
538  return NULL;
539  }
540  if (act < 2)
541  {
542  if (act != 0) res = m_act[1]->m;
543  omFreeBin((void *)m_act[1], smprec_bin);
544  return res;
545  }
546  normalize = 0;
547  this->smInitPerm();
548  this->smPivot();
549  this->smSign();
550  this->smSelectPR();
551  this->sm1Elim();
552  crd++;
553  m_res[crd] = piv;
554  this->smColDel();
555  act--;
556  this->smZeroElim();
557  if (sign == 0)
558  {
559  this->smActDel();
560  return NULL;
561  }
562  if (act < 2)
563  {
564  this->smFinalMult();
565  this->smPivDel();
566  if (act != 0) res = m_act[1]->m;
567  omFreeBin((void *)m_act[1], smprec_bin);
568  return res;
569  }
570  loop
571  {
572  this->smNewPivot();
573  this->smSign();
574  this->smSelectPR();
575  this->smMultCol();
576  this->smHElim();
577  crd++;
578  m_res[crd] = piv;
579  this->smColDel();
580  act--;
581  this->smZeroElim();
582  if (sign == 0)
583  {
584  this->smPivDel();
585  this->smActDel();
586  return NULL;
587  }
588  if (act < 2)
589  {
590  if (TEST_OPT_PROT) PrintS(".\n");
591  this->smFinalMult();
592  this->smPivDel();
593  if (act != 0) res = m_act[1]->m;
594  omFreeBin((void *)m_act[1], smprec_bin);
595  return res;
596  }
597  }
598 }
599 
600 /*
601 * the new Bareiss elimination:
602 * - with x unreduced last rows, pivots from here are not allowed
603 * - the method will finish for number of unreduced columns < y
604 */
606 {
607  if ((x > 0) && (x < nrows))
608  {
609  tored -= x;
610  this->smToredElim();
611  }
612  if (y < 1) y = 1;
613  if (act <= y)
614  {
615  this->smCopToRes();
616  return;
617  }
618  normalize = this->smCheckNormalize();
619  if (normalize) this->smNormalize();
620  this->smPivot();
621  this->smSelectPR();
622  this->sm1Elim();
623  crd++;
624  this->smColToRow();
625  act--;
626  this->smRowToCol();
627  this->smZeroElim();
628  if (tored != nrows)
629  this->smToredElim();
630  if (act <= y)
631  {
632  this->smFinalMult();
633  this->smCopToRes();
634  return;
635  }
636  loop
637  {
638  if (normalize) this->smNormalize();
639  this->smNewPivot();
640  this->smSelectPR();
641  this->smMultCol();
642  this->smHElim();
643  crd++;
644  this->smColToRow();
645  act--;
646  this->smRowToCol();
647  this->smZeroElim();
648  if (tored != nrows)
649  this->smToredElim();
650  if (act <= y)
651  {
652  if (TEST_OPT_PROT) PrintS(".\n");
653  this->smFinalMult();
654  this->smCopToRes();
655  return;
656  }
657  }
658 }
659 
660 /* ----------------- pivot method ------------------ */
661 
662 /*
663 * prepare smPivot, compute weights for rows and columns
664 * and the weight for all points
665 */
667 {
668  float wc, wp, w;
669  smpoly a;
670  int i;
671 
672  wp = 0.0;
673  for (i=tored; i; i--) wrw[i] = 0.0; // ???
674  for (i=act; i; i--)
675  {
676  wc = 0.0;
677  a = m_act[i];
678  loop
679  {
680  if (a->pos > tored)
681  break;
682  w = a->f = sm_PolyWeight(a,_R);
683  wc += w;
684  wrw[a->pos] += w;
685  a = a->n;
686  if (a == NULL)
687  break;
688  }
689  wp += wc;
690  wcl[i] = wc;
691  }
692  wpoints = wp;
693 }
694 
695 /*
696 * compute pivot
697 */
699 {
700  float wopt = 1.0e30;
701  float wc, wr, wp, w;
702  smpoly a;
703  int i, copt, ropt;
704 
705  this->smWeights();
706  for (i=act; i; i--)
707  {
708  a = m_act[i];
709  loop
710  {
711  if (a->pos > tored)
712  break;
713  w = a->f;
714  wc = wcl[i]-w;
715  wr = wrw[a->pos]-w;
716  if ((wr<0.25) || (wc<0.25)) // row or column with only one point
717  {
718  if (w<wopt)
719  {
720  wopt = w;
721  copt = i;
722  ropt = a->pos;
723  }
724  }
725  else // elimination
726  {
727  wp = w*(wpoints-wcl[i]-wr);
728  wp += wr*wc;
729  if (wp < wopt)
730  {
731  wopt = wp;
732  copt = i;
733  ropt = a->pos;
734  }
735  }
736  a = a->n;
737  if (a == NULL)
738  break;
739  }
740  }
741  rpiv = ropt;
742  cpiv = copt;
743  if (cpiv != act)
744  {
745  a = m_act[act];
746  m_act[act] = m_act[cpiv];
747  m_act[cpiv] = a;
748  }
749 }
750 
751 /*
752 * prepare smPivot, compute weights for rows and columns
753 * and the weight for all points
754 */
756 {
757  float wc, wp, w, hp = piv->f;
758  smpoly a;
759  int i, f, e = crd;
760 
761  wp = 0.0;
762  for (i=tored; i; i--) wrw[i] = 0.0; // ???
763  for (i=act; i; i--)
764  {
765  wc = 0.0;
766  a = m_act[i];
767  loop
768  {
769  if (a->pos > tored)
770  break;
771  w = a->f;
772  f = a->e;
773  if (f < e)
774  {
775  w *= hp;
776  if (f) w /= m_res[f]->f;
777  }
778  wc += w;
779  wrw[a->pos] += w;
780  a = a->n;
781  if (a == NULL)
782  break;
783  }
784  wp += wc;
785  wcl[i] = wc;
786  }
787  wpoints = wp;
788 }
789 
790 /*
791 * compute pivot
792 */
794 {
795  float wopt = 1.0e30, hp = piv->f;
796  float wc, wr, wp, w;
797  smpoly a;
798  int i, copt, ropt, f, e = crd;
799 
800  this->smNewWeights();
801  for (i=act; i; i--)
802  {
803  a = m_act[i];
804  loop
805  {
806  if (a->pos > tored)
807  break;
808  w = a->f;
809  f = a->e;
810  if (f < e)
811  {
812  w *= hp;
813  if (f) w /= m_res[f]->f;
814  }
815  wc = wcl[i]-w;
816  wr = wrw[a->pos]-w;
817  if ((wr<0.25) || (wc<0.25)) // row or column with only one point
818  {
819  if (w<wopt)
820  {
821  wopt = w;
822  copt = i;
823  ropt = a->pos;
824  }
825  }
826  else // elimination
827  {
828  wp = w*(wpoints-wcl[i]-wr);
829  wp += wr*wc;
830  if (wp < wopt)
831  {
832  wopt = wp;
833  copt = i;
834  ropt = a->pos;
835  }
836  }
837  a = a->n;
838  if (a == NULL)
839  break;
840  }
841  }
842  rpiv = ropt;
843  cpiv = copt;
844  if (cpiv != act)
845  {
846  a = m_act[act];
847  m_act[act] = m_act[cpiv];
848  m_act[cpiv] = a;
849  }
850 }
851 
852 /* ----------------- elimination ------------------ */
853 
854 /* first step of elimination */
856 {
857  poly p = piv->m; // pivotelement
858  smpoly c = m_act[act]; // pivotcolumn
859  smpoly r = red; // row to reduce
860  smpoly res, a, b;
861  poly w, ha, hb;
862 
863  if ((c == NULL) || (r == NULL))
864  {
865  while (r!=NULL) sm_ElemDelete(&r,_R);
866  return;
867  }
868  do
869  {
870  a = m_act[r->pos];
871  res = dumm;
872  res->n = NULL;
873  b = c;
874  w = r->m;
875  loop // combine the chains a and b: p*a + w*b
876  {
877  if (a == NULL)
878  {
879  do
880  {
881  res = res->n = smElemCopy(b);
882  res->m = pp_Mult_qq(b->m, w,_R);
883  res->e = 1;
884  res->f = sm_PolyWeight(res,_R);
885  b = b->n;
886  } while (b != NULL);
887  break;
888  }
889  if (a->pos < b->pos)
890  {
891  res = res->n = a;
892  a = a->n;
893  }
894  else if (a->pos > b->pos)
895  {
896  res = res->n = smElemCopy(b);
897  res->m = pp_Mult_qq(b->m, w,_R);
898  res->e = 1;
899  res->f = sm_PolyWeight(res,_R);
900  b = b->n;
901  }
902  else
903  {
904  ha = pp_Mult_qq(a->m, p,_R);
905  p_Delete(&a->m,_R);
906  hb = pp_Mult_qq(b->m, w,_R);
907  ha = p_Add_q(ha, hb,_R);
908  if (ha != NULL)
909  {
910  a->m = ha;
911  a->e = 1;
912  a->f = sm_PolyWeight(a,_R);
913  res = res->n = a;
914  a = a->n;
915  }
916  else
917  {
918  sm_ElemDelete(&a,_R);
919  }
920  b = b->n;
921  }
922  if (b == NULL)
923  {
924  res->n = a;
925  break;
926  }
927  }
928  m_act[r->pos] = dumm->n;
929  sm_ElemDelete(&r,_R);
930  } while (r != NULL);
931 }
932 
933 /* higher steps of elimination */
935 {
936  poly hp = this->smMultPoly(piv);
937  poly gp = piv->m; // pivotelement
938  smpoly c = m_act[act]; // pivotcolumn
939  smpoly r = red; // row to reduce
940  smpoly res, a, b;
941  poly ha, hr, x, y;
942  int e, ip, ir, ia;
943 
944  if ((c == NULL) || (r == NULL))
945  {
946  while(r!=NULL) sm_ElemDelete(&r,_R);
947  p_Delete(&hp,_R);
948  return;
949  }
950  e = crd+1;
951  ip = piv->e;
952  do
953  {
954  a = m_act[r->pos];
955  res = dumm;
956  res->n = NULL;
957  b = c;
958  hr = r->m;
959  ir = r->e;
960  loop // combine the chains a and b: (hp,gp)*a(l) + hr*b(h)
961  {
962  if (a == NULL)
963  {
964  do
965  {
966  res = res->n = smElemCopy(b);
967  x = SM_MULT(b->m, hr, m_res[ir]->m,_R);
968  b = b->n;
969  if(ir) SM_DIV(x, m_res[ir]->m,_R);
970  res->m = x;
971  res->e = e;
972  res->f = sm_PolyWeight(res,_R);
973  } while (b != NULL);
974  break;
975  }
976  if (a->pos < b->pos)
977  {
978  res = res->n = a;
979  a = a->n;
980  }
981  else if (a->pos > b->pos)
982  {
983  res = res->n = smElemCopy(b);
984  x = SM_MULT(b->m, hr, m_res[ir]->m,_R);
985  b = b->n;
986  if(ir) SM_DIV(x, m_res[ir]->m,_R);
987  res->m = x;
988  res->e = e;
989  res->f = sm_PolyWeight(res,_R);
990  }
991  else
992  {
993  ha = a->m;
994  ia = a->e;
995  if (ir >= ia)
996  {
997  if (ir > ia)
998  {
999  x = SM_MULT(ha, m_res[ir]->m, m_res[ia]->m,_R);
1000  p_Delete(&ha,_R);
1001  ha = x;
1002  if (ia) SM_DIV(ha, m_res[ia]->m,_R);
1003  ia = ir;
1004  }
1005  x = SM_MULT(ha, gp, m_res[ia]->m,_R);
1006  p_Delete(&ha,_R);
1007  y = SM_MULT(b->m, hr, m_res[ia]->m,_R);
1008  }
1009  else if (ir >= ip)
1010  {
1011  if (ia < crd)
1012  {
1013  x = SM_MULT(ha, m_res[crd]->m, m_res[ia]->m,_R);
1014  p_Delete(&ha,_R);
1015  ha = x;
1016  SM_DIV(ha, m_res[ia]->m,_R);
1017  }
1018  y = hp;
1019  if(ir > ip)
1020  {
1021  y = SM_MULT(y, m_res[ir]->m, m_res[ip]->m,_R);
1022  if (ip) SM_DIV(y, m_res[ip]->m,_R);
1023  }
1024  ia = ir;
1025  x = SM_MULT(ha, y, m_res[ia]->m,_R);
1026  if (y != hp) p_Delete(&y,_R);
1027  p_Delete(&ha,_R);
1028  y = SM_MULT(b->m, hr, m_res[ia]->m,_R);
1029  }
1030  else
1031  {
1032  x = SM_MULT(hr, m_res[ia]->m, m_res[ir]->m,_R);
1033  if (ir) SM_DIV(x, m_res[ir]->m,_R);
1034  y = SM_MULT(b->m, x, m_res[ia]->m,_R);
1035  p_Delete(&x,_R);
1036  x = SM_MULT(ha, gp, m_res[ia]->m,_R);
1037  p_Delete(&ha,_R);
1038  }
1039  ha = p_Add_q(x, y,_R);
1040  if (ha != NULL)
1041  {
1042  if (ia) SM_DIV(ha, m_res[ia]->m,_R);
1043  a->m = ha;
1044  a->e = e;
1045  a->f = sm_PolyWeight(a,_R);
1046  res = res->n = a;
1047  a = a->n;
1048  }
1049  else
1050  {
1051  a->m = NULL;
1052  sm_ElemDelete(&a,_R);
1053  }
1054  b = b->n;
1055  }
1056  if (b == NULL)
1057  {
1058  res->n = a;
1059  break;
1060  }
1061  }
1062  m_act[r->pos] = dumm->n;
1063  sm_ElemDelete(&r,_R);
1064  } while (r != NULL);
1065  p_Delete(&hp,_R);
1066 }
1067 
1068 /* ----------------- transfer ------------------ */
1069 
1070 /*
1071 * select the pivotrow and store it to red and piv
1072 */
1074 {
1075  smpoly b = dumm;
1076  smpoly a, ap;
1077  int i;
1078 
1079  if (TEST_OPT_PROT)
1080  {
1081  if ((crd+1)%10)
1082  PrintS(".");
1083  else
1084  PrintS(".\n");
1085  }
1086  a = m_act[act];
1087  if (a->pos < rpiv)
1088  {
1089  do
1090  {
1091  ap = a;
1092  a = a->n;
1093  } while (a->pos < rpiv);
1094  ap->n = a->n;
1095  }
1096  else
1097  m_act[act] = a->n;
1098  piv = a;
1099  a->n = NULL;
1100  for (i=1; i<act; i++)
1101  {
1102  a = m_act[i];
1103  if (a->pos < rpiv)
1104  {
1105  loop
1106  {
1107  ap = a;
1108  a = a->n;
1109  if ((a == NULL) || (a->pos > rpiv))
1110  break;
1111  if (a->pos == rpiv)
1112  {
1113  ap->n = a->n;
1114  a->m = p_Neg(a->m,_R);
1115  b = b->n = a;
1116  b->pos = i;
1117  break;
1118  }
1119  }
1120  }
1121  else if (a->pos == rpiv)
1122  {
1123  m_act[i] = a->n;
1124  a->m = p_Neg(a->m,_R);
1125  b = b->n = a;
1126  b->pos = i;
1127  }
1128  }
1129  b->n = NULL;
1130  red = dumm->n;
1131 }
1132 
1133 /*
1134 * store the pivotcol in m_row
1135 * m_act[cols][pos(rows)] => m_row[rows][pos(cols)]
1136 */
1138 {
1139  smpoly c = m_act[act];
1140  smpoly h;
1141 
1142  while (c != NULL)
1143  {
1144  h = c;
1145  c = c->n;
1146  h->n = m_row[h->pos];
1147  m_row[h->pos] = h;
1148  h->pos = crd;
1149  }
1150 }
1151 
1152 /*
1153 * store the pivot and the assosiated row in m_row
1154 * to m_res (result):
1155 * piv + m_row[rows][pos(cols)] => m_res[cols][pos(rows)]
1156 */
1158 {
1159  smpoly r = m_row[rpiv];
1160  smpoly a, ap, h;
1161 
1162  m_row[rpiv] = NULL;
1163  perm[crd] = rpiv;
1164  piv->pos = crd;
1165  m_res[crd] = piv;
1166  while (r != NULL)
1167  {
1168  ap = m_res[r->pos];
1169  loop
1170  {
1171  a = ap->n;
1172  if (a == NULL)
1173  {
1174  ap->n = h = r;
1175  r = r->n;
1176  h->n = a;
1177  h->pos = crd;
1178  break;
1179  }
1180  ap = a;
1181  }
1182  }
1183 }
1184 
1185 /* ----------------- C++ stuff ------------------ */
1186 
1187 /*
1188 * clean m_act from zeros (after elim)
1189 */
1191 {
1192  int i = 0;
1193  int j;
1194 
1195  loop
1196  {
1197  i++;
1198  if (i > act) return;
1199  if (m_act[i] == NULL) break;
1200  }
1201  j = i;
1202  loop
1203  {
1204  j++;
1205  if (j > act) break;
1206  if (m_act[j] != NULL)
1207  {
1208  m_act[i] = m_act[j];
1209  i++;
1210  }
1211  }
1212  act -= (j-i);
1213  sign = 0;
1214 }
1215 
1216 /*
1217 * clean m_act from cols not to reduced (after elim)
1218 * put them to m_res
1219 */
1221 {
1222  int i = 0;
1223  int j;
1224 
1225  loop
1226  {
1227  i++;
1228  if (i > act) return;
1229  if (m_act[i]->pos > tored)
1230  {
1231  m_res[inred] = m_act[i];
1232  inred--;
1233  break;
1234  }
1235  }
1236  j = i;
1237  loop
1238  {
1239  j++;
1240  if (j > act) break;
1241  if (m_act[j]->pos > tored)
1242  {
1243  m_res[inred] = m_act[j];
1244  inred--;
1245  }
1246  else
1247  {
1248  m_act[i] = m_act[j];
1249  i++;
1250  }
1251  }
1252  act -= (j-i);
1253  sign = 0;
1254 }
1255 
1256 /*
1257 * copy m_act to m_res
1258 */
1260 {
1261  smpoly a,ap,r,h;
1262  int i,j,k,l;
1263 
1264  i = 0;
1265  if (act)
1266  {
1267  a = m_act[act]; // init perm
1268  do
1269  {
1270  i++;
1271  perm[crd+i] = a->pos;
1272  a = a->n;
1273  } while ((a != NULL) && (a->pos <= tored));
1274  for (j=act-1;j;j--) // load all positions of perm
1275  {
1276  a = m_act[j];
1277  k = 1;
1278  loop
1279  {
1280  if (perm[crd+k] >= a->pos)
1281  {
1282  if (perm[crd+k] > a->pos)
1283  {
1284  for (l=i;l>=k;l--) perm[crd+l+1] = perm[crd+l];
1285  perm[crd+k] = a->pos;
1286  i++;
1287  }
1288  a = a->n;
1289  if ((a == NULL) || (a->pos > tored)) break;
1290  }
1291  k++;
1292  if ((k > i) && (a->pos <= tored))
1293  {
1294  do
1295  {
1296  i++;
1297  perm[crd+i] = a->pos;
1298  a = a->n;
1299  } while ((a != NULL) && (a->pos <= tored));
1300  break;
1301  }
1302  }
1303  }
1304  }
1305  for (j=act;j;j--) // renumber m_act
1306  {
1307  k = 1;
1308  a = m_act[j];
1309  while ((a != NULL) && (a->pos <= tored))
1310  {
1311  if (perm[crd+k] == a->pos)
1312  {
1313  a->pos = crd+k;
1314  a = a->n;
1315  }
1316  k++;
1317  }
1318  }
1319  tored = crd+i;
1320  for(k=1;k<=i;k++) // clean this from m_row
1321  {
1322  j = perm[crd+k];
1323  if (m_row[j] != NULL)
1324  {
1325  r = m_row[j];
1326  m_row[j] = NULL;
1327  do
1328  {
1329  ap = m_res[r->pos];
1330  loop
1331  {
1332  a = ap->n;
1333  if (a == NULL)
1334  {
1335  h = ap->n = r;
1336  r = r->n;
1337  h->n = NULL;
1338  h->pos = crd+k;
1339  break;
1340  }
1341  ap = a;
1342  }
1343  } while (r!=NULL);
1344  }
1345  }
1346  while(act) // clean m_act
1347  {
1348  crd++;
1349  m_res[crd] = m_act[act];
1350  act--;
1351  }
1352  for (i=1;i<=tored;i++) // take the rest of m_row
1353  {
1354  if(m_row[i] != NULL)
1355  {
1356  tored++;
1357  r = m_row[i];
1358  m_row[i] = NULL;
1359  perm[tored] = i;
1360  do
1361  {
1362  ap = m_res[r->pos];
1363  loop
1364  {
1365  a = ap->n;
1366  if (a == NULL)
1367  {
1368  h = ap->n = r;
1369  r = r->n;
1370  h->n = NULL;
1371  h->pos = tored;
1372  break;
1373  }
1374  ap = a;
1375  }
1376  } while (r!=NULL);
1377  }
1378  }
1379  for (i=tored+1;i<=nrows;i++) // take the rest of m_row
1380  {
1381  if(m_row[i] != NULL)
1382  {
1383  r = m_row[i];
1384  m_row[i] = NULL;
1385  do
1386  {
1387  ap = m_res[r->pos];
1388  loop
1389  {
1390  a = ap->n;
1391  if (a == NULL)
1392  {
1393  h = ap->n = r;
1394  r = r->n;
1395  h->n = NULL;
1396  h->pos = i;
1397  break;
1398  }
1399  ap = a;
1400  }
1401  } while (r!=NULL);
1402  }
1403  }
1404  while (inred < ncols) // take unreducable
1405  {
1406  crd++;
1407  inred++;
1408  m_res[crd] = m_res[inred];
1409  }
1410 }
1411 
1412 /*
1413 * multiply and divide the column, that goes to result
1414 */
1416 {
1417  smpoly a = m_act[act];
1418  int e = crd;
1419  poly ha;
1420  int f;
1421 
1422  while (a != NULL)
1423  {
1424  f = a->e;
1425  if (f < e)
1426  {
1427  ha = SM_MULT(a->m, m_res[e]->m, m_res[f]->m,_R);
1428  p_Delete(&a->m,_R);
1429  if (f) SM_DIV(ha, m_res[f]->m,_R);
1430  a->m = ha;
1431  if (normalize) p_Normalize(a->m,_R);
1432  }
1433  a = a->n;
1434  }
1435 }
1436 
1437 /*
1438 * multiply and divide the m_act finaly
1439 */
1441 {
1442  smpoly a;
1443  poly ha;
1444  int i, f;
1445  int e = crd;
1446 
1447  for (i=act; i; i--)
1448  {
1449  a = m_act[i];
1450  do
1451  {
1452  f = a->e;
1453  if (f < e)
1454  {
1455  ha = SM_MULT(a->m, m_res[e]->m, m_res[f]->m, _R);
1456  p_Delete(&a->m,_R);
1457  if (f) SM_DIV(ha, m_res[f]->m, _R);
1458  a->m = ha;
1459  }
1460  if (normalize) p_Normalize(a->m, _R);
1461  a = a->n;
1462  } while (a != NULL);
1463  }
1464 }
1465 
1466 /*
1467 * check for denominators
1468 */
1470 {
1471  int i;
1472  smpoly a;
1473 
1474  for (i=act; i; i--)
1475  {
1476  a = m_act[i];
1477  do
1478  {
1479  if(sm_HaveDenom(a->m,_R)) return 1;
1480  a = a->n;
1481  } while (a != NULL);
1482  }
1483  return 0;
1484 }
1485 
1486 /*
1487 * normalize
1488 */
1490 {
1491  smpoly a;
1492  int i;
1493  int e = crd;
1494 
1495  for (i=act; i; i--)
1496  {
1497  a = m_act[i];
1498  do
1499  {
1500  if (e == a->e) p_Normalize(a->m,_R);
1501  a = a->n;
1502  } while (a != NULL);
1503  }
1504 }
1505 
1506 /*
1507 * multiply and divide the element, save poly
1508 */
1510 {
1511  int f = a->e;
1512  poly r, h;
1513 
1514  if (f < crd)
1515  {
1516  h = r = a->m;
1517  h = SM_MULT(h, m_res[crd]->m, m_res[f]->m, _R);
1518  if (f) SM_DIV(h, m_res[f]->m, _R);
1519  a->m = h;
1520  if (normalize) p_Normalize(a->m,_R);
1521  a->f = sm_PolyWeight(a,_R);
1522  return r;
1523  }
1524  else
1525  return NULL;
1526 }
1527 
1528 /*
1529 * delete the m_act finaly
1530 */
1532 {
1533  smpoly a;
1534  int i;
1535 
1536  for (i=act; i; i--)
1537  {
1538  a = m_act[i];
1539  do
1540  {
1541  sm_ElemDelete(&a,_R);
1542  } while (a != NULL);
1543  }
1544 }
1545 
1546 /*
1547 * delete the pivotcol
1548 */
1550 {
1551  smpoly a = m_act[act];
1552 
1553  while (a != NULL)
1554  {
1555  sm_ElemDelete(&a,_R);
1556  }
1557 }
1558 
1559 /*
1560 * delete pivot elements
1561 */
1563 {
1564  int i=crd;
1565 
1566  while (i != 0)
1567  {
1568  sm_ElemDelete(&m_res[i],_R);
1569  i--;
1570  }
1571 }
1572 
1573 /*
1574 * the sign of the determinant
1575 */
1577 {
1578  int j,i;
1579  if (act > 2)
1580  {
1581  if (cpiv!=act) sign=-sign;
1582  if ((act%2)==0) sign=-sign;
1583  i=1;
1584  j=perm[1];
1585  while(j<rpiv)
1586  {
1587  sign=-sign;
1588  i++;
1589  j=perm[i];
1590  }
1591  while(perm[i]!=0)
1592  {
1593  perm[i]=perm[i+1];
1594  i++;
1595  }
1596  }
1597  else
1598  {
1599  if (cpiv!=1) sign=-sign;
1600  if (rpiv!=perm[1]) sign=-sign;
1601  }
1602 }
1603 
1605 {
1606  int i;
1607  for (i=act;i;i--) perm[i]=i;
1608 }
1609 
1610 /* ----------------- arithmetic ------------------ */
1611 /*2
1612 * exact division a/b
1613 * a destroyed, b NOT destroyed
1614 */
1615 void sm_PolyDiv(poly a, poly b, const ring R)
1616 {
1617  const number x = pGetCoeff(b);
1618  number y, yn;
1619  poly t, h, dummy;
1620  int i;
1621 
1622  if (pNext(b) == NULL)
1623  {
1624  do
1625  {
1626  if (!p_LmIsConstantComp(b,R))
1627  {
1628  for (i=rVar(R); i; i--)
1629  p_SubExp(a,i,p_GetExp(b,i,R),R);
1630  p_Setm(a,R);
1631  }
1632  y = n_Div(pGetCoeff(a),x,R->cf);
1633  n_Normalize(y,R->cf);
1634  p_SetCoeff(a,y,R);
1635  pIter(a);
1636  } while (a != NULL);
1637  return;
1638  }
1639  dummy = p_Init(R);
1640  do
1641  {
1642  for (i=rVar(R); i; i--)
1643  p_SubExp(a,i,p_GetExp(b,i,R),R);
1644  p_Setm(a,R);
1645  y = n_Div(pGetCoeff(a),x,R->cf);
1646  n_Normalize(y,R->cf);
1647  p_SetCoeff(a,y,R);
1648  yn = n_InpNeg(n_Copy(y,R->cf),R->cf);
1649  t = pNext(b);
1650  h = dummy;
1651  do
1652  {
1653  h = pNext(h) = p_Init(R);
1654  //pSetComp(h,0);
1655  for (i=rVar(R); i; i--)
1656  p_SetExp(h,i,p_GetExp(a,i,R)+p_GetExp(t,i,R),R);
1657  p_Setm(h,R);
1658  pSetCoeff0(h,n_Mult(yn, pGetCoeff(t),R->cf));
1659  pIter(t);
1660  } while (t != NULL);
1661  n_Delete(&yn,R->cf);
1662  pNext(h) = NULL;
1663  a = pNext(a) = p_Add_q(pNext(a), pNext(dummy),R);
1664  } while (a!=NULL);
1665  p_LmFree(dummy, R);
1666 }
1667 
1668 
1669 //disable that, as it fails with coef buckets
1670 //#define X_MAS
1671 #ifdef X_MAS
1672 // Note: the following in not addapted to SW :(
1673 /*
1674 /// returns the part of (a*b)/exp(lead(c)) with nonegative exponents
1675 poly smMultDiv(poly a, poly b, const poly c)
1676 {
1677  poly pa, e, res, r;
1678  BOOLEAN lead;
1679  int la, lb, lr;
1680 
1681  if ((c == NULL) || pLmIsConstantComp(c))
1682  {
1683  return pp_Mult_qq(a, b);
1684  }
1685 
1686  pqLength(a, b, la, lb, SM_MIN_LENGTH_BUCKET);
1687 
1688  // we iter over b, so, make sure b is the shortest
1689  // such that we minimize our iterations
1690  if (lb > la)
1691  {
1692  r = a;
1693  a = b;
1694  b = r;
1695  lr = la;
1696  la = lb;
1697  lb = lr;
1698  }
1699  res = NULL;
1700  e = p_Init();
1701  lead = FALSE;
1702  while (!lead)
1703  {
1704  pSetCoeff0(e,pGetCoeff(b));
1705  if (smIsNegQuot(e, b, c))
1706  {
1707  lead = pLmDivisibleByNoComp(e, a);
1708  r = smSelectCopy_ExpMultDiv(a, e, b, c);
1709  }
1710  else
1711  {
1712  lead = TRUE;
1713  r = pp_Mult__mm(a, e);
1714  }
1715  if (lead)
1716  {
1717  if (res != NULL)
1718  {
1719  smFindRef(&pa, &res, r);
1720  if (pa == NULL)
1721  lead = FALSE;
1722  }
1723  else
1724  {
1725  pa = res = r;
1726  }
1727  }
1728  else
1729  res = p_Add_q(res, r);
1730  pIter(b);
1731  if (b == NULL)
1732  {
1733  pLmFree(e);
1734  return res;
1735  }
1736  }
1737 
1738  if (!TEST_OPT_NOT_BUCKETS && lb >= SM_MIN_LENGTH_BUCKET)
1739  {
1740  // use buckets if minimum length is smaller than threshold
1741  spolyrec rp;
1742  poly append;
1743  // find the last monomial before pa
1744  if (res == pa)
1745  {
1746  append = &rp;
1747  pNext(append) = res;
1748  }
1749  else
1750  {
1751  append = res;
1752  while (pNext(append) != pa)
1753  {
1754  assume(pNext(append) != NULL);
1755  pIter(append);
1756  }
1757  }
1758  kBucket_pt bucket = kBucketCreate(currRing);
1759  kBucketInit(bucket, pNext(append), 0);
1760  do
1761  {
1762  pSetCoeff0(e,pGetCoeff(b));
1763  if (smIsNegQuot(e, b, c))
1764  {
1765  lr = la;
1766  r = pp_Mult_Coeff_mm_DivSelect_MultDiv(a, lr, e, b, c);
1767  if (pLmDivisibleByNoComp(e, a))
1768  append = kBucket_ExtractLarger_Add_q(bucket, append, r, &lr);
1769  else
1770  kBucket_Add_q(bucket, r, &lr);
1771  }
1772  else
1773  {
1774  r = pp_Mult__mm(a, e);
1775  append = kBucket_ExtractLarger_Add_q(bucket, append, r, &la);
1776  }
1777  pIter(b);
1778  } while (b != NULL);
1779  pNext(append) = kBucketClear(bucket);
1780  kBucketDestroy(&bucket);
1781  pTest(append);
1782  }
1783  else
1784  {
1785  // use old sm stuff
1786  do
1787  {
1788  pSetCoeff0(e,pGetCoeff(b));
1789  if (smIsNegQuot(e, b, c))
1790  {
1791  r = smSelectCopy_ExpMultDiv(a, e, b, c);
1792  if (pLmDivisibleByNoComp(e, a))
1793  smCombineChain(&pa, r);
1794  else
1795  pa = p_Add_q(pa,r);
1796  }
1797  else
1798  {
1799  r = pp_Mult__mm(a, e);
1800  smCombineChain(&pa, r);
1801  }
1802  pIter(b);
1803  } while (b != NULL);
1804  }
1805  pLmFree(e);
1806  return res;
1807 }
1808 */
1809 #else
1810 
1811 /*
1812 * returns the part of (a*b)/exp(lead(c)) with nonegative exponents
1813 */
1814 poly sm_MultDiv(poly a, poly b, const poly c, const ring R)
1815 {
1816  poly pa, e, res, r;
1817  BOOLEAN lead;
1818 
1819  if ((c == NULL) || p_LmIsConstantComp(c,R))
1820  {
1821  return pp_Mult_qq(a, b, R);
1822  }
1823  if (smSmaller(a, b))
1824  {
1825  r = a;
1826  a = b;
1827  b = r;
1828  }
1829  res = NULL;
1830  e = p_Init(R);
1831  lead = FALSE;
1832  while (!lead)
1833  {
1834  pSetCoeff0(e,pGetCoeff(b));
1835  if (sm_IsNegQuot(e, b, c, R))
1836  {
1837  lead = p_LmDivisibleByNoComp(e, a, R);
1838  r = sm_SelectCopy_ExpMultDiv(a, e, b, c, R);
1839  }
1840  else
1841  {
1842  lead = TRUE;
1843  r = pp_Mult_mm(a, e,R);
1844  }
1845  if (lead)
1846  {
1847  if (res != NULL)
1848  {
1849  sm_FindRef(&pa, &res, r, R);
1850  if (pa == NULL)
1851  lead = FALSE;
1852  }
1853  else
1854  {
1855  pa = res = r;
1856  }
1857  }
1858  else
1859  res = p_Add_q(res, r, R);
1860  pIter(b);
1861  if (b == NULL)
1862  {
1863  p_LmFree(e, R);
1864  return res;
1865  }
1866  }
1867  do
1868  {
1869  pSetCoeff0(e,pGetCoeff(b));
1870  if (sm_IsNegQuot(e, b, c, R))
1871  {
1872  r = sm_SelectCopy_ExpMultDiv(a, e, b, c, R);
1873  if (p_LmDivisibleByNoComp(e, a,R))
1874  sm_CombineChain(&pa, r, R);
1875  else
1876  pa = p_Add_q(pa,r,R);
1877  }
1878  else
1879  {
1880  r = pp_Mult_mm(a, e, R);
1881  sm_CombineChain(&pa, r, R);
1882  }
1883  pIter(b);
1884  } while (b != NULL);
1885  p_LmFree(e, R);
1886  return res;
1887 }
1888 #endif
1889 /*n
1890 * exact division a/b
1891 * a is a result of smMultDiv
1892 * a destroyed, b NOT destroyed
1893 */
1894 void sm_SpecialPolyDiv(poly a, poly b, const ring R)
1895 {
1896  if (pNext(b) == NULL)
1897  {
1898  sm_PolyDivN(a, pGetCoeff(b),R);
1899  return;
1900  }
1901  sm_ExactPolyDiv(a, b, R);
1902 }
1903 
1904 
1905 /* ------------ internals arithmetic ------------- */
1906 static void sm_ExactPolyDiv(poly a, poly b, const ring R)
1907 {
1908  const number x = pGetCoeff(b);
1909  poly tail = pNext(b), e = p_Init(R);
1910  poly h;
1911  number y, yn;
1912  int lt = pLength(tail);
1913 
1914  if (lt + 1 >= SM_MIN_LENGTH_BUCKET && !TEST_OPT_NOT_BUCKETS)
1915  {
1917  kBucketInit(bucket, pNext(a), 0);
1918  int lh = 0;
1919  do
1920  {
1921  y = n_Div(pGetCoeff(a), x, R->cf);
1922  n_Normalize(y, R->cf);
1923  p_SetCoeff(a,y, R);
1924  yn = n_InpNeg(n_Copy(y, R->cf), R->cf);
1925  pSetCoeff0(e,yn);
1926  lh = lt;
1927  if (sm_IsNegQuot(e, a, b, R))
1928  {
1929  h = pp_Mult_Coeff_mm_DivSelect_MultDiv(tail, lh, e, a, b, R);
1930  }
1931  else
1932  h = pp_Mult_mm(tail, e, R);
1933  n_Delete(&yn, R->cf);
1934  kBucket_Add_q(bucket, h, &lh);
1935 
1936  a = pNext(a) = kBucketExtractLm(bucket);
1937  } while (a!=NULL);
1938  kBucketDestroy(&bucket);
1939  }
1940  else
1941  {
1942  do
1943  {
1944  y = n_Div(pGetCoeff(a), x, R->cf);
1945  n_Normalize(y, R->cf);
1946  p_SetCoeff(a,y, R);
1947  yn = n_InpNeg(n_Copy(y, R->cf), R->cf);
1948  pSetCoeff0(e,yn);
1949  if (sm_IsNegQuot(e, a, b, R))
1950  h = sm_SelectCopy_ExpMultDiv(tail, e, a, b, R);
1951  else
1952  h = pp_Mult_mm(tail, e, R);
1953  n_Delete(&yn, R->cf);
1954  a = pNext(a) = p_Add_q(pNext(a), h, R);
1955  } while (a!=NULL);
1956  }
1957  p_LmFree(e, R);
1958 }
1959 
1960 // obachman --> Wilfried: check the following
1961 static BOOLEAN sm_IsNegQuot(poly a, const poly b, const poly c, const ring R)
1962 {
1963  if (p_LmDivisibleByNoComp(c, b,R))
1964  {
1965  p_ExpVectorDiff(a, b, c,R);
1966  // Hmm: here used to be a p_Setm(a): but it is unnecessary,
1967  // if b and c are correct
1968  return FALSE;
1969  }
1970  else
1971  {
1972  int i;
1973  for (i=rVar(R); i>0; i--)
1974  {
1975  if(p_GetExp(c,i,R) > p_GetExp(b,i,R))
1976  p_SetExp(a,i,p_GetExp(c,i,R)-p_GetExp(b,i,R),R);
1977  else
1978  p_SetExp(a,i,0,R);
1979  }
1980  // here we actually might need a p_Setm, if a is to be used in
1981  // comparisons
1982  return TRUE;
1983  }
1984 }
1985 
1986 static void sm_ExpMultDiv(poly t, const poly b, const poly c, const ring R)
1987 {
1988  p_Test(t,R);
1989  p_LmTest(b,R);
1990  p_LmTest(c,R);
1991  poly bc = p_New(R);
1992 
1993  p_ExpVectorDiff(bc, b, c, R);
1994 
1995  while(t!=NULL)
1996  {
1997  p_ExpVectorAdd(t, bc, R);
1998  pIter(t);
1999  }
2000  p_LmFree(bc, R);
2001 }
2002 
2003 
2004 static void sm_PolyDivN(poly a, const number x, const ring R)
2005 {
2006  number y;
2007 
2008  do
2009  {
2010  y = n_Div(pGetCoeff(a),x, R->cf);
2011  n_Normalize(y, R->cf);
2012  p_SetCoeff(a,y, R);
2013  pIter(a);
2014  } while (a != NULL);
2015 }
2016 
2018 {
2019  loop
2020  {
2021  pIter(b);
2022  if (b == NULL) return TRUE;
2023  pIter(a);
2024  if (a == NULL) return FALSE;
2025  }
2026 }
2027 
2028 static void sm_CombineChain(poly *px, poly r, const ring R)
2029 {
2030  poly pa = *px, pb;
2031  number x;
2032  int i;
2033 
2034  loop
2035  {
2036  pb = pNext(pa);
2037  if (pb == NULL)
2038  {
2039  pa = pNext(pa) = r;
2040  break;
2041  }
2042  i = p_LmCmp(pb, r,R);
2043  if (i > 0)
2044  pa = pb;
2045  else
2046  {
2047  if (i == 0)
2048  {
2049  x = n_Add(pGetCoeff(pb), pGetCoeff(r),R->cf);
2050  p_LmDelete(&r,R);
2051  if (n_IsZero(x,R->cf))
2052  {
2053  p_LmDelete(&pb,R);
2054  pNext(pa) = p_Add_q(pb,r,R);
2055  }
2056  else
2057  {
2058  pa = pb;
2059  p_SetCoeff(pa,x,R);
2060  pNext(pa) = p_Add_q(pNext(pa), r, R);
2061  }
2062  }
2063  else
2064  {
2065  pa = pNext(pa) = r;
2066  pNext(pa) = p_Add_q(pb, pNext(pa),R);
2067  }
2068  break;
2069  }
2070  }
2071  *px = pa;
2072 }
2073 
2074 
2075 static void sm_FindRef(poly *ref, poly *px, poly r, const ring R)
2076 {
2077  number x;
2078  int i;
2079  poly pa = *px, pp = NULL;
2080 
2081  loop
2082  {
2083  i = p_LmCmp(pa, r,R);
2084  if (i > 0)
2085  {
2086  pp = pa;
2087  pIter(pa);
2088  if (pa==NULL)
2089  {
2090  pNext(pp) = r;
2091  break;
2092  }
2093  }
2094  else
2095  {
2096  if (i == 0)
2097  {
2098  x = n_Add(pGetCoeff(pa), pGetCoeff(r),R->cf);
2099  p_LmDelete(&r,R);
2100  if (n_IsZero(x,R->cf))
2101  {
2102  p_LmDelete(&pa,R);
2103  if (pp!=NULL)
2104  pNext(pp) = p_Add_q(pa,r,R);
2105  else
2106  *px = p_Add_q(pa,r,R);
2107  }
2108  else
2109  {
2110  pp = pa;
2111  p_SetCoeff(pp,x,R);
2112  pNext(pp) = p_Add_q(pNext(pp), r, R);
2113  }
2114  }
2115  else
2116  {
2117  if (pp!=NULL)
2118  pp = pNext(pp) = r;
2119  else
2120  *px = pp = r;
2121  pNext(pp) = p_Add_q(pa, pNext(r),R);
2122  }
2123  break;
2124  }
2125  }
2126  *ref = pp;
2127 }
2128 
2129 /* ----------------- internal 'C' stuff ------------------ */
2130 
2131 static void sm_ElemDelete(smpoly *r, const ring R)
2132 {
2133  smpoly a = *r, b = a->n;
2134 
2135  p_Delete(&a->m, R);
2136  omFreeBin((void *)a, smprec_bin);
2137  *r = b;
2138 }
2139 
2141 {
2143  memcpy(r, a, sizeof(smprec));
2144 /* r->m = pCopy(r->m); */
2145  return r;
2146 }
2147 
2148 /*
2149 * from poly to smpoly
2150 * do not destroy p
2151 */
2152 static smpoly sm_Poly2Smpoly(poly q, const ring R)
2153 {
2154  poly pp;
2155  smpoly res, a;
2156  long x;
2157 
2158  if (q == NULL)
2159  return NULL;
2160  a = res = (smpoly)omAllocBin(smprec_bin);
2161  a->pos = x = p_GetComp(q,R);
2162  a->m = q;
2163  a->e = 0;
2164  loop
2165  {
2166  p_SetComp(q,0,R);
2167  pp = q;
2168  pIter(q);
2169  if (q == NULL)
2170  {
2171  a->n = NULL;
2172  return res;
2173  }
2174  if (p_GetComp(q,R) != x)
2175  {
2176  a = a->n = (smpoly)omAllocBin(smprec_bin);
2177  pNext(pp) = NULL;
2178  a->pos = x = p_GetComp(q,R);
2179  a->m = q;
2180  a->e = 0;
2181  }
2182  }
2183 }
2184 
2185 /*
2186 * from smpoly to poly
2187 * destroy a
2188 */
2189 static poly sm_Smpoly2Poly(smpoly a, const ring R)
2190 {
2191  smpoly b;
2192  poly res, pp, q;
2193  long x;
2194 
2195  if (a == NULL)
2196  return NULL;
2197  x = a->pos;
2198  q = res = a->m;
2199  loop
2200  {
2201  p_SetComp(q,x,R);
2202  pp = q;
2203  pIter(q);
2204  if (q == NULL)
2205  break;
2206  }
2207  loop
2208  {
2209  b = a;
2210  a = a->n;
2211  omFreeBin((void *)b, smprec_bin);
2212  if (a == NULL)
2213  return res;
2214  x = a->pos;
2215  q = pNext(pp) = a->m;
2216  loop
2217  {
2218  p_SetComp(q,x,R);
2219  pp = q;
2220  pIter(q);
2221  if (q == NULL)
2222  break;
2223  }
2224  }
2225 }
2226 
2227 /*
2228 * weigth of a polynomial, for pivot strategy
2229 */
2230 static float sm_PolyWeight(smpoly a, const ring R)
2231 {
2232  poly p = a->m;
2233  int i;
2234  float res = (float)n_Size(pGetCoeff(p),R->cf);
2235 
2236  if (pNext(p) == NULL)
2237  {
2238  for(i=rVar(R); i>0; i--)
2239  {
2240  if (p_GetExp(p,i,R) != 0) return res+1.0;
2241  }
2242  return res;
2243  }
2244  else
2245  {
2246  i = 0;
2247  res = 0.0;
2248  do
2249  {
2250  i++;
2251  res += (float)n_Size(pGetCoeff(p),R->cf);
2252  pIter(p);
2253  }
2254  while (p);
2255  return res+(float)i;
2256  }
2257 }
2258 
2259 static BOOLEAN sm_HaveDenom(poly a, const ring R)
2260 {
2261  BOOLEAN sw;
2262  number x;
2263 
2264  while (a != NULL)
2265  {
2266  x = n_GetDenom(pGetCoeff(a),R->cf);
2267  sw = n_IsOne(x,R->cf);
2268  n_Delete(&x,R->cf);
2269  if (!sw)
2270  {
2271  return TRUE;
2272  }
2273  pIter(a);
2274  }
2275  return FALSE;
2276 }
2277 
2278 static number sm_Cleardenom(ideal id, const ring R)
2279 {
2280  poly a;
2281  number x,y,res=n_Init(1,R->cf);
2282  BOOLEAN sw=FALSE;
2283 
2284  for (int i=0; i<IDELEMS(id); i++)
2285  {
2286  a = id->m[i];
2287  sw = sm_HaveDenom(a,R);
2288  if (sw) break;
2289  }
2290  if (!sw) return res;
2291  for (int i=0; i<IDELEMS(id); i++)
2292  {
2293  a = id->m[i];
2294  if (a!=NULL)
2295  {
2296  x = n_Copy(pGetCoeff(a),R->cf);
2297  p_Cleardenom(a, R);
2298  y = n_Div(x,pGetCoeff(a),R->cf);
2299  n_Delete(&x,R->cf);
2300  x = n_Mult(res,y,R->cf);
2301  n_Normalize(x,R->cf);
2302  n_Delete(&res,R->cf);
2303  res = x;
2304  }
2305  }
2306  return res;
2307 }
2308 
2309 /* ----------------- gauss elimination ------------------ */
2310 /* in structs.h */
2311 typedef struct smnrec sm_nrec;
2312 typedef sm_nrec * smnumber;
2313 struct smnrec{
2314  smnumber n; // the next element
2315  int pos; // position
2316  number m; // the element
2317 };
2318 
2320 
2321 /* declare internal 'C' stuff */
2322 static void sm_NumberDelete(smnumber *, const ring R);
2323 static smnumber smNumberCopy(smnumber);
2324 static smnumber sm_Poly2Smnumber(poly, const ring);
2325 static poly sm_Smnumber2Poly(number,const ring);
2326 static BOOLEAN smCheckSolv(ideal);
2327 
2328 /* class for sparse number matrix:
2329 */
2331 private:
2332  int nrows, ncols; // dimension of the problem
2333  int act; // number of unreduced columns (start: ncols)
2334  int crd; // number of reduced columns (start: 0)
2335  int tored; // border for rows to reduce
2336  int sing; // indicator for singular problem
2337  int rpiv; // row-position of the pivot
2338  int *perm; // permutation of rows
2339  number *sol; // field for solution
2340  int *wrw, *wcl; // weights of rows and columns
2341  smnumber * m_act; // unreduced columns
2342  smnumber * m_res; // reduced columns (result)
2343  smnumber * m_row; // reduced part of rows
2344  smnumber red; // row to reduce
2345  smnumber piv; // pivot
2346  smnumber dumm; // allocated dummy
2347  ring _R;
2348  void smColToRow();
2349  void smRowToCol();
2350  void smSelectPR();
2351  void smRealPivot();
2352  void smZeroToredElim();
2353  void smGElim();
2354  void smAllDel();
2355 public:
2356  sparse_number_mat(ideal, const ring);
2357  ~sparse_number_mat();
2358  int smIsSing() { return sing; }
2359  void smTriangular();
2360  void smSolv();
2361  ideal smRes2Ideal();
2362 };
2363 
2364 /* ----------------- basics (used from 'C') ------------------ */
2365 /*2
2366 * returns the solution of a linear equation
2367 * solution of M*x = r (M has dimension nXn) =>
2368 * I = module(transprose(M)) + r*gen(n+1)
2369 * uses Gauss-elimination
2370 */
2371 ideal sm_CallSolv(ideal I, const ring R)
2372 {
2373  sparse_number_mat *linsolv;
2374  ring tmpR;
2375  ideal rr;
2376 
2377  if (id_IsConstant(I,R)==FALSE)
2378  {
2379  WerrorS("symbol in equation");
2380  return NULL;
2381  }
2382  I->rank = id_RankFreeModule(I,R);
2383  if (smCheckSolv(I)) return NULL;
2384  tmpR=sm_RingChange(R,1);
2385  rr=idrCopyR(I,R, tmpR);
2386  linsolv = new sparse_number_mat(rr,tmpR);
2387  rr=NULL;
2388  linsolv->smTriangular();
2389  if (linsolv->smIsSing() == 0)
2390  {
2391  linsolv->smSolv();
2392  rr = linsolv->smRes2Ideal();
2393  }
2394  else
2395  WerrorS("singular problem for linsolv");
2396  delete linsolv;
2397  if (rr!=NULL)
2398  rr = idrMoveR(rr,tmpR,R);
2399  sm_KillModifiedRing(tmpR);
2400  return rr;
2401 }
2402 
2403 /*
2404 * constructor, destroy smat
2405 */
2407 {
2408  int i;
2409  poly* pmat;
2410  _R=R;
2411 
2412  crd = sing = 0;
2413  act = ncols = smat->ncols;
2414  tored = nrows = smat->rank;
2415  i = tored+1;
2416  perm = (int *)omAlloc(sizeof(int)*i);
2417  m_row = (smnumber *)omAlloc0(sizeof(smnumber)*i);
2418  wrw = (int *)omAlloc(sizeof(int)*i);
2419  i = ncols+1;
2420  wcl = (int *)omAlloc(sizeof(int)*i);
2421  m_act = (smnumber *)omAlloc(sizeof(smnumber)*i);
2422  m_res = (smnumber *)omAlloc0(sizeof(smnumber)*i);
2423  dumm = (smnumber)omAllocBin(smnrec_bin);
2424  pmat = smat->m;
2425  for(i=ncols; i; i--)
2426  {
2427  m_act[i] = sm_Poly2Smnumber(pmat[i-1],_R);
2428  }
2429  omFreeSize((ADDRESS)pmat,smat->ncols*sizeof(poly));
2431 }
2432 
2433 /*
2434 * destructor
2435 */
2437 {
2438  int i;
2439  omFreeBin((ADDRESS)dumm, smnrec_bin);
2440  i = ncols+1;
2441  omFreeSize((ADDRESS)m_res, sizeof(smnumber)*i);
2442  omFreeSize((ADDRESS)m_act, sizeof(smnumber)*i);
2443  omFreeSize((ADDRESS)wcl, sizeof(int)*i);
2444  i = nrows+1;
2445  omFreeSize((ADDRESS)wrw, sizeof(int)*i);
2446  omFreeSize((ADDRESS)m_row, sizeof(smnumber)*i);
2447  omFreeSize((ADDRESS)perm, sizeof(int)*i);
2448 }
2449 
2450 /*
2451 * triangularization by Gauss-elimination
2452 */
2454 {
2455  tored--;
2456  this->smZeroToredElim();
2457  if (sing != 0) return;
2458  while (act > 1)
2459  {
2460  this->smRealPivot();
2461  this->smSelectPR();
2462  this->smGElim();
2463  crd++;
2464  this->smColToRow();
2465  act--;
2466  this->smRowToCol();
2467  this->smZeroToredElim();
2468  if (sing != 0) return;
2469  }
2470  if (TEST_OPT_PROT) PrintS(".\n");
2471  piv = m_act[1];
2472  rpiv = piv->pos;
2473  m_act[1] = piv->n;
2474  piv->n = NULL;
2475  crd++;
2476  this->smColToRow();
2477  act--;
2478  this->smRowToCol();
2479 }
2480 
2481 /*
2482 * solve the triangular system
2483 */
2485 {
2486  int i, j;
2487  number x, y, z;
2488  smnumber s, d, r = m_row[nrows];
2489 
2490  m_row[nrows] = NULL;
2491  sol = (number *)omAlloc0(sizeof(number)*(crd+1));
2492  while (r != NULL) // expand the rigth hand side
2493  {
2494  sol[r->pos] = r->m;
2495  s = r;
2496  r = r->n;
2497  omFreeBin((ADDRESS)s, smnrec_bin);
2498  }
2499  i = crd; // solve triangular system
2500  if (sol[i] != NULL)
2501  {
2502  x = sol[i];
2503  sol[i] = n_Div(x, m_res[i]->m,_R->cf);
2504  n_Delete(&x,_R->cf);
2505  }
2506  i--;
2507  while (i > 0)
2508  {
2509  x = NULL;
2510  d = m_res[i];
2511  s = d->n;
2512  while (s != NULL)
2513  {
2514  j = s->pos;
2515  if (sol[j] != NULL)
2516  {
2517  z = n_Mult(sol[j], s->m,_R->cf);
2518  if (x != NULL)
2519  {
2520  y = x;
2521  x = n_Sub(y, z,_R->cf);
2522  n_Delete(&y,_R->cf);
2523  n_Delete(&z,_R->cf);
2524  }
2525  else
2526  x = n_InpNeg(z,_R->cf);
2527  }
2528  s = s->n;
2529  }
2530  if (sol[i] != NULL)
2531  {
2532  if (x != NULL)
2533  {
2534  y = n_Add(x, sol[i],_R->cf);
2535  n_Delete(&x,_R->cf);
2536  if (n_IsZero(y,_R->cf))
2537  {
2538  n_Delete(&sol[i],_R->cf);
2539  sol[i] = NULL;
2540  }
2541  else
2542  sol[i] = y;
2543  }
2544  }
2545  else
2546  sol[i] = x;
2547  if (sol[i] != NULL)
2548  {
2549  x = sol[i];
2550  sol[i] = n_Div(x, d->m,_R->cf);
2551  n_Delete(&x,_R->cf);
2552  }
2553  i--;
2554  }
2555  this->smAllDel();
2556 }
2557 
2558 /*
2559 * transform the result to an ideal
2560 */
2562 {
2563  int i, j;
2564  ideal res = idInit(crd, 1);
2565 
2566  for (i=crd; i; i--)
2567  {
2568  j = perm[i]-1;
2569  res->m[j] = sm_Smnumber2Poly(sol[i],_R);
2570  }
2571  omFreeSize((ADDRESS)sol, sizeof(number)*(crd+1));
2572  return res;
2573 }
2574 
2575 /* ----------------- pivot method ------------------ */
2576 
2577 /*
2578 * compute pivot
2579 */
2581 {
2582  smnumber a;
2583  number x, xo;
2584  int i, copt, ropt;
2585 
2586  xo=n_Init(0,_R->cf);
2587  for (i=act; i; i--)
2588  {
2589  a = m_act[i];
2590  while ((a!=NULL) && (a->pos<=tored))
2591  {
2592  x = a->m;
2593  if (n_GreaterZero(x,_R->cf))
2594  {
2595  if (n_Greater(x,xo,_R->cf))
2596  {
2597  n_Delete(&xo,_R->cf);
2598  xo = n_Copy(x,_R->cf);
2599  copt = i;
2600  ropt = a->pos;
2601  }
2602  }
2603  else
2604  {
2605  xo = n_InpNeg(xo,_R->cf);
2606  if (n_Greater(xo,x,_R->cf))
2607  {
2608  n_Delete(&xo,_R->cf);
2609  xo = n_Copy(x,_R->cf);
2610  copt = i;
2611  ropt = a->pos;
2612  }
2613  xo = n_InpNeg(xo,_R->cf);
2614  }
2615  a = a->n;
2616  }
2617  }
2618  rpiv = ropt;
2619  if (copt != act)
2620  {
2621  a = m_act[act];
2622  m_act[act] = m_act[copt];
2623  m_act[copt] = a;
2624  }
2625  n_Delete(&xo,_R->cf);
2626 }
2627 
2628 /* ----------------- elimination ------------------ */
2629 
2630 /* one step of Gauss-elimination */
2632 {
2633  number p = n_Invers(piv->m,_R->cf); // pivotelement
2634  smnumber c = m_act[act]; // pivotcolumn
2635  smnumber r = red; // row to reduce
2636  smnumber res, a, b;
2637  number w, ha, hb;
2638 
2639  if ((c == NULL) || (r == NULL))
2640  {
2641  while (r!=NULL) sm_NumberDelete(&r,_R);
2642  return;
2643  }
2644  do
2645  {
2646  a = m_act[r->pos];
2647  res = dumm;
2648  res->n = NULL;
2649  b = c;
2650  w = n_Mult(r->m, p,_R->cf);
2651  n_Delete(&r->m,_R->cf);
2652  r->m = w;
2653  loop // combine the chains a and b: a + w*b
2654  {
2655  if (a == NULL)
2656  {
2657  do
2658  {
2659  res = res->n = smNumberCopy(b);
2660  res->m = n_Mult(b->m, w,_R->cf);
2661  b = b->n;
2662  } while (b != NULL);
2663  break;
2664  }
2665  if (a->pos < b->pos)
2666  {
2667  res = res->n = a;
2668  a = a->n;
2669  }
2670  else if (a->pos > b->pos)
2671  {
2672  res = res->n = smNumberCopy(b);
2673  res->m = n_Mult(b->m, w,_R->cf);
2674  b = b->n;
2675  }
2676  else
2677  {
2678  hb = n_Mult(b->m, w,_R->cf);
2679  ha = n_Add(a->m, hb,_R->cf);
2680  n_Delete(&hb,_R->cf);
2681  n_Delete(&a->m,_R->cf);
2682  if (n_IsZero(ha,_R->cf))
2683  {
2684  sm_NumberDelete(&a,_R);
2685  }
2686  else
2687  {
2688  a->m = ha;
2689  res = res->n = a;
2690  a = a->n;
2691  }
2692  b = b->n;
2693  }
2694  if (b == NULL)
2695  {
2696  res->n = a;
2697  break;
2698  }
2699  }
2700  m_act[r->pos] = dumm->n;
2701  sm_NumberDelete(&r,_R);
2702  } while (r != NULL);
2703  n_Delete(&p,_R->cf);
2704 }
2705 
2706 /* ----------------- transfer ------------------ */
2707 
2708 /*
2709 * select the pivotrow and store it to red and piv
2710 */
2712 {
2713  smnumber b = dumm;
2714  smnumber a, ap;
2715  int i;
2716 
2717  if (TEST_OPT_PROT)
2718  {
2719  if ((crd+1)%10)
2720  PrintS(".");
2721  else
2722  PrintS(".\n");
2723  }
2724  a = m_act[act];
2725  if (a->pos < rpiv)
2726  {
2727  do
2728  {
2729  ap = a;
2730  a = a->n;
2731  } while (a->pos < rpiv);
2732  ap->n = a->n;
2733  }
2734  else
2735  m_act[act] = a->n;
2736  piv = a;
2737  a->n = NULL;
2738  for (i=1; i<act; i++)
2739  {
2740  a = m_act[i];
2741  if (a->pos < rpiv)
2742  {
2743  loop
2744  {
2745  ap = a;
2746  a = a->n;
2747  if ((a == NULL) || (a->pos > rpiv))
2748  break;
2749  if (a->pos == rpiv)
2750  {
2751  ap->n = a->n;
2752  a->m = n_InpNeg(a->m,_R->cf);
2753  b = b->n = a;
2754  b->pos = i;
2755  break;
2756  }
2757  }
2758  }
2759  else if (a->pos == rpiv)
2760  {
2761  m_act[i] = a->n;
2762  a->m = n_InpNeg(a->m,_R->cf);
2763  b = b->n = a;
2764  b->pos = i;
2765  }
2766  }
2767  b->n = NULL;
2768  red = dumm->n;
2769 }
2770 
2771 /*
2772 * store the pivotcol in m_row
2773 * m_act[cols][pos(rows)] => m_row[rows][pos(cols)]
2774 */
2776 {
2777  smnumber c = m_act[act];
2778  smnumber h;
2779 
2780  while (c != NULL)
2781  {
2782  h = c;
2783  c = c->n;
2784  h->n = m_row[h->pos];
2785  m_row[h->pos] = h;
2786  h->pos = crd;
2787  }
2788 }
2789 
2790 /*
2791 * store the pivot and the assosiated row in m_row
2792 * to m_res (result):
2793 * piv + m_row[rows][pos(cols)] => m_res[cols][pos(rows)]
2794 */
2796 {
2797  smnumber r = m_row[rpiv];
2798  smnumber a, ap, h;
2799 
2800  m_row[rpiv] = NULL;
2801  perm[crd] = rpiv;
2802  piv->pos = crd;
2803  m_res[crd] = piv;
2804  while (r != NULL)
2805  {
2806  ap = m_res[r->pos];
2807  loop
2808  {
2809  a = ap->n;
2810  if (a == NULL)
2811  {
2812  ap->n = h = r;
2813  r = r->n;
2814  h->n = a;
2815  h->pos = crd;
2816  break;
2817  }
2818  ap = a;
2819  }
2820  }
2821 }
2822 
2823 /* ----------------- C++ stuff ------------------ */
2824 
2825 /*
2826 * check singularity
2827 */
2829 {
2830  smnumber a;
2831  int i = act;
2832 
2833  loop
2834  {
2835  if (i == 0) return;
2836  a = m_act[i];
2837  if ((a==NULL) || (a->pos>tored))
2838  {
2839  sing = 1;
2840  this->smAllDel();
2841  return;
2842  }
2843  i--;
2844  }
2845 }
2846 
2847 /*
2848 * delete all smnumber's
2849 */
2851 {
2852  smnumber a;
2853  int i;
2854 
2855  for (i=act; i; i--)
2856  {
2857  a = m_act[i];
2858  while (a != NULL)
2859  sm_NumberDelete(&a,_R);
2860  }
2861  for (i=crd; i; i--)
2862  {
2863  a = m_res[i];
2864  while (a != NULL)
2865  sm_NumberDelete(&a,_R);
2866  }
2867  if (act)
2868  {
2869  for (i=nrows; i; i--)
2870  {
2871  a = m_row[i];
2872  while (a != NULL)
2873  sm_NumberDelete(&a,_R);
2874  }
2875  }
2876 }
2877 
2878 /* ----------------- internal 'C' stuff ------------------ */
2879 
2880 static void sm_NumberDelete(smnumber *r, const ring R)
2881 {
2882  smnumber a = *r, b = a->n;
2883 
2884  n_Delete(&a->m,R->cf);
2885  omFreeBin((ADDRESS)a, smnrec_bin);
2886  *r = b;
2887 }
2888 
2889 static smnumber smNumberCopy(smnumber a)
2890 {
2891  smnumber r = (smnumber)omAllocBin(smnrec_bin);
2892  memcpy(r, a, sizeof(smnrec));
2893  return r;
2894 }
2895 
2896 /*
2897 * from poly to smnumber
2898 * do not destroy p
2899 */
2900 static smnumber sm_Poly2Smnumber(poly q, const ring R)
2901 {
2902  smnumber a, res;
2903  poly p = q;
2904 
2905  if (p == NULL)
2906  return NULL;
2907  a = res = (smnumber)omAllocBin(smnrec_bin);
2908  a->pos = p_GetComp(p,R);
2909  a->m = pGetCoeff(p);
2910  n_New(&pGetCoeff(p),R->cf);
2911  loop
2912  {
2913  pIter(p);
2914  if (p == NULL)
2915  {
2916  p_Delete(&q,R);
2917  a->n = NULL;
2918  return res;
2919  }
2920  a = a->n = (smnumber)omAllocBin(smnrec_bin);
2921  a->pos = p_GetComp(p,R);
2922  a->m = pGetCoeff(p);
2923  n_New(&pGetCoeff(p),R->cf);
2924  }
2925 }
2926 
2927 /*
2928 * from smnumber to poly
2929 * destroy a
2930 */
2931 static poly sm_Smnumber2Poly(number a, const ring R)
2932 {
2933  poly res;
2934 
2935  if (a == NULL) return NULL;
2936  res = p_Init(R);
2937  pSetCoeff0(res, a);
2938  return res;
2939 }
2940 
2941 /*2
2942 * check the input
2943 */
2944 static BOOLEAN smCheckSolv(ideal I)
2945 { int i = I->ncols;
2946  if ((i == 0) || (i != I->rank-1))
2947  {
2948  WerrorS("wrong dimensions for linsolv");
2949  return TRUE;
2950  }
2951  for(;i;i--)
2952  {
2953  if(I->m[i-1] == NULL)
2954  {
2955  WerrorS("singular input for linsolv");
2956  return TRUE;
2957  }
2958  }
2959  return FALSE;
2960 }
#define n_New(n, r)
Definition: coeffs.h:444
static FORCE_INLINE number n_Sub(number a, number b, const coeffs r)
return the difference of &#39;a&#39; and &#39;b&#39;, i.e., a-b
Definition: coeffs.h:673
static poly sm_Smpoly2Poly(smpoly, const ring)
Definition: sparsmat.cc:2189
float f
Definition: sparsmat.cc:58
#define omAllocBin(bin)
Definition: omAllocDecl.h:205
static FORCE_INLINE BOOLEAN n_Greater(number a, number b, const coeffs r)
ordered fields: TRUE iff &#39;a&#39; is larger than &#39;b&#39;; in Z/pZ: TRUE iff la > lb, where la and lb are the l...
Definition: coeffs.h:515
int smGetRed()
Definition: sparsmat.cc:178
sparse_number_mat(ideal, const ring)
Definition: sparsmat.cc:2406
const CanonicalForm int s
Definition: facAbsFact.cc:55
ring sm_RingChange(const ring origR, long bound)
Definition: sparsmat.cc:262
smnumber * m_res
Definition: sparsmat.cc:2342
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:57
void sm_SpecialPolyDiv(poly a, poly b, const ring R)
Definition: sparsmat.cc:1894
static void sm_ExpMultDiv(poly, const poly, const poly, const ring)
Definition: sparsmat.cc:1986
static BOOLEAN p_LmIsConstantComp(const poly p, const ring r)
Definition: p_polys.h:932
static poly normalize(poly next_p, ideal add_generators, syStrategy syzstr, int *g_l, int *p_l, int crit_comp)
Definition: syz3.cc:1024
void kBucketInit(kBucket_pt bucket, poly lm, int length)
Definition: kbuckets.cc:467
const poly a
Definition: syzextra.cc:212
omBin_t * omBin
Definition: omStructs.h:12
static CanonicalForm bound(const CFMatrix &M)
Definition: cf_linsys.cc:460
#define Print
Definition: emacs.cc:83
void sm_PolyDiv(poly a, poly b, const ring R)
Definition: sparsmat.cc:1615
void smWeights()
Definition: sparsmat.cc:666
smpoly * m_res
Definition: sparsmat.cc:143
smpoly dumm
Definition: sparsmat.cc:147
static poly pp_Mult_Coeff_mm_DivSelect_MultDiv(poly p, int &lp, poly m, poly a, poly b, const ring currRing)
Definition: sparsmat.cc:81
#define TEST_OPT_PROT
Definition: options.h:98
loop
Definition: myNF.cc:98
#define FALSE
Definition: auxiliary.h:95
return P p
Definition: myNF.cc:203
static BOOLEAN p_LmDivisibleByNoComp(poly a, poly b, const ring r)
Definition: p_polys.h:1754
static FORCE_INLINE BOOLEAN n_IsOne(number n, const coeffs r)
TRUE iff &#39;n&#39; represents the one element.
Definition: coeffs.h:472
void smRowToCol()
Definition: sparsmat.cc:1157
void smNormalize()
Definition: sparsmat.cc:1489
omBin sip_sideal_bin
Definition: simpleideals.cc:30
static unsigned long p_SetComp(poly p, unsigned long c, ring r)
Definition: p_polys.h:242
void smPivot()
Definition: sparsmat.cc:698
poly sm_CallDet(ideal I, const ring R)
Definition: sparsmat.cc:358
int smGetSign()
Definition: sparsmat.cc:176
#define p_GetComp(p, r)
Definition: monomials.h:72
poly prMoveR(poly &p, ring src_r, ring dest_r)
Definition: prCopy.cc:91
ideal smRes2Ideal()
Definition: sparsmat.cc:2561
static FORCE_INLINE number n_Init(long i, const coeffs r)
a number representing i in the given coeff field/ring r
Definition: coeffs.h:542
int rows() const
Definition: intvec.h:88
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
float wpoints
Definition: sparsmat.cc:140
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:580
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
int smCheckNormalize()
Definition: sparsmat.cc:1469
omBin smprec_bin
Definition: sparsmat.cc:79
static BOOLEAN smSmaller(poly, poly)
Definition: sparsmat.cc:2017
static poly pp_Mult_mm(poly p, poly m, const ring r)
Definition: p_polys.h:957
#define TRUE
Definition: auxiliary.h:99
BOOLEAN id_IsConstant(ideal id, const ring r)
test if the ideal has only constant polynomials NOTE: zero ideal/module is also constant ...
void smActDel()
Definition: sparsmat.cc:1531
number m
Definition: sparsmat.cc:2316
static FORCE_INLINE void n_Normalize(number &n, const coeffs r)
inplace-normalization of n; produces some canonical representation of n;
Definition: coeffs.h:582
void * ADDRESS
Definition: auxiliary.h:116
sm_prec * smpoly
Definition: sparsmat.cc:51
void WerrorS(const char *s)
Definition: feFopen.cc:24
int k
Definition: cfEzgcd.cc:93
static number sm_Cleardenom(ideal, const ring)
Definition: sparsmat.cc:2278
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 omAlloc(size)
Definition: omAllocDecl.h:210
static void sm_NumberDelete(smnumber *, const ring R)
Definition: sparsmat.cc:2880
long sm_ExpBound(ideal m, int di, int ra, int t, const ring currRing)
Definition: sparsmat.cc:193
smnumber * m_act
Definition: sparsmat.cc:2341
static number p_SetCoeff(poly p, number n, ring r)
Definition: p_polys.h:407
static void p_LmFree(poly p, ring)
Definition: p_polys.h:678
int normalize
Definition: sparsmat.cc:138
poly kBucketExtractLm(kBucket_pt bucket)
Definition: kbuckets.cc:485
void smZeroToredElim()
Definition: sparsmat.cc:2828
poly pp
Definition: myNF.cc:296
static long p_SubExp(poly p, int v, long ee, ring r)
Definition: p_polys.h:608
#define TEST_OPT_NOT_BUCKETS
Definition: options.h:100
int e
Definition: sparsmat.cc:56
smpoly red
Definition: sparsmat.cc:145
#define pIter(p)
Definition: monomials.h:44
poly res
Definition: myNF.cc:322
static FORCE_INLINE number n_Mult(number a, number b, const coeffs r)
return the product of &#39;a&#39; and &#39;b&#39;, i.e., a*b
Definition: coeffs.h:640
#define M
Definition: sirandom.c:24
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:10
static poly sm_Smnumber2Poly(number, const ring)
Definition: sparsmat.cc:2931
void smSelectPR()
Definition: sparsmat.cc:1073
void smNewBareiss(int, int)
Definition: sparsmat.cc:605
sm_nrec * smnumber
Definition: sparsmat.cc:2312
const ring r
Definition: syzextra.cc:208
smpoly n
Definition: sparsmat.cc:54
void kBucketDestroy(kBucket_pt *bucket_pt)
Definition: kbuckets.cc:200
Definition: intvec.h:14
long id_RankFreeModule(ideal s, ring lmRing, ring tailRing)
return the maximal component number found in any polynomial in s
static smnumber smNumberCopy(smnumber)
Definition: sparsmat.cc:2889
BOOLEAN rComplete(ring r, int force)
this needs to be called whenever a new ring is created: new fields in ring are created (like VarOffse...
Definition: ring.cc:3351
void smNewPivot()
Definition: sparsmat.cc:793
static long p_GetExp(const poly p, const unsigned long iBitmask, const int VarOffset)
get a single variable exponent : the integer VarOffset encodes:
Definition: p_polys.h:464
int j
Definition: myNF.cc:70
static BOOLEAN p_IsConstant(const poly p, const ring r)
Definition: p_polys.h:1876
static FORCE_INLINE number n_Add(number a, number b, const coeffs r)
return the sum of &#39;a&#39; and &#39;b&#39;, i.e., a+b
Definition: coeffs.h:660
int pos
Definition: sparsmat.cc:55
static void sm_FindRef(poly *, poly *, poly, const ring)
Definition: sparsmat.cc:2075
static smpoly sm_Poly2Smpoly(poly, const ring)
Definition: sparsmat.cc:2152
ring rCopy0(const ring r, BOOLEAN copy_qideal, BOOLEAN copy_ordering)
Definition: ring.cc:1321
static BOOLEAN smCheckSolv(ideal)
Definition: sparsmat.cc:2944
static poly pp_Mult_qq(poly p, poly q, const ring r)
Definition: p_polys.h:1070
int nrows
Definition: cf_linsys.cc:32
const ring R
Definition: DebugPrint.cc:36
BOOLEAN sm_CheckDet(ideal I, int d, BOOLEAN sw, const ring r)
Definition: sparsmat.cc:307
static BOOLEAN sm_HaveDenom(poly, const ring)
Definition: sparsmat.cc:2259
void smColDel()
Definition: sparsmat.cc:1549
float * wrw
Definition: sparsmat.cc:141
P bucket
Definition: myNF.cc:79
ideal idrMoveR(ideal &id, ring src_r, ring dest_r)
Definition: prCopy.cc:249
smnumber n
Definition: sparsmat.cc:2314
All the auxiliary stuff.
static int p_LmCmp(poly p, poly q, const ring r)
Definition: p_polys.h:1467
static FORCE_INLINE number n_Invers(number a, const coeffs r)
return the multiplicative inverse of &#39;a&#39;; raise an error if &#39;a&#39; is not invertible ...
Definition: coeffs.h:568
#define SM_MIN_LENGTH_BUCKET
Definition: sparsmat.cc:45
static FORCE_INLINE number n_InpNeg(number n, const coeffs r)
in-place negation of n MUST BE USED: n = n_InpNeg(n) (no copy is returned)
Definition: coeffs.h:561
smpoly piv
Definition: sparsmat.cc:146
static int si_max(const int a, const int b)
Definition: auxiliary.h:121
int i
Definition: cfEzgcd.cc:123
void PrintS(const char *s)
Definition: reporter.cc:284
static poly p_Mult_nn(poly p, number n, const ring r)
Definition: p_polys.h:895
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:501
static void p_ExpVectorAdd(poly p1, poly p2, const ring r)
Definition: p_polys.h:1334
static unsigned pLength(poly a)
Definition: p_polys.h:189
#define IDELEMS(i)
Definition: simpleideals.h:24
static FORCE_INLINE BOOLEAN n_IsZero(number n, const coeffs r)
TRUE iff &#39;n&#39; represents the zero element.
Definition: coeffs.h:468
int * perm
Definition: sparsmat.cc:139
#define p_LmTest(p, r)
Definition: p_polys.h:161
static void sm_PolyDivN(poly, const number, const ring)
Definition: sparsmat.cc:2004
void smPivDel()
Definition: sparsmat.cc:1562
CanonicalForm gp
Definition: cfModGcd.cc:4043
#define p_Test(p, r)
Definition: p_polys.h:160
smnumber * m_row
Definition: sparsmat.cc:2343
smpoly * m_act
Definition: sparsmat.cc:142
smpoly * m_row
Definition: sparsmat.cc:144
void p_Normalize(poly p, const ring r)
Definition: p_polys.cc:3633
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:843
#define omGetSpecBin(size)
Definition: omBin.h:11
void smHElim()
Definition: sparsmat.cc:934
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:38
void sm1Elim()
Definition: sparsmat.cc:855
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
static void p_ExpVectorDiff(poly pr, poly p1, poly p2, const ring r)
Definition: p_polys.h:1397
static unsigned long p_SetExp(poly p, const unsigned long e, const unsigned long iBitmask, const int VarOffset)
set a single variable exponent : VarOffset encodes the position in p->exp
Definition: p_polys.h:483
#define NULL
Definition: omList.c:10
ideal smRes2Mod()
Definition: sparsmat.cc:504
static FORCE_INLINE number n_Copy(number n, const coeffs r)
return a copy of &#39;n&#39;
Definition: coeffs.h:455
void sm_KillModifiedRing(ring r)
Definition: sparsmat.cc:293
ideal sm_CallSolv(ideal I, const ring R)
Definition: sparsmat.cc:2371
static void sm_CombineChain(poly *, poly, const ring)
Definition: sparsmat.cc:2028
static FORCE_INLINE number n_Div(number a, number b, const coeffs r)
return the quotient of &#39;a&#39; and &#39;b&#39;, i.e., a/b; raises an error if &#39;b&#39; is not invertible in r exceptio...
Definition: coeffs.h:619
static smnumber sm_Poly2Smnumber(poly, const ring)
Definition: sparsmat.cc:2900
void smInitPerm()
Definition: sparsmat.cc:1604
int int ncols
Definition: cf_linsys.cc:32
poly smDet()
Definition: sparsmat.cc:531
void smColToRow()
Definition: sparsmat.cc:1137
const CanonicalForm & w
Definition: facAbsFact.cc:55
static smpoly smElemCopy(smpoly)
Definition: sparsmat.cc:2140
#define SM_MULT
Definition: sparsmat.h:23
sparse_mat(ideal, const ring)
Definition: sparsmat.cc:442
Variable x
Definition: cfModGcd.cc:4023
static BOOLEAN p_IsConstantPoly(const poly p, const ring r)
Definition: p_polys.h:1890
static FORCE_INLINE number n_GetDenom(number &n, const coeffs r)
return the denominator of n (if elements of r are by nature not fractional, result is 1) ...
Definition: coeffs.h:607
#define pNext(p)
Definition: monomials.h:43
static FORCE_INLINE BOOLEAN n_Equal(number a, number b, const coeffs r)
TRUE iff &#39;a&#39; and &#39;b&#39; represent the same number; they may have different representations.
Definition: coeffs.h:464
void rKillModifiedRing(ring r)
Definition: ring.cc:2962
ideal idrCopyR(ideal id, ring src_r, ring dest_r)
Definition: prCopy.cc:193
static void p_Setm(poly p, const ring r)
Definition: p_polys.h:228
void smSign()
Definition: sparsmat.cc:1576
#define pSetCoeff0(p, n)
Definition: monomials.h:67
poly sm_MultDiv(poly a, poly b, const poly c, const ring R)
Definition: sparsmat.cc:1814
static void smMinSelect(long *, int, int)
Definition: sparsmat.cc:240
#define SM_DIV
Definition: sparsmat.h:24
void smMultCol()
Definition: sparsmat.cc:1415
void sm_CallBareiss(ideal I, int x, int y, ideal &M, intvec **iv, const ring R)
Definition: sparsmat.cc:403
static void sm_ExactPolyDiv(poly, poly, const ring)
Definition: sparsmat.cc:1906
void smCopToRes()
Definition: sparsmat.cc:1259
static poly p_Neg(poly p, const ring r)
Definition: p_polys.h:1013
static void p_LmDelete(poly p, const ring r)
Definition: p_polys.h:706
static omBin smnrec_bin
Definition: sparsmat.cc:2319
static float sm_PolyWeight(smpoly, const ring)
Definition: sparsmat.cc:2230
static scmon act
Definition: hdegree.cc:1078
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete &#39;p&#39;
Definition: coeffs.h:459
void smToredElim()
Definition: sparsmat.cc:1220
static FORCE_INLINE BOOLEAN n_GreaterZero(number n, const coeffs r)
ordered fields: TRUE iff &#39;n&#39; is positive; in Z/pZ: TRUE iff 0 < m <= roundedBelow(p/2), where m is the long representing n in C: TRUE iff (Im(n) != 0 and Im(n) >= 0) or (Im(n) == 0 and Re(n) >= 0) in K(a)/<p(a)>: TRUE iff (n != 0 and (LC(n) > 0 or deg(n) > 0)) in K(t_1, ..., t_n): TRUE iff (LC(numerator(n) is a constant and > 0) or (LC(numerator(n) is not a constant) in Z/2^kZ: TRUE iff 0 < n <= 2^(k-1) in Z/mZ: TRUE iff the internal mpz is greater than zero in Z: TRUE iff n > 0
Definition: coeffs.h:498
void smZeroElim()
Definition: sparsmat.cc:1190
polyrec * poly
Definition: hilb.h:10
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:877
static poly sm_SelectCopy_ExpMultDiv(poly p, poly m, poly a, poly b, const ring currRing)
Definition: sparsmat.cc:103
#define omFreeBin(addr, bin)
Definition: omAllocDecl.h:259
static BOOLEAN rOrd_is_Comp_dp(const ring r)
Definition: ring.h:765
void smFinalMult()
Definition: sparsmat.cc:1440
kBucket_pt kBucketCreate(const ring bucket_ring)
Creation/Destruction of buckets.
Definition: kbuckets.cc:193
static poly p_New(const ring, omBin bin)
Definition: p_polys.h:659
int perm[100]
poly m
Definition: sparsmat.cc:57
static Poly * h
Definition: janet.cc:978
int BOOLEAN
Definition: auxiliary.h:86
static poly p_Init(const ring r, omBin bin)
Definition: p_polys.h:1243
const poly b
Definition: syzextra.cc:213
poly smMultPoly(smpoly)
Definition: sparsmat.cc:1509
poly p_Cleardenom(poly p, const ring r)
Definition: p_polys.cc:2715
static void sm_ElemDelete(smpoly *, const ring)
Definition: sparsmat.cc:2131
void smNewWeights()
Definition: sparsmat.cc:755
ideal idrCopyR_NoSort(ideal id, ring src_r, ring dest_r)
Definition: prCopy.cc:206
static FORCE_INLINE int n_Size(number n, const coeffs r)
return a non-negative measure for the complexity of n; return 0 only when n represents zero; (used fo...
Definition: coeffs.h:574
void smToIntvec(intvec *)
Definition: sparsmat.cc:520
static int sign(int x)
Definition: ring.cc:3328
void Werror(const char *fmt,...)
Definition: reporter.cc:189
int pos
Definition: sparsmat.cc:2315
static poly pp_Mult_Coeff_mm_DivSelect(poly p, const poly m, const ring r)
Definition: p_polys.h:996
#define omAlloc0(size)
Definition: omAllocDecl.h:211
int l
Definition: cfEzgcd.cc:94
smpoly * smGetAct()
Definition: sparsmat.cc:177
static long p_MaxComp(poly p, ring lmRing, ring tailRing)
Definition: p_polys.h:287
void kBucket_Add_q(kBucket_pt bucket, poly q, int *l)
Add to Bucket a poly ,i.e. Bpoly == q+Bpoly.
Definition: kbuckets.cc:628
static BOOLEAN sm_IsNegQuot(poly, const poly, const poly, const ring)
Definition: sparsmat.cc:1961