ap.h
Go to the documentation of this file.
1 /********************************************************************
2 AP library.
3 See www.alglib.net or alglib.sources.ru for details.
4 ********************************************************************/
5 
6 #ifndef AP_H
7 #define AP_H
8 
9 #include <stdlib.h>
10 #include <math.h>
11 
12 /********************************************************************
13 Checking of the array boundaries mode.
14 ********************************************************************/
15 //#define NO_AP_ASSERT
16 #define AP_ASSERT
17 
18 #ifndef AP_ASSERT
19 #define NO_AP_ASSERT
20 #endif
21 
22 #ifdef NO_AP_ASSERT
23 #ifdef AP_ASSERT
24 #undef NO_AP_ASSERT
25 #endif
26 #endif
27 
28 
29 /********************************************************************
30 This symbol is used for debugging. Do not define it and do not remove
31 comments.
32 ********************************************************************/
33 //#define UNSAFE_MEM_COPY
34 
35 
36 /********************************************************************
37 Namespace of a standard library AlgoPascal.
38 ********************************************************************/
39 namespace ap
40 {
41 
42 
43 /********************************************************************
44 Exception class.
45 ********************************************************************/
46 class ap_error
47 {
48 public:
49  static void make_assertion(bool bClause)
50  { if(!bClause) throw ap_error(); };
51 private:
52 };
53 
54 /********************************************************************
55 Class defining a complex number with double precision.
56 ********************************************************************/
57 class complex;
58 
59 class complex
60 {
61 public:
62  complex():x(0.0),y(0.0){};
63  complex(const double &_x):x(_x),y(0.0){};
64  complex(const double &_x, const double &_y):x(_x),y(_y){};
65  complex(const complex &z):x(z.x),y(z.y){};
66 
67  complex& operator= (const double& v){ x = v; y = 0.0; return *this; };
68  complex& operator+=(const double& v){ x += v; return *this; };
69  complex& operator-=(const double& v){ x -= v; return *this; };
70  complex& operator*=(const double& v){ x *= v; y *= v; return *this; };
71  complex& operator/=(const double& v){ x /= v; y /= v; return *this; };
72 
73  complex& operator= (const complex& z){ x = z.x; y = z.y; return *this; };
74  complex& operator+=(const complex& z){ x += z.x; y += z.y; return *this; };
75  complex& operator-=(const complex& z){ x -= z.x; y -= z.y; return *this; };
76  complex& operator*=(const complex& z){ double t = x*z.x-y*z.y; y = x*z.y+y*z.x; x = t; return *this; };
77  complex& operator/=(const complex& z)
78  {
79  ap::complex result;
80  double e;
81  double f;
82  if( fabs(z.y)<fabs(z.x) )
83  {
84  e = z.y/z.x;
85  f = z.x+z.y*e;
86  result.x = (z.x+z.y*e)/f;
87  result.y = (z.y-z.x*e)/f;
88  }
89  else
90  {
91  e = z.x/z.y;
92  f = z.y+z.x*e;
93  result.x = (z.y+z.x*e)/f;
94  result.y = (-z.x+z.y*e)/f;
95  }
96  *this = result;
97  return *this;
98  };
99 
100  double x, y;
101 };
102 
103 const complex operator/(const complex& lhs, const complex& rhs);
104 const bool operator==(const complex& lhs, const complex& rhs);
105 const bool operator!=(const complex& lhs, const complex& rhs);
106 const complex operator+(const complex& lhs);
107 const complex operator-(const complex& lhs);
108 const complex operator+(const complex& lhs, const complex& rhs);
109 const complex operator+(const complex& lhs, const double& rhs);
110 const complex operator+(const double& lhs, const complex& rhs);
111 const complex operator-(const complex& lhs, const complex& rhs);
112 const complex operator-(const complex& lhs, const double& rhs);
113 const complex operator-(const double& lhs, const complex& rhs);
114 const complex operator*(const complex& lhs, const complex& rhs);
115 const complex operator*(const complex& lhs, const double& rhs);
116 const complex operator*(const double& lhs, const complex& rhs);
117 const complex operator/(const complex& lhs, const complex& rhs);
118 const complex operator/(const double& lhs, const complex& rhs);
119 const complex operator/(const complex& lhs, const double& rhs);
120 const double abscomplex(const complex &z);
121 const complex conj(const complex &z);
122 const complex csqr(const complex &z);
123 
124 
125 /********************************************************************
126 Template defining vector in memory. It is used by the basic
127 subroutines of linear algebra.
128 
129 Vector consists of Length elements of type T, starting from an element,
130 which Data is pointed to. Interval between adjacent elements equals
131 the value of Step.
132 
133 The class provides an access for reading only.
134 ********************************************************************/
135 template<class T>
136 class const_raw_vector
137 {
138 public:
139  const_raw_vector(const T *Data, int Length, int Step):
140  pData(const_cast<T*>(Data)),iLength(Length),iStep(Step){};
141 
142  const T* GetData() const
143  { return pData; };
144 
145  int GetLength() const
146  { return iLength; };
147 
148  int GetStep() const
149  { return iStep; };
150 protected:
151  T *pData;
152  int iLength, iStep;
153 };
154 
155 
156 /********************************************************************
157 Template defining vector in memory, derived from const_raw_vector.
158 It is used by the basic subroutines of linear algebra.
159 
160 Vector consists of Length elements of type T, starting from an element,
161 which Data is pointed to. Interval between adjacent elements equals
162 the value of Step.
163 
164 The class provides an access both for reading and writing.
165 ********************************************************************/
166 template<class T>
167 class raw_vector : public const_raw_vector<T>
168 {
169 public:
170  raw_vector(T *Data, int Length, int Step):const_raw_vector<T>(Data, Length, Step){};
171 
172  T* GetData()
173  { return const_raw_vector<T>::pData; };
174 };
175 
176 
177 /********************************************************************
178 Scalar product
179 ********************************************************************/
180 template<class T>
181 T vdotproduct(const_raw_vector<T> v1, const_raw_vector<T> v2)
182 {
183  ap_error::make_assertion(v1.GetLength()==v2.GetLength());
184  if( v1.GetStep()==1 && v2.GetStep()==1 )
185  {
186  //
187  // fast
188  //
189  T r = 0;
190  const T *p1 = v1.GetData();
191  const T *p2 = v2.GetData();
192  int imax = v1.GetLength()/4;
193  int i;
194  for(i=imax; i!=0; i--)
195  {
196  r += (*p1)*(*p2) + p1[1]*p2[1] + p1[2]*p2[2] + p1[3]*p2[3];
197  p1+=4;
198  p2+=4;
199  }
200  for(i=0; i<v1.GetLength()%4; i++)
201  r += (*(p1++))*(*(p2++));
202  return r;
203  }
204  else
205  {
206  //
207  // general
208  //
209  int offset11 = v1.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
210  int offset21 = v2.GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
211  T r = 0;
212  const T *p1 = v1.GetData();
213  const T *p2 = v2.GetData();
214  int imax = v1.GetLength()/4;
215  int i;
216  for(i=0; i<imax; i++)
217  {
218  r += (*p1)*(*p2) + p1[offset11]*p2[offset21] + p1[offset12]*p2[offset22] + p1[offset13]*p2[offset23];
219  p1+=offset14;
220  p2+=offset24;
221  }
222  for(i=0; i<v1.GetLength()%4; i++)
223  {
224  r += (*p1)*(*p2);
225  p1+=offset11;
226  p2+=offset21;
227  }
228  return r;
229  }
230 }
231 
232 
233 /********************************************************************
234 Copy one vector into another
235 ********************************************************************/
236 template<class T>
237 void vmove(raw_vector<T> vdst, const_raw_vector<T> vsrc)
238 {
239  ap_error::make_assertion(vdst.GetLength()==vsrc.GetLength());
240  if( vdst.GetStep()==1 && vsrc.GetStep()==1 )
241  {
242  //
243  // fast
244  //
245  T *p1 = vdst.GetData();
246  const T *p2 = vsrc.GetData();
247  int imax = vdst.GetLength()/2;
248  int i;
249  for(i=imax; i!=0; i--)
250  {
251  *p1 = *p2;
252  p1[1] = p2[1];
253  p1 += 2;
254  p2 += 2;
255  }
256  if(vdst.GetLength()%2 != 0)
257  *p1 = *p2;
258  return;
259  }
260  else
261  {
262  //
263  // general
264  //
265  int offset11 = vdst.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
266  int offset21 = vsrc.GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
267  T *p1 = vdst.GetData();
268  const T *p2 = vsrc.GetData();
269  int imax = vdst.GetLength()/4;
270  int i;
271  for(i=0; i<imax; i++)
272  {
273  *p1 = *p2;
274  p1[offset11] = p2[offset21];
275  p1[offset12] = p2[offset22];
276  p1[offset13] = p2[offset23];
277  p1 += offset14;
278  p2 += offset24;
279  }
280  for(i=0; i<vdst.GetLength()%4; i++)
281  {
282  *p1 = *p2;
283  p1 += vdst.GetStep();
284  p2 += vsrc.GetStep();
285  }
286  return;
287  }
288 }
289 
290 
291 /********************************************************************
292 Copy one vector multiplied by -1 into another.
293 ********************************************************************/
294 template<class T>
295 void vmoveneg(raw_vector<T> vdst, const_raw_vector<T> vsrc)
296 {
297  ap_error::make_assertion(vdst.GetLength()==vsrc.GetLength());
298  if( vdst.GetStep()==1 && vsrc.GetStep()==1 )
299  {
300  //
301  // fast
302  //
303  T *p1 = vdst.GetData();
304  const T *p2 = vsrc.GetData();
305  int imax = vdst.GetLength()/2;
306  int i;
307  for(i=0; i<imax; i++)
308  {
309  *p1 = -*p2;
310  p1[1] = -p2[1];
311  p1 += 2;
312  p2 += 2;
313  }
314  if(vdst.GetLength()%2 != 0)
315  *p1 = -*p2;
316  return;
317  }
318  else
319  {
320  //
321  // general
322  //
323  int offset11 = vdst.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
324  int offset21 = vsrc.GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
325  T *p1 = vdst.GetData();
326  const T *p2 = vsrc.GetData();
327  int imax = vdst.GetLength()/4;
328  int i;
329  for(i=imax; i!=0; i--)
330  {
331  *p1 = -*p2;
332  p1[offset11] = -p2[offset21];
333  p1[offset12] = -p2[offset22];
334  p1[offset13] = -p2[offset23];
335  p1 += offset14;
336  p2 += offset24;
337  }
338  for(i=0; i<vdst.GetLength()%4; i++)
339  {
340  *p1 = -*p2;
341  p1 += vdst.GetStep();
342  p2 += vsrc.GetStep();
343  }
344  return;
345  }
346 }
347 
348 
349 /********************************************************************
350 Copy one vector multiplied by a number into another vector.
351 ********************************************************************/
352 template<class T, class T2>
353 void vmove(raw_vector<T> vdst, const_raw_vector<T> vsrc, T2 alpha)
354 {
355  ap_error::make_assertion(vdst.GetLength()==vsrc.GetLength());
356  if( vdst.GetStep()==1 && vsrc.GetStep()==1 )
357  {
358  //
359  // fast
360  //
361  T *p1 = vdst.GetData();
362  const T *p2 = vsrc.GetData();
363  int imax = vdst.GetLength()/4;
364  int i;
365  for(i=imax; i!=0; i--)
366  {
367  *p1 = alpha*(*p2);
368  p1[1] = alpha*p2[1];
369  p1[2] = alpha*p2[2];
370  p1[3] = alpha*p2[3];
371  p1 += 4;
372  p2 += 4;
373  }
374  for(i=0; i<vdst.GetLength()%4; i++)
375  *(p1++) = alpha*(*(p2++));
376  return;
377  }
378  else
379  {
380  //
381  // general
382  //
383  int offset11 = vdst.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
384  int offset21 = vsrc.GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
385  T *p1 = vdst.GetData();
386  const T *p2 = vsrc.GetData();
387  int imax = vdst.GetLength()/4;
388  int i;
389  for(i=0; i<imax; i++)
390  {
391  *p1 = alpha*(*p2);
392  p1[offset11] = alpha*p2[offset21];
393  p1[offset12] = alpha*p2[offset22];
394  p1[offset13] = alpha*p2[offset23];
395  p1 += offset14;
396  p2 += offset24;
397  }
398  for(i=0; i<vdst.GetLength()%4; i++)
399  {
400  *p1 = alpha*(*p2);
401  p1 += vdst.GetStep();
402  p2 += vsrc.GetStep();
403  }
404  return;
405  }
406 }
407 
408 
409 /********************************************************************
410 Vector addition
411 ********************************************************************/
412 template<class T>
413 void vadd(raw_vector<T> vdst, const_raw_vector<T> vsrc)
414 {
415  ap_error::make_assertion(vdst.GetLength()==vsrc.GetLength());
416  if( vdst.GetStep()==1 && vsrc.GetStep()==1 )
417  {
418  //
419  // fast
420  //
421  T *p1 = vdst.GetData();
422  const T *p2 = vsrc.GetData();
423  int imax = vdst.GetLength()/4;
424  int i;
425  for(i=imax; i!=0; i--)
426  {
427  *p1 += *p2;
428  p1[1] += p2[1];
429  p1[2] += p2[2];
430  p1[3] += p2[3];
431  p1 += 4;
432  p2 += 4;
433  }
434  for(i=0; i<vdst.GetLength()%4; i++)
435  *(p1++) += *(p2++);
436  return;
437  }
438  else
439  {
440  //
441  // general
442  //
443  int offset11 = vdst.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
444  int offset21 = vsrc.GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
445  T *p1 = vdst.GetData();
446  const T *p2 = vsrc.GetData();
447  int imax = vdst.GetLength()/4;
448  int i;
449  for(i=0; i<imax; i++)
450  {
451  *p1 += *p2;
452  p1[offset11] += p2[offset21];
453  p1[offset12] += p2[offset22];
454  p1[offset13] += p2[offset23];
455  p1 += offset14;
456  p2 += offset24;
457  }
458  for(i=0; i<vdst.GetLength()%4; i++)
459  {
460  *p1 += *p2;
461  p1 += vdst.GetStep();
462  p2 += vsrc.GetStep();
463  }
464  return;
465  }
466 }
467 
468 
469 /********************************************************************
470 Add one vector multiplied by a number to another vector.
471 ********************************************************************/
472 template<class T, class T2>
473 void vadd(raw_vector<T> vdst, const_raw_vector<T> vsrc, T2 alpha)
474 {
475  ap_error::make_assertion(vdst.GetLength()==vsrc.GetLength());
476  if( vdst.GetStep()==1 && vsrc.GetStep()==1 )
477  {
478  //
479  // fast
480  //
481  T *p1 = vdst.GetData();
482  const T *p2 = vsrc.GetData();
483  int imax = vdst.GetLength()/4;
484  int i;
485  for(i=imax; i!=0; i--)
486  {
487  *p1 += alpha*(*p2);
488  p1[1] += alpha*p2[1];
489  p1[2] += alpha*p2[2];
490  p1[3] += alpha*p2[3];
491  p1 += 4;
492  p2 += 4;
493  }
494  for(i=0; i<vdst.GetLength()%4; i++)
495  *(p1++) += alpha*(*(p2++));
496  return;
497  }
498  else
499  {
500  //
501  // general
502  //
503  int offset11 = vdst.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
504  int offset21 = vsrc.GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
505  T *p1 = vdst.GetData();
506  const T *p2 = vsrc.GetData();
507  int imax = vdst.GetLength()/4;
508  int i;
509  for(i=0; i<imax; i++)
510  {
511  *p1 += alpha*(*p2);
512  p1[offset11] += alpha*p2[offset21];
513  p1[offset12] += alpha*p2[offset22];
514  p1[offset13] += alpha*p2[offset23];
515  p1 += offset14;
516  p2 += offset24;
517  }
518  for(i=0; i<vdst.GetLength()%4; i++)
519  {
520  *p1 += alpha*(*p2);
521  p1 += vdst.GetStep();
522  p2 += vsrc.GetStep();
523  }
524  return;
525  }
526 }
527 
528 
529 /********************************************************************
530 Vector subtraction
531 ********************************************************************/
532 template<class T>
533 void vsub(raw_vector<T> vdst, const_raw_vector<T> vsrc)
534 {
535  ap_error::make_assertion(vdst.GetLength()==vsrc.GetLength());
536  if( vdst.GetStep()==1 && vsrc.GetStep()==1 )
537  {
538  //
539  // fast
540  //
541  T *p1 = vdst.GetData();
542  const T *p2 = vsrc.GetData();
543  int imax = vdst.GetLength()/4;
544  int i;
545  for(i=imax; i!=0; i--)
546  {
547  *p1 -= *p2;
548  p1[1] -= p2[1];
549  p1[2] -= p2[2];
550  p1[3] -= p2[3];
551  p1 += 4;
552  p2 += 4;
553  }
554  for(i=0; i<vdst.GetLength()%4; i++)
555  *(p1++) -= *(p2++);
556  return;
557  }
558  else
559  {
560  //
561  // general
562  //
563  int offset11 = vdst.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
564  int offset21 = vsrc.GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
565  T *p1 = vdst.GetData();
566  const T *p2 = vsrc.GetData();
567  int imax = vdst.GetLength()/4;
568  int i;
569  for(i=0; i<imax; i++)
570  {
571  *p1 -= *p2;
572  p1[offset11] -= p2[offset21];
573  p1[offset12] -= p2[offset22];
574  p1[offset13] -= p2[offset23];
575  p1 += offset14;
576  p2 += offset24;
577  }
578  for(i=0; i<vdst.GetLength()%4; i++)
579  {
580  *p1 -= *p2;
581  p1 += vdst.GetStep();
582  p2 += vsrc.GetStep();
583  }
584  return;
585  }
586 }
587 
588 
589 /********************************************************************
590 Subtract one vector multiplied by a number from another vector.
591 ********************************************************************/
592 template<class T, class T2>
593 void vsub(raw_vector<T> vdst, const_raw_vector<T> vsrc, T2 alpha)
594 {
595  vadd(vdst, vsrc, -alpha);
596 }
597 
598 
599 /********************************************************************
600 In-place vector multiplication
601 ********************************************************************/
602 template<class T, class T2>
603 void vmul(raw_vector<T> vdst, T2 alpha)
604 {
605  if( vdst.GetStep()==1 )
606  {
607  //
608  // fast
609  //
610  T *p1 = vdst.GetData();
611  int imax = vdst.GetLength()/4;
612  int i;
613  for(i=imax; i!=0; i--)
614  {
615  *p1 *= alpha;
616  p1[1] *= alpha;
617  p1[2] *= alpha;
618  p1[3] *= alpha;
619  p1 += 4;
620  }
621  for(i=0; i<vdst.GetLength()%4; i++)
622  *(p1++) *= alpha;
623  return;
624  }
625  else
626  {
627  //
628  // general
629  //
630  int offset11 = vdst.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
631  T *p1 = vdst.GetData();
632  int imax = vdst.GetLength()/4;
633  int i;
634  for(i=0; i<imax; i++)
635  {
636  *p1 *= alpha;
637  p1[offset11] *= alpha;
638  p1[offset12] *= alpha;
639  p1[offset13] *= alpha;
640  p1 += offset14;
641  }
642  for(i=0; i<vdst.GetLength()%4; i++)
643  {
644  *p1 *= alpha;
645  p1 += vdst.GetStep();
646  }
647  return;
648  }
649 }
650 
651 
652 /********************************************************************
653 Template of a dynamical one-dimensional array
654 ********************************************************************/
655 template<class T>
656 class template_1d_array
657 {
658 public:
659  template_1d_array()
660  {
661  m_Vec=0;
662  m_iVecSize = 0;
663  };
664 
665  ~template_1d_array()
666  {
667  if(m_Vec)
668  delete[] m_Vec;
669  };
670 
671  template_1d_array(const template_1d_array &rhs)
672  {
673  m_iVecSize = rhs.m_iVecSize;
674  m_iLow = rhs.m_iLow;
675  m_iHigh = rhs.m_iHigh;
676  if(rhs.m_Vec)
677  {
678  m_Vec = new T[m_iVecSize];
679  #ifndef UNSAFE_MEM_COPY
680  for(int i=0; i<m_iVecSize; i++)
681  m_Vec[i] = rhs.m_Vec[i];
682  #else
683  memcpy(m_Vec, rhs.m_Vec, m_iVecSize*sizeof(T));
684  #endif
685  }
686  else
687  m_Vec=0;
688  };
689 
690 
691  const template_1d_array& operator=(const template_1d_array &rhs)
692  {
693  if( this==&rhs )
694  return *this;
695 
696  m_iLow = rhs.m_iLow;
697  m_iHigh = rhs.m_iHigh;
698  m_iVecSize = rhs.m_iVecSize;
699  if(m_Vec)
700  delete[] m_Vec;
701  if(rhs.m_Vec)
702  {
703  m_Vec = new T[m_iVecSize];
704  #ifndef UNSAFE_MEM_COPY
705  for(int i=0; i<m_iVecSize; i++)
706  m_Vec[i] = rhs.m_Vec[i];
707  #else
708  memcpy(m_Vec, rhs.m_Vec, m_iVecSize*sizeof(T));
709  #endif
710  }
711  else
712  m_Vec=0;
713  return *this;
714  };
715 
716 
717  const T& operator()(int i) const
718  {
719  #ifndef NO_AP_ASSERT
720  ap_error::make_assertion(i>=m_iLow && i<=m_iHigh);
721  #endif
722  return m_Vec[ i-m_iLow ];
723  };
724 
725 
726  T& operator()(int i)
727  {
728  #ifndef NO_AP_ASSERT
729  ap_error::make_assertion(i>=m_iLow && i<=m_iHigh);
730  #endif
731  return m_Vec[ i-m_iLow ];
732  };
733 
734 
735  void setbounds( int iLow, int iHigh )
736  {
737  if(m_Vec)
738  delete[] m_Vec;
739  m_iLow = iLow;
740  m_iHigh = iHigh;
741  m_iVecSize = iHigh-iLow+1;
742  m_Vec = new T[m_iVecSize];
743  };
744 
745 
746  void setcontent( int iLow, int iHigh, const T *pContent )
747  {
748  setbounds(iLow, iHigh);
749  for(int i=iLow; i<=iHigh; i++)
750  (*this)(i) = pContent[i-iLow];
751  };
752 
753 
754  T* getcontent()
755  {
756  return m_Vec;
757  };
758 
759  const T* getcontent() const
760  {
761  return m_Vec;
762  };
763 
764 
765  int getlowbound(int iBoundNum = 0) const
766  {
767  return m_iLow;
768  };
769 
770 
771  int gethighbound(int iBoundNum = 0) const
772  {
773  return m_iHigh;
774  };
775 
776  raw_vector<T> getvector(int iStart, int iEnd)
777  {
778  if( iStart>iEnd || wrongIdx(iStart) || wrongIdx(iEnd) )
779  return raw_vector<T>(0, 0, 1);
780  else
781  return raw_vector<T>(m_Vec+iStart-m_iLow, iEnd-iStart+1, 1);
782  };
783 
784 
785  const_raw_vector<T> getvector(int iStart, int iEnd) const
786  {
787  if( iStart>iEnd || wrongIdx(iStart) || wrongIdx(iEnd) )
788  return const_raw_vector<T>(0, 0, 1);
789  else
790  return const_raw_vector<T>(m_Vec+iStart-m_iLow, iEnd-iStart+1, 1);
791  };
792 private:
793  bool wrongIdx(int i) const { return i<m_iLow || i>m_iHigh; };
794 
795  T *m_Vec;
796  long m_iVecSize;
797  long m_iLow, m_iHigh;
798 };
799 
800 
801 
802 /********************************************************************
803 Template of a dynamical two-dimensional array
804 ********************************************************************/
805 template<class T>
806 class template_2d_array
807 {
808 public:
809  template_2d_array()
810  {
811  m_Vec=0;
812  m_iVecSize=0;
813  };
814 
815  ~template_2d_array()
816  {
817  if(m_Vec)
818  delete[] m_Vec;
819  };
820 
821  template_2d_array(const template_2d_array &rhs)
822  {
823  m_iVecSize = rhs.m_iVecSize;
824  m_iLow1 = rhs.m_iLow1;
825  m_iLow2 = rhs.m_iLow2;
826  m_iHigh1 = rhs.m_iHigh1;
827  m_iHigh2 = rhs.m_iHigh2;
828  m_iConstOffset = rhs.m_iConstOffset;
829  m_iLinearMember = rhs.m_iLinearMember;
830  if(rhs.m_Vec)
831  {
832  m_Vec = new T[m_iVecSize];
833  #ifndef UNSAFE_MEM_COPY
834  for(int i=0; i<m_iVecSize; i++)
835  m_Vec[i] = rhs.m_Vec[i];
836  #else
837  memcpy(m_Vec, rhs.m_Vec, m_iVecSize*sizeof(T));
838  #endif
839  }
840  else
841  m_Vec=0;
842  };
843  const template_2d_array& operator=(const template_2d_array &rhs)
844  {
845  if( this==&rhs )
846  return *this;
847 
848  m_iLow1 = rhs.m_iLow1;
849  m_iLow2 = rhs.m_iLow2;
850  m_iHigh1 = rhs.m_iHigh1;
851  m_iHigh2 = rhs.m_iHigh2;
852  m_iConstOffset = rhs.m_iConstOffset;
853  m_iLinearMember = rhs.m_iLinearMember;
854  m_iVecSize = rhs.m_iVecSize;
855  if(m_Vec)
856  delete[] m_Vec;
857  if(rhs.m_Vec)
858  {
859  m_Vec = new T[m_iVecSize];
860  #ifndef UNSAFE_MEM_COPY
861  for(int i=0; i<m_iVecSize; i++)
862  m_Vec[i] = rhs.m_Vec[i];
863  #else
864  memcpy(m_Vec, rhs.m_Vec, m_iVecSize*sizeof(T));
865  #endif
866  }
867  else
868  m_Vec=0;
869  return *this;
870  };
871 
872  const T& operator()(int i1, int i2) const
873  {
874  #ifndef NO_AP_ASSERT
875  ap_error::make_assertion(i1>=m_iLow1 && i1<=m_iHigh1);
876  ap_error::make_assertion(i2>=m_iLow2 && i2<=m_iHigh2);
877  #endif
878  return m_Vec[ m_iConstOffset + i2 +i1*m_iLinearMember];
879  };
880 
881  T& operator()(int i1, int i2)
882  {
883  #ifndef NO_AP_ASSERT
884  ap_error::make_assertion(i1>=m_iLow1 && i1<=m_iHigh1);
885  ap_error::make_assertion(i2>=m_iLow2 && i2<=m_iHigh2);
886  #endif
887  return m_Vec[ m_iConstOffset + i2 +i1*m_iLinearMember];
888  };
889 
890  void setbounds( int iLow1, int iHigh1, int iLow2, int iHigh2 )
891  {
892  if(m_Vec)
893  delete[] m_Vec;
894  m_iVecSize = (iHigh1-iLow1+1)*(iHigh2-iLow2+1);
895  m_Vec = new T[m_iVecSize];
896  m_iLow1 = iLow1;
897  m_iHigh1 = iHigh1;
898  m_iLow2 = iLow2;
899  m_iHigh2 = iHigh2;
900  m_iConstOffset = -m_iLow2-m_iLow1*(m_iHigh2-m_iLow2+1);
901  m_iLinearMember = (m_iHigh2-m_iLow2+1);
902  };
903 
904  void setcontent( int iLow1, int iHigh1, int iLow2, int iHigh2, const T *pContent )
905  {
906  setbounds(iLow1, iHigh1, iLow2, iHigh2);
907  for(int i=0; i<m_iVecSize; i++)
908  m_Vec[i]=pContent[i];
909  };
910 
911  T* getcontent()
912  {
913  return m_Vec;
914  };
915 
916  const T* getcontent() const
917  {
918  return m_Vec;
919  };
920 
921  int getlowbound(int iBoundNum) const
922  {
923  return iBoundNum==1 ? m_iLow1 : m_iLow2;
924  };
925 
926  int gethighbound(int iBoundNum) const
927  {
928  return iBoundNum==1 ? m_iHigh1 : m_iHigh2;
929  };
930 
931  raw_vector<T> getcolumn(int iColumn, int iRowStart, int iRowEnd)
932  {
933  if( (iRowStart>iRowEnd) || wrongColumn(iColumn) || wrongRow(iRowStart) ||wrongRow(iRowEnd) )
934  return raw_vector<T>(0, 0, 1);
935  else
936  return raw_vector<T>(&((*this)(iRowStart, iColumn)), iRowEnd-iRowStart+1, m_iLinearMember);
937  };
938 
939  raw_vector<T> getrow(int iRow, int iColumnStart, int iColumnEnd)
940  {
941  if( (iColumnStart>iColumnEnd) || wrongRow(iRow) || wrongColumn(iColumnStart) || wrongColumn(iColumnEnd))
942  return raw_vector<T>(0, 0, 1);
943  else
944  return raw_vector<T>(&((*this)(iRow, iColumnStart)), iColumnEnd-iColumnStart+1, 1);
945  };
946 
947  const_raw_vector<T> getcolumn(int iColumn, int iRowStart, int iRowEnd) const
948  {
949  if( (iRowStart>iRowEnd) || wrongColumn(iColumn) || w