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