source: git/libpolys/misc/intvec.cc @ 6ce030f

fieker-DuValspielwiese
Last change on this file since 6ce030f was 6ce030f, checked in by Oleksandr Motsak <motsak@…>, 12 years ago
removal of the $Id$ svn tag from everywhere NOTE: the git SHA1 may be used instead (only on special places) NOTE: the libraries Singular/LIB/*.lib still contain the marker due to our current use of svn
  • Property mode set to 100644
File size: 14.5 KB
Line 
1/*****************************************
2*  Computer Algebra System SINGULAR      *
3*****************************************/
4/*
5* ABSTRACT: class intvec: lists/vectors of integers
6*/
7#ifndef INTVEC_CC
8#define INTVEC_CC
9
10#include "config.h"
11#include <misc/auxiliary.h>
12
13// #include <findexec/feFopen.h>
14#include <misc/intvec.h>
15#include <omalloc/omalloc.h>
16
17/*0 implementation*/
18
19// omBin intvec_bin = omGetSpecBin(sizeof(intvec));
20#if 0
21intvec::intvec(intvec* iv)
22{
23  row = iv->rows();
24  col = iv->cols();
25  v   = (int *)omAlloc(sizeof(int)*row*col);
26  for (int i=0; i<row*col; i++)
27  {
28    v[i] = (*iv)[i];
29  }
30}
31#endif
32
33intvec::intvec(int s, int e)
34{
35  int inc;
36  col = 1;
37  if (s<e)
38  {
39    row =  e-s+1;
40    inc =  1;
41  }
42  else
43  {
44    row = s-e+1;
45    inc = -1;
46  }
47  v = (int *)omAlloc(sizeof(int)*row);
48  for (int i=0; i<row; i++)
49  {
50    v[i] = s;
51    s+=inc;
52  }
53}
54
55intvec::intvec(int r, int c, int init)
56{
57  row = r;
58  col = c;
59  int l = r*c;
60  if (l>0) /*(r>0) && (c>0) */
61    v = (int *)omAlloc(sizeof(int)*l);
62  else
63    v = NULL;
64  for (int i=0; i<l; i++)
65  {
66    v[i] = init;
67  }
68}
69
70char * intvec::ivString(int not_mat,int spaces, int dim) const
71{
72  //Print("ivString:this=%x,v=%x,row=%d\n",this,v,row);
73#ifndef OM_NDEBUG
74  omCheckAddr((void *)this);
75  if (v!=NULL) omCheckAddr((void *)v);
76#endif
77  StringSetS("");
78  if ((col == 1)&&(not_mat))
79  {
80    int i=0;
81    for (; i<row-1; i++)
82    {
83      StringAppend("%d,",v[i]);
84    }
85    if (i<row)
86    {
87      StringAppend("%d",v[i]);
88    }
89  }
90  else
91  {
92    for (int j=0; j<row; j++)
93    {
94      if (j<row-1)
95      {
96        for (int i=0; i<col; i++)
97        {
98          StringAppend("%d%c",v[j*col+i],',');
99        }
100      }
101      else
102      {
103        for (int i=0; i<col; i++)
104        {
105          StringAppend("%d%c",v[j*col+i],i<col-1 ? ',' : ' ');
106        }
107      }
108      if (j+1<row)
109      {
110        if (dim > 1) StringAppendS("\n");
111        if (spaces>0) StringAppend("%-*.*s",spaces,spaces," ");
112      }
113    }
114  }
115  return StringAppendS("");
116}
117
118void intvec::resize(int new_length)
119{
120  assume(new_length > 0 && col == 1);
121  v = (int*) omRealloc0Size(v, row*sizeof(int), new_length*sizeof(int));
122  row = new_length;
123}
124
125char * intvec::String(int dim) const
126{
127  return omStrDup(ivString(1, 0, dim));
128}
129
130#ifndef NDEBUG
131// debug only
132void intvec::view () const
133{
134  Print ("intvec: {rows: %d, cols: %d, length: %d, Values: \n", rows(), cols(), length());
135
136  for (int i = 0; i < rows(); i++)
137  {
138    Print ("Row[%3d]:", i);
139    for (int j = 0; j < cols(); j++)
140      Print (" %5d", this->operator[]((i)*cols()+j) );
141    PrintLn ();
142  }
143  PrintS ("}\n");
144}
145#endif
146
147void intvec::show(int notmat,int spaces) const
148{
149  if (spaces>0)
150  {
151    PrintNSpaces(spaces);
152    PrintS(ivString(notmat,spaces));
153  }
154  else
155    PrintS(ivString(notmat,0));
156}
157
158void intvec::operator+=(int intop)
159{
160  for (int i=0; i<row*col; i++) { v[i] += intop; }
161}
162
163void intvec::operator-=(int intop)
164{
165  for (int i=0; i<row*col; i++) { v[i] -= intop; }
166}
167
168void intvec::operator*=(int intop)
169{
170  for (int i=0; i<row*col; i++) { v[i] *= intop; }
171}
172
173void intvec::operator/=(int intop)
174{
175  if (intop == 0) return;
176  int bb=ABS(intop);
177  for (int i=0; i<row*col; i++)
178  {
179    int r=v[i];
180    int c=r%bb;
181    if (c<0) c+=bb;
182    r=(r-c)/intop;
183    v[i]=r;
184  }
185}
186
187void intvec::operator%=(int intop)
188{
189  if (intop == 0) return;
190  int bb=ABS(intop);
191  for (int i=0; i<row*col; i++)
192  {
193    int r=v[i];
194    int c=r%bb;
195    if (c<0) c+=bb;
196    v[i]=c;
197  }
198}
199
200int intvec::compare(const intvec* op) const
201{
202  if ((col!=1) ||(op->cols()!=1))
203  {
204    if((col!=op->cols())
205    || (row!=op->rows()))
206      return -2;
207  }
208  int i;
209  for (i=0; i<si_min(length(),op->length()); i++)
210  {
211    if (v[i] > (*op)[i])
212      return 1;
213    if (v[i] < (*op)[i])
214      return -1;
215  }
216  // this can only happen for intvec: (i.e. col==1)
217  for (; i<row; i++)
218  {
219    if (v[i] > 0)
220      return 1;
221    if (v[i] < 0)
222      return -1;
223  }
224  for (; i<op->rows(); i++)
225  {
226    if (0 > (*op)[i])
227      return 1;
228    if (0 < (*op)[i])
229      return -1;
230  }
231  return 0;
232}
233int intvec::compare(int o) const
234{
235  for (int i=0; i<row*col; i++)
236  {
237    if (v[i] <o) return -1;
238    if (v[i] >o) return 1;
239  }
240  return 0;
241}
242
243#if 0
244intvec * ivCopy(intvec * o)
245{
246  intvec * iv=new intvec(o);
247  return iv;
248}
249#endif
250
251intvec * ivAdd(intvec * a, intvec * b)
252{
253  intvec * iv;
254  int mn, ma, i;
255  if (a->cols() != b->cols()) return NULL;
256  mn = si_min(a->rows(),b->rows());
257  ma = si_max(a->rows(),b->rows());
258  if (a->cols() == 1)
259  {
260    iv = new intvec(ma);
261    for (i=0; i<mn; i++) (*iv)[i] = (*a)[i] + (*b)[i];
262    if (ma > mn)
263    {
264      if (ma == a->rows())
265      {
266        for(i=mn; i<ma; i++) (*iv)[i] = (*a)[i];
267      }
268      else
269      {
270        for(i=mn; i<ma; i++) (*iv)[i] = (*b)[i];
271      }
272    }
273    return iv;
274  }
275  if (mn != ma) return NULL;
276  iv = new intvec(a);
277  for (i=0; i<mn*a->cols(); i++) { (*iv)[i] += (*b)[i]; }
278  return iv;
279}
280
281intvec * ivSub(intvec * a, intvec * b)
282{
283  intvec * iv;
284  int mn, ma, i;
285  if (a->cols() != b->cols()) return NULL;
286  mn = si_min(a->rows(),b->rows());
287  ma = si_max(a->rows(),b->rows());
288  if (a->cols() == 1)
289  {
290    iv = new intvec(ma);
291    for (i=0; i<mn; i++) (*iv)[i] = (*a)[i] - (*b)[i];
292    if (ma > mn)
293    {
294      if (ma == a->rows())
295      {
296        for(i=mn; i<ma; i++) (*iv)[i] = (*a)[i];
297      }
298      else
299      {
300        for(i=mn; i<ma; i++) (*iv)[i] = -(*b)[i];
301      }
302    }
303    return iv;
304  }
305  if (mn != ma) return NULL;
306  iv = new intvec(a);
307  for (i=0; i<mn*a->cols(); i++) { (*iv)[i] -= (*b)[i]; }
308  return iv;
309}
310
311intvec * ivTranp(intvec * o)
312{
313  int i, j, r = o->rows(), c = o->cols();
314  intvec * iv= new intvec(c, r, 0);
315  for (i=0; i<r; i++)
316  {
317    for (j=0; j<c; j++)
318      (*iv)[j*r+i] = (*o)[i*c+j];
319  }
320  return iv;
321}
322
323int ivTrace(intvec * o)
324{
325  int i, s = 0, m = si_min(o->rows(),o->cols()), c = o->cols();
326  for (i=0; i<m; i++)
327  {
328    s += (*o)[i*c+i];
329  }
330  return s;
331}
332
333intvec * ivMult(intvec * a, intvec * b)
334{
335  int i, j, k, sum,
336      ra = a->rows(), ca = a->cols(),
337      rb = b->rows(), cb = b->cols();
338  intvec * iv;
339  if (ca != rb) return NULL;
340  iv = new intvec(ra, cb, 0);
341  for (i=0; i<ra; i++)
342  {
343    for (j=0; j<cb; j++)
344    {
345      sum = 0;
346      for (k=0; k<ca; k++)
347        sum += (*a)[i*ca+k]*(*b)[k*cb+j];
348      (*iv)[i*cb+j] = sum;
349    }
350  }
351  return iv;
352}
353
354/*2
355*computes a triangular matrix
356*/
357//void ivTriangMat(intvec * imat)
358//{
359//  int i=0,j=imat->rows(),k=j*imat->cols()-1;
360//
361//  ivTriangIntern(imat,i,j);
362//  i *= imat->cols();
363//  for(j=k;j>=i;j--)
364//    (*imat)[j] = 0;
365//}
366
367/* def. internals */
368static int ivColPivot(intvec *, int, int, int, int);
369static void ivNegRow(intvec *, int);
370static void ivSaveRow(intvec *, int);
371static void ivSetRow(intvec *, int, int);
372static void ivFreeRow(intvec *, int, int);
373static void ivReduce(intvec *, int, int, int, int);
374static void ivZeroElim(intvec *,int, int, int &);
375static void ivRowContent(intvec *, int, int);
376static void ivKernFromRow(intvec *, intvec *, intvec *,
377                          int, int, int);
378static intvec * ivOptimizeKern(intvec *);
379static int ivGcd(int, int);
380static void ivOptRecursive(intvec *, intvec *, intvec *,
381                          int &, int &, int);
382static void ivOptSolve(intvec *, intvec *, int &, int &);
383static void ivContent(intvec *);
384static int ivL1Norm(intvec *);
385static int ivCondNumber(intvec *, int); 
386
387/* Triangulierung in intmat.cc */
388void ivTriangIntern(intvec *imat, int &ready, int &all)
389{
390  int rpiv, colpos=0, rowpos=0;
391  int ia=ready, ie=all;
392
393  do
394  {
395    rowpos++;
396    do
397    {
398      colpos++;
399      rpiv = ivColPivot(imat, colpos, rowpos, ia, ie);
400    } while (rpiv==0);
401    if (rpiv>ia)
402    {
403      if (rowpos!=rpiv)
404      {
405        ivSaveRow(imat, rpiv);
406        ivFreeRow(imat, rowpos, rpiv);
407        ivSetRow(imat, rowpos, colpos);
408        rpiv = rowpos;
409      }
410      ia++;
411      if (ia==imat->cols())
412      {
413        ready = ia;
414        all = ie;
415        return;
416      }
417    }
418    ivReduce(imat, rpiv, colpos, ia, ie);
419    ivZeroElim(imat, colpos, ia, ie);
420  } while (ie>ia);
421  ready = ia;
422  all = ie;
423}
424
425/* Kernberechnung in intmat.cc */
426intvec * ivSolveKern(intvec *imat, int dimtr)
427{
428  int d=imat->cols();
429  int kdim=d-dimtr;
430  intvec *perm = new intvec(dimtr+1);
431  intvec *kern = new intvec(kdim,d,0);
432  intvec *res;
433  int c, cp, r, t;
434
435  t = kdim;
436  c = 1;
437  for (r=1;r<=dimtr;r++)
438  {
439    while (IMATELEM(*imat,r,c)==0) c++;
440    (*perm)[r] = c;
441    c++;
442  }
443  c = d;
444  for (r=dimtr;r>0;r--)
445  {
446    cp = (*perm)[r];
447    if (cp!=c)
448    {
449      ivKernFromRow(kern, imat, perm, t, r, c);
450      t -= (c-cp);
451      if (t==0)
452        break;
453      c = cp-1;
454    }
455    else
456      c--;
457  }
458  if (kdim>1)
459    res = ivOptimizeKern(kern);
460  else
461    res = ivTranp(kern);
462  delete kern;
463  delete perm;
464  return res;
465}
466
467/* internals */
468static int ivColPivot(intvec *imat, int colpos, int rowpos, int ready, int all)
469{
470  int rpiv;
471
472  if (IMATELEM(*imat,rowpos,colpos)!=0)
473    return rowpos;
474  for (rpiv=ready+1;rpiv<=all;rpiv++)
475  {
476    if (IMATELEM(*imat,rpiv,colpos)!=0)
477      return rpiv;
478  }
479  return 0;
480}
481
482static void ivNegRow(intvec *imat, int rpiv)
483{
484  int i;
485  for (i=imat->cols();i!=0;i--)
486    IMATELEM(*imat,rpiv,i) = -(IMATELEM(*imat,rpiv,i));
487}
488
489static void ivSaveRow(intvec *imat, int rpiv)
490{
491  int i, j=imat->rows();
492
493  for (i=imat->cols();i!=0;i--)
494    IMATELEM(*imat,j,i) = IMATELEM(*imat,rpiv,i);
495}
496
497static void ivSetRow(intvec *imat, int rowpos, int colpos)
498{
499  int i, j=imat->rows();
500
501  for (i=imat->cols();i!=0;i--)
502    IMATELEM(*imat,rowpos,i) = IMATELEM(*imat,j,i);
503  ivRowContent(imat, rowpos, colpos);
504}
505
506static void ivFreeRow(intvec *imat, int rowpos, int rpiv)
507{
508  int i, j;
509
510  for (j=rpiv-1;j>=rowpos;j--)
511  {
512    for (i=imat->cols();i!=0;i--)
513      IMATELEM(*imat,j+1,i) = IMATELEM(*imat,j,i);
514  }
515}
516
517static void ivReduce(intvec *imat, int rpiv, int colpos,
518                     int ready, int all)
519{
520  int tgcd, ce, m1, m2, j, i;
521  int piv = IMATELEM(*imat,rpiv,colpos);
522
523  for (j=all;j>ready;j--)
524  {
525    ce = IMATELEM(*imat,j,colpos);
526    if (ce!=0)
527    {
528      IMATELEM(*imat,j,colpos) = 0;
529      m1 = piv;
530      m2 = ce;
531      tgcd = ivGcd(m1, m2);
532      if (tgcd != 1)
533      {
534        m1 /= tgcd;
535        m2 /= tgcd;
536      }
537      for (i=imat->cols();i>colpos;i--)
538      {
539        IMATELEM(*imat,j,i) = IMATELEM(*imat,j,i)*m1-
540                              IMATELEM(*imat,rpiv,i)*m2;
541      }
542      ivRowContent(imat, j, colpos+1);
543    }
544  }
545}
546
547static void ivZeroElim(intvec *imat, int colpos,
548                     int ready, int &all)
549{
550  int j, i, k, l;
551
552  k = ready;
553  for (j=ready+1;j<=all;j++)
554  {
555    for (i=imat->cols();i>colpos;i--)
556    {
557      if (IMATELEM(*imat,j,i)!=0)
558      {
559        k++;
560        if (k<j)
561        {
562          for (l=imat->cols();l>colpos;l--)
563            IMATELEM(*imat,k,l) = IMATELEM(*imat,j,l);
564        }
565        break;
566      }
567    }
568  }
569  all = k;
570}
571
572static void ivRowContent(intvec *imat, int rowpos, int colpos)
573{
574  int tgcd, m;
575  int i=imat->cols();
576
577  loop
578  {
579    tgcd = IMATELEM(*imat,rowpos,i--);
580    if (tgcd!=0) break;
581    if (i<colpos) return;
582  }
583  if (tgcd<0) tgcd = -tgcd;
584  if (tgcd==1) return;
585  loop
586  {
587    m = IMATELEM(*imat,rowpos,i--);
588    if (m!=0) tgcd= ivGcd(tgcd, m);
589    if (tgcd==1) return;
590    if (i<colpos) break;
591  }
592  for (i=imat->cols();i>=colpos;i--)
593    IMATELEM(*imat,rowpos,i) /= tgcd;
594}
595
596static void ivKernFromRow(intvec *kern, intvec *imat,
597                          intvec *perm, int pos, int r, int c)
598{
599  int piv, cp, g, i, j, k, s;
600
601  for (i=c;i>(*perm)[r];i--)
602  {
603    IMATELEM(*kern,pos,i) = 1;
604    for (j=r;j!=0;j--)
605    {
606      cp = (*perm)[j];
607      s=0;
608      for(k=c;k>cp;k--)
609        s += IMATELEM(*imat,j,k)*IMATELEM(*kern,pos,k);
610      if (s!=0)
611      {
612        piv = IMATELEM(*imat,j,cp);
613        g = ivGcd(piv,s);
614        if (g!=1)
615        {
616          s /= g;
617          piv /= g;
618        }
619        for(k=c;k>cp;k--)
620          IMATELEM(*kern,pos,k) *= piv;
621        IMATELEM(*kern,pos,cp) = -s;
622        ivRowContent(kern,pos,cp);
623      }
624    }
625    if (IMATELEM(*kern,pos,i)<0)
626      ivNegRow(kern,pos);
627    pos--;
628  }
629}
630
631static int ivGcd(int a,int b)
632{
633  int x;
634
635  if (a<0) a=-a;
636  if (b<0) b=-b;
637  if (b>a)
638  {
639    x=b;
640    b=a;
641    a=x;
642  }
643  while (b!=0)
644  {
645    x = a % b;
646    a = b;
647    b = x;
648  }
649  return a;
650}
651
652static intvec * ivOptimizeKern(intvec *kern)
653{
654  int i,l,j,c=kern->cols(),r=kern->rows();
655  intvec *res=new intvec(c);
656
657  //if (TEST_OPT_PROT)
658  //  Warn(" %d linear independent solutions\n",r);
659  for (i=r;i>1;i--)
660  {
661    for (j=c;j>0;j--)
662    {
663      (*res)[j-1] += IMATELEM(*kern,i,j);
664    }
665  }
666  ivContent(res);
667  if (r<11)
668  {
669    l = ivCondNumber(res,-c);
670    j = ivL1Norm(res);
671    ivOptRecursive(res, NULL, kern, l, j, r);
672  }
673  return res;
674}
675
676static void ivOptRecursive(intvec *res, intvec *w, intvec *kern,
677                          int &l, int &j, int pos)
678{
679  int m, k, d;
680  intvec *h;
681
682  d=kern->rows();
683  d=96/(d*d);
684  if (d<3) d=3;
685  if (w!=0)
686    h = new intvec(w);
687  else
688    h = new intvec(res->rows());
689  for (m=d;m>0;m--)
690  {
691    for(k=h->rows()-1;k>=0;k--)
692      (*h)[k] += IMATELEM(*kern,pos,k+1);
693    if(pos>1)
694      ivOptRecursive(res, h, kern, l, j, pos-1);
695    else
696      ivOptSolve(res, h, l, j);
697  }
698  delete h;
699  if (pos>1)
700    ivOptRecursive(res, w, kern, l, j, pos-1);
701  else if (w!=NULL)
702    ivOptSolve(res, w, l, j);
703}
704
705static void ivOptSolve(intvec *res, intvec *w, int &l, int &j)
706{
707  int l0, j0, k;
708
709  l0 = ivCondNumber(w, l);
710  if (l0==l)
711  {
712    ivContent(w);
713    j0 = ivL1Norm(w);
714    if(j0<j)
715    {
716      j = j0;
717      for(k=w->rows()-1;k>=0;k--)
718        (*res)[k] = (*w)[k];
719    }
720    return;
721  }
722  if(l0>l)
723  {
724    l = l0;
725    ivContent(w);
726    j = ivL1Norm(w);
727    for(k=w->rows()-1;k>=0;k--)
728      (*res)[k] = (*w)[k];
729  } 
730}
731
732static int ivL1Norm(intvec *w)
733{
734  int i, j, s = 0;
735
736  for (i=w->rows()-1;i>=0;i--)
737  {
738    j = (*w)[i];
739    if (j>0)
740      s += j;
741    else
742      s -= j;
743  }
744  return s;
745}
746
747static int ivCondNumber(intvec *w, int l)
748{
749  int l0=0, i;
750
751  if (l<0)
752  {
753    for (i=w->rows()-1;i>=0;i--)
754    {
755      if ((*w)[i]<0) l0--;
756    }
757    if (l0==0)
758    {
759      for (i=w->rows()-1;i>=0;i--)
760      {
761        if ((*w)[i]>0) l0++;
762      }
763    }
764    return l0;
765  }
766  else
767  {
768    for (i=w->rows()-1;i>=0;i--)
769    {
770      if ((*w)[i]<0) return -1;
771    }
772    for (i=w->rows()-1;i>=0;i--)
773    {
774      if ((*w)[i]>0) l0++;
775    }
776    return l0;
777  }
778}
779
780static void ivContent(intvec *w)
781{
782  int tgcd, m;
783  int i=w->rows()-1;
784
785  loop
786  {
787    tgcd = (*w)[i--];
788    if (tgcd!=0) break;
789    if (i<0) return;
790  }
791  if (tgcd<0) tgcd = -tgcd;
792  if (tgcd==1) return;
793  loop
794  {
795    m = (*w)[i--];
796    if (m!=0) tgcd= ivGcd(tgcd, m);
797    if (tgcd==1) return;
798    if (i<0) break;
799  }
800  for (i=w->rows()-1;i>=0;i--)
801    (*w)[i] /= tgcd;
802}
803
804#endif
Note: See TracBrowser for help on using the repository browser.