source: git/kernel/intvec.cc @ fec53d

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