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