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