source: git/kernel/intvec.cc @ e9c3b2

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