source: git/libpolys/polys/sparsmat.cc @ b78996

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