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