source: git/libpolys/polys/sparsmat.cc @ 9ccaaf

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