source: git/Singular/sparsmat.cc @ a5ceaf

spielwiese
Last change on this file since a5ceaf was a5ceaf, checked in by Wilfred Pohl <pohl@…>, 24 years ago
det for sparse git-svn-id: file:///usr/local/Singular/svn/trunk@3088 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 37.3 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: sparsmat.cc,v 1.6 1999-06-02 14:09:34 pohl Exp $ */
5
6/*
7* ABSTRACT: operations with sparse matrices (bareiss, ...)
8*/
9
10#include "mod2.h"
11#include "tok.h"
12#include "febase.h"
13#include "structs.h"
14#include "intvec.h"
15#include "lists.h"
16#include "ring.h"
17#include "polys.h"
18#include "ipid.h"
19#include "ideals.h"
20#include "numbers.h"
21#include "sparsmat.h"
22
23/* ----------------- macros ------------------ */
24/* #define OLD_DIV */
25#ifdef OLD_DIV
26#define SM_MULT(A,B,C) smMult(A,B)
27#define SM_DIV smPolyDiv
28#else
29#define SM_MULT smMultDiv
30#define SM_DIV smSpecialPolyDiv
31#endif
32/* ----------------- general definitions ------------------ */
33
34/* elements of a 'sparse' matrix in SINGULAR */
35typedef struct smprec sm_prec;
36typedef sm_prec * smpoly;
37struct smprec{
38  smpoly n;            // the next element
39  int pos;             // position
40  int e;               // level
41  poly m;              // the element
42  float f;             // complexity of the element
43};
44
45/* declare internal 'C' stuff */
46static void smDebug() { int k=0; }
47static void smExactPolyDiv(poly, poly);
48static BOOLEAN smIsScalar(const poly);
49static BOOLEAN smIsNegQuot(poly, const poly, const poly);
50static poly smEMult(poly, const poly);
51static BOOLEAN smCheckLead(const poly, const poly);
52static poly smDMult(poly, const poly);
53static void smPolyDivN(poly a, const number);
54static BOOLEAN smSmaller(poly, poly);
55static void smCombineChain(poly *, poly);
56static void smFindRef(poly *, poly *, poly);
57
58static void smElemDelete(smpoly *);
59static smpoly smElemCopy(smpoly);
60static float smPolyWeight(smpoly);
61static smpoly smPoly2Smpoly(poly p);
62static poly smSmpoly2Poly(smpoly a);
63
64/* class for sparse matrix:
65*      3 parts of matrix during the algorithm
66*      m_act[cols][pos(rows)] => m_row[rows][pos(cols)] => m_res[cols][pos(rows)]
67*      input                     pivotcols as rows         result
68*      pivot                     like a stack              from pivot and pivotcols
69*      elimination                                         rows reordered according
70*                                                          to pivot choise
71*                                                          stored in perm
72*      a step is as follows
73*      - search of pivot (smPivot)
74*      - swap pivot column and last column (smSwapC)
75*        select pivotrow to piv and red (smSelectPR)
76*        consider sign
77*      - elimination (smInitElim, sm0Elim, sm1Elim)
78*        clear zero column as result of elimination (smZeroElim)
79*      - tranfer from
80*        piv and m_row to m_res (smRowToCol)
81*        m_act to m_row (smColToRow)
82*/
83class sparse_mat{
84private:
85  int nrows, ncols;    // dimension of the problem
86  int sign;            // for determinant (start: 1)
87  int act;             // number of unreduced columns (start: ncols)
88  int crd;             // number of reduced columns (start: 0)
89  int tored;           // border for rows to reduce
90  int rpiv, cpiv;      // position of the pivot
91  unsigned short *perm;// permutation of rows
92  float wpoints;       // weight of all points
93  float *wrw, *wcl;    // weights of rows and columns
94  smpoly * m_act;      // unreduced columns
95  smpoly * m_res;      // reduced columns (result)
96  smpoly * m_row;      // reduced part of rows
97  smpoly red;          // row to reduce
98  smpoly piv, oldpiv;  // pivot and previous pivot
99  smpoly dumm;         // allocated dummy
100  void smColToRow();
101  void smRowToCol();
102  void smFinalMult();
103  void smSparseHomog();
104  void smWeights();
105  void smPivot();
106  void smNewWeights();
107  void smNewPivot();
108  void smZeroElim();
109  void smToredElim();
110  void smCopToRes();
111  void smSelectPR();
112  void smElim();
113  void sm1Elim();
114  void smHElim();
115  void smMultCol();
116  poly smMultPoly(smpoly);
117  void smActDel();
118  void smColDel();
119  void smPivDel();
120  void smSign();
121public:
122  sparse_mat(ideal);
123  ~sparse_mat();
124  int smGetSign() { return sign; }
125  smpoly * smGetAct() { return m_act; }
126  int smGetRed() { return tored; }
127  ideal smRes2Mod();
128  poly smDet();
129  void smBareiss(int, int);
130  void smNewBareiss(int, int);
131  void smToIntvec(intvec *);
132};
133
134/* ----------------- ops with rings ------------------ */
135ideal smRingCopy(ideal I, ring *ri, sip_sring &tmpR)
136{
137  ring origR =NULL;
138  ideal II;
139  if (currRing->order[0]!=ringorder_c)
140  {
141    origR =currRing;
142    tmpR=*origR;
143    int *ord=(int*)Alloc0(3*sizeof(int));
144    int *block0=(int*)Alloc(3*sizeof(int));
145    int *block1=(int*)Alloc(3*sizeof(int));
146    ord[0]=ringorder_c;
147    ord[1]=ringorder_dp;
148    tmpR.order=ord;
149    block0[1]=1;
150    tmpR.block0=block0;
151    block1[1]=tmpR.N;
152    tmpR.block1=block1;
153    rComplete(&tmpR,1);
154    rChangeCurrRing(&tmpR,TRUE);
155    // fetch data from the old ring
156    II=idInit(IDELEMS(I),I->rank);
157    int k;
158    for (k=0;k<IDELEMS(I);k++) II->m[k] = pFetchCopy(origR, I->m[k]);
159  }
160  else
161  {
162    II=idCopy(I);
163  }
164  *ri = origR;
165  return II;
166}
167
168void smRingClean(ring origR, ip_sring &tmpR)
169{
170  rChangeCurrRing(origR,TRUE);
171  rUnComplete(&tmpR);
172  Free((ADDRESS)tmpR.order,3*sizeof(int));
173  Free((ADDRESS)tmpR.block0,3*sizeof(int));
174  Free((ADDRESS)tmpR.block1,3*sizeof(int));
175}
176
177/* ----------------- basics (used from 'C') ------------------ */
178/*2
179*returns the determinant of the module I;
180*uses Bareiss algorithm
181*/
182poly smCallDet(ideal I)
183{
184  int r=idRankFreeModule(I);
185  if (I->ncols != r)
186  {
187    Werror("det of %d x %d module (matrix)",r,I->ncols);
188    return NULL;
189  }
190  poly res,save;
191  ring origR;
192  sip_sring tmpR;
193  ideal II=smRingCopy(I,&origR,tmpR);
194  sparse_mat *det = new sparse_mat(II);
195
196  idDelete(&II);
197  if (det->smGetAct() == NULL)
198  {
199    delete det;
200    if (origR!=NULL) smRingClean(origR,tmpR);
201    return NULL;
202  }
203  res=det->smDet();
204  if(det->smGetSign()<0) res=pNeg(res);
205  delete det;
206  if (origR!=NULL)
207  {
208    rChangeCurrRing(origR,TRUE);
209    save = res;
210    res = pFetchCopy(&tmpR, save);
211    rChangeCurrRing(&tmpR,FALSE);
212    pDelete(&save);
213    smRingClean(origR,tmpR);
214  }
215  return res;
216}
217
218lists smCallBareiss(ideal I, int x, int y)
219{
220  lists res=(lists)Alloc(sizeof(slists));
221  ideal II = idCopy(I);
222  sparse_mat *bareiss = new sparse_mat(II);
223  intvec *v;
224  ideal m;
225
226  idDelete(&II);
227  if (bareiss->smGetAct() == NULL)
228  {
229    Free((ADDRESS)res, sizeof(slists));
230    return NULL;
231  }
232  bareiss->smBareiss(x, y);
233  m = bareiss->smRes2Mod();
234  v = new intvec(bareiss->smGetRed());
235  bareiss->smToIntvec(v);
236  delete bareiss;
237  res->Init(2);
238  res->m[0].rtyp=MODUL_CMD;
239  res->m[0].data=(void *)m;
240  res->m[1].rtyp=INTVEC_CMD;
241  res->m[1].data=(void *)v;
242  return res;
243}
244
245lists smCallNewBareiss(ideal I, int x, int y)
246{
247  ring origR;
248  sip_sring tmpR;
249  lists res=(lists)Alloc(sizeof(slists));
250  ideal II=smRingCopy(I,&origR,tmpR);
251  sparse_mat *bareiss = new sparse_mat(II);
252  ideal mm=II;
253  intvec *v;
254  ideal m;
255
256  if (bareiss->smGetAct() == NULL)
257  {
258    delete bareiss;
259    if (origR!=NULL) smRingClean(origR,tmpR);
260    v=new intvec(1,pVariables);
261  }
262  else
263  {
264    idDelete(&II);
265    bareiss->smNewBareiss(x, y);
266    m = bareiss->smRes2Mod();
267    v = new intvec(bareiss->smGetRed());
268    bareiss->smToIntvec(v);
269    delete bareiss;
270    if (origR!=NULL)
271    {
272      rChangeCurrRing(origR,TRUE);
273      mm=idInit(IDELEMS(m),m->rank);
274      int k;
275      for (k=0;k<IDELEMS(m);k++) mm->m[k] = pFetchCopy(&tmpR, m->m[k]);
276      rChangeCurrRing(&tmpR,FALSE);
277      idDelete(&m);
278      smRingClean(origR,tmpR);
279    }
280    else
281    {
282      mm=m;
283    }
284  }
285  res->Init(2);
286  res->m[0].rtyp=MODUL_CMD;
287  res->m[0].data=(void *)mm;
288  res->m[1].rtyp=INTVEC_CMD;
289  res->m[1].data=(void *)v;
290  return res;
291}
292
293/*
294* constructor
295*/
296sparse_mat::sparse_mat(ideal smat)
297{
298  int i;
299  polyset pmat;
300
301  ncols = smat->ncols;
302  nrows = idRankFreeModule(smat);
303  if (nrows <= 0)
304  {
305    m_act = NULL;
306    return;
307  }
308  sign = 1;
309  act = ncols;
310  crd = 0;
311  tored = nrows; // without border
312  i = tored+1;
313  perm = (unsigned short *)Alloc(sizeof(unsigned short)*(i+1));
314  perm[i] = 0;
315  m_row = (smpoly *)Alloc0(sizeof(smpoly)*i);
316  wrw = (float *)Alloc(sizeof(float)*i);
317  i = ncols+1;
318  wcl = (float *)Alloc(sizeof(float)*i);
319  m_act = (smpoly *)Alloc(sizeof(smpoly)*i);
320  m_res = (smpoly *)Alloc0(sizeof(smpoly)*i);
321  dumm = (smpoly)Alloc(sizeof(smprec));
322  m_res[0] = (smpoly)Alloc(sizeof(smprec));
323  m_res[0]->m = NULL;
324  pmat = smat->m;
325  for(i=ncols; i; i--)
326  {
327    m_act[i] = smPoly2Smpoly(pmat[i-1]);
328    pmat[i-1] = NULL;
329  }
330  this->smZeroElim();
331  oldpiv = NULL;
332}
333
334/*
335* destructor
336*/
337sparse_mat::~sparse_mat()
338{
339  int i;
340  if (m_act == NULL) return;
341  Free((ADDRESS)m_res[0], sizeof(smprec));
342  Free((ADDRESS)dumm, sizeof(smprec));
343  i = ncols+1;
344  Free((ADDRESS)m_res, sizeof(smpoly)*i);
345  Free((ADDRESS)m_act, sizeof(smpoly)*i);
346  Free((ADDRESS)wcl, sizeof(float)*i);
347  i = nrows+1;
348  Free((ADDRESS)wrw, sizeof(float)*i);
349  Free((ADDRESS)m_row, sizeof(smpoly)*i);
350  Free((ADDRESS)perm, sizeof(unsigned short)*(i+1));
351}
352
353/*
354* transform the result to a module
355*/
356ideal sparse_mat::smRes2Mod()
357{
358  ideal res = idInit(crd, crd);
359  int i;
360
361  for (i=crd; i; i--) res->m[i-1] = smSmpoly2Poly(m_res[i]);
362  res->rank = idRankFreeModule(res);
363  return res;
364}
365
366/*
367* permutation of rows
368*/
369void sparse_mat::smToIntvec(intvec *v)
370{
371  int i;
372
373  for (i=v->rows()-1; i>=0; i--)
374    (*v)[i] = perm[i+1];
375}
376
377/* ---------------- the algorithm's ------------------ */
378/*
379* the determinant (up to sign), uses new Bareiss elimination
380*/
381poly sparse_mat::smDet()
382{
383  int y;
384  poly res = NULL;
385
386  if (sign == 0)
387  {
388    this->smActDel();
389    return NULL;
390  }
391  if (act < 2)
392  {
393    if (act != 0) res = m_act[1]->m;
394    Free((void *)m_act[1], sizeof(smprec));
395    return res;
396  }
397  this->smPivot();
398  this->smSign();
399  this->smSelectPR();
400  this->sm1Elim();
401  crd++;
402  this->smColDel();
403  act--;
404  this->smZeroElim();
405  if (sign == 0)
406  {
407    this->smActDel();
408    return NULL;
409  }
410  if (act < 2)
411  {
412    if (act != 0) res = m_act[1]->m;
413    Free((void *)m_act[1], sizeof(smprec));
414    return res;
415  }
416  loop
417  {
418    m_res[crd] = piv;
419    this->smNewPivot();
420    this->smSign();
421    this->smSelectPR();
422    this->smMultCol();
423    this->smHElim();
424    crd++;
425    this->smColDel();
426    act--;
427    this->smZeroElim();
428    if (sign == 0)
429    {
430      this->smPivDel();
431      this->smActDel();
432      return NULL;
433    }
434    if (act < 2)
435    {
436      this->smPivDel();
437      if (act != 0) res = m_act[1]->m;
438      Free((void *)m_act[1], sizeof(smprec));   
439      return res;
440    }
441  }
442}
443
444/*
445* the Bareiss elimination:
446*   - with x unreduced last rows, pivots from here are not allowed
447*   - the method will finish for number of unreduced columns < y
448*/
449void sparse_mat::smBareiss(int x, int y)
450{
451  if ((x > 0) && (x < nrows))
452  {
453    tored -= x;
454    this->smToredElim();
455  }
456  if (y < 2) y = 2;
457  if (act < y)
458  {
459    if (act != 0)
460    {
461      this->smCopToRes();
462    }
463    return;
464  }
465  this->smPivot();
466  this->smSelectPR();
467  this->smElim();
468  crd++;
469  this->smColToRow();
470  act--;
471  this->smRowToCol();
472  this->smZeroElim();
473  if (tored != nrows)
474    this->smToredElim();
475  if (act < y)
476  {
477    this->smCopToRes();
478    return;
479  }
480  loop
481  {
482    this->smPivot();
483    oldpiv = piv;
484    this->smSelectPR();
485    this->smElim();
486    crd++;
487    this->smColToRow();
488    act--;
489    this->smRowToCol();
490    this->smZeroElim();
491    if (tored != nrows)
492      this->smToredElim();
493    if (act < y)
494    {
495      if (act != 0)
496      {
497        this->smCopToRes();
498      }
499      else
500        tored = crd;
501      return;
502    }
503  }
504}
505
506/*
507* the new Bareiss elimination:
508*   - with x unreduced last rows, pivots from here are not allowed
509*   - the method will finish for number of unreduced columns < y
510*/
511void sparse_mat::smNewBareiss(int x, int y)
512{
513  if ((x > 0) && (x < nrows))
514  {
515    tored -= x;
516    this->smToredElim();
517  }
518  if (y < 2) y = 2;
519  if (act < y)
520  {
521    if (act != 0)
522    {
523      this->smCopToRes();
524    }
525    return;
526  }
527  this->smPivot();
528  this->smSelectPR();
529  this->sm1Elim();
530  crd++;
531  this->smColToRow();
532  act--;
533  this->smRowToCol();
534  this->smZeroElim();
535  if (tored != nrows)
536    this->smToredElim();
537  if (act < y)
538  {
539    this->smFinalMult();
540    this->smCopToRes();
541    return;
542  }
543  loop
544  {
545    this->smNewPivot();
546    this->smSelectPR();
547    this->smMultCol();
548    this->smHElim();
549    crd++;
550    this->smColToRow();
551    act--;
552    this->smRowToCol();
553    this->smZeroElim();
554    if (tored != nrows)
555      this->smToredElim();
556    if (act < y)
557    {
558      if (act != 0)
559      {
560        this->smFinalMult();
561        this->smCopToRes();
562      }
563      else
564        tored = crd;
565      return;
566    }
567  }
568}
569
570/* ----------------- pivot method ------------------ */
571
572/*
573* prepare smPivot, compute weights for rows and columns
574* and the weight for all points
575*/
576void sparse_mat::smWeights()
577{
578  float wc, wp, w;
579  smpoly a;
580  int i;
581
582  wp = 0.0;
583  for (i=tored; i; i--) wrw[i] = 0.0; // ???
584  for (i=act; i; i--)
585  {
586    wc = 0.0;
587    a = m_act[i];
588    loop
589    {
590      if (a->pos > tored)
591        break;
592      w = a->f = smPolyWeight(a);
593      wc += w;
594      wrw[a->pos] += w;
595      a = a->n;
596      if (a == NULL)
597        break;
598    }
599    wp += wc;
600    wcl[i] = wc;
601  }
602  wpoints = wp;
603}
604
605/*
606* compute pivot
607*/
608void sparse_mat::smPivot()
609{
610  float wopt = 1.0e30;
611  float wc, wr, wp, w;
612  smpoly a;
613  int i, copt, ropt;
614
615  this->smWeights();
616  for (i=act; i; i--)
617  {
618    a = m_act[i];
619    loop
620    {
621      if (a->pos > tored)
622        break;
623      w = a->f;
624      wc = wcl[i]-w;
625      wr = wrw[a->pos]-w;
626      if ((wr<0.25) || (wc<0.25)) // row or column with only one point
627      {
628        if (w<wopt)
629        {
630          wopt = w;
631          copt = i;
632          ropt = a->pos;
633        }
634      }
635      else // elimination
636      {
637        wp = w*(wpoints-wcl[i]-wr);
638        wp += wr*wc;
639        if (wp < wopt)
640        {
641          wopt = wp;
642          copt = i;
643          ropt = a->pos;
644        } 
645      }
646      a = a->n;
647      if (a == NULL)
648        break;
649    }
650  }
651  rpiv = ropt;
652  cpiv = copt;
653  if (cpiv != act)
654  {
655    a = m_act[act];
656    m_act[act] = m_act[cpiv];
657    m_act[cpiv] = a;
658  }
659}
660
661/*
662* prepare smPivot, compute weights for rows and columns
663* and the weight for all points
664*/
665void sparse_mat::smNewWeights()
666{
667  float wc, wp, w, hp = piv->f;
668  smpoly a;
669  int i, f, e = crd;
670
671  wp = 0.0;
672  for (i=tored; i; i--) wrw[i] = 0.0; // ???
673  for (i=act; i; i--)
674  {
675    wc = 0.0;
676    a = m_act[i];
677    loop
678    {
679      if (a->pos > tored)
680        break;
681      w = a->f;
682      f = a->e;
683      if (f < e)
684      {
685        w *= hp;
686        if (f) w /= m_res[f]->f;
687      }
688      wc += w;
689      wrw[a->pos] += w;
690      a = a->n;
691      if (a == NULL)
692        break;
693    }
694    wp += wc;
695    wcl[i] = wc;
696  }
697  wpoints = wp;
698}
699
700/*
701* compute pivot
702*/
703void sparse_mat::smNewPivot()
704{
705  float wopt = 1.0e30, hp = piv->f;
706  float wc, wr, wp, w;
707  smpoly a;
708  int i, copt, ropt, f, e = crd;
709
710  this->smNewWeights();
711  for (i=act; i; i--)
712  {
713    a = m_act[i];
714    loop
715    {
716      if (a->pos > tored)
717        break;
718      w = a->f;
719      f = a->e;
720      if (f < e)
721      {
722        w *= hp;
723        if (f) w /= m_res[f]->f;
724      }
725      wc = wcl[i]-w;
726      wr = wrw[a->pos]-w;
727      if ((wr<0.25) || (wc<0.25)) // row or column with only one point
728      {
729        if (w<wopt)
730        {
731          wopt = w;
732          copt = i;
733          ropt = a->pos;
734        }
735      }
736      else // elimination
737      {
738        wp = w*(wpoints-wcl[i]-wr);
739        wp += wr*wc;
740        if (wp < wopt)
741        {
742          wopt = wp;
743          copt = i;
744          ropt = a->pos;
745        } 
746      }
747      a = a->n;
748      if (a == NULL)
749        break;
750    }
751  }
752  rpiv = ropt;
753  cpiv = copt;
754  if (cpiv != act)
755  {
756    a = m_act[act];
757    m_act[act] = m_act[cpiv];
758    m_act[cpiv] = a;
759  }
760}
761
762/* ----------------- elimination ------------------ */
763
764/* steps of elimination */
765void sparse_mat::smElim()
766{
767  poly p = piv->m;        // pivotelement
768  smpoly c = m_act[act];  // pivotcolumn
769  smpoly r = red;         // row to reduce
770  poly q;
771  smpoly res, a, b;
772  poly w, ha, hb;
773  int i;
774
775  if (oldpiv != NULL) q = oldpiv->m; // previous pivot
776  else q = NULL; 
777  if ((c == NULL) || (r == NULL))
778  {
779    while (r) smElemDelete(&r);
780    for (i=1; i<act; i++)
781    {
782      a = m_act[i];
783      do
784      {
785        ha = SM_MULT(a->m, p, q);
786        pDelete(&a->m);
787        if (q) SM_DIV(ha, q);
788        a->m = ha;
789        a = a->n;
790      } while (a != NULL);
791    }
792    return;
793  }
794  for (i=1; i<act; i++)
795  {
796    a = m_act[i];
797    if ((r == NULL) || (i != r->pos))  // cols without elimination
798    {
799      do
800      {
801        ha = SM_MULT(a->m, p, q);
802        pDelete(&a->m);
803        if (q) SM_DIV(ha, q);
804        a->m = ha;
805        a = a->n;
806      } while (a != NULL);
807    }
808    else                    // cols with elimination
809    {
810      res = dumm;
811      res->n = NULL;
812      b = c;
813      w = r->m;
814      loop                  // combine the chains a and b: p*a + w*b
815      {
816        if (a->pos < b->pos)
817        {
818          ha = SM_MULT(a->m, p, q);
819          pDelete(&a->m);
820          if (q) SM_DIV(ha, q);
821          a->m = ha;
822          res = res->n = a;
823          a = a->n;
824        }
825        else if (a->pos > b->pos)
826        {
827          res = res->n = smElemCopy(b);
828          hb = SM_MULT(b->m, w, q);
829          b = b->n;
830          if (q) SM_DIV(hb, q);
831          res->m = hb;
832        }     
833        else
834        {
835          ha = SM_MULT(a->m, p, q);
836          pDelete(&a->m);
837          hb = SM_MULT(b->m, w, q);
838          ha = pAdd(ha, hb);
839          if (ha != NULL)
840          {
841            if (q) SM_DIV(ha, q);
842            a->m = ha;
843            res = res->n = a;
844            a = a->n;
845          }
846          else
847          {
848            smElemDelete(&a);
849          }
850          b = b->n;
851        }
852        if (a == NULL)
853        {
854          if (b != NULL)
855          {
856            do
857            {
858              res = res->n = smElemCopy(b);
859              hb = SM_MULT(b->m, w, q);
860              if (q) SM_DIV(hb, q);
861              res->m = hb;
862              b = b->n;
863            } while (b != NULL);
864          }
865          else
866            res->n = NULL;
867          break;
868        }
869        if (b == NULL)
870        {
871          do
872          {
873            ha = SM_MULT(a->m, p, q);
874            pDelete(&a->m);
875            if (q) SM_DIV(ha, q);
876            a->m = ha;
877            res = res->n = a;
878            a = a->n;
879          } while (a != NULL);
880          break;
881        }
882      }
883      m_act[i] = dumm->n;
884      if (r) smElemDelete(&r);
885    }
886  }
887}
888
889/* first step of elimination */
890void sparse_mat::sm1Elim()
891{
892  poly p = piv->m;        // pivotelement
893  smpoly c = m_act[act];  // pivotcolumn
894  smpoly r = red;         // row to reduce
895  smpoly res, a, b;
896  poly w, ha, hb;
897
898  if ((c == NULL) || (r == NULL))
899    return; 
900  do
901  {
902    a = m_act[r->pos];
903    res = dumm;
904    res->n = NULL;
905    b = c;
906    w = r->m;
907    loop                  // combine the chains a and b: p*a + w*b
908    {
909      if (a->pos < b->pos)
910      {
911        res = res->n = a;
912        a = a->n;
913      }
914      else if (a->pos > b->pos)
915      {
916        res = res->n = smElemCopy(b);
917        res->m = smMult(b->m, w);
918        res->e = 1;
919        res->f = smPolyWeight(res);
920        b = b->n;
921      }     
922      else
923      {
924        ha = smMult(a->m, p);
925        pDelete(&a->m);
926        hb = smMult(b->m, w);
927        ha = pAdd(ha, hb);
928        if (ha != NULL)
929        {
930          a->m = ha;
931          a->e = 1;
932          a->f = smPolyWeight(a);
933          res = res->n = a;
934          a = a->n;
935        }
936        else
937        {
938          smElemDelete(&a);
939        }
940        b = b->n;
941      }
942      if (b == NULL)
943      {
944        res->n = a;
945        break;
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    }
960    m_act[r->pos] = dumm->n;
961    smElemDelete(&r);
962  } while (r != NULL);
963}
964
965/* higher steps of elimination */
966void sparse_mat::smHElim()
967{
968  poly hp = this->smMultPoly(piv);
969  poly gp = piv->m;       // pivotelement
970  smpoly c = m_act[act];  // pivotcolumn
971  smpoly r = red;         // row to reduce
972  smpoly res, a, b;
973  poly ha, hr, gc, x, y;
974  int e, ip, ir, ia, lev;
975
976  if ((c == NULL) || (r == NULL))
977  {
978    pDelete(&hp);
979    return;
980  }
981  e = crd+1;
982  ip = piv->e;
983  do
984  {
985    a = m_act[r->pos];
986    res = dumm;
987    res->n = NULL;
988    b = c;
989    hr = r->m;
990    ir = r->e;
991    loop                  // combine the chains a and b: (hp,gp)*a(l) + hr*b(h)
992    {
993      if (a->pos < b->pos)
994      {
995        res = res->n = a;
996        a = a->n;
997      }
998      else if (a->pos > b->pos)
999      {
1000        res = res->n = smElemCopy(b);
1001        x = SM_MULT(b->m, hr, m_res[ir]->m);
1002        b = b->n;
1003        if(ir) SM_DIV(x, m_res[ir]->m);
1004        res->m = x;
1005        res->e = e;
1006        res->f = smPolyWeight(res);
1007      }     
1008      else
1009      {
1010        y = NULL;
1011        x = hr;
1012        ha = a->m;
1013        ia = a->e;
1014        if (ir > ia)
1015        {
1016          gc = SM_MULT(ha, m_res[ir]->m, m_res[ia]->m);
1017          pDelete(&ha);
1018          ha = gc;
1019          if (ia) SM_DIV(ha, m_res[ia]->m);
1020          ia = ir;
1021        }
1022        else if (ir < ia)
1023        {
1024          if (ip < ia)
1025          {
1026            gc = SM_MULT(ha, m_res[crd]->m, m_res[ia]->m);
1027            pDelete(&ha);
1028            ha = gc;
1029            if (ia) SM_DIV(ha, m_res[ia]->m);
1030            y = hp;
1031            ia = ip;
1032            if (ir > ia)
1033            {
1034              y = SM_MULT(y, m_res[ir]->m, m_res[ia]->m);
1035              if (ia) SM_DIV(y, m_res[ia]->m);
1036              ia = ir;
1037            }
1038            else if (ir < ia)
1039            {
1040              x = SM_MULT(x, m_res[ia]->m, m_res[ir]->m);
1041              if (ir) SM_DIV(x, m_res[ir]->m);
1042            }
1043          }
1044          else
1045          {
1046            x = SM_MULT(x, m_res[ia]->m, m_res[ir]->m);
1047            if (ir) SM_DIV(x, m_res[ir]->m);
1048          }
1049        }
1050        if (y == NULL) y = gp;
1051        gc = SM_MULT(ha, y, m_res[ia]->m);
1052        pDelete(&ha);
1053        x = SM_MULT(x, b->m, m_res[ia]->m);
1054        x = pAdd(x, gc);
1055        if (x != NULL)
1056        {
1057          if (ia) SM_DIV(x, m_res[ia]->m);
1058          a->m = x;
1059          a->e = e;
1060          a->f = smPolyWeight(a);
1061          res = res->n = a;
1062          a = a->n;
1063        }
1064        else
1065        {
1066          a->m = NULL;
1067          smElemDelete(&a);
1068        }
1069        b = b->n;
1070      }
1071      if (b == NULL)
1072      {
1073        res->n = a;
1074        break;
1075      }
1076      if (a == NULL)
1077      {
1078        do
1079        {
1080          res = res->n = smElemCopy(b);
1081          x = SM_MULT(b->m, hr, m_res[ir]->m);
1082          b = b->n;
1083          if(ir) SM_DIV(x, m_res[ir]->m);
1084          res->m = x;
1085          res->e = e;
1086          res->f = smPolyWeight(res);
1087        } while (b != NULL);
1088        break;
1089      }
1090    }
1091    m_act[r->pos] = dumm->n;
1092    smElemDelete(&r);
1093  } while (r != NULL);
1094  pDelete(&hp);
1095}
1096
1097/* ----------------- transfer ------------------ */
1098
1099/*
1100* select the pivotrow and store it to red and piv
1101*/
1102void sparse_mat::smSelectPR()
1103{
1104  smpoly b = dumm;
1105  smpoly a, ap;
1106  int i;
1107
1108  a = m_act[act];
1109  if (a->pos < rpiv)
1110  {
1111    do
1112    {
1113      ap = a;
1114      a = a->n;
1115    } while (a->pos < rpiv);
1116    ap->n = a->n;
1117  }
1118  else
1119    m_act[act] = a->n;
1120  piv = a;
1121  a->n = NULL;
1122  for (i=1; i<act; i++)
1123  {
1124    a = m_act[i];
1125    if (a->pos < rpiv)
1126    {
1127      loop
1128      {
1129        ap = a;
1130        a = a->n;
1131        if ((a == NULL) || (a->pos > rpiv))
1132          break;
1133        if (a->pos == rpiv)
1134        {
1135          ap->n = a->n;
1136          a->m = pNeg(a->m);
1137          b = b->n = a;
1138          b->pos = i;
1139          break;
1140        }
1141      }
1142    }
1143    else if (a->pos == rpiv)
1144    {
1145      m_act[i] = a->n;
1146      a->m = pNeg(a->m);
1147      b = b->n = a;
1148      b->pos = i;
1149    }
1150  }
1151  b->n = NULL;
1152  red = dumm->n;
1153}
1154
1155/*
1156* store the pivotcol in m_row
1157*   m_act[cols][pos(rows)] => m_row[rows][pos(cols)]
1158*/
1159void sparse_mat::smColToRow()
1160{
1161  smpoly c = m_act[act];
1162  smpoly h;
1163
1164  while (c != NULL)
1165  {
1166    h = c;
1167    c = c->n;
1168    h->n = m_row[h->pos];
1169    m_row[h->pos] = h;
1170    h->pos = crd;
1171  }
1172}
1173
1174/*
1175* store the pivot and the assosiated row in m_row
1176* to m_res (result):
1177*   piv + m_row[rows][pos(cols)] => m_res[cols][pos(rows)]
1178*/
1179void sparse_mat::smRowToCol()
1180{
1181  smpoly r = m_row[rpiv];
1182  smpoly a, ap, h;
1183
1184  m_row[rpiv] = NULL;
1185  perm[crd] = rpiv;
1186  piv->pos = crd;
1187  m_res[crd] = piv;
1188  while (r != NULL)
1189  {
1190    ap = m_res[r->pos];
1191    loop
1192    {
1193      a = ap->n;
1194      if (a == NULL)
1195      {
1196        ap->n = h = r;
1197        r = r->n;
1198        h->n = a;
1199        h->pos = crd;
1200        break;
1201      }
1202      ap = a;
1203    }
1204  }
1205}
1206
1207/* ----------------- C++ stuff ------------------ */
1208
1209/*
1210*  clean m_act from zeros (after elim)
1211*/
1212void sparse_mat::smZeroElim()
1213{
1214  int i = 0;
1215  int j;
1216
1217  loop
1218  {
1219    i++;
1220    if (i > act) return;
1221    if (m_act[i] == NULL) break;
1222  }
1223  j = i;
1224  loop
1225  {
1226    j++;
1227    if (j > act) break;
1228    if (m_act[j] != NULL)
1229    {
1230      m_act[i] = m_act[j];
1231      i++;
1232    }
1233  }
1234  act -= (j-i);
1235  sign = 0;
1236}
1237
1238/*
1239*  clean m_act from cols not to reduced (after elim)
1240*  put them to m_res
1241*/
1242void sparse_mat::smToredElim()
1243{
1244  int i = 0;
1245  int j;
1246
1247  loop
1248  {
1249    i++;
1250    if (i > act) return;
1251    if (m_act[i]->pos > tored)
1252    {
1253      m_res[crd] = m_act[i];
1254      crd++;
1255      break;
1256    }
1257  }
1258  j = i;
1259  loop
1260  {
1261    j++;
1262    if (j > act) break;
1263    if (m_act[i]->pos > tored)
1264    {
1265      m_res[crd] = m_act[i];
1266      crd++;
1267    }
1268    else
1269    {
1270      m_act[i] = m_act[j];
1271      i++;
1272    }
1273  }
1274  act -= (j-i);
1275  sign = 0;
1276}
1277
1278/*
1279*  copy m_act to m_res
1280*/
1281void sparse_mat::smCopToRes()
1282{
1283  smpoly p, r, ap, a, h;
1284  int c, i;
1285
1286  if (act == 0) return;
1287  if (act > 1)
1288  {
1289    Werror("Not implemented");
1290  }
1291  crd++;
1292  c = crd;
1293  p = m_res[c] = m_act[act];
1294  do
1295  {
1296    r = m_row[p->pos];
1297    m_row[p->pos] = NULL;
1298    perm[c] = p->pos;
1299    p->pos = c;
1300    while (r != NULL)
1301    {
1302      ap = m_res[r->pos];
1303      loop
1304      {
1305        a = ap->n;
1306        if (a == NULL)
1307        {
1308          ap->n = h = r;
1309          r = r->n;
1310          h->n = a;
1311          h->pos = c;
1312          break;
1313        }
1314        ap = a;
1315      }
1316    }
1317    c++;
1318    p = p->n;
1319  } while (p != NULL);
1320  for (i=1;i<=nrows;i++)
1321  {
1322    if(m_row[i] != NULL)
1323    {
1324      r = m_row[i];
1325      m_row[i] = NULL;
1326      do
1327      {
1328        ap = m_res[r->pos];
1329        loop
1330        {
1331          a = ap->n;
1332          if (a == NULL)
1333          {
1334            ap->n = h = r;
1335            r = r->n;
1336            h->n = a;
1337            h->pos = c;
1338            break;
1339          }
1340          ap = a;
1341        }
1342      } while (r!=NULL);
1343      c++;
1344    }
1345  }
1346  tored = c-1;
1347}
1348
1349/*
1350* multiply and divide the column, that goes to result
1351*/
1352void sparse_mat::smMultCol()
1353{
1354  smpoly a = m_act[act];
1355  int e = crd;
1356  poly ha;
1357  int f;
1358
1359  while (a != NULL)
1360  {
1361    f = a->e;
1362    if (f < e)
1363    {
1364      ha = SM_MULT(a->m, m_res[e]->m, m_res[f]->m);
1365      pDelete(&a->m);
1366      if (f) SM_DIV(ha, m_res[f]->m);
1367      a->m = ha;
1368    }
1369    a = a->n;
1370  }
1371}
1372
1373/*
1374* multiply and divide the m_act finaly
1375*/
1376void sparse_mat::smFinalMult()
1377{
1378  smpoly a = m_act[act];
1379  int e = crd;
1380  poly ha;
1381  int i, f;
1382
1383  for (i=act; i; i--)
1384  {
1385    a = m_act[i];
1386    do
1387    {
1388      f = a->e;
1389      if (f < e)
1390      {
1391        ha = SM_MULT(a->m, m_res[e]->m, m_res[f]->m);
1392        pDelete(&a->m);
1393        if (f) SM_DIV(ha, m_res[f]->m);
1394        a->m = ha;
1395      }
1396      a = a->n;
1397    } while (a != NULL);
1398  }
1399}
1400
1401/*
1402* multiply and divide the element, save poly
1403*/
1404poly sparse_mat::smMultPoly(smpoly a)
1405{
1406  int f = a->e;
1407  poly r, h;
1408
1409  if (f < crd)
1410  {
1411    h = r = a->m;
1412    h = SM_MULT(h, m_res[crd]->m, m_res[f]->m);
1413    if (f) SM_DIV(h, m_res[f]->m);
1414    a->m = h;
1415    a->f = smPolyWeight(a);
1416    return r;
1417  }
1418  else
1419    return NULL;
1420}
1421
1422/*
1423* delete the m_act finaly
1424*/
1425void sparse_mat::smActDel()
1426{
1427  smpoly a;
1428  int i;
1429
1430  for (i=act; i; i--)
1431  {
1432    a = m_act[i];
1433    do
1434    {
1435      smElemDelete(&a);
1436    } while (a != NULL);
1437  }
1438}
1439
1440/*
1441* delete the pivotcol
1442*/
1443void sparse_mat::smColDel()
1444{
1445  smpoly a = m_act[act];
1446
1447  while (a != NULL)
1448  {
1449    smElemDelete(&a);
1450  }
1451}
1452
1453/*
1454* delete pivot elements
1455*/
1456void sparse_mat::smPivDel()
1457{
1458  int i=crd-1;
1459
1460  while (i != 0)
1461  {
1462    smElemDelete(&m_res[i]);
1463    i--;
1464  }
1465}
1466
1467/*
1468* the sign of the determinant
1469*/
1470void sparse_mat::smSign()
1471{
1472  int j,i=1;
1473  if (cpiv!=act) sign=-sign;
1474  if ((act%2)==0) sign=-sign;
1475  j=perm[1];
1476  while(j<rpiv)
1477  {
1478    sign=-sign;
1479    i++;
1480    j=perm[i];
1481  }
1482  while(perm[i]!=0)
1483  {
1484    perm[i]=perm[i+1];
1485    i++;
1486  }
1487}
1488
1489/* ----------------- arithmetic ------------------ */
1490
1491/*
1492*  returns a*b
1493*  a,b NOT destroyed
1494*/
1495poly smMult(poly a, poly b)
1496{
1497  poly pa, res, r;
1498
1499  if (smSmaller(a, b))
1500  {
1501    r = a;
1502    a = b;
1503    b = r;
1504  }
1505  if (pNext(b) == NULL)
1506  {
1507    if (smIsScalar(b))
1508      return pMultCopyN(a, pGetCoeff(b));
1509    else
1510      return smEMult(a, b);
1511  }
1512  pa = res = smEMult(a, b);
1513  pIter(b);
1514  do
1515  {
1516    r = smEMult(a, b);
1517    smCombineChain(&pa, r);
1518    pIter(b);
1519  } while (b != NULL);
1520  return res;
1521}
1522
1523/*2
1524* exact division a/b
1525* a destroyed, b NOT destroyed
1526*/
1527void smPolyDiv(poly a, poly b)
1528{
1529  const number x = pGetCoeff(b);
1530  number y, yn;
1531  poly t, h, dummy;
1532  int i;
1533
1534  if (pNext(b) == NULL)
1535  {
1536    do
1537    {
1538      if (!smIsScalar(b))
1539      {
1540        for (i=pVariables; i; i--)
1541          pSubExp(a,i,pGetExp(b,i));
1542        pSetm(a);
1543      }
1544      y = nDiv(pGetCoeff(a),x);
1545      nNormalize(y);
1546      pSetCoeff(a,y);
1547      pIter(a);
1548    } while (a != NULL);
1549    return;
1550  }
1551  dummy = pInit();
1552  do
1553  {
1554    for (i=pVariables; i; i--)
1555      pSubExp(a,i,pGetExp(b,i));
1556    pSetm(a);
1557    y = nDiv(pGetCoeff(a),x);
1558    nNormalize(y);
1559    pSetCoeff(a,y);
1560    yn = nNeg(nCopy(y));
1561    t = pNext(b);
1562    h = dummy;
1563    do
1564    {
1565      h = pNext(h) = pInit();
1566      //pSetComp(h,0);
1567      for (i=pVariables; i; i--)
1568        pSetExp(h,i,pGetExp(a,i)+pGetExp(t,i));
1569      pSetm(h);
1570      pSetCoeff0(h,nMult(yn, pGetCoeff(t)));
1571      pIter(t);
1572    } while (t != NULL);
1573    nDelete(&yn);
1574    pNext(h) = NULL;
1575    a = pNext(a) = pAdd(pNext(a), pNext(dummy));
1576  } while (a!=NULL);
1577  pFree1(dummy);
1578}
1579
1580/*
1581*  returns the part of (a*b)/exp(lead(c)) with nonegative exponents
1582*/
1583poly smMultDiv(poly a, poly b, const poly c)
1584{
1585  poly pa, e, res, r;
1586  BOOLEAN lead;
1587
1588  if (smSmaller(a, b))
1589  {
1590    r = a;
1591    a = b;
1592    b = r;
1593  }
1594  if ((c == NULL) || smIsScalar(c))
1595  {
1596    if (pNext(b) == NULL)
1597    {
1598      if (smIsScalar(b))
1599        return pMultCopyN(a, pGetCoeff(b));
1600      else
1601        return smEMult(a, b);
1602    }
1603    pa = res = smEMult(a, b);
1604    pIter(b);
1605    do
1606    {
1607      r = smEMult(a, b);
1608      smCombineChain(&pa, r);
1609      pIter(b);
1610    } while (b != NULL);
1611    return res;
1612  }
1613  res = NULL;
1614  e = pInit();
1615  lead = FALSE;
1616  while (!lead)
1617  {
1618    pSetCoeff0(e,pGetCoeff(b));
1619    if (smIsNegQuot(e, b, c))
1620    {
1621      lead = smCheckLead(a, e);
1622      r = smDMult(a, e);
1623    }
1624    else
1625    {
1626      lead = TRUE;
1627      r = smEMult(a, e);
1628    }
1629    if (lead)
1630    {
1631      if (res != NULL)
1632      {
1633        smFindRef(&pa, &res, r);
1634        if (pa == NULL)
1635          lead = FALSE;
1636      }
1637      else
1638      {
1639        pa = res = r;
1640      }
1641    }
1642    else
1643      res = pAdd(res, r);
1644    pIter(b);
1645    if (b == NULL)
1646    {
1647      pFree1(e);
1648      return res;
1649    }
1650  }
1651  do
1652  {
1653    pSetCoeff0(e,pGetCoeff(b));
1654    if (smIsNegQuot(e, b, c))
1655    {
1656      r = smDMult(a, e);
1657      if (smCheckLead(a, e))
1658        smCombineChain(&pa, r);
1659      else
1660        pa = pAdd(pa,r); 
1661    }
1662    else
1663    {
1664      r = smEMult(a, e);
1665      smCombineChain(&pa, r);
1666    }
1667    pIter(b);
1668  } while (b != NULL);
1669  pFree1(e);
1670  return res;
1671}
1672
1673/*n
1674* exact division a/b
1675* a is a result of smMultDiv
1676* a destroyed, b NOT destroyed
1677*/
1678void smSpecialPolyDiv(poly a, poly b)
1679{
1680  if (pNext(b) == NULL)
1681  {
1682    smPolyDivN(a, pGetCoeff(b));
1683    return;
1684  }
1685  smExactPolyDiv(a, b);
1686}
1687
1688/* ------------ internals arithmetic ------------- */
1689static void smExactPolyDiv(poly a, poly b)
1690{
1691  const number x = pGetCoeff(b);
1692  poly tail = pNext(b), e = pInit();
1693  poly h;
1694  number y, yn;
1695
1696  do
1697  {
1698    y = nDiv(pGetCoeff(a), x);
1699    nNormalize(y);
1700    pSetCoeff(a,y);
1701    yn = nNeg(nCopy(y));
1702    pSetCoeff0(e,yn);
1703    if (smIsNegQuot(e, a, b)) 
1704      h = smDMult(tail, e);
1705    else
1706      h = smEMult(tail, e);
1707    nDelete(&yn);
1708    a = pNext(a) = pAdd(pNext(a), h);
1709  } while (a!=NULL);
1710  pFree1(e);
1711}
1712
1713static BOOLEAN smIsScalar(const poly Term)
1714{
1715  int i;
1716
1717  for (i=pVariables; i; i--)
1718  {
1719    if (pGetExp(Term,i)) return FALSE;
1720  }
1721  return TRUE;
1722}
1723
1724static BOOLEAN smIsNegQuot(poly a, const poly b, const poly c)
1725{
1726  int i;
1727
1728  for (i=pVariables; i; i--)
1729  {
1730    pSetExp(a,i,pGetExp(b,i)-pGetExp(c,i));
1731    if (pGetExp(a,i) < 0)
1732    {
1733      while(--i) pSetExp(a,i,pGetExp(b,i)-pGetExp(c,i));
1734      return TRUE;
1735    }
1736  }
1737  return FALSE;
1738}
1739
1740static poly smEMult(poly t, const poly e)
1741{
1742  const number y = pGetCoeff(e);
1743  poly res, h;
1744  int i;
1745
1746  h = res = pInit();
1747  loop
1748  {
1749    //pSetComp(h,0);
1750    for (i=pVariables; i; i--)
1751      pSetExp(h,i,pGetExp(e,i)+pGetExp(t,i));
1752    pSetm(h);
1753    pSetCoeff0(h,nMult(y,pGetCoeff(t)));
1754    pIter(t);
1755    if (t == NULL)
1756    {
1757      pNext(h) = NULL;
1758      return res;
1759    }
1760    h = pNext(h) = pInit();
1761  }
1762}
1763
1764static BOOLEAN smCheckLead(const poly t, const poly e)
1765{
1766  int i;
1767  for (i=pVariables; i; i--)
1768  {
1769    if ((pGetExp(e,i)+pGetExp(t,i)) < 0)
1770      return FALSE;
1771  }
1772  return TRUE;
1773}
1774
1775static poly smDMult(poly t, const poly e)
1776{
1777  const number y = pGetCoeff(e);
1778  number x;
1779  poly res, h, r;
1780  int i;
1781  EXPONENT_TYPE w;
1782
1783  r = h = res = pInit();
1784  loop
1785  {
1786    x = pGetCoeff(t);
1787    for (i=pVariables; i; i--)
1788    {
1789      w = pGetExp(e,i)+pGetExp(t,i);
1790      if (w < 0) break;
1791      pSetExp(h,i,w);
1792    }
1793    pIter(t);
1794    if (i == 0)
1795    {
1796      //pSetComp(h,0);
1797      pSetm(h);
1798      pSetCoeff0(h,nMult(y,x));
1799      if (t == NULL)
1800      {
1801        pNext(h) = NULL;
1802        return res;
1803      }
1804      r = h;
1805      h = pNext(h) = pInit();
1806    }
1807    else if (t == NULL)
1808    {
1809      if (r != h)
1810        pNext(r) = NULL;
1811      else
1812        res = NULL;
1813      pFree1(h);
1814      return res;
1815    }
1816  }
1817}
1818
1819static void smPolyDivN(poly a, const number x)
1820{
1821  number y;
1822
1823  do
1824  {
1825    y = nDiv(pGetCoeff(a),x);
1826    nNormalize(y);
1827    pSetCoeff(a,y);
1828    pIter(a);
1829  } while (a != NULL);
1830}
1831
1832static BOOLEAN smSmaller(poly a, poly b)
1833{
1834  loop
1835  {
1836    pIter(b);
1837    if (b == NULL) return TRUE;
1838    pIter(a);
1839    if (a == NULL) return FALSE;
1840  }
1841}
1842
1843static void smCombineChain(poly *px, poly r)
1844{
1845  poly pa = *px, pb;
1846  number x;
1847  int i;
1848
1849  loop
1850  {
1851    pb = pNext(pa);
1852    if (pb == NULL)
1853    {
1854      pa = pNext(pa) = r;
1855      break;
1856    }
1857    i = pComp0(pb, r);
1858    if (i > 0)
1859      pa = pb;
1860    else
1861    {
1862      if (i == 0)
1863      {
1864        x = nAdd(pGetCoeff(pb), pGetCoeff(r));
1865        pDelete1(&r);
1866        if (nIsZero(x))
1867        {
1868          pDelete1(&pb);
1869          pNext(pa) = pAdd(pb,r);
1870        }
1871        else
1872        {
1873          pa = pb;
1874          pSetCoeff(pa,x);
1875          pNext(pa) = pAdd(pNext(pa), r);
1876        }
1877      }
1878      else
1879      {
1880        pa = pNext(pa) = r;
1881        pNext(pa) = pAdd(pb, pNext(pa));
1882      }
1883      break;
1884    }
1885  }
1886  *px = pa;
1887}
1888
1889static void smFindRef(poly *ref, poly *px, poly r)
1890{
1891  number x;
1892  int i;
1893  poly pa = *px, pp = NULL;
1894
1895  loop
1896  {
1897    i = pComp0(pa, r);
1898    if (i > 0)
1899    {
1900      pp = pa;
1901      pIter(pa);
1902      if (pa==NULL)
1903      {
1904        pNext(pp) = r;
1905        break;
1906      }
1907    }
1908    else
1909    {
1910      if (i == 0)
1911      {
1912        x = nAdd(pGetCoeff(pa), pGetCoeff(r));
1913        pDelete1(&r);
1914        if (nIsZero(x))
1915        {
1916          pDelete1(&pa);
1917          if (pp!=NULL)
1918            pNext(pp) = pAdd(pa,r);
1919          else
1920            *px = pAdd(pa,r);
1921        }
1922        else
1923        {
1924          pp = pa;
1925          pSetCoeff(pp,x);
1926          pNext(pp) = pAdd(pNext(pp), r);
1927        }
1928      }
1929      else
1930      {
1931        if (pp!=NULL)
1932          pp = pNext(pp) = r;
1933        else
1934          *px = pp = r;
1935        pNext(pp) = pAdd(pa, pNext(r));
1936      }
1937      break;
1938    }
1939  }
1940  *ref = pp;
1941}
1942
1943/* ----------------- internal 'C' stuff ------------------ */
1944
1945static void smElemDelete(smpoly *r)
1946{
1947  smpoly a = *r, b = a->n;
1948
1949  pDelete(&a->m);
1950  Free((void *)a, sizeof(smprec));
1951  *r = b;
1952}
1953
1954static smpoly smElemCopy(smpoly a)
1955{
1956  smpoly r = (smpoly)Alloc(sizeof(smprec));
1957
1958  memcpy(r, a, sizeof(smprec));
1959/*  r->m = pCopy(r->m); */
1960  return r;
1961}
1962
1963/*
1964* from poly to smpoly
1965* do not destroy p
1966*/
1967static smpoly smPoly2Smpoly(poly q)
1968{
1969  poly pp;
1970  smpoly res, a;
1971  EXPONENT_TYPE x;
1972
1973  if (q == NULL)
1974    return NULL;
1975  a = res = (smpoly)Alloc(sizeof(smprec));
1976  a->pos = x = pGetComp(q);
1977  a->m = q;
1978  a->e = 0;
1979  loop
1980  {
1981    pSetComp(q,0);
1982    pp = q;
1983    pIter(q);
1984    if (q == NULL)
1985    {
1986      a->n = NULL;
1987      return res;
1988    }
1989    if (pGetComp(q) != x)
1990    {
1991      a = a->n = (smpoly)Alloc(sizeof(smprec));
1992      pNext(pp) = NULL;
1993      a->pos = x = pGetComp(q);
1994      a->m = q;
1995      a->e = 0;
1996    }
1997  }
1998}
1999
2000/*
2001* from smpoly to poly
2002* destroy a
2003*/
2004static poly smSmpoly2Poly(smpoly a)
2005{
2006  smpoly b;
2007  poly res, pp, q;
2008  EXPONENT_TYPE x;
2009
2010  if (a == NULL)
2011    return NULL;
2012  x = a->pos;
2013  q = res = a->m;
2014  loop
2015  {
2016    pSetComp(q,x);
2017    pp = q;
2018    pIter(q);
2019    if (q == NULL)
2020      break;
2021  }
2022  loop
2023  {
2024    b = a;
2025    a = a->n;
2026    Free((void *)b, sizeof(smprec));
2027    if (a == NULL)
2028      return res;
2029    x = a->pos;
2030    q = pNext(pp) = a->m;
2031    loop
2032    {
2033      pSetComp(q,x);
2034      pp = q;
2035      pIter(q);
2036      if (q == NULL)
2037        break;
2038    }
2039  }
2040}
2041
2042/*
2043* weigth of a polynomial, for pivot strategy
2044*/
2045static float smPolyWeight(smpoly a)
2046{
2047  poly p = a->m;
2048  int i;
2049  float res = (float)nSize(pGetCoeff(p));
2050
2051  if (pNext(p) == NULL)
2052  {
2053    for(i=pVariables; i>0; i--)
2054    {
2055      if (pGetExp(p,i) != 0) return res+1.0;
2056    }
2057    return res;
2058  }
2059  else
2060  {
2061    i = 0;
2062    res = 0.0;
2063    do
2064    {
2065      i++;
2066      res += (float)nSize(pGetCoeff(p));
2067      pIter(p);
2068    }
2069    while (p);
2070    return res+(float)i;
2071  }
2072}
Note: See TracBrowser for help on using the repository browser.