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 #include "coeffrings.h"
24 
25 #include "monomials/ring.h"
26 #include "monomials/p_polys.h"
27 
28 #include "simpleideals.h"
29 
30 
31 #include "sparsmat.h"
32 #include "prCopy.h"
33 
34 #include "templates/p_Procs.h"
35 
36 #include "kbuckets.h"
37 #include "operations/p_Mult_q.h"
38 
39 // define SM_NO_BUCKETS, if sparsemat stuff should not use buckets
40 // #define SM_NO_BUCKETS
41 
42 // this is also influenced by TEST_OPT_NOTBUCKETS
43 #ifndef SM_NO_BUCKETS
44 // buckets do carry a small additional overhead: only use them if
45 // min-length of polys is >= SM_MIN_LENGTH_BUCKET
46 #define SM_MIN_LENGTH_BUCKET MIN_LENGTH_BUCKET - 5
47 #else
48 #define SM_MIN_LENGTH_BUCKET MAX_INT_VAL
49 #endif
50 
51 typedef struct smprec sm_prec;
52 typedef sm_prec * smpoly;
53 struct smprec
54 {
55  smpoly n; // the next element
56  int pos; // position
57  int e; // level
58  poly m; // the element
59  float f; // complexity of the element
60 };
61 
62 
63 /* declare internal 'C' stuff */
64 static void sm_ExactPolyDiv(poly, poly, const ring);
65 static BOOLEAN sm_IsNegQuot(poly, const poly, const poly, const ring);
66 static void sm_ExpMultDiv(poly, const poly, const poly, const ring);
67 static void sm_PolyDivN(poly, const number, const ring);
68 static BOOLEAN smSmaller(poly, poly);
69 static void sm_CombineChain(poly *, poly, const ring);
70 static void sm_FindRef(poly *, poly *, poly, const ring);
71 
72 static void sm_ElemDelete(smpoly *, const ring);
73 static smpoly smElemCopy(smpoly);
74 static float sm_PolyWeight(smpoly, const ring);
75 static smpoly sm_Poly2Smpoly(poly, const ring);
76 static poly sm_Smpoly2Poly(smpoly, const ring);
77 static BOOLEAN sm_HaveDenom(poly, const ring);
78 static number sm_Cleardenom(ideal, const ring);
79 
81 
83  poly a, poly b, const ring currRing)
84 {
85  if (rOrd_is_Comp_dp(currRing) && currRing->ExpL_Size > 2)
86  {
87  // pp_Mult_Coeff_mm_DivSelectMult only works for (c/C,dp) and
88  // ExpL_Size > 2
89  // should be generalized, at least to dp with ExpL_Size == 2
90  // (is the case for 1 variable)
91  int shorter;
92  p = currRing->p_Procs->pp_Mult_Coeff_mm_DivSelectMult(p, m, a, b,
93  shorter, currRing);
94  lp -= shorter;
95  }
96  else
97  {
98  p = pp_Mult_Coeff_mm_DivSelect(p, lp, m, currRing);
99  sm_ExpMultDiv(p, a, b, currRing);
100  }
101  return p;
102 }
103 
105 {
106  int lp = 0;
107  return pp_Mult_Coeff_mm_DivSelect_MultDiv(p, lp, m, a, b, currRing);
108 }
109 
110 
111 /* class for sparse matrix:
112 * 3 parts of matrix during the algorithm
113 * m_act[cols][pos(rows)] => m_row[rows][pos(cols)] => m_res[cols][pos(rows)]
114 * input pivotcols as rows result
115 * pivot like a stack from pivot and pivotcols
116 * elimination rows reordered according
117 * to pivot choise
118 * stored in perm
119 * a step is as follows
120 * - search of pivot (smPivot)
121 * - swap pivot column and last column (smSwapC)
122 * select pivotrow to piv and red (smSelectPR)
123 * consider sign
124 * - elimination (smInitElim, sm0Elim, sm1Elim)
125 * clear zero column as result of elimination (smZeroElim)
126 * - tranfer from
127 * piv and m_row to m_res (smRowToCol)
128 * m_act to m_row (smColToRow)
129 */
131 private:
132  int nrows, ncols; // dimension of the problem
133  int sign; // for determinant (start: 1)
134  int act; // number of unreduced columns (start: ncols)
135  int crd; // number of reduced columns (start: 0)
136  int tored; // border for rows to reduce
137  int inred; // unreducable part
138  int rpiv, cpiv; // position of the pivot
139  int normalize; // Normalization flag
140  int *perm; // permutation of rows
141  float wpoints; // weight of all points
142  float *wrw, *wcl; // weights of rows and columns
143  smpoly * m_act; // unreduced columns
144  smpoly * m_res; // reduced columns (result)
145  smpoly * m_row; // reduced part of rows
146  smpoly red; // row to reduce
147  smpoly piv, oldpiv; // pivot and previous pivot
148  smpoly dumm; // allocated dummy
149  ring _R;
150 
151  void smColToRow();
152  void smRowToCol();
153  void smFinalMult();
154  void smSparseHomog();
155  void smWeights();
156  void smPivot();
157  void smNewWeights();
158  void smNewPivot();
159  void smZeroElim();
160  void smToredElim();
161  void smCopToRes();
162  void smSelectPR();
163  void sm1Elim();
164  void smHElim();
165  void smMultCol();
167  void smActDel();
168  void smColDel();
169  void smPivDel();
170  void smSign();
171  void smInitPerm();
172  int smCheckNormalize();
173  void smNormalize();
174 public:
175  sparse_mat(ideal, const ring);
176  ~sparse_mat();
177  int smGetSign() { return sign; }
178  smpoly * smGetAct() { return m_act; }
179  int smGetRed() { return tored; }
180  ideal smRes2Mod();
181  poly smDet();
182  void smNewBareiss(int, int);
183  void smToIntvec(intvec *);
184 };
185 
186 /*
187 * estimate maximal exponent for det or minors,
188 * the module m has di vectors and maximal rank ra,
189 * estimate yields for the tXt minors,
190 * we have di,ra >= t
191 */
192 static void smMinSelect(long *, int, int);
193 
194 long sm_ExpBound( ideal m, int di, int ra, int t, const ring currRing)
195 {
196  poly p;
197  long kr, kc;
198  long *r, *c;
199  int al, bl, i, j, k;
200 
201  if (ra==0) ra=1;
202  al = di*sizeof(long);
203  c = (long *)omAlloc(al);
204  bl = ra*sizeof(long);
205  r = (long *)omAlloc0(bl);
206  for (i=di-1;i>=0;i--)
207  {
208  kc = 0;
209  p = m->m[i];
210  while(p!=NULL)
211  {
212  k = p_GetComp(p, currRing)-1;
213  kr = r[k];
214  for (j=rVar(currRing);j>0;j--)
215  {
216  if(p_GetExp(p,j, currRing)>kc)
217  kc=p_GetExp(p,j, currRing);
218  if(p_GetExp(p,j, currRing)>kr)
219  kr=p_GetExp(p,j, currRing);
220  }
221  r[k] = kr;
222  pIter(p);
223  }
224  c[i] = kc;
225  }
226  if (t<di) smMinSelect(c, t, di);
227  if (t<ra) smMinSelect(r, t, ra);
228  kr = kc = 0;
229  for (j=t-1;j>=0;j--)
230  {
231  kr += r[j];
232  kc += c[j];
233  }
234  omFreeSize((ADDRESS)c, al);
235  omFreeSize((ADDRESS)r, bl);
236  if (kr<kc) kc = kr;
237  if (kr<1) kr = 1;
238  return kr;
239 }
240 
241 static void smMinSelect(long *c, int t, int d)
242 {
243  long m;
244  int pos, i;
245  do
246  {
247  d--;
248  pos = d;
249  m = c[pos];
250  for (i=d-1;i>=0;i--)
251  {
252  if(c[i]<m)
253  {
254  pos = i;
255  m = c[i];
256  }
257  }
258  for (i=pos;i<d;i++) c[i] = c[i+1];
259  } while (d>t);
260 }
261 
262 /* ----------------- ops with rings ------------------ */
263 ring sm_RingChange(const ring origR, long bound)
264 {
265 // *origR =currRing;
266  ring tmpR=rCopy0(origR,FALSE,FALSE);
267  int *ord=(int*)omAlloc0(3*sizeof(int));
268  int *block0=(int*)omAlloc(3*sizeof(int));
269  int *block1=(int*)omAlloc(3*sizeof(int));
270  ord[0]=ringorder_c;
271  ord[1]=ringorder_dp;
272  tmpR->order=ord;
273  tmpR->OrdSgn=1;
274  block0[1]=1;
275  tmpR->block0=block0;
276  block1[1]=tmpR->N;
277  tmpR->block1=block1;
278  tmpR->bitmask = 2*bound;
279  tmpR->wvhdl = (int **)omAlloc0((3) * sizeof(int*));
280 
281 // ???
282 // if (tmpR->bitmask > currRing->bitmask) tmpR->bitmask = currRing->bitmask;
283 
284  rComplete(tmpR,1);
285  if (origR->qideal!=NULL)
286  {
287  tmpR->qideal= idrCopyR_NoSort(origR->qideal, origR, tmpR);
288  }
289  if (TEST_OPT_PROT)
290  Print("[%ld:%d]", (long) tmpR->bitmask, tmpR->ExpL_Size);
291  return tmpR;
292 }
293 
295 {
296  if (r->qideal!=NULL) id_Delete(&(r->qideal),r);
298 }
299 
300 /*2
301 * Bareiss or Chinese remainder ?
302 * I is dXd
303 * sw = TRUE -> I Matrix
304 * FALSE -> I Module
305 * return True -> change Type
306 * FALSE -> same Type
307 */
308 BOOLEAN sm_CheckDet(ideal I, int d, BOOLEAN sw, const ring r)
309 {
310  int s,t,i;
311  poly p;
312 
313  if (d>100)
314  return sw;
315  if (!rField_is_Q(r))
316  return sw;
317  s = t = 0;
318  if (sw)
319  {
320  for(i=IDELEMS(I)-1;i>=0;i--)
321  {
322  p=I->m[i];
323  if (p!=NULL)
324  {
325  if(!p_IsConstant(p,r))
326  return sw;
327  s++;
328  t+=n_Size(pGetCoeff(p),r->cf);
329  }
330  }
331  }
332  else
333  {
334  for(i=IDELEMS(I)-1;i>=0;i--)
335  {
336  p=I->m[i];
337  if (!p_IsConstantPoly(p,r))
338  return sw;
339  while (p!=NULL)
340  {
341  s++;
342  t+=n_Size(pGetCoeff(p),r->cf);
343  pIter(p);
344  }
345  }
346  }
347  s*=15;
348  if (t>s)
349  return !sw;
350  else
351  return sw;
352 }
353 
354 /* ----------------- basics (used from 'C') ------------------ */
355 /*2
356 *returns the determinant of the module I;
357 *uses Bareiss algorithm
358 */
359 poly sm_CallDet(ideal I,const ring R)
360 {
361  if (I->ncols != I->rank)
362  {
363  Werror("det of %ld x %d module (matrix)",I->rank,I->ncols);
364  return NULL;
365  }
366  int r=id_RankFreeModule(I,R);
367  if (I->ncols != r) // some 0-lines at the end
368  {
369  return NULL;
370  }
371  long bound=sm_ExpBound(I,r,r,r,R);
372  number diag,h=n_Init(1,R->cf);
373  poly res;
374  ring tmpR;
375  sparse_mat *det;
376  ideal II;
377 
378  tmpR=sm_RingChange(R,bound);
379  II = idrCopyR(I, R, tmpR);
380  diag = sm_Cleardenom(II,tmpR);
381  det = new sparse_mat(II,tmpR);
382  id_Delete(&II,tmpR);
383  if (det->smGetAct() == NULL)
384  {
385  delete det;
386  sm_KillModifiedRing(tmpR);
387  return NULL;
388  }
389  res=det->smDet();
390  if(det->smGetSign()<0) res=p_Neg(res,tmpR);
391  delete det;
392  res = prMoveR(res, tmpR, R);
393  sm_KillModifiedRing(tmpR);
394  if (!n_Equal(diag,h,R->cf))
395  {
396  p_Mult_nn(res,diag,R);
397  p_Normalize(res,R);
398  }
399  n_Delete(&diag,R->cf);
400  n_Delete(&h,R->cf);
401  return res;
402 }
403 
404 void sm_CallBareiss(ideal I, int x, int y, ideal & M, intvec **iv, const ring R)
405 {
406  int r=id_RankFreeModule(I,R),t=r;
407  int c=IDELEMS(I),s=c;
408  long bound;
409  ring tmpR;
410  sparse_mat *bareiss;
411 
412  if ((x>0) && (x<t))
413  t-=x;
414  if ((y>1) && (y<s))
415  s-=y;
416  if (t>s) t=s;
417  bound=sm_ExpBound(I,c,r,t,R);
418  tmpR=sm_RingChange(R,bound);
419  ideal II = idrCopyR(I, R, tmpR);
420  bareiss = new sparse_mat(II,tmpR);
421  if (bareiss->smGetAct() == NULL)
422  {
423  delete bareiss;
424  *iv=new intvec(1,rVar(tmpR));
425  }
426  else
427  {
428  id_Delete(&II,tmpR);
429  bareiss->smNewBareiss(x, y);
430  II = bareiss->smRes2Mod();
431  *iv = new intvec(bareiss->smGetRed());
432  bareiss->smToIntvec(*iv);
433  delete bareiss;
434  II = idrMoveR(II,tmpR,R);
435  }
436  sm_KillModifiedRing(tmpR);
437  M=II;
438 }
439 
440 /*
441 * constructor
442 */
443 sparse_mat::sparse_mat(ideal smat, const ring RR)
444 {
445  int i;
446  poly* pmat;
447  _R=RR;
448 
449  ncols = smat->ncols;
450  nrows = id_RankFreeModule(smat,RR);
451  if (nrows <= 0)
452  {
453  m_act = NULL;
454  return;
455  }
456  sign = 1;
457  inred = act = ncols;
458  crd = 0;
459  tored = nrows; // without border
460  i = tored+1;
461  perm = (int *)omAlloc(sizeof(int)*(i+1));
462  perm[i] = 0;
463  m_row = (smpoly *)omAlloc0(sizeof(smpoly)*i);
464  wrw = (float *)omAlloc(sizeof(float)*i);
465  i = ncols+1;
466  wcl = (float *)omAlloc(sizeof(float)*i);
467  m_act = (smpoly *)omAlloc(sizeof(smpoly)*i);
468  m_res = (smpoly *)omAlloc0(sizeof(smpoly)*i);
471  m_res[0]->m = NULL;
472  pmat = smat->m;
473  for(i=ncols; i; i--)
474  {
475  m_act[i] = sm_Poly2Smpoly(pmat[i-1], RR);
476  pmat[i-1] = NULL;
477  }
478  this->smZeroElim();
479  oldpiv = NULL;
480 }
481 
482 /*
483 * destructor
484 */
486 {
487  int i;
488  if (m_act == NULL) return;
491  i = ncols+1;
492  omFreeSize((ADDRESS)m_res, sizeof(smpoly)*i);
493  omFreeSize((ADDRESS)m_act, sizeof(smpoly)*i);
494  omFreeSize((ADDRESS)wcl, sizeof(float)*i);
495  i = nrows+1;
496  omFreeSize((ADDRESS)wrw, sizeof(float)*i);
497  omFreeSize((ADDRESS)m_row, sizeof(smpoly)*i);
498  omFreeSize((ADDRESS)perm, sizeof(int)*(i+1));
499 }
500 
501 /*
502 * transform the result to a module
503 */
504 
506 {
507  ideal res = idInit(crd, crd);
508  int i;
509 
510  for (i=crd; i; i--) res->m[i-1] = sm_Smpoly2Poly(m_res[i],_R);
511  res->rank = id_RankFreeModule(res,_R);
512  return res;
513 }
514 
515 /*
516 * permutation of rows
517 */
519 {
520  int i;
521 
522  for (i=v->rows()-1; i>=0; i--)
523  (*v)[i] = perm[i+1];
524 }
525 /* ---------------- the algorithm's ------------------ */
526 /*
527 * the determinant (up to sign), uses new Bareiss elimination
528 */
530 {
531  poly res = NULL;
532 
533  if (sign == 0)
534  {
535  this->smActDel();
536  return NULL;
537  }
538  if (act < 2)
539  {
540  if (act != 0) res = m_act[1]->m;
541  omFreeBin((void *)m_act[1], smprec_bin);
542  return res;
543  }
544  normalize = 0;
545  this->smInitPerm();
546  this->smPivot();
547  this->smSign();
548  this->smSelectPR();
549  this->sm1Elim();
550  crd++;
551  m_res[crd] = piv;
552  this->smColDel();
553  act--;
554  this->smZeroElim();
555  if (sign == 0)
556  {
557  this->smActDel();
558  return NULL;
559  }
560  if (act < 2)
561  {
562  this->smFinalMult();
563  this->smPivDel();
564  if (act != 0) res = m_act[1]->m;
565  omFreeBin((void *)m_act[1], smprec_bin);
566  return res;
567  }
568  loop
569  {
570  this->smNewPivot();
571  this->smSign();
572  this->smSelectPR();
573  this->smMultCol();
574  this->smHElim();
575  crd++;
576  m_res[crd] = piv;
577  this->smColDel();
578  act--;
579  this->smZeroElim();
580  if (sign == 0)
581  {
582  this->smPivDel();
583  this->smActDel();
584  return NULL;
585  }
586  if (act < 2)
587  {
588  if (TEST_OPT_PROT) PrintS(".\n");
589  this->smFinalMult();
590  this->smPivDel();
591  if (act != 0) res = m_act[1]->m;
592  omFreeBin((void *)m_act[1], smprec_bin);
593  return res;
594  }
595  }
596 }
597 
598 /*
599 * the new Bareiss elimination:
600 * - with x unreduced last rows, pivots from here are not allowed
601 * - the method will finish for number of unreduced columns < y
602 */
604 {
605  if ((x > 0) && (x < nrows))
606  {
607  tored -= x;
608  this->smToredElim();
609  }
610  if (y < 1) y = 1;
611  if (act <= y)
612  {
613  this->smCopToRes();
614  return;
615  }
616  normalize = this->smCheckNormalize();
617  if (normalize) this->smNormalize();
618  this->smPivot();
619  this->smSelectPR();
620  this->sm1Elim();
621  crd++;
622  this->smColToRow();
623  act--;
624  this->smRowToCol();
625  this->smZeroElim();
626  if (tored != nrows)
627  this->smToredElim();
628  if (act <= y)
629  {
630  this->smFinalMult();
631  this->smCopToRes();
632  return;
633  }
634  loop
635  {
636  if (normalize) this->smNormalize();
637  this->smNewPivot();
638  this->smSelectPR();
639  this->smMultCol();
640  this->smHElim();
641  crd++;
642  this->smColToRow();
643  act--;
644  this->smRowToCol();
645  this->smZeroElim();
646  if (tored != nrows)
647  this->smToredElim();
648  if (act <= y)
649  {
650  if (TEST_OPT_PROT) PrintS(".\n");
651  this->smFinalMult();
652  this->smCopToRes();
653  return;
654  }
655  }
656 }
657 
658 /* ----------------- pivot method ------------------ */
659 
660 /*
661 * prepare smPivot, compute weights for rows and columns
662 * and the weight for all points
663 */
665 {
666  float wc, wp, w;
667  smpoly a;
668  int i;
669 
670  wp = 0.0;
671  for (i=tored; i; i--) wrw[i] = 0.0; // ???
672  for (i=act; i; i--)
673  {
674  wc = 0.0;
675  a = m_act[i];
676  loop
677  {
678  if (a->pos > tored)
679  break;
680  w = a->f = sm_PolyWeight(a,_R);
681  wc += w;
682  wrw[a->pos] += w;
683  a = a->n;
684  if (a == NULL)
685  break;
686  }
687  wp += wc;
688  wcl[i] = wc;
689  }
690  wpoints = wp;
691 }
692 
693 /*
694 * compute pivot
695 */
697 {
698  float wopt = 1.0e30;
699  float wc, wr, wp, w;
700  smpoly a;
701  int i, copt, ropt;
702 
703  this->smWeights();
704  for (i=act; i; i--)
705  {
706  a = m_act[i];
707  loop
708  {
709  if (a->pos > tored)
710  break;
711  w = a->f;
712  wc = wcl[i]-w;
713  wr = wrw[a->pos]-w;
714  if ((wr<0.25) || (wc<0.25)) // row or column with only one point
715  {
716  if (w<wopt)
717  {
718  wopt = w;
719  copt = i;
720  ropt = a->pos;
721  }
722  }
723  else // elimination
724  {
725  wp = w*(wpoints-wcl[i]-wr);
726  wp += wr*wc;
727  if (wp < wopt)
728  {
729  wopt = wp;
730  copt = i;
731  ropt = a->pos;
732  }
733  }
734  a = a->n;
735  if (a == NULL)
736  break;
737  }
738  }
739  rpiv = ropt;
740  cpiv = copt;
741  if (cpiv != act)
742  {
743  a = m_act[act];
744  m_act[act] = m_act[cpiv];
745  m_act[cpiv] = a;
746  }
747 }
748 
749 /*
750 * prepare smPivot, compute weights for rows and columns
751 * and the weight for all points
752 */
754 {
755  float wc, wp, w, hp = piv->f;
756  smpoly a;
757  int i, f, e = crd;
758 
759  wp = 0.0;
760  for (i=tored; i; i--) wrw[i] = 0.0; // ???
761  for (i=act; i; i--)
762  {
763  wc = 0.0;
764  a = m_act[i];
765  loop
766  {
767  if (a->pos > tored)
768  break;
769  w = a->f;
770  f = a->e;
771  if (f < e)
772  {
773  w *= hp;
774  if (f) w /= m_res[f]->f;
775  }
776  wc += w;
777  wrw[a->pos] += w;
778  a = a->n;
779  if (a == NULL)
780  break;
781  }
782  wp += wc;
783  wcl[i] = wc;
784  }
785  wpoints = wp;
786 }
787 
788 /*
789 * compute pivot
790 */
792 {
793  float wopt = 1.0e30, hp = piv->f;
794  float wc, wr, wp, w;
795  smpoly a;
796  int i, copt, ropt, f, e = crd;
797 
798  this->smNewWeights();
799  for (i=act; i; i--)
800  {
801  a = m_act[i];
802  loop
803  {
804  if (a->pos > tored)
805  break;
806  w = a->f;
807  f = a->e;
808  if (f < e)
809  {
810  w *= hp;
811  if (f) w /= m_res[f]->f;
812  }
813  wc = wcl[i]-w;
814  wr = wrw[a->pos]-w;
815  if ((wr<0.25) || (wc<0.25)) // row or column with only one point
816  {
817  if (w<wopt)
818  {
819  wopt = w;
820  copt = i;
821  ropt = a->pos;
822  }
823  }
824  else // elimination
825  {
826  wp = w*(wpoints-wcl[i]-wr);
827  wp += wr*wc;
828  if (wp < wopt)
829  {
830  wopt = wp;
831  copt = i;
832  ropt = a->pos;
833  }
834  }
835  a = a->n;
836  if (a == NULL)
837  break;
838  }
839  }
840  rpiv = ropt;
841  cpiv = copt;
842  if (cpiv != act)
843  {
844  a = m_act[act];
845  m_act[act] = m_act[cpiv];
846  m_act[cpiv] = a;
847  }
848 }
849 
850 /* ----------------- elimination ------------------ */
851 
852 /* first step of elimination */
854 {
855  poly p = piv->m; // pivotelement
856  smpoly c = m_act[act]; // pivotcolumn
857  smpoly r = red; // row to reduce
858  smpoly res, a, b;
859  poly w, ha, hb;
860 
861  if ((c == NULL) || (r == NULL))
862  {
863  while (r!=NULL) sm_ElemDelete(&r,_R);
864  return;
865  }
866  do
867  {
868  a = m_act[r->pos];
869  res = dumm;
870  res->n = NULL;
871  b = c;
872  w = r->m;
873  loop // combine the chains a and b: p*a + w*b
874  {
875  if (a == NULL)
876  {
877  do
878  {
879  res = res->n = smElemCopy(b);
880  res->m = pp_Mult_qq(b->m, w,_R);
881  res->e = 1;
882  res->f = sm_PolyWeight(res,_R);
883  b = b->n;
884  } while (b != NULL);
885  break;
886  }
887  if (a->pos < b->pos)
888  {
889  res = res->n = a;
890  a = a->n;
891  }
892  else if (a->pos > b->pos)
893  {
894  res = res->n = smElemCopy(b);
895  res->m = pp_Mult_qq(b->m, w,_R);
896  res->e = 1;
897  res->f = sm_PolyWeight(res,_R);
898  b = b->n;
899  }
900  else
901  {
902  ha = pp_Mult_qq(a->m, p,_R);
903  p_Delete(&a->m,_R);
904  hb = pp_Mult_qq(b->m, w,_R);
905  ha = p_Add_q(ha, hb,_R);
906  if (ha != NULL)
907  {
908  a->m = ha;
909  a->e = 1;
910  a->f = sm_PolyWeight(a,_R);
911  res = res->n = a;
912  a = a->n;
913  }
914  else
915  {
916  sm_ElemDelete(&a,_R);
917  }
918  b = b->n;
919  }
920  if (b == NULL)
921  {
922  res->n = a;
923  break;
924  }
925  }
926  m_act[r->pos] = dumm->n;
927  sm_ElemDelete(&r,_R);
928  } while (r != NULL);
929 }
930 
931 /* higher steps of elimination */
933 {
934  poly hp = this->smMultPoly(piv);
935  poly gp = piv->m; // pivotelement
936  smpoly c = m_act[act]; // pivotcolumn
937  smpoly r = red; // row to reduce
938  smpoly res, a, b;
939  poly ha, hr, x, y;
940  int e, ip, ir, ia;
941 
942  if ((c == NULL) || (r == NULL))
943  {
944  while(r!=NULL) sm_ElemDelete(&r,_R);
945  p_Delete(&hp,_R);
946  return;
947  }
948  e = crd+1;
949  ip = piv->e;
950  do
951  {
952  a = m_act[r->pos];
953  res = dumm;
954  res->n = NULL;
955  b = c;
956  hr = r->m;
957  ir = r->e;
958  loop // combine the chains a and b: (hp,gp)*a(l) + hr*b(h)
959  {
960  if (a == NULL)
961  {
962  do
963  {
964  res = res->n = smElemCopy(b);
965  x = SM_MULT(b->m, hr, m_res[ir]->m,_R);
966  b = b->n;
967  if(ir) SM_DIV(x, m_res[ir]->m,_R);
968  res->m = x;
969  res->e = e;
970  res->f = sm_PolyWeight(res,_R);
971  } while (b != NULL);
972  break;
973  }
974  if (a->pos < b->pos)
975  {
976  res = res->n = a;
977  a = a->n;
978  }
979  else if (a->pos > b->pos)
980  {
981  res = res->n = smElemCopy(b);
982  x = SM_MULT(b->m, hr, m_res[ir]->m,_R);
983  b = b->n;
984  if(ir) SM_DIV(x, m_res[ir]->m,_R);
985  res->m = x;
986  res->e = e;
987  res->f = sm_PolyWeight(res,_R);
988  }
989  else
990  {
991  ha = a->m;
992  ia = a->e;
993  if (ir >= ia)
994  {
995  if (ir > ia)
996  {
997  x = SM_MULT(ha, m_res[ir]->m, m_res[ia]->m,_R);
998  p_Delete(&ha,_R);
999  ha = x;
1000  if (ia) SM_DIV(ha, m_res[ia]->m,_R);
1001  ia = ir;
1002  }
1003  x = SM_MULT(ha, gp, m_res[ia]->m,_R);
1004  p_Delete(&ha,_R);
1005  y = SM_MULT(b->m, hr, m_res[ia]->m,_R);
1006  }
1007  else if (ir >= ip)
1008  {
1009  if (ia < crd)
1010  {
1011  x = SM_MULT(ha, m_res[crd]->m, m_res[ia]->m,_R);
1012  p_Delete(&ha,_R);
1013  ha = x;
1014  SM_DIV(ha, m_res[ia]->m,_R);
1015  }
1016  y = hp;
1017  if(ir > ip)
1018  {
1019  y = SM_MULT(y, m_res[ir]->m, m_res[ip]->m,_R);
1020  if (ip) SM_DIV(y, m_res[ip]->m,_R);
1021  }
1022  ia = ir;
1023  x = SM_MULT(ha, y, m_res[ia]->m,_R);
1024  if (y != hp) p_Delete(&y,_R);
1025  p_Delete(&ha,_R);
1026  y = SM_MULT(b->m, hr, m_res[ia]->m,_R);
1027  }
1028  else
1029  {
1030  x = SM_MULT(hr, m_res[ia]->m, m_res[ir]->m,_R);
1031  if (ir) SM_DIV(x, m_res[ir]->m,_R);
1032  y = SM_MULT(b->m, x, m_res[ia]->m,_R);
1033  p_Delete(&x,_R);
1034  x = SM_MULT(ha, gp, m_res[ia]->m,_R);
1035  p_Delete(&ha,_R);
1036  }
1037  ha = p_Add_q(x, y,_R);
1038  if (ha != NULL)
1039  {
1040  if (ia) SM_DIV(ha, m_res[ia]->m,_R);
1041  a->m = ha;
1042  a->e = e;
1043  a->f = sm_PolyWeight(a,_R);
1044  res = res->n = a;
1045  a = a->n;
1046  }
1047  else
1048  {
1049  a->m = NULL;
1050  sm_ElemDelete(&a,_R);
1051  }
1052  b = b->n;
1053  }
1054  if (b == NULL)
1055  {
1056  res->n = a;
1057  break;
1058  }
1059  }
1060  m_act[r->pos] = dumm->n;
1061  sm_ElemDelete(&r,_R);
1062  } while (r != NULL);
1063  p_Delete(&hp,_R);
1064 }
1065 
1066 /* ----------------- transfer ------------------ */
1067 
1068 /*
1069 * select the pivotrow and store it to red and piv
1070 */
1072 {
1073  smpoly b = dumm;
1074  smpoly a, ap;
1075  int i;
1076 
1077  if (TEST_OPT_PROT)
1078  {
1079  if ((crd+1)%10)
1080  PrintS(".");
1081  else
1082  PrintS(".\n");
1083  }
1084  a = m_act[act];
1085  if (a->pos < rpiv)
1086  {
1087  do
1088  {
1089  ap = a;
1090  a = a->n;
1091  } while (a->pos < rpiv);
1092  ap->n = a->n;
1093  }
1094  else
1095  m_act[act] = a->n;
1096  piv = a;
1097  a->n = NULL;
1098  for (i=1; i<act; i++)
1099  {
1100  a = m_act[i];
1101  if (a->pos < rpiv)
1102  {
1103  loop
1104  {
1105  ap = a;
1106  a = a->n;
1107  if ((a == NULL) || (a->pos > rpiv))
1108  break;
1109  if (a->pos == rpiv)
1110  {
1111  ap->n = a->n;
1112  a->m = p_Neg(a->m,_R);
1113  b = b->n = a;
1114  b->pos = i;
1115  break;
1116  }
1117  }
1118  }
1119  else if (a->pos == rpiv)
1120  {
1121  m_act[i] = a->n;
1122  a->m = p_Neg(a->m,_R);
1123  b = b->n = a;
1124  b->pos = i;
1125  }
1126  }
1127  b->n = NULL;
1128  red = dumm->n;
1129 }
1130 
1131 /*
1132 * store the pivotcol in m_row
1133 * m_act[cols][pos(rows)] => m_row[rows][pos(cols)]
1134 */
1136 {
1137  smpoly c = m_act[act];
1138  smpoly h;
1139 
1140  while (c != NULL)
1141  {
1142  h = c;
1143  c = c->n;
1144  h->n = m_row[h->pos];
1145  m_row[h->pos] = h;
1146  h->pos = crd;
1147  }
1148 }
1149 
1150 /*
1151 * store the pivot and the assosiated row in m_row
1152 * to m_res (result):
1153 * piv + m_row[rows][pos(cols)] => m_res[cols][pos(rows)]
1154 */
1156 {
1157  smpoly r = m_row[rpiv];
1158  smpoly a, ap, h;
1159 
1160  m_row[rpiv] = NULL;
1161  perm[crd] = rpiv;
1162  piv->pos = crd;
1163  m_res[crd] = piv;
1164  while (r != NULL)
1165  {
1166  ap = m_res[r->pos];
1167  loop
1168  {
1169  a = ap->n;
1170  if (a == NULL)
1171  {
1172  ap->n = h = r;
1173  r = r->n;
1174  h->n = a;
1175  h->pos = crd;
1176  break;
1177  }
1178  ap = a;
1179  }
1180  }
1181 }
1182 
1183 /* ----------------- C++ stuff ------------------ */
1184 
1185 /*
1186 * clean m_act from zeros (after elim)
1187 */
1189 {
1190  int i = 0;
1191  int j;
1192 
1193  loop
1194  {
1195  i++;
1196  if (i > act) return;
1197  if (m_act[i] == NULL) break;
1198  }
1199  j = i;
1200  loop
1201  {
1202  j++;
1203  if (j > act) break;
1204  if (m_act[j] != NULL)
1205  {
1206  m_act[i] = m_act[j];
1207  i++;
1208  }
1209  }
1210  act -= (j-i);
1211  sign = 0;
1212 }
1213 
1214 /*
1215 * clean m_act from cols not to reduced (after elim)
1216 * put them to m_res
1217 */
1219 {
1220  int i = 0;
1221  int j;
1222 
1223  loop
1224  {
1225  i++;
1226  if (i > act) return;
1227  if (m_act[i]->pos > tored)
1228  {
1229  m_res[inred] = m_act[i];
1230  inred--;
1231  break;
1232  }
1233  }
1234  j = i;
1235  loop
1236  {
1237  j++;
1238  if (j > act) break;
1239  if (m_act[j]->pos > tored)
1240  {
1241  m_res[inred] = m_act[j];
1242  inred--;
1243  }
1244  else
1245  {
1246  m_act[i] = m_act[j];
1247  i++;
1248  }
1249  }
1250  act -= (j-i);
1251  sign = 0;
1252 }
1253 
1254 /*
1255 * copy m_act to m_res
1256 */
1258 {
1259  smpoly a,ap,r,h;
1260  int i,j,k,l;
1261 
1262  i = 0;
1263  if (act)
1264  {
1265  a = m_act[act]; // init perm
1266  do
1267  {
1268  i++;
1269  perm[crd+i] = a->pos;
1270  a = a->n;
1271  } while ((a != NULL) && (a->pos <= tored));
1272  for (j=act-1;j;j--) // load all positions of perm
1273  {
1274  a = m_act[j];
1275  k = 1;
1276  loop
1277  {
1278  if (perm[crd+k] >= a->pos)
1279  {
1280  if (perm[crd+k] > a->pos)
1281  {
1282  for (l=i;l>=k;l--) perm[crd+l+1] = perm[crd+l];
1283  perm[crd+k] = a->pos;
1284  i++;
1285  }
1286  a = a->n;
1287  if ((a == NULL) || (a->pos > tored)) break;
1288  }
1289  k++;
1290  if ((k > i) && (a->pos <= tored))
1291  {
1292  do
1293  {
1294  i++;
1295  perm[crd+i] = a->pos;
1296  a = a->n;
1297  } while ((a != NULL) && (a->pos <= tored));
1298  break;
1299  }
1300  }
1301  }
1302  }
1303  for (j=act;j;j--) // renumber m_act
1304  {
1305  k = 1;
1306  a = m_act[j];
1307  while ((a != NULL) && (a->pos <= tored))
1308  {
1309  if (perm[crd+k] == a->pos)
1310  {
1311  a->pos = crd+k;
1312  a = a->n;
1313  }
1314  k++;
1315  }
1316  }
1317  tored = crd+i;
1318  for(k=1;k<=i;k++) // clean this from m_row
1319  {
1320  j = perm[crd+k];
1321  if (m_row[j] != NULL)
1322  {
1323  r = m_row[j];
1324  m_row[j] = NULL;
1325  do
1326  {
1327  ap = m_res[r->pos];
1328  loop
1329  {
1330  a = ap->n;
1331  if (a == NULL)
1332  {
1333  h = ap->n = r;
1334  r = r->n;
1335  h->n = NULL;
1336  h->pos = crd+k;
1337  break;
1338  }
1339  ap = a;
1340  }
1341  } while (r!=NULL);
1342  }
1343  }
1344  while(act) // clean m_act
1345  {
1346  crd++;
1347  m_res[crd] = m_act[act];
1348  act--;
1349  }
1350  for (i=1;i<=tored;i++) // take the rest of m_row
1351  {
1352  if(m_row[i] != NULL)
1353  {
1354  tored++;
1355  r = m_row[i];
1356  m_row[i] = NULL;
1357  perm[tored] = i;
1358  do
1359  {
1360  ap = m_res[r->pos];
1361  loop
1362  {
1363  a = ap->n;
1364  if (a == NULL)
1365  {
1366  h = ap->n = r;
1367  r = r->n;
1368  h->n = NULL;
1369  h->pos = tored;
1370  break;
1371  }
1372  ap = a;
1373  }
1374  } while (r!=NULL);
1375  }
1376  }
1377  for (i=tored+1;i<=nrows;i++) // take the rest of m_row
1378  {
1379  if(m_row[i] != NULL)
1380  {
1381  r = m_row[i];
1382  m_row[i] = NULL;
1383  do
1384  {
1385  ap = m_res[r->pos];
1386  loop
1387  {
1388  a = ap->n;
1389  if (a == NULL)
1390  {
1391  h = ap->n = r;
1392  r = r->n;
1393  h->n = NULL;
1394  h->pos = i;
1395  break;
1396  }
1397  ap = a;
1398  }
1399  } while (r!=NULL);
1400  }
1401  }
1402  while (inred < ncols) // take unreducable
1403  {
1404  crd++;
1405  inred++;
1406  m_res[crd] = m_res[inred];
1407  }
1408 }
1409 
1410 /*
1411 * multiply and divide the column, that goes to result
1412 */
1414 {
1415  smpoly a = m_act[act];
1416  int e = crd;
1417  poly ha;
1418  int f;
1419 
1420  while (a != NULL)
1421  {
1422  f = a->e;
1423  if (f < e)
1424  {
1425  ha = SM_MULT(a->m, m_res[e]->m, m_res[f]->m,_R);
1426  p_Delete(&a->m,_R);
1427  if (f) SM_DIV(ha, m_res[f]->m,_R);
1428  a->m = ha;
1429  if (normalize) p_Normalize(a->m,_R);
1430  }
1431  a = a->n;
1432  }
1433 }
1434 
1435 /*
1436 * multiply and divide the m_act finaly
1437 */
1439 {
1440  smpoly a;
1441  poly ha;
1442  int i, f;
1443  int e = crd;
1444 
1445  for (i=act; i; i--)
1446  {
1447  a = m_act[i];
1448  do
1449  {
1450  f = a->e;
1451  if (f < e)
1452  {
1453  ha = SM_MULT(a->m, m_res[e]->m, m_res[f]->m, _R);
1454  p_Delete(&a->m,_R);
1455  if (f) SM_DIV(ha, m_res[f]->m, _R);
1456  a->m = ha;
1457  }
1458  if (normalize) p_Normalize(a->m, _R);
1459  a = a->n;
1460  } while (a != NULL);
1461  }
1462 }
1463 
1464 /*
1465 * check for denominators
1466 */
1468 {
1469  int i;
1470  smpoly a;
1471 
1472  for (i=act; i; i--)
1473  {
1474  a = m_act[i];
1475  do
1476  {
1477  if(sm_HaveDenom(a->m,_R)) return 1;
1478  a = a->n;
1479  } while (a != NULL);
1480  }
1481  return 0;
1482 }
1483 
1484 /*
1485 * normalize
1486 */
1488 {
1489  smpoly a;
1490  int i;
1491  int e = crd;
1492 
1493  for (i=act; i; i--)
1494  {
1495  a = m_act[i];
1496  do
1497  {
1498  if (e == a->e) p_Normalize(a->m,_R);
1499  a = a->n;
1500  } while (a != NULL);
1501  }
1502 }
1503 
1504 /*
1505 * multiply and divide the element, save poly
1506 */
1508 {
1509  int f = a->e;
1510  poly r, h;
1511 
1512  if (f < crd)
1513  {
1514  h = r = a->m;
1515  h = SM_MULT(h, m_res[crd]->m, m_res[f]->m, _R);
1516  if (f) SM_DIV(h, m_res[f]->m, _R);
1517  a->m = h;
1518  if (normalize) p_Normalize(a->m,_R);
1519  a->f = sm_PolyWeight(a,_R);
1520  return r;
1521  }
1522  else
1523  return NULL;
1524 }
1525 
1526 /*
1527 * delete the m_act finaly
1528 */
1530 {
1531  smpoly a;
1532  int i;
1533 
1534  for (i=act; i; i--)
1535  {
1536  a = m_act[i];
1537  do
1538  {
1539  sm_ElemDelete(&a,_R);
1540  } while (a != NULL);
1541  }
1542 }
1543 
1544 /*
1545 * delete the pivotcol
1546 */
1548 {
1549  smpoly a = m_act[act];
1550 
1551  while (a != NULL)
1552  {
1553  sm_ElemDelete(&a,_R);
1554  }
1555 }
1556 
1557 /*
1558 * delete pivot elements
1559 */
1561 {
1562  int i=crd;
1563 
1564  while (i != 0)
1565  {
1566  sm_ElemDelete(&m_res[i],_R);
1567  i--;
1568  }
1569 }
1570 
1571 /*
1572 * the sign of the determinant
1573 */
1575 {
1576  int j,i;
1577  if (act > 2)
1578  {
1579  if (cpiv!=act) sign=-sign;
1580  if ((act%2)==0) sign=-sign;
1581  i=1;
1582  j=perm[1];
1583  while(j<rpiv)
1584  {
1585  sign=-sign;
1586  i++;
1587  j=perm[i];
1588  }
1589  while(perm[i]!=0)
1590  {
1591  perm[i]=perm[i+1];
1592  i++;
1593  }
1594  }
1595  else
1596  {
1597  if (cpiv!=1) sign=-sign;
1598  if (rpiv!=perm[1]) sign=-sign;
1599  }
1600 }
1601 
1603 {
1604  int i;
1605  for (i=act;i;i--) perm[i]=i;
1606 }
1607 
1608 /* ----------------- arithmetic ------------------ */
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 
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
1887 /*n
1888 * exact division a/b
1889 * a is a result of smMultDiv
1890 * a destroyed, b NOT destroyed
1891 */
1892 void sm_SpecialPolyDiv(poly a, poly b, const ring R)
1893 {
1894  if (pNext(b) == NULL)
1895  {
1896  sm_PolyDivN(a, pGetCoeff(b),R);
1897  return;
1898  }
1899  sm_ExactPolyDiv(a, b, R);
1900 }
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);
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 */
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);
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);
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:441
static FORCE_INLINE number n_Sub(number a, number b, const coeffs r)
return the difference of 'a' and 'b', i.e., a-b
Definition: coeffs.h:668
static poly sm_Smpoly2Poly(smpoly, const ring)
Definition: sparsmat.cc:2187
float f
Definition: sparsmat.cc:59
#define omAllocBin(bin)
Definition: omAllocDecl.h:205
static FORCE_INLINE BOOLEAN n_Greater(number a, number b, const coeffs r)
ordered fields: TRUE iff 'a' is larger than 'b'; in Z/pZ: TRUE iff la > lb, where la and lb are the l...
Definition: coeffs.h:512
int smGetRed()
Definition: sparsmat.cc:179
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:263
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:1892
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:937
void kBucketInit(kBucket_pt bucket, poly lm, int length)
Definition: kbuckets.cc:471
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:1613
void smWeights()
Definition: sparsmat.cc:664
smpoly * m_res
Definition: sparsmat.cc:144
smpoly dumm
Definition: sparsmat.cc:148
static poly pp_Mult_Coeff_mm_DivSelect_MultDiv(poly p, int &lp, poly m, poly a, poly b, const ring currRing)
Definition: sparsmat.cc:82
#define TEST_OPT_PROT
Definition: options.h:98
loop
Definition: myNF.cc:98
#define FALSE
Definition: auxiliary.h:140
return P p
Definition: myNF.cc:203
static BOOLEAN p_LmDivisibleByNoComp(poly a, poly b, const ring r)
Definition: p_polys.h:1662
static FORCE_INLINE BOOLEAN n_IsOne(number n, const coeffs r)
TRUE iff 'n' represents the one element.
Definition: coeffs.h:469
void smRowToCol()
Definition: sparsmat.cc:1155
void smNormalize()
Definition: sparsmat.cc:1487
omBin sip_sideal_bin
Definition: simpleideals.cc:30
static unsigned long p_SetComp(poly p, unsigned long c, ring r)
Definition: p_polys.h:236
void smPivot()
Definition: sparsmat.cc:696
poly sm_CallDet(ideal I, const ring R)
Definition: sparsmat.cc:359
int smGetSign()
Definition: sparsmat.cc:177
#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:539
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
float wpoints
Definition: sparsmat.cc:141
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:540
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
int smCheckNormalize()
Definition: sparsmat.cc:1467
omBin smprec_bin
Definition: sparsmat.cc:80
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:962
#define TRUE
Definition: auxiliary.h:144
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:1529
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:579
void * ADDRESS
Definition: auxiliary.h:161
sm_prec * smpoly
Definition: sparsmat.cc:52
void WerrorS(const char *s)
Definition: feFopen.cc:23
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:194
smnumber * m_act
Definition: sparsmat.cc:2339
kBucket_pt kBucketCreate(ring bucket_ring)
Creation/Destruction of buckets.
Definition: kbuckets.cc:197
static number p_SetCoeff(poly p, number n, ring r)
Definition: p_polys.h:401
static void p_LmFree(poly p, ring)
Definition: p_polys.h:679
int normalize
Definition: sparsmat.cc:139
static int pLength(poly a)
Definition: p_polys.h:189
poly kBucketExtractLm(kBucket_pt bucket)
Definition: kbuckets.cc:489
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:609
#define TEST_OPT_NOT_BUCKETS
Definition: options.h:100
int e
Definition: sparsmat.cc:57
smpoly red
Definition: sparsmat.cc:146
#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 'a' and 'b', i.e., a*b
Definition: coeffs.h:635
#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:12
static poly sm_Smnumber2Poly(number, const ring)
Definition: sparsmat.cc:2929
void smSelectPR()
Definition: sparsmat.cc:1071
void smNewBareiss(int, int)
Definition: sparsmat.cc:603
sm_nrec * smnumber
Definition: sparsmat.cc:2310
const ring r
Definition: syzextra.cc:208
smpoly n
Definition: sparsmat.cc:55
void kBucketDestroy(kBucket_pt *bucket_pt)
Definition: kbuckets.cc:204
Definition: intvec.h:16
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:3435
void smNewPivot()
Definition: sparsmat.cc:791
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:465
int j
Definition: myNF.cc:70
static BOOLEAN p_IsConstant(const poly p, const ring r)
Definition: p_polys.h:1784
static FORCE_INLINE number n_Add(number a, number b, const coeffs r)
return the sum of 'a' and 'b', i.e., a+b
Definition: coeffs.h:655
int pos
Definition: sparsmat.cc:56
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:1318
static BOOLEAN smCheckSolv(ideal)
Definition: sparsmat.cc:2942
static poly pp_Mult_qq(poly p, poly q, const ring r)
Definition: p_polys.h:1075
BOOLEAN sm_CheckDet(ideal I, int d, BOOLEAN sw, const ring r)
Definition: sparsmat.cc:308
static BOOLEAN sm_HaveDenom(poly, const ring)
Definition: sparsmat.cc:2257
void smColDel()
Definition: sparsmat.cc:1547
float * wrw
Definition: sparsmat.cc:142
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:1472
static FORCE_INLINE number n_Invers(number a, const coeffs r)
return the multiplicative inverse of 'a'; raise an error if 'a' is not invertible ...
Definition: coeffs.h:565
#define SM_MIN_LENGTH_BUCKET
Definition: sparsmat.cc:46
int m
Definition: cfEzgcd.cc:119
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:558
smpoly piv
Definition: sparsmat.cc:147
FILE * f
Definition: checklibs.c:7
int i
Definition: cfEzgcd.cc:123
void PrintS(const char *s)
Definition: reporter.cc:294
static poly p_Mult_nn(poly p, number n, const ring r)
Definition: p_polys.h:902
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:461
static void p_ExpVectorAdd(poly p1, poly p2, const ring r)
Definition: p_polys.h:1339
#define IDELEMS(i)
Definition: simpleideals.h:24
static FORCE_INLINE BOOLEAN n_IsZero(number n, const coeffs r)
TRUE iff 'n' represents the zero element.
Definition: coeffs.h:465
int * perm
Definition: sparsmat.cc:140
#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:1560
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:143
smpoly * m_row
Definition: sparsmat.cc:145
void p_Normalize(poly p, const ring r)
Definition: p_polys.cc:3621
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:850
#define omGetSpecBin(size)
Definition: omBin.h:11
void smHElim()
Definition: sparsmat.cc:932
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:38
void sm1Elim()
Definition: sparsmat.cc:853
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:1402
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:484
float * wcl
Definition: sparsmat.cc:142
#define NULL
Definition: omList.c:10
ideal smRes2Mod()
Definition: sparsmat.cc:505
static FORCE_INLINE number n_Copy(number n, const coeffs r)
return a copy of 'n'
Definition: coeffs.h:452
void sm_KillModifiedRing(ring r)
Definition: sparsmat.cc:294
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 'a' and 'b', i.e., a/b; raises an error if 'b' is not invertible in r exceptio...
Definition: coeffs.h:615
int rows() const
Definition: intvec.h:88
static smnumber sm_Poly2Smnumber(poly, const ring)
Definition: sparsmat.cc:2898
#define R
Definition: sirandom.c:26
void smInitPerm()
Definition: sparsmat.cc:1602
poly smDet()
Definition: sparsmat.cc:529
void smColToRow()
Definition: sparsmat.cc:1135
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:443
Variable x
Definition: cfModGcd.cc:4023
static BOOLEAN p_IsConstantPoly(const poly p, const ring r)
Definition: p_polys.h:1798
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:604
#define pNext(p)
Definition: monomials.h:43
static FORCE_INLINE BOOLEAN n_Equal(number a, number b, const coeffs r)
TRUE iff 'a' and 'b' represent the same number; they may have different representations.
Definition: coeffs.h:461
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:436
void smSign()
Definition: sparsmat.cc:1574
#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:241
#define SM_DIV
Definition: sparsmat.h:24
void smMultCol()
Definition: sparsmat.cc:1413
void sm_CallBareiss(ideal I, int x, int y, ideal &M, intvec **iv, const ring R)
Definition: sparsmat.cc:404
static void sm_ExactPolyDiv(poly, poly, const ring)
Definition: sparsmat.cc:1904
void smCopToRes()
Definition: sparsmat.cc:1257
static poly p_Neg(poly p, const ring r)
Definition: p_polys.h:1018
static void p_LmDelete(poly p, const ring r)
Definition: p_polys.h:707
static omBin smnrec_bin
Definition: sparsmat.cc:2317
static float sm_PolyWeight(smpoly, const ring)
Definition: sparsmat.cc:2228
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
Definition: coeffs.h:456
void smToredElim()
Definition: sparsmat.cc:1218
static FORCE_INLINE BOOLEAN n_GreaterZero(number n, const coeffs r)
ordered fields: TRUE iff 'n' 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)/: 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:495
void smZeroElim()
Definition: sparsmat.cc:1188
polyrec * poly
Definition: hilb.h:10
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:884
static poly sm_SelectCopy_ExpMultDiv(poly p, poly m, poly a, poly b, const ring currRing)
Definition: sparsmat.cc:104
#define omFreeBin(addr, bin)
Definition: omAllocDecl.h:259
static BOOLEAN rOrd_is_Comp_dp(const ring r)
Definition: ring.h:771
void smFinalMult()
Definition: sparsmat.cc:1438
static poly p_New(const ring, omBin bin)
Definition: p_polys.h:660
poly m
Definition: sparsmat.cc:58
static Poly * h
Definition: janet.cc:978
int BOOLEAN
Definition: auxiliary.h:131
void smSparseHomog()
static poly p_Init(const ring r, omBin bin)
Definition: p_polys.h:1248
const poly b
Definition: syzextra.cc:213
poly smMultPoly(smpoly)
Definition: sparsmat.cc:1507
poly p_Cleardenom(poly p, const ring r)
Definition: p_polys.cc:2682
static void sm_ElemDelete(smpoly *, const ring)
Definition: sparsmat.cc:2129
void smNewWeights()
Definition: sparsmat.cc:753
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:571
void smToIntvec(intvec *)
Definition: sparsmat.cc:518
void Werror(const char *fmt,...)
Definition: reporter.cc:199
int pos
Definition: sparsmat.cc:2313
smpoly oldpiv
Definition: sparsmat.cc:147
static poly pp_Mult_Coeff_mm_DivSelect(poly p, const poly m, const ring r)
Definition: p_polys.h:1001
#define omAlloc0(size)
Definition: omAllocDecl.h:211
int l
Definition: cfEzgcd.cc:94
smpoly * smGetAct()
Definition: sparsmat.cc:178
void kBucket_Add_q(kBucket_pt bucket, poly q, int *l)
Add to Bucket a poly ,i.e. Bpoly == q+Bpoly.
Definition: kbuckets.cc:636
static BOOLEAN sm_IsNegQuot(poly, const poly, const poly, const ring)
Definition: sparsmat.cc:1959