source: git/libpolys/misc/intvec.cc @ 71f6dbd

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