source: git/Singular/sparsmat.cc @ e03e67

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