source: git/gfanlib/gfanlib_vector.h @ 26b713

jengelh-datetimespielwiese
Last change on this file since 26b713 was 26b713, checked in by Yue Ren <ren@…>, 10 years ago
chg: updated gfanlib package to newest version in master
  • Property mode set to 100644
File size: 11.1 KB
Line 
1/*
2 * lib_zvector.h
3 *
4 *  Created on: Sep 6, 2010
5 *      Author: anders
6 */
7
8#ifndef GFANLIB_ZVECTOR_H_INCLUDED
9#define GFANLIB_ZVECTOR_H_INCLUDED
10
11#include <vector>
12#include <list>
13#include <assert.h>
14#include <algorithm>
15#include <iostream>
16
17#include "gfanlib_z.h"
18#include "gfanlib_q.h"
19
20namespace gfan{
21
22inline void outOfRange(int i, int n)
23{
24  std::cerr<<"Index out of range. i="<<i<<" n="<<n<<std::endl;
25  assert(0);
26}
27
28template <class typ> class Vector{
29private:
30  std::vector<typ> v;
31public:
32  //--------------
33  // Constructors
34  //--------------
35  Vector(const Vector &a):v(a.v){
36  }
37  Vector(int n):v(n){
38    assert(n>=0);
39  };
40/*  Vektor(const typ *p, int n):v(n){
41    assert(n>=0);
42    support=0;for(int i=0;i<n;i++)v[i]=p[i];
43  };*/
44  ~Vector(){
45  };
46  Vector(){
47  };
48  static Vector standardVector(int n, int i)
49    {
50      Vector v(n);
51      v[i]=typ(1);
52      return v;
53    }
54  static Vector allOnes(int n)
55    {
56      Vector v(n);
57      for(int i=0;i<n;i++)
58        v[i]=typ(1);
59      return v;
60    }
61  //-------------
62  // Conversions
63  //-------------
64/*  template<class T> Vektor(const Vektor<T>& c)
65          :v(c.size())
66        {
67          for(int i=0;i<size();i++)v[i]=typ(c[i]);
68        }
69*/
70  //--------
71  // Access
72  //--------
73  typ& operator[](int n)
74    {
75      if(!(n>=0 && n<(int)v.size()))outOfRange(n,v.size());
76      return (v[n]);
77    }
78  const typ& operator[](int n)const{assert(n>=0 && n<(int)v.size());return (v[n]);}
79  const typ& UNCHECKEDACCESS(int n)const{return (v[n]);}
80
81  //-------------
82  // STL related
83  //-------------
84  unsigned int size()const{return v.size();};
85  void resize(int n){v.resize(n,typ());};
86  void grow(int i){if(size()<i)resize(i);}
87  void push_back(typ a)
88  {
89    v.push_back(a);
90  }
91  void sort()
92    {
93      std::sort(v.begin(),v.end());
94    }
95  bool nextPermutation()
96    {
97      return std::next_permutation(v.begin(),v.end());
98    }
99  bool operator<(const Vector & b)const
100    {
101      if(size()<b.size())return true;
102      if(size()>b.size())return false;
103      for(int i=0;i<size();i++)
104        {
105          if(v[i]<b[i])return true;
106          if(b[i]<v[i])return false;
107        }
108      return false;
109    }
110
111  //-----------------
112  // Arithmetic fast
113  //-----------------
114  typ sum()const{typ f;for(typename std::vector<typ>::const_iterator i=v.begin();i!=v.end();i++)f+=*i;return f;};
115  Vector& operator+=(const Vector& q){assert(size()==q.size());typename std::vector<typ>::const_iterator j=q.v.begin();for(typename std::vector<typ>::iterator i=v.begin();i!=v.end();++i,++j)(*i)+=(*j);return *this;}
116  Vector& operator-=(const Vector& q){assert(size()==q.size());typename std::vector<typ>::const_iterator j=q.v.begin();for(typename std::vector<typ>::iterator i=v.begin();i!=v.end();++i,++j)(*i)-=(*j);return *this;}
117  Vector& operator/=(const Vector& q){assert(size()==q.size());typename std::vector<typ>::const_iterator j=q.v.begin();for(typename std::vector<typ>::iterator i=v.begin();i!=v.end();++i,++j)(*i)/=(*j);return *this;}
118  inline friend typ dot(const Vector& p, const Vector& q){assert(p.size()==q.size());typ s;typename std::vector<typ>::const_iterator j=q.v.begin();for(typename std::vector<typ>::const_iterator i=p.v.begin();i!=p.v.end();++i,++j)s+=(*i)*(*j);return s;}
119//  inline friend int64 dotLong(const Vektor& p, const Vektor& q){assert(p.size()==q.size());int64 s=0;for(int i=0;i<p.size();i++)s+=(int64)p[i]*(int64)q[i];return s;}
120  bool operator==(const Vector & q)const{if(size()!=q.size())return false;typename std::vector<typ>::const_iterator j=q.v.begin();for(typename std::vector<typ>::const_iterator i=v.begin();i!=v.end();++i,++j)if((*i)!=(*j))return false;return true;}
121  bool operator!=(const Vector & q)const {return !(operator==(q));}
122  bool isZero() const
123    {
124      for(typename std::vector<typ>::const_iterator i=v.begin();i!=v.end();++i)if(!(*i).isZero())return false;
125      return true;
126    }
127  bool isPositive() const
128    {
129      for(typename std::vector<typ>::const_iterator i=v.begin();i!=v.end();++i)if((*i).sign()<=0)return false;
130      return true;
131    }
132  bool isNonNegative() const
133    {
134      for(typename std::vector<typ>::const_iterator i=v.begin();i!=v.end();++i)if(i->sign()<0)return false;
135      return true;
136    }
137/*  int max()const
138  {
139    int ret=-0x7fffffff; //not completely correct, but kind of works for 64bit
140    for(int i=0;i<v.size();i++)if(ret<v[i])ret=v[i];
141    return ret;
142  }
143  int min()const
144  {
145    int ret=0x7fffffff;
146    for(int i=0;i<v.size();i++)if(ret>v[i])ret=v[i];
147    return ret;
148  }
149*/
150  friend bool dependent(const Vector& p, const Vector& q)
151    {
152          /*
153      typ pp=dot(p,p);
154      typ qq=dot(q,q);
155      typ pq=dot(p,q);
156      return pq*pq==pp*qq;
157*/
158          int n=p.size();
159          assert(n==q.size());
160          int i;
161          for(i=0;i<n;i++)
162          {
163            if(!p.v[i].isZero())break;
164          }
165          if(i==n)return true;
166          if(q.v[i].isZero())return q.isZero();
167          typ a=p.v[i];
168          typ b=q.v[i];
169          for(int j=0;j<n;j++)
170            if(a*q.v[j]!=b*p.v[j])return false;
171          return true;
172    }
173
174  //-----------------
175  // Arithmetic slow
176  //-----------------
177  inline friend Vector operator*(typ s, const Vector& q){Vector p=q;for(int i=0;i<q.size();i++)p[i]*=s;return p;}
178//  inline friend Vektor operator/(const Vektor& q, typ s){Vektor p=q;for(int i=0;i<q.size();i++)p[i]/=s;return p;}
179/*  inline friend Vector operator*(const Vektor& p, const Vektor& q){assert(p.size()==q.size());Vektor p1=p;for(int i=0;i<p.size();i++)p1.v[i]*=q.v[i];return p1;}
180  inline friend Vektor operator+(const Vektor& p, const Vektor& q){assert(p.size()==q.size());Vektor p1=p;for(int i=0;i<p.size();i++)p1[i]+=q[i];return p1;}
181*/
182  inline friend Vector operator/(const Vector& p, typ const &s){Vector ret(p.size());for(unsigned i=0;i<p.size();i++)ret[i]=p[i]/s;return ret;}
183  inline friend Vector operator+(const Vector& p, const Vector& q){assert(p.size()==q.size());Vector p1=p;for(int i=0;i<p.size();i++)p1[i]+=q[i];return p1;}
184  inline friend Vector operator-(const Vector& p, const Vector& q){assert(p.size()==q.size());Vector p1=p;for(int i=0;i<p.size();i++)p1[i]-=q[i];return p1;}
185  inline friend Vector max(const Vector& p, const Vector& q){assert(p.size()==q.size());Vector p1=p;for(int i=0;i<p.size();i++)if(p1[i]<q[i])p1[i]=q[i];return p1;}
186  inline friend Vector min(const Vector& p, const Vector& q){assert(p.size()==q.size());Vector p1=p;for(int i=0;i<p.size();i++)if(p1[i]>q[i])p1[i]=q[i];return p1;}
187
188  friend Vector operator-(const Vector &b)
189  {
190    Vector ret(b.size());
191    for(int i=0;i<b.size();i++)ret[i]=-b[i];
192    return ret;
193  }
194  //------------------
195  // Monomial related
196  //------------------
197/*  int divides(const Vektor& q) const
198    {
199      assert(size()==q.size());
200      int n=v.size();
201      for(int i=0;i<n;i++)
202        {
203          if(v[i]>0)if(q.v[i]<v[i])return 0;
204        }
205      return 1;
206    }
207  inline friend bool relativelyPrime(const Vektor& p, const Vektor& q)
208    {
209      assert(p.size()==q.size());
210      int n=p.size();
211      for(int t=0;t<n;t++)if((p[t]>0)&&(q[t]>0)) return false;
212      return true;
213    }
214  Vektor supportVector()const
215    {
216      Vektor r(v.size());
217      for(int i=0;i<size();i++)
218        r[i]=(v[i]!=0);
219      return r;
220    }
221*/
222  //------------------------------
223  // Subvectors and concatenation
224  //------------------------------
225  Vector subvector(int begin, int end)const
226    {
227      assert(begin>=0);
228      assert(end<=size());
229      assert(end>=begin);
230      Vector ret(end-begin);
231      for(int i=0;i<end-begin;i++)
232        ret[i]=v[begin+i];
233      return ret;
234    }
235/*  Vektor subvector(list<int> const &chosenIndices)const
236  {
237    Vektor ret(chosenIndices.size());
238    int I=0;
239    for(list<int>::const_iterator i=chosenIndices.begin();i!=chosenIndices.end();i++,I++)ret[I]=v[*i];
240    return ret;
241  }
242*/
243  friend Vector concatenation(Vector const &a, Vector const &b)
244  {
245    Vector ret(a.size()+b.size());
246    for(int i=0;i<a.size();i++)ret[i]=a[i];
247    for(int i=0;i<b.size();i++)ret[i+a.size()]=b[i];
248    return ret;
249  }
250  //----------------------
251  // Zero/nonZero entries
252  //----------------------
253/*  int indexOfLargestNonzeroEntry()const
254  {
255    int ret=-1;
256    for(int i=0;i<v.size();i++)
257      {
258        if(v[i])ret=i;
259      }
260    return ret;
261  }
262  Vektor supportIndices()const
263  {
264    Vektor ret(0);
265    for(int i=0;i<v.size();i++)
266      if(v[i]!=0)ret.push_back(i);
267    return ret;
268  }
269  Vektor supportAsZeroOneVector()const
270  {
271    Vektor ret(v.size());
272    for(int i=0;i<v.size();i++)ret[i]=bool(v[i]);
273    return ret;
274  }
275  void calcsupport(void)
276    {
277      support=0;
278      for(int i=0;i<v.size();i++)support=(support<<1)|(((v[i]>0)==true)&1);
279    }
280*/
281  friend std::ostream &operator<<(std::ostream &f, Vector const &a)
282  {
283    f<<"(";
284    for(typename std::vector<typ>::const_iterator i=a.v.begin();i!=a.v.end();i++)
285    {
286      if(i!=a.v.begin()) f<<",";
287      f<<*i;
288    }
289    return f<<")";
290  }
291  typ gcd()const
292  {
293    typ temp1,temp2;
294    typ ret(1);
295    for(int i=0;i<size();i++)
296      ret=typ::gcd(ret,v[i],temp1,temp2);
297    return ret;
298  }
299  Vector normalized()const
300  {
301    assert(!typ::isField());
302    return (*this)/gcd();
303  }
304};
305
306typedef Vector<Integer> ZVector;
307typedef Vector<Rational> QVector;
308typedef Vector<int> IntVector;
309
310inline QVector ZToQVector(ZVector const &v)
311{
312  QVector ret(v.size());
313  for(unsigned i=0;i<v.size();i++)ret[i]=Rational(v[i]);
314  return ret;
315}
316
317
318inline IntVector ZToIntVector(ZVector const &v)
319{
320  IntVector ret(v.size());
321  for(unsigned i=0;i<v.size();i++)ret[i]=v[i].toInt();
322  return ret;
323}
324
325
326inline ZVector IntToZVector(IntVector const &v)
327{
328  ZVector ret(v.size());
329  for(unsigned i=0;i<v.size();i++)ret[i]=Integer(v[i]);
330  return ret;
331}
332
333inline ZVector QToZVectorPrimitive(QVector const &v)
334{
335    int n=v.size();
336    ZVector ret(n);
337
338    mpz_t lcm;
339    mpz_t gcd;
340    mpz_init_set_ui(lcm, 1);
341    mpz_init_set_ui(gcd, 0);
342
343    mpq_t a;
344    mpq_init(a);
345    for(int j=0;j<n;j++)
346      {
347        v[j].setGmp(a);
348        if(mpz_cmp_si(mpq_denref(a),1)!=0)
349          mpz_lcm(lcm,lcm,mpq_denref(a));
350        if(mpz_sgn(mpq_numref(a))!=0)
351          mpz_gcd(gcd,gcd,mpq_numref(a));
352      }
353    mpq_clear(a);
354    if(mpz_sgn(gcd)!=0)//v is non-zero
355      {
356        if((mpz_cmp_si(lcm,1)==0)&&(mpz_cmp_si(gcd,1)==0)) //gcd=lcm=1
357          {
358            mpq_t a;
359            mpq_init(a);
360
361            for(int i=0;i<n;i++)
362              {
363                v[i].setGmp(a);
364                ret[i]=Integer(mpq_numref(a));
365              }
366            mpq_clear(a);
367          }
368        else
369          {
370            mpq_t a;
371            mpq_init(a);
372            mpz_t tempA;
373            mpz_t tempB;
374            mpz_init(tempA);
375            mpz_init(tempB);
376            for(int i=0;i<n;i++)
377              {
378                v[i].setGmp(a);
379                mpz_set(tempA,mpq_denref(a));
380                mpz_set(tempB,mpq_numref(a));
381                mpz_mul(tempA,gcd,tempA);
382                mpz_mul(tempB,lcm,tempB);
383                mpz_divexact(tempA,tempB,tempA);
384                ret[i]=Integer(tempA);
385              }
386            mpz_clear(tempB);
387            mpz_clear(tempA);
388            mpq_clear(a);
389          }
390      }
391    mpz_clear(gcd);
392    mpz_clear(lcm);
393
394    return ret;
395}
396
397}
398
399
400#endif /* LIB_ZVECTOR_H_ */
Note: See TracBrowser for help on using the repository browser.