source: git/Singular/sparsmat.cc @ 48aa42

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