source: git/Singular/sparsmat.cc @ 512a2b

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