source: git/kernel/sparsmat.cc @ ca371d

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