vmmlib  1.8.0
Templatized C++ vector and matrix math library
 All Classes Namespaces Files Functions Variables Typedefs Pages
vector.hpp
1 /*
2  * Copyright (c) 2006-2015, Visualization and Multimedia Lab,
3  * University of Zurich <http://vmml.ifi.uzh.ch>,
4  * Eyescale Software GmbH,
5  * Blue Brain Project, EPFL
6  *
7  * This file is part of VMMLib <https://github.com/VMML/vmmlib/>
8  *
9  * Redistribution and use in source and binary forms, with or without
10  * modification, are permitted provided that the following conditions are met:
11  *
12  * Redistributions of source code must retain the above copyright notice, this
13  * list of conditions and the following disclaimer. Redistributions in binary
14  * form must reproduce the above copyright notice, this list of conditions and
15  * the following disclaimer in the documentation and/or other materials provided
16  * with the distribution. Neither the name of the Visualization and Multimedia
17  * Lab, University of Zurich nor the names of its contributors may be used to
18  * endorse or promote products derived from this software without specific prior
19  * written permission.
20  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
21  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
22  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
23  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
24  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
25  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
26  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
27  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
28  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
29  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
30  * POSSIBILITY OF SUCH DAMAGE.
31  */
32 
33 #ifndef __VMML__VECTOR__HPP__
34 #define __VMML__VECTOR__HPP__
35 
36 #include <vmmlib/vmmlib_config.hpp>
37 #include <vmmlib/math.hpp>
38 #include <vmmlib/enable_if.hpp>
39 #include <vmmlib/exception.hpp>
40 
41 #include <iostream>
42 #include <iomanip>
43 #include <vector>
44 #include <string>
45 #include <cstring>
46 #include <limits>
47 #include <algorithm>
48 
49 namespace vmml
50 {
51 
52 template< size_t M, typename T = float >
53 class vector
54 {
55 public:
56  typedef T value_type;
57  typedef T* iterator;
58  typedef const T* const_iterator;
59  typedef std::reverse_iterator< iterator > reverse_iterator;
60  typedef std::reverse_iterator< const_iterator > const_reverse_iterator;
61 
62  static const size_t DIMENSION = M;
63 
64  // constructors
65  vector() : array() {} // http://stackoverflow.com/questions/5602030
66  explicit vector( const T& a ); // sets all components to a;
67  vector( const T& x, const T& y );
68  vector( const T& x, const T& y, const T& z );
69  vector( const T& x, const T& y, const T& z, const T& w );
70 
71 #ifndef SWIG
72  // initializes the first M-1 values from vector_, the last from last_
73  vector( const vector< M-1, T >& vector_, T last_ );
74 #endif
75 
76  vector( const T* values );
77 
78 #ifdef __OSG_MATH
79  template< typename OSGVEC3 >
80  explicit vector( const OSGVEC3& from,
81  typename enable_if< M == 3, OSGVEC3 >::type* = 0 );
82 #endif
83 
84  // vec< M > with homogeneous coordinates <-> vec< M-1 > conversion ctor
85  // to-homogenous-coordinates ctor
86  template< size_t N >
87  vector( const vector< N, T >& source_,
88  typename enable_if< N == M - 1 >::type* = 0 );
89 
90  // from-homogenous-coordinates vector
91  template< size_t N >
92  vector( const vector< N, T >& source_,
93  typename enable_if< N == M + 1 >::type* = 0 );
94 
95  template< typename U > vector( const vector< M, U >& source_ );
96 
97  // iterators
98  inline iterator begin();
99  inline iterator end();
100  inline const_iterator begin() const;
101  inline const_iterator end() const;
102  inline reverse_iterator rbegin();
103  inline reverse_iterator rend();
104  inline const_reverse_iterator rbegin() const;
105  inline const_reverse_iterator rend() const;
106 
107 # ifndef VMMLIB_NO_CONVERSION_OPERATORS
108  // conversion operators
109  inline operator T*();
110  inline operator const T*() const;
111 # else
112  inline T& operator[]( size_t index );
113  inline const T& operator[]( size_t index ) const;
114 # endif
115 
116  // accessors
117  inline T& operator()( size_t index );
118  inline const T& operator()( size_t index ) const;
119 
120  inline T& at( size_t index );
121  inline const T& at( size_t index ) const;
122 
123  // element accessors for M <= 4;
124  inline T& x();
125  inline T& y();
126  inline T& z();
127  inline T& w();
128  inline const T& x() const;
129  inline const T& y() const;
130  inline const T& z() const;
131  inline const T& w() const;
132 
133  // pixel color element accessors for M<= 4
134  inline T& r();
135  inline T& g();
136  inline T& b();
137  inline T& a();
138  inline const T& r() const;
139  inline const T& g() const;
140  inline const T& b() const;
141  inline const T& a() const;
142 
143  bool operator==( const vector& other ) const;
144  bool operator!=( const vector& other ) const;
145  bool equals( const vector& other,
146  T tolerance = std::numeric_limits< T >::epsilon() ) const;
147  bool operator<( const vector& other ) const;
148 
149  // remember kids: c_arrays are dangerous and evil!
150  vector& operator=( const T* c_array );
151  T operator=( T filler );
152 
153  vector& operator=( const vector& other );
154 
155  // returns void to avoid 'silent' loss of precision when chaining
156  template< typename U > void operator=( const vector< M, U >& other );
157 
158  // to-homogenous-coordinates assignment operator
159  // non-chainable because of sfinae
160  template< size_t N >
161  typename enable_if< N == M - 1 >::type*
162  operator=( const vector< N, T >& source_ );
163 
164  // from-homogenous-coordinates assignment operator
165  // non-chainable because of sfinae
166  template< size_t N >
167  typename enable_if< N == M + 1 >::type*
168  operator=( const vector< N, T >& source_ );
169 
170  vector operator*( const vector& other ) const;
171  vector operator/( const vector& other ) const;
172  vector operator+( const vector& other ) const;
173  vector operator-( const vector& other ) const;
174 
175  void operator*=( const vector& other );
176  void operator/=( const vector& other );
177  void operator+=( const vector& other );
178  void operator-=( const vector& other );
179 
180  vector operator*( const T other ) const;
181  vector operator/( const T other ) const;
182  vector operator+( const T other ) const;
183  vector operator-( const T other ) const;
184 
185  void operator*=( const T other );
186  void operator/=( const T other );
187  void operator+=( const T other );
188  void operator-=( const T other );
189 
190  vector operator-() const;
191 
192  const vector& negate();
193 
194  void set( T a ); // sets all components to a;
195 #ifndef SWIG
196  void set( const vector< M-1, T >& v, T a );
197 #endif
198  template< size_t N >
199  void set( const vector< N, T >& v );
200 
201  // sets the first few components to a certain value
202  void set( T x, T y );
203  void set( T x, T y, T z );
204  void set( T x, T y, T z, T w );
205 
206  template< typename input_iterator_t >
207  void iter_set( input_iterator_t begin_, input_iterator_t end_ );
208 
209  // compute the cross product of two vectors
210  // note: there's also a free function:
211  // vector<> cross( const vector<>, const vector<> )
212 
213  // result = vec1.cross( vec2 ) => retval result = vec1 x vec2
214  template< typename TT >
215  vector cross( const vector< M, TT >& rhs,
216  typename enable_if< M == 3, TT >::type* = 0 ) const;
217 
218  // result.cross( vec1, vec2 ) => (this) = vec1 x vec2
219  template< typename TT >
220  void cross( const vector< M, TT >& a, const vector< M, TT >& b,
221  typename enable_if< M == 3, TT >::type* = 0 );
222 
223 
224  // compute the dot product of two vectors
225  // note: there's also a free function:
226  // T dot( const vector<>, const vector<> );
227  inline T dot( const vector& other ) const;
228 
229 
230  // normalize the vector
231  // note: there's also a free function:
232  // vector<> normalize( const vector<> );
233  inline T normalize();
234 
235  //sets all vector components to random values
236  //remember to set srand( seed );
237  //if seed is set to -1, srand( seed ) was set outside set_random
238  //otherwise srand( seed ) will be called with the given seed
239  void set_random( int seed = -1 );
240 
241  inline T length() const;
242  inline T squared_length() const;
243 
244  inline T distance( const vector& other ) const;
245  inline T squared_distance( const vector& other ) const;
246 
248  T product() const;
249 
250  template< typename TT >
251  vector< 3, T > rotate( const T theta, vector< M, TT > axis,
252  typename enable_if< M == 3, TT >::type* = 0 ) const;
253 
254  // right hand system, CCW triangle
255  // (*this) = normal of v0,v1,v2
256  void compute_normal( const vector& v0, const vector& v1, const vector& v2 );
257  // retval = normal of (this), v1, v2
258  vector compute_normal( const vector& v1, const vector& v2 ) const;
259 
261  template< size_t N >
262  vector< N, T >& get_sub_vector( size_t offset = 0,
263  typename enable_if< M >= N >::type* = 0 );
264 
266  template< size_t N >
267  const vector< N, T >& get_sub_vector( size_t offset = 0,
268  typename enable_if< M >= N >::type* = 0 ) const;
269 
270  // sphere functions - sphere layout: center xyz, radius w
271  template< typename TT >
272  inline vector< 3, T > project_point_onto_sphere(
273  const vector< 3, TT >& point,
274  typename enable_if< M == 4, TT >::type* = 0 ) const;
275 
276  // returns a negative distance if the point lies in the sphere
277  template< typename TT >
278  inline T distance_to_sphere( const vector< 3, TT >& point,
279  typename enable_if< M == 4, TT >::type* = 0 ) const;
280 
281  // plane functions - plane layout; normal xyz, distance w
282  template< typename TT >
283  inline T distance_to_plane( const vector< 3, TT >& point,
284  typename enable_if< M == 4, TT >::type* = 0 ) const;
285 
286  template< typename TT >
287  inline vector< 3, T > project_point_onto_plane(
288  const vector< 3, TT >& point,
289  typename enable_if< M == 4, TT >::type* = 0 ) const;
290 
291  // returns the index of the minimal resp. maximal value in the vector
292  size_t find_min_index() const;
293  size_t find_max_index() const;
294 
295  // returns the index of the minimal resp. maximal value in the vector
296  size_t find_abs_min_index() const;
297  size_t find_abs_max_index() const;
298 
299  // returns minimal resp. maximal value in the vector
300  T& find_min();
301  T& find_max();
302  const T& find_min() const;
303  const T& find_max() const;
304 
305  void clamp( const T& min = 0.0, const T& max = 1.0 );
306 
307  template< typename TT >
308  void scale_to( vector< M, TT >& scaled_vector,
309  T min_value = -1.0, T max_value = 1.0 ) const;
310 
311  inline static size_t size(); // returns M
312 
313  bool is_unit_vector() const;
314 
315  // perturbs each component by randomly + or - the perturbation parameter
316  void perturb( T perturbation = 0.0001 );
317 
318  void sqrt_elementwise();
319  double norm() const; //l2 norm
320 
321  // computes the reciprocal value for each component, x = 1/x;
322  // WARNING: might result in nans if division by 0!
323  void reciprocal();
324  // computes the reciprocal value for each component, x = 1/x;
325  // checks every component for 0, sets to max value if zero.
326  void reciprocal_safe();
327 
328  template< typename TT >
329  void cast_from( const vector< M, TT >& other );
330 
331  size_t nnz() const;
332 
333  // test each component of the vector for isnan and isinf
334  // inline bool is_valid() const; -> moved to class validator
335 
336  friend std::ostream& operator<< ( std::ostream& os, const vector& vector_ )
337  {
338  const std::ios::fmtflags flags = os.flags();
339  const int prec = os.precision();
340 
341  os.setf( std::ios::right, std::ios::adjustfield );
342  os.precision( 5 );
343  os << "[ ";
344  for( size_t index = 0; index < M; ++index )
345  os << std::setw(10) << vector_.at( index ) << " ";
346  os << "]";
347  os.precision( prec );
348  os.setf( flags );
349  return os;
350  }
351 
352  T array[ M ];
353 
354 #ifndef SWIG
355  // Vector3 defaults
356  static const vector FORWARD;
357  static const vector BACKWARD;
358  static const vector UP;
359  static const vector DOWN;
360  static const vector LEFT;
361  static const vector RIGHT;
362 
363  static const vector ONE;
364  static const vector ZERO;
365 
366  // Unit vectors
367  static const vector UNIT_X;
368  static const vector UNIT_Y;
369  static const vector UNIT_Z;
370 #endif
371 
372 }; // class vector
373 
374 //
375 // typedefs and statics
376 //
377 #ifndef SWIG
378 template< size_t M, typename T >
379 const vector< M, T > vector< M, T >::FORWARD( 0, 0, -1 );
380 template< size_t M, typename T >
381 const vector< M, T > vector< M, T >::BACKWARD( 0, 0, 1 );
382 template< size_t M, typename T >
383 const vector< M, T > vector< M, T >::UP( 0, 1, 0 );
384 template< size_t M, typename T >
385 const vector< M, T > vector< M, T >::DOWN( 0, -1, 0 );
386 template< size_t M, typename T >
387 const vector< M, T > vector< M, T >::LEFT( -1, 0, 0 );
388 template< size_t M, typename T >
389 const vector< M, T > vector< M, T >::RIGHT( 1, 0, 0 );
390 template< size_t M, typename T >
391 
392 const vector< M, T > vector< M, T >::ONE( static_cast< T >( 1 ));
393 template< size_t M, typename T >
394 const vector< M, T > vector< M, T >::ZERO( static_cast< T >( 0 ));
395 template< size_t M, typename T >
396 
397 const vector< M, T > vector< M, T >::UNIT_X( 1, 0, 0 );
398 template< size_t M, typename T >
399 const vector< M, T > vector< M, T >::UNIT_Y( 0, 1, 0 );
400 template< size_t M, typename T >
401 const vector< M, T > vector< M, T >::UNIT_Z( 0, 0, 1 );
402 #endif
403 
404 #ifndef VMMLIB_NO_TYPEDEFS
405 # ifdef _MSC_VER
406  typedef UINT8 uint8_t;
407 # endif
408 typedef vmml::vector< 2, int > Vector2i;
409 typedef vmml::vector< 3, int > Vector3i;
410 typedef vmml::vector< 4, int > Vector4i;
411 typedef vmml::vector< 2, unsigned > Vector2ui;
412 typedef vmml::vector< 3, unsigned > Vector3ui;
413 typedef vmml::vector< 4, unsigned > Vector4ui;
414 typedef vmml::vector< 3, double > Vector3d;
415 typedef vmml::vector< 4, double > Vector4d;
416 typedef vmml::vector< 2, float > Vector2f;
417 typedef vmml::vector< 3, float > Vector3f;
418 typedef vmml::vector< 4, float > Vector4f;
419 typedef vmml::vector< 3, uint8_t > Vector3ub;
420 typedef vmml::vector< 4, uint8_t > Vector4ub;
421 #endif
422 
423 //
424 // some free functions for convenience
425 //
426 
427 template< size_t M, typename T >
428 bool equals( const vector< M, T >& a, const vector< M, T >& b )
429 {
430  return a.equals( b );
431 }
432 
433 
434 // allows float * vector, not only vector * float
435 template< size_t M, typename T >
436 static vector< M, T > operator* ( T factor, const vector< M, T >& vector_ )
437 {
438  return vector_ * factor;
439 }
440 
441 
442 template< size_t M, typename T >
443 inline T dot( const vector< M, T >& first, const vector< M, T >& second )
444 {
445  return first.dot( second );
446 }
447 
448 
449 template< size_t M, typename T >
450 inline vector< M, T > cross( const vector< 3, T >& a, const vector< 3, T >& b )
451 {
452  return a.cross( b );
453 }
454 
455 
456 template< size_t M, typename T >
457 inline vector< M, T > normalize( const vector< M, T >& vector_ )
458 {
459  vector< M, T > v( vector_ );
460  v.normalize();
461  return v;
462 }
463 
464 template< typename T >
465 inline vector< 4, T > compute_plane( const vector< 3, T >& a,
466  const vector< 3, T >& b,
467  const vector< 3, T >& c )
468 {
469  const vector< 3, T > cb = b - c;
470  const vector< 3, T > ba = a - b;
471 
472  vector< 4, T > plane = cb.cross( ba );
473  plane.normalize();
474  plane.w() = -plane.x() * a.x() - plane.y() * a.y() - plane.z() * a.z();
475  return plane;
476 }
477 
478 template< size_t M, typename T >
479 vector< M, T >::vector( const T& _a )
480 {
481  for( iterator it = begin(), it_end = end(); it != it_end; ++it )
482  {
483  *it = _a;
484  }
485 }
486 
487 template< size_t M, typename T >
488 vector< M, T >::vector( const T& _x, const T& _y )
489 {
490  array[ 0 ] = _x;
491  array[ 1 ] = _y;
492 }
493 
494 
495 template< size_t M, typename T >
496 vector< M, T >::vector( const T& _x, const T& _y, const T& _z )
497 {
498  array[ 0 ] = _x;
499  array[ 1 ] = _y;
500  array[ 2 ] = _z;
501 }
502 
503 
504 
505 template< size_t M, typename T >
506 vector< M, T >::vector( const T& _x, const T& _y, const T& _z, const T& _w )
507 {
508  array[ 0 ] = _x;
509  array[ 1 ] = _y;
510  array[ 2 ] = _z;
511  array[ 3 ] = _w;
512 }
513 
514 
515 template< size_t M, typename T >
516 vector< M, T >::vector( const T* values )
517 {
518  memcpy( array, values, M * sizeof( T ));
519 }
520 
521 #ifdef __OSG_MATH
522 template< size_t M, typename T >
523 template< typename OSGVEC3 >
524 vector< M, T >::vector( const OSGVEC3& from,
525  typename enable_if< M == 3, OSGVEC3 >::type* )
526 {
527  array[ 0 ] = from.x();
528  array[ 1 ] = from.y();
529  array[ 2 ] = from.z();
530 }
531 #endif
532 
533 #ifndef SWIG
534 template< size_t M, typename T >
535 // initializes the first M-1 values from vector_, the last from last_
536 vector< M, T >::vector( const vector< M-1, T >& vector_, T last_ )
537 {
538  typename vector< M-1, T >::const_iterator
539  it = vector_.begin(), it_end = vector_.end();
540 
541  iterator my_it = begin();
542 
543  for( ; it != it_end; ++it, ++my_it )
544  {
545  (*my_it) = *it;
546  }
547  (*my_it) = last_;
548 }
549 #endif
550 
551 
552 
553 // to-homogenous-coordinates ctor
554 template< size_t M, typename T >
555 template< size_t N >
556 vector< M, T >::vector( const vector< N, T >& source_,
557  typename enable_if< N == M - 1 >::type* )
558 {
559  (*this) = source_;
560 }
561 
562 
563 
564 
565 // from-homogenous-coordinates ctor
566 template< size_t M, typename T >
567 template< size_t N >
568 vector< M, T >::vector( const vector< N, T >& source_,
569  typename enable_if< N == M + 1 >::type* )
570 {
571  (*this) = source_;
572 }
573 
574 
575 template< size_t M, typename T >
576 template< typename U >
577 vector< M, T >::vector( const vector< M, U >& source_ )
578 {
579  (*this) = source_;
580 }
581 
582 
583 
584 template< size_t M, typename T > void vector< M, T >::set( T _a )
585 {
586  for( iterator it = begin(), it_end = end(); it != it_end; ++it )
587  *it = _a;
588 }
589 
590 
591 #ifndef SWIG
592 template< size_t M, typename T >
593 void vector< M, T >::set( const vector< M-1, T >& v, T _a )
594 {
595  memcpy( array, v.array, sizeof( T ) * (M-1) );
596  at( M-1 ) = _a;
597 }
598 #endif
599 
600 template< size_t M, typename T > template< size_t N >
601 void vector< M, T >::set( const vector< N, T >& v )
602 {
603  size_t minimum = M;
604  if (N < M) minimum = N;
605  memcpy( array, v.array, sizeof( T ) * minimum );
606 }
607 
608 template< size_t M, typename T >
609 void vector< M, T >::set( T _x, T _y )
610 {
611  array[ 0 ] = _x;
612  array[ 1 ] = _y;
613 }
614 
615 
616 template< size_t M, typename T >
617 void vector< M, T >::set( T _x, T _y, T _z )
618 {
619  array[ 0 ] = _x;
620  array[ 1 ] = _y;
621  array[ 2 ] = _z;
622 }
623 
624 
625 
626 template< size_t M, typename T >
627 void vector< M, T >::set( T _x, T _y, T _z, T _w )
628 {
629  array[ 0 ] = _x;
630  array[ 1 ] = _y;
631  array[ 2 ] = _z;
632  array[ 3 ] = _w;
633 }
634 
635 
636 template< size_t M, typename T >
637 inline T&
638 vector< M, T >::operator()( size_t index )
639 {
640  return at( index );
641 }
642 
643 
644 
645 template< size_t M, typename T >
646 inline const T&
647 vector< M, T >::operator()( size_t index ) const
648 {
649  return at( index );
650 }
651 
652 
653 
654 template< size_t M, typename T >
655 inline T&
656 vector< M, T >::at( size_t index )
657 {
658  #ifdef VMMLIB_SAFE_ACCESSORS
659  if ( index >= M )
660  {
661  VMMLIB_ERROR( "at() - index out of bounds", VMMLIB_HERE );
662  }
663  #endif
664  return array[ index ];
665 }
666 
667 
668 
669 template< size_t M, typename T >
670 inline const T&
671 vector< M, T >::at( size_t index ) const
672 {
673  #ifdef VMMLIB_SAFE_ACCESSORS
674  if ( index >= M )
675  {
676  VMMLIB_ERROR( "at() - index out of bounds", VMMLIB_HERE );
677  }
678  #endif
679  return array[ index ];
680 }
681 
682 
683 #ifndef VMMLIB_NO_CONVERSION_OPERATORS
684 
685 template< size_t M, typename T >
686 vector< M, T >::operator T*()
687 {
688  return array;
689 }
690 
691 
692 
693 template< size_t M, typename T >
694 vector< M, T >::operator const T*() const
695 {
696  return array;
697 }
698 #else
699 
700 template< size_t M, typename T >
701 T&
702 vector< M, T >::operator[]( size_t index )
703 {
704  return at( index );
705 }
706 
707 template< size_t M, typename T >
708 const T&
709 vector< M, T >::operator[]( size_t index ) const
710 {
711  return at( index );
712 }
713 
714 
715 #endif
716 
717 
718 #if 0
719 template< size_t M, typename T >
720 inline T&
721 vector< M, T >::operator[]( size_t index )
722 {
723  return at( index );
724 }
725 
726 
727 
728 template< size_t M, typename T >
729 inline const T&
730 vector< M, T >::operator[]( size_t index ) const
731 {
732  return at( index );
733 }
734 #endif
735 
736 
737 template< size_t M, typename T >
738 vector< M, T >
739 vector< M, T >::operator*( const vector< M, T >& other ) const
740 {
741  vector< M, T > result;
742  for( size_t index = 0; index < M; ++index )
743  result.at( index ) = at( index ) * other.at( index );
744  return result;
745 }
746 
747 
748 
749 template< size_t M, typename T >
750 vector< M, T >
751 vector< M, T >::operator/( const vector< M, T >& other ) const
752 {
753  vector< M, T > result;
754  for( size_t index = 0; index < M; ++index )
755  result.at( index ) = at( index ) / other.at( index );
756  return result;
757 }
758 
759 
760 
761 template< size_t M, typename T >
762 vector< M, T >
763 vector< M, T >::operator+( const vector< M, T >& other ) const
764 {
765  vector< M, T > result;
766  for( size_t index = 0; index < M; ++index )
767  result.at( index ) = at( index ) + other.at( index );
768  return result;
769 }
770 
771 
772 
773 template< size_t M, typename T >
774 vector< M, T >
775 vector< M, T >::operator-( const vector< M, T >& other ) const
776 {
777  vector< M, T > result;
778  for( size_t index = 0; index < M; ++index )
779  result.at( index ) = at( index ) - other.at( index );
780  return result;
781 }
782 
783 
784 
785 
786 template< size_t M, typename T >
787 void
788 vector< M, T >::operator*=( const vector< M, T >& other )
789 {
790  for( size_t index = 0; index < M; ++index )
791  at( index ) *= other.at( index );
792 }
793 
794 
795 
796 template< size_t M, typename T >
797 void
798 vector< M, T >::operator/=( const vector< M, T >& other )
799 {
800  for( size_t index = 0; index < M; ++index )
801  at( index ) /= other.at( index );
802 }
803 
804 
805 
806 template< size_t M, typename T >
807 void
808 vector< M, T >::operator+=( const vector< M, T >& other )
809 {
810  for( size_t index = 0; index < M; ++index )
811  at( index ) += other.at( index );
812 }
813 
814 
815 
816 template< size_t M, typename T >
817 void
818 vector< M, T >::operator-=( const vector< M, T >& other )
819 {
820  for( size_t index = 0; index < M; ++index )
821  at( index ) -= other.at( index );
822 }
823 
824 
825 
826 template< size_t M, typename T >
827 vector< M, T >
828 vector< M, T >::operator*( const T other ) const
829 {
830  vector< M, T > result;
831  for( size_t index = 0; index < M; ++index )
832  result.at( index ) = at( index ) * other;
833  return result;
834 }
835 
836 
837 
838 template< size_t M, typename T >
839 vector< M, T >
840 vector< M, T >::operator/( const T other ) const
841 {
842  vector< M, T > result;
843  for( size_t index = 0; index < M; ++index )
844  result.at( index ) = at( index ) / other;
845  return result;
846 }
847 
848 
849 
850 template< size_t M, typename T >
851 vector< M, T >
852 vector< M, T >::operator+( const T other ) const
853 {
854  vector< M, T > result;
855  for( size_t index = 0; index < M; ++index )
856  result.at( index ) = at( index ) + other;
857  return result;
858 }
859 
860 
861 
862 template< size_t M, typename T >
863 vector< M, T >
864 vector< M, T >::operator-( const T other ) const
865 {
866  vector< M, T > result;
867  for( size_t index = 0; index < M; ++index )
868  result.at( index ) = at( index ) - other;
869  return result;
870 }
871 
872 
873 
874 
875 template< size_t M, typename T >
876 void
877 vector< M, T >::operator*=( const T other )
878 {
879  for( size_t index = 0; index < M; ++index )
880  at( index ) *= other;
881 }
882 
883 
884 
885 template< size_t M, typename T >
886 void
887 vector< M, T >::operator/=( const T other )
888 {
889  for( size_t index = 0; index < M; ++index )
890  at( index ) /= other;
891 }
892 
893 
894 
895 template< size_t M, typename T >
896 void
897 vector< M, T >::operator+=( const T other )
898 {
899  for( size_t index = 0; index < M; ++index )
900  at( index ) += other;
901 }
902 
903 
904 
905 template< size_t M, typename T >
906 void
907 vector< M, T >::operator-=( const T other )
908 {
909  for( size_t index = 0; index < M; ++index )
910  at( index ) -= other;
911 }
912 
913 
914 
915 template< size_t M, typename T >
916 vector< M, T >
917 vector< M, T >::operator-() const
918 {
919  vector< M, T > v( *this );
920  return v.negate();
921 }
922 
923 
924 
925 template< size_t M, typename T >
926 const vector< M, T >&
927 vector< M, T >::negate()
928 {
929  for( size_t index = 0; index < M; ++index )
930  array[ index ] = -array[ index ];
931  return *this;
932 }
933 
934 
935 
936 template< size_t M, typename T >
937 inline T&
938 vector< M, T >::x()
939 {
940  return array[ 0 ];
941 }
942 
943 
944 
945 template< size_t M, typename T >
946 inline T&
947 vector< M, T >::y()
948 {
949  return array[ 1 ];
950 }
951 
952 
953 
954 template< size_t M, typename T >
955 inline T&
956 vector< M, T >::z()
957 {
958  return array[ 2 ];
959 }
960 
961 
962 
963 template< size_t M, typename T >
964 inline T&
965 vector< M, T >::w()
966 {
967  return array[ 3 ];
968 }
969 
970 
971 
972 template< size_t M, typename T >
973 inline const T&
974 vector< M, T >::x() const
975 {
976  return array[ 0 ];
977 }
978 
979 
980 
981 template< size_t M, typename T >
982 inline const T&
983 vector< M, T >::y() const
984 {
985  return array[ 1 ];
986 }
987 
988 
989 
990 template< size_t M, typename T >
991 inline const T&
992 vector< M, T >::z() const
993 {
994  return array[ 2 ];
995 }
996 
997 
998 
999 template< size_t M, typename T >
1000 inline const T&
1001 vector< M, T >::w() const
1002 {
1003  return array[ 3 ];
1004 }
1005 
1006 
1007 template< size_t M, typename T >
1008 inline T&
1009 vector< M, T >::r()
1010 {
1011  return array[ 0 ];
1012 }
1013 
1014 
1015 
1016 template< size_t M, typename T >
1017 inline T&
1018 vector< M, T >::g()
1019 {
1020  return array[ 1 ];
1021 }
1022 
1023 
1024 
1025 template< size_t M, typename T >
1026 inline T&
1027 vector< M, T >::b()
1028 {
1029  return array[ 2 ];
1030 }
1031 
1032 
1033 
1034 template< size_t M, typename T >
1035 inline T&
1036 vector< M, T >::a()
1037 {
1038  return array[ 3 ];
1039 }
1040 
1041 
1042 
1043 template< size_t M, typename T >
1044 inline const T&
1045 vector< M, T >::r() const
1046 {
1047  return array[ 0 ];
1048 }
1049 
1050 
1051 
1052 template< size_t M, typename T >
1053 inline const T&
1054 vector< M, T >::g() const
1055 {
1056  return array[ 1 ];
1057 }
1058 
1059 
1060 
1061 template< size_t M, typename T >
1062 inline const T&
1063 vector< M, T >::b() const
1064 {
1065  return array[ 2 ];
1066 }
1067 
1068 
1069 
1070 template< size_t M, typename T >
1071 inline const T&
1072 vector< M, T >::a() const
1073 {
1074  return array[ 3 ];
1075 }
1076 
1077 // result = vec1.cross( vec2 ) => result = vec1 x vec2
1078 template< size_t M, typename T >
1079 template< typename TT >
1080 inline vector< M, T > vector< M, T >::cross( const vector< M, TT >& rhs,
1081  typename enable_if< M == 3, TT >::type* ) const
1082 {
1083  vector< M, T > result;
1084  result.cross( *this, rhs );
1085  return result;
1086 }
1087 
1088 
1089 
1090 // result.cross( vec1, vec2 ) => (this) = vec1 x vec2
1091 template< size_t M, typename T >
1092 template< typename TT >
1093 void vector< M, T >::cross( const vector< M, TT >& aa,
1094  const vector< M, TT >& bb,
1095  typename enable_if< M == 3, TT >::type* )
1096 {
1097  array[ 0 ] = aa.y() * bb.z() - aa.z() * bb.y();
1098  array[ 1 ] = aa.z() * bb.x() - aa.x() * bb.z();
1099  array[ 2 ] = aa.x() * bb.y() - aa.y() * bb.x();
1100 }
1101 
1102 
1103 
1104 template< size_t M, typename T >
1105 inline T vector< M, T >::dot( const vector< M, T >& other ) const
1106 {
1107  T tmp = 0.0;
1108  for( size_t index = 0; index < M; ++index )
1109  tmp += at( index ) * other.at( index );
1110 
1111  return tmp;
1112 }
1113 
1114 
1115 template< size_t M, typename T >
1116 inline T vector< M, T >::normalize()
1117 {
1118  const T len = length();
1119 
1120  if ( len == 0 )
1121  return 0;
1122 
1123  const T tmp = 1.0 / len;
1124  (*this) *= tmp;
1125  return len;
1126 }
1127 
1128 template< size_t M, typename T >
1129 inline T vector< M, T >::length() const
1130 {
1131  return std::sqrt( squared_length() );
1132 }
1133 
1134 template< size_t M, typename T >
1135 inline T vector< M, T >::squared_length() const
1136 {
1137  T _squared_length = 0.0;
1138  for( const_iterator it = begin(), it_end = end(); it != it_end; ++it )
1139  _squared_length += (*it) * (*it);
1140 
1141  return _squared_length;
1142 }
1143 
1144 
1145 
1146 template< size_t M, typename T >
1147 inline T
1148 vector< M, T >::distance( const vector< M, T >& other ) const
1149 {
1150  return std::sqrt( squared_distance( other ) );
1151 }
1152 
1153 
1154 
1155 template< size_t M, typename T >
1156 inline T vector< M, T >::squared_distance( const vector< M, T >& other ) const
1157 {
1158  vector< M, T > tmp( *this );
1159  tmp -= other;
1160  return tmp.squared_length();
1161 }
1162 
1163 template< size_t M, typename T > inline T vector< M, T >::product() const
1164 {
1165  T result = at( 0 );
1166  for( size_t i = 1; i < M; ++i )
1167  result *= at( i );
1168  return result;
1169 }
1170 
1171 template< size_t M, typename T >
1172 void vector< M, T >::compute_normal( const vector< M, T >& aa,
1173  const vector< M, T >& bb,
1174  const vector< M, T >& cc )
1175 {
1176  vector< M, T > u,v;
1177  // right hand system, CCW triangle
1178  u = bb - aa;
1179  v = cc - aa;
1180  cross( u, v );
1181  normalize();
1182 }
1183 
1184 
1185 
1186 template< size_t M, typename T >
1187 vector< M, T >
1188 vector< M, T >::compute_normal( const vector< M, T >& bb,
1189  const vector< M, T >& cc ) const
1190 {
1191  vector< M, T > tmp;
1192  tmp.compute_normal( *this, bb, cc);
1193  return tmp;
1194 }
1195 
1196 template< size_t M, typename T >
1197 template< typename TT >
1198 vector< 3, T > vector< M, T >::rotate( const T theta, vector< M, TT > axis,
1199  typename enable_if< M == 3, TT >::type* ) const
1200 {
1201  axis.normalize();
1202 
1203  const T costheta = std::cos( theta );
1204  const T sintheta = std::sin( theta );
1205 
1206  return vector< 3, T >(
1207  (costheta + ( 1.0f - costheta ) * axis.x() * axis.x() ) * x() +
1208  (( 1 - costheta ) * axis.x() * axis.y() - axis.z() * sintheta ) * y() +
1209  (( 1 - costheta ) * axis.x() * axis.z() + axis.y() * sintheta ) * z(),
1210 
1211  (( 1 - costheta ) * axis.x() * axis.y() + axis.z() * sintheta ) * x() +
1212  ( costheta + ( 1 - costheta ) * axis.y() * axis.y() ) * y() +
1213  (( 1 - costheta ) * axis.y() * axis.z() - axis.x() * sintheta ) * z(),
1214 
1215  (( 1 - costheta ) * axis.x() * axis.z() - axis.y() * sintheta ) * x() +
1216  (( 1 - costheta ) * axis.y() * axis.z() + axis.x() * sintheta ) * y() +
1217  ( costheta + ( 1 - costheta ) * axis.z() * axis.z() ) * z() );
1218 }
1219 
1220 
1221 // sphere layout: center xyz, radius w
1222 template< size_t M, typename T >
1223 template< typename TT >
1224 inline vector< 3, T >
1225 vector< M, T >::
1226 project_point_onto_sphere( const vector< 3, TT >& point,
1227  typename enable_if< M == 4, TT >::type* ) const
1228 {
1229  const vector< 3, T >& _center = get_sub_vector< 3 >( 0 );
1230 
1231  vector< 3, T > projected_point( point );
1232  projected_point -= _center;
1233  projected_point.normalize();
1234  projected_point *= w();
1235  return _center + projected_point;
1236 }
1237 
1238 
1239 
1240 // sphere layout: center xyz, radius w
1241 template< size_t M, typename T >
1242 template< typename TT >
1243 inline T
1244 vector< M, T >::
1245 distance_to_sphere( const vector< 3, TT >& point,
1246  typename enable_if< M == 4, TT >::type* ) const
1247 {
1248  const vector< 3, T >& center_ = get_sub_vector< 3 >( 0 );
1249  return ( point - center_ ).length() - w();
1250 }
1251 
1252 template< size_t M, typename T > template< size_t N > inline
1253 vector< N, T >& vector< M, T >::get_sub_vector( size_t offset,
1254  typename enable_if< M >= N >::type* )
1255 {
1256  assert( offset <= M - N );
1257  return reinterpret_cast< vector< N, T >& >( *( begin() + offset ) );
1258 }
1259 
1260 template< size_t M, typename T > template< size_t N > inline
1261 const vector< N, T >& vector< M, T >::get_sub_vector( size_t offset,
1262  typename enable_if< M >= N >::type* ) const
1263 {
1264  assert( offset <= M - N );
1265  return reinterpret_cast< const vector< N, T >& >( *( begin() + offset ) );
1266 }
1267 
1268 
1269 // plane: normal xyz, distance w
1270 template< size_t M, typename T > template< typename TT >
1271 inline T vector< M, T >::distance_to_plane( const vector< 3, TT >& point,
1272  typename enable_if< M == 4, TT >::type* ) const
1273 {
1274  const vector< 3, T >& normal = get_sub_vector< 3 >( 0 );
1275  return normal.dot( point ) + w();
1276 }
1277 
1278 
1279 
1280 // plane: normal xyz, distance w
1281 template< size_t M, typename T >
1282 template< typename TT >
1283 vector< 3, T >
1284 vector< M, T >::project_point_onto_plane( const vector< 3, TT >& point,
1285  typename enable_if< M == 4, TT >::type* ) const
1286 {
1287  const vector< 3, T >& normal = get_sub_vector< 3 >( 0 );
1288  return point - ( normal * distance_to_plane( point ) );
1289 }
1290 
1291 
1292 
1293 template< size_t M, typename T >
1294 bool vector< M, T >::operator==( const vector< M, T >& other ) const
1295 {
1296  return memcmp( array, other.array, sizeof( array )) == 0;
1297 }
1298 
1299 
1300 template< size_t M, typename T >
1301 bool vector< M, T >::operator!=( const vector< M, T >& other ) const
1302 {
1303  return ! this->operator==( other );
1304 }
1305 
1306 
1307 template< size_t M, typename T >
1308 bool
1309 vector< M, T >::equals( const vector< M, T >& other, T tolerance ) const
1310 {
1311  for( size_t index = 0; index < M; ++index )
1312  if( fabs( at( index ) - other( index ) ) >= tolerance )
1313  return false;
1314  return true;
1315 
1316 }
1317 
1318 
1319 template< size_t M, typename T >
1320 bool
1321 vector< M, T >::operator<( const vector< M, T >& other ) const
1322 {
1323  for(size_t index = 0; index < M; ++index )
1324  {
1325  if (at( index ) < other.at( index )) return true;
1326  if (other.at( index ) < at( index )) return false;
1327  }
1328  return false;
1329 }
1330 
1331 
1332 // to-homogenous-coordinates assignment operator
1333 // non-chainable because of sfinae
1334 template< size_t M, typename T > template< size_t N >
1335 typename enable_if< N == M - 1 >::type*
1336 vector< M, T >::operator=( const vector< N, T >& source_ )
1337 {
1338  std::copy( source_.begin(), source_.end(), begin() );
1339  at( M - 1 ) = static_cast< T >( 1.0 );
1340  return 0;
1341 }
1342 
1343 
1344 // from-homogenous-coordinates assignment operator
1345 // non-chainable because of sfinae
1346 template< size_t M, typename T > template< size_t N >
1347 typename enable_if< N == M + 1 >::type*
1348 vector< M, T >::operator=( const vector< N, T >& source_ )
1349 {
1350  const T w_reci = static_cast< T >( 1.0 ) / source_( M );
1351  iterator it = begin(), it_end = end();
1352  for( size_t index = 0; it != it_end; ++it, ++index )
1353  *it = source_( index ) * w_reci;
1354  return 0;
1355 }
1356 
1357 
1358 template< size_t M, typename T >
1359 vector< M, T >& vector< M, T >::operator=( const T* c_array )
1360 {
1361  iter_set( c_array, c_array + M );
1362  return *this;
1363 }
1364 
1365 
1366 
1367 template< size_t M, typename T >
1368 T vector< M, T >::operator=( T filler_value )
1369 {
1370  for( size_t index = 0; index < M; ++index )
1371  at( index ) = filler_value;
1372  return filler_value;
1373 }
1374 
1375 
1376 
1377 
1378 template< size_t M, typename T >
1379 vector< M, T >& vector< M, T >::operator=( const vector< M, T >& other )
1380 {
1381  if( this != &other )
1382  memcpy( array, other.array, M * sizeof( T ) );
1383  return *this;
1384 }
1385 
1386 
1387 
1388 // returns void to avoid 'silent' loss of precision when chaining
1389 template< size_t M, typename T > template< typename U >
1390 void vector< M, T >::operator=( const vector< M, U >& source_ )
1391 {
1392  typedef typename vector< M, U >::const_iterator u_c_iter;
1393  u_c_iter it = source_.begin(), it_end = source_.end();
1394  for( iterator my_it = begin(); it != it_end; ++it, ++my_it )
1395  *my_it = static_cast< T >( *it );
1396 }
1397 
1398 
1399 
1400 template< size_t M, typename T >
1401 template< typename input_iterator_t >
1402 void
1403 vector< M, T >::iter_set( input_iterator_t begin_, input_iterator_t end_ )
1404 {
1405  input_iterator_t in_it = begin_;
1406  iterator it = begin(), it_end = end();
1407  for( ; it != it_end && in_it != end_; ++it, ++in_it )
1408  (*it) = static_cast< T >( *in_it );
1409 }
1410 
1411 template< size_t M, typename T >
1412 void vector< M, T >::clamp( const T& min, const T& max )
1413 {
1414  for( size_t i = 0; i < M; ++i )
1415  {
1416  if( array[i] < min )
1417  array[i] = min;
1418  if( array[i] > max )
1419  array[i] = max;
1420  }
1421 }
1422 
1423 
1424 
1425 template< size_t M, typename T >
1426 template< typename TT >
1427 void
1428 vector< M, T >::scale_to( vector< M, TT >& result_,
1429  T min_value, T max_value ) const
1430 {
1431  T range = max_value-min_value;
1432  T half_range = range * 0.5;
1433  T scale = ( 1.0 / range ) * static_cast< T >( std::numeric_limits< TT >::max() );
1434 
1435  for( size_t index = 0; index < M; ++index )
1436  {
1437  result_.at( index )
1438  = static_cast< TT >( ( at( index ) + half_range ) * scale );
1439  }
1440 
1441 }
1442 
1443 
1444 
1445 template< size_t M, typename T >
1446 inline size_t
1447 vector< M, T >::size()
1448 {
1449  return M;
1450 }
1451 
1452 
1453 
1454 template< size_t M, typename T >
1455 size_t
1456 vector< M, T >::find_min_index() const
1457 {
1458  return std::min_element( begin(), end() ) - begin();
1459 }
1460 
1461 
1462 
1463 template< size_t M, typename T >
1464 size_t
1465 vector< M, T >::find_max_index() const
1466 {
1467  return std::max_element( begin(), end() ) - begin();
1468 }
1469 
1470 
1471 
1472 template< size_t M, typename T >
1473 size_t
1474 vector< M, T >::find_abs_min_index() const
1475 {
1476  return std::min_element( begin(), end(), vmml::math::abs_less< T >() ) - begin();
1477 }
1478 
1479 
1480 
1481 template< size_t M, typename T >
1482 size_t
1483 vector< M, T >::find_abs_max_index() const
1484 {
1485  return std::max_element( begin(), end(), vmml::math::abs_greater< T >() ) - begin();
1486 }
1487 
1488 
1489 
1490 template< size_t M, typename T >
1491 T&
1492 vector< M, T >::find_min()
1493 {
1494  return *std::min_element( begin(), end() );
1495 }
1496 
1497 
1498 
1499 template< size_t M, typename T >
1500 const T&
1501 vector< M, T >::find_min() const
1502 {
1503  return *std::min_element( begin(), end() );
1504 }
1505 
1506 
1507 
1508 template< size_t M, typename T >
1509 T&
1510 vector< M, T >::find_max()
1511 {
1512  return *std::max_element( begin(), end() );
1513 }
1514 
1515 
1516 
1517 template< size_t M, typename T >
1518 const T&
1519 vector< M, T >::find_max() const
1520 {
1521  return *std::max_element( begin(), end() );
1522 }
1523 
1524 
1525 template< size_t M, typename T >
1526 inline typename vector< M, T >::iterator
1527 vector< M, T >::begin()
1528 {
1529  return array;
1530 }
1531 
1532 
1533 template< size_t M, typename T >
1534 inline typename vector< M, T >::iterator
1535 vector< M, T >::end()
1536 {
1537  return array + M; ;
1538 }
1539 
1540 
1541 template< size_t M, typename T >
1542 inline typename vector< M, T >::const_iterator
1543 vector< M, T >::begin() const
1544 {
1545  return array;
1546 }
1547 
1548 
1549 template< size_t M, typename T >
1550 inline typename vector< M, T >::const_iterator
1551 vector< M, T >::end() const
1552 {
1553  return array + M; ;
1554 }
1555 
1556 
1557 
1558 template< size_t M, typename T >
1559 inline typename vector< M, T >::reverse_iterator
1560 vector< M, T >::rbegin()
1561 {
1562  return array + M - 1;
1563 }
1564 
1565 
1566 template< size_t M, typename T >
1567 inline typename vector< M, T >::reverse_iterator
1568 vector< M, T >::rend()
1569 {
1570  return array - 1;
1571 }
1572 
1573 
1574 template< size_t M, typename T >
1575 inline typename vector< M, T >::const_reverse_iterator
1576 vector< M, T >::rbegin() const
1577 {
1578  return array + M - 1;
1579 }
1580 
1581 
1582 template< size_t M, typename T >
1583 inline typename vector< M, T >::const_reverse_iterator
1584 vector< M, T >::rend() const
1585 {
1586  return array - 1;
1587 }
1588 
1589 
1590 
1591 template< size_t M, typename T >
1592 bool
1593 vector< M, T >::is_unit_vector() const
1594 {
1595  const_iterator it = begin(), it_end = end();
1596  bool one = false;
1597  for( ; it != it_end; ++it )
1598  {
1599  if ( *it == 1.0 )
1600  {
1601  if ( one )
1602  return false;
1603  one = true;
1604  }
1605  else if ( *it != 0.0 )
1606  {
1607  return false;
1608  }
1609  }
1610  return one;
1611 }
1612 
1613 
1614 
1615 
1616 template< size_t M, typename T >
1617 void
1618 vector< M, T >::perturb( T perturbation )
1619 {
1620  for( iterator it = begin(), it_end = end(); it != it_end; ++it )
1621  {
1622  (*it) += ( rand() & 1u ) ? perturbation : -perturbation;
1623  }
1624 
1625 }
1626 
1627 template< size_t M, typename T >
1628 void
1629 vector< M, T >::sqrt_elementwise()
1630 {
1631  for( iterator it = begin(), it_end = end(); it != it_end; ++it )
1632  {
1633  (*it) = std::sqrt(*it);
1634  }
1635 }
1636 
1637 
1638 
1639 template< size_t M, typename T >
1640 void
1641 vector< M, T >::reciprocal()
1642 {
1643  for( iterator it = begin(), it_end = end(); it != it_end; ++it )
1644  {
1645  (*it) = static_cast< T >( 1.0 ) / (*it);
1646  }
1647 }
1648 
1649 
1650 
1651 template< size_t M, typename T >
1652 void
1653 vector< M, T >::reciprocal_safe()
1654 {
1655  for( iterator it = begin(), it_end = end(); it != it_end; ++it )
1656  {
1657  T& v = *it;
1658 
1659  if ( v == static_cast< T >( 0 ) )
1660  v = std::numeric_limits< T >::max();
1661  else
1662  v = static_cast< T >( 1.0 ) / v;
1663  }
1664 }
1665 
1666 
1667 
1668 template< size_t M, typename T >
1669 template< typename TT >
1670 void
1671 vector< M, T >::cast_from( const vector< M, TT >& other )
1672 {
1673  typedef vmml::vector< M, TT > vector_tt_type ;
1674  typedef typename vector_tt_type::const_iterator tt_const_iterator;
1675 
1676  iterator it = begin(), it_end = end();
1677  tt_const_iterator other_it = other.begin();
1678  for( ; it != it_end; ++it, ++other_it )
1679  {
1680  *it = static_cast< T >( *other_it );
1681  }
1682 }
1683 
1684 template< size_t M, typename T >
1685 size_t
1686 vector< M, T >::nnz() const
1687 {
1688  size_t counter = 0;
1689 
1690  const_iterator it = begin(),
1691  it_end = end();
1692  for( ; it != it_end; ++it)
1693  {
1694  if ( *it != 0 ) {
1695  ++counter;
1696  }
1697  }
1698 
1699  return counter;
1700 }
1701 
1702 
1703 template< size_t M, typename T >
1704 double
1705 vector< M, T >::norm( ) const
1706 {
1707  double norm_v = 0.0;
1708 
1709  const_iterator it = begin(), it_end = end();
1710  for( ; it != it_end; ++it )
1711  {
1712  norm_v += *it * *it;
1713  }
1714 
1715  return std::sqrt(norm_v);
1716 }
1717 
1718 template< size_t M, typename T >
1719 void
1720 vector< M, T >::set_random( int seed )
1721 {
1722  if ( seed >= 0 )
1723  srand( seed );
1724 
1725  for( size_t i = 0; i < M; ++i )
1726  {
1727  const double fillValue = double( rand( )) / double( RAND_MAX );
1728  at( i ) = -1.0 + 2.0 * fillValue;
1729  }
1730 }
1731 
1732 
1733 } // namespace vmml
1734 
1735 #endif
vmml::enable_if
Definition: enable_if.hpp:46
vmml::math::abs_greater
Definition: math.hpp:92
vmml::vector::array
T array[M]
storage
Definition: vector.hpp:352
vmml::vector::get_sub_vector
vector< N, T > & get_sub_vector(size_t offset=0, typename enable_if< M >=N >::type *=0)
Definition: vector.hpp:1253
vmml::math::abs_less
Definition: math.hpp:82
vmml::vector::product
T product() const
Definition: vector.hpp:1163