source: git/libpolys/misc/intvec.cc @ 16f511

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