source: git/kernel/intvec.cc @ a41f3aa

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