source: git/Singular/sparsmat.cc @ c232af

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