source: git/libpolys/misc/intvec.cc @ 9c83f2

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