vmmlib  1.7.0
 All Classes Namespaces Functions Pages
vector.hpp
1 /*
2  * Copyright (c) 2006-2014, 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_vector_ ) const;
245  inline T squared_distance( const vector& other_vector_ ) const;
246 
247  template< typename TT >
248  vector< 3, T > rotate( const T theta, vector< M, TT > axis,
249  typename enable_if< M == 3, TT >::type* = 0 ) const;
250 
251  // right hand system, CCW triangle
252  // (*this) = normal of v0,v1,v2
253  void compute_normal( const vector& v0, const vector& v1, const vector& v2 );
254  // retval = normal of (this), v1, v2
255  vector compute_normal( const vector& v1, const vector& v2 ) const;
256 
257  template< size_t N >
258  void get_sub_vector( vector< N, T >& sub_v_, size_t offset = 0,
259  typename enable_if< M >= N >::type* = 0 );
260 
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 
265  template< size_t N >
266  const vector< N, T >& get_sub_vector( size_t offset = 0,
267  typename enable_if< M >= N >::type* = 0 ) const;
268 
269  // sphere functions - sphere layout: center xyz, radius w
270  template< typename TT >
271  inline vector< 3, T > project_point_onto_sphere(
272  const vector< 3, TT >& point,
273  typename enable_if< M == 4, TT >::type* = 0 ) const;
274 
275  // returns a negative distance if the point lies in the sphere
276  template< typename TT >
277  inline T distance_to_sphere( const vector< 3, TT >& point,
278  typename enable_if< M == 4, TT >::type* = 0 ) const;
279 
280  // plane functions - plane layout; normal xyz, distance w
281  template< typename TT >
282  inline T distance_to_plane( const vector< 3, TT >& point,
283  typename enable_if< M == 4, TT >::type* = 0 ) const;
284 
285  template< typename TT >
286  inline vector< 3, T > project_point_onto_plane(
287  const vector< 3, TT >& point,
288  typename enable_if< M == 4, TT >::type* = 0 ) const;
289 
290  // returns the index of the minimal resp. maximal value in the vector
291  size_t find_min_index() const;
292  size_t find_max_index() const;
293 
294  // returns the index of the minimal resp. maximal value in the vector
295  size_t find_abs_min_index() const;
296  size_t find_abs_max_index() const;
297 
298  // returns minimal resp. maximal value in the vector
299  T& find_min();
300  T& find_max();
301  const T& find_min() const;
302  const T& find_max() const;
303 
304  void clamp( const T& min = 0.0, const T& max = 1.0 );
305 
306  template< typename TT >
307  void scale_to( vector< M, TT >& scaled_vector,
308  T min_value = -1.0, T max_value = 1.0 ) const;
309 
310  inline static size_t size(); // returns M
311 
312  bool is_unit_vector() const;
313 
314  // perturbs each component by randomly + or - the perturbation parameter
315  void perturb( T perturbation = 0.0001 );
316 
317  void sqrt_elementwise();
318  double norm() const; //l2 norm
319 
320  // computes the reciprocal value for each component, x = 1/x;
321  // WARNING: might result in nans if division by 0!
322  void reciprocal();
323  // computes the reciprocal value for each component, x = 1/x;
324  // checks every component for 0, sets to max value if zero.
325  void reciprocal_safe();
326 
327  template< typename TT >
328  void cast_from( const vector< M, TT >& other );
329 
330  size_t nnz() const;
331 
332  // test each component of the vector for isnan and isinf
333  // inline bool is_valid() const; -> moved to class validator
334 
335  friend std::ostream& operator<< ( std::ostream& os, const vector& vector_ )
336  {
337 #ifdef EQFABRIC_API_H
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 #else
350  os << "(";
351  size_t index = 0;
352  for( ; index < M - 1; ++index )
353  {
354  os << vector_.at( index ) << ", ";
355  }
356  os << vector_.at( index ) << ")";
357 #endif
358  return os;
359  }
360 
361 
362 
363  // storage
364  VMMLIB_ALIGN( T array[ M ] );
365 
366 #ifndef SWIG
367  // Vector3 defaults
368  static const vector FORWARD;
369  static const vector BACKWARD;
370  static const vector UP;
371  static const vector DOWN;
372  static const vector LEFT;
373  static const vector RIGHT;
374 
375  static const vector ONE;
376  static const vector ZERO;
377 
378  // Unit vectors
379  static const vector UNIT_X;
380  static const vector UNIT_Y;
381  static const vector UNIT_Z;
382 #endif
383 
384 }; // class vector
385 
386 //
387 // typedefs and statics
388 //
389 
390 #ifndef VMMLIB_NO_TYPEDEFS
391 typedef vector< 2, float > vec2f;
392 typedef vector< 2, double > vec2d;
393 typedef vector< 3, float > vec3f;
394 typedef vector< 3, double > vec3d;
395 typedef vector< 4, float > vec4f;
396 typedef vector< 4, double > vec4d;
397 #endif
398 
399 #ifndef SWIG
400 template< size_t M, typename T >
401 const vector< M, T > vector< M, T >::FORWARD( 0, 0, -1 );
402 template< size_t M, typename T >
403 const vector< M, T > vector< M, T >::BACKWARD( 0, 0, 1 );
404 template< size_t M, typename T >
405 const vector< M, T > vector< M, T >::UP( 0, 1, 0 );
406 template< size_t M, typename T >
407 const vector< M, T > vector< M, T >::DOWN( 0, -1, 0 );
408 template< size_t M, typename T >
409 const vector< M, T > vector< M, T >::LEFT( -1, 0, 0 );
410 template< size_t M, typename T >
411 const vector< M, T > vector< M, T >::RIGHT( 1, 0, 0 );
412 template< size_t M, typename T >
413 
414 const vector< M, T > vector< M, T >::ONE( static_cast< T >( 1 ));
415 template< size_t M, typename T >
416 const vector< M, T > vector< M, T >::ZERO( static_cast< T >( 0 ));
417 template< size_t M, typename T >
418 
419 const vector< M, T > vector< M, T >::UNIT_X( 1, 0, 0 );
420 template< size_t M, typename T >
421 const vector< M, T > vector< M, T >::UNIT_Y( 0, 1, 0 );
422 template< size_t M, typename T >
423 const vector< M, T > vector< M, T >::UNIT_Z( 0, 0, 1 );
424 #endif
425 
426 //
427 // some free functions for convenience
428 //
429 
430 template< size_t M, typename T >
431 bool equals( const vector< M, T >& a, const vector< M, T >& b )
432 {
433  return a.equals( b );
434 }
435 
436 
437 // allows float * vector, not only vector * float
438 template< size_t M, typename T >
439 static vector< M, T >
440 operator* ( T factor, const vector< M, T >& vector_ )
441 {
442  return vector_ * factor;
443 }
444 
445 
446 template< size_t M, typename T >
447 inline T dot( const vector< M, T >& first, const vector< M, T >& second )
448 {
449  return first.dot( second );
450 }
451 
452 
453 template< size_t M, typename T >
454 inline vector< M, T > cross( const vector< 3, T >& a, const vector< 3, T >& b )
455 {
456  return a.cross( b );
457 }
458 
459 
460 template< size_t M, typename T >
461 inline vector< M, T > normalize( const vector< M, T >& vector_ )
462 {
463  vector< M, T > v( vector_ );
464  v.normalize();
465  return v;
466 }
467 
468 template< typename T >
469 inline vector< 4, T > compute_plane( const vector< 3, T >& a,
470  const vector< 3, T >& b,
471  const vector< 3, T >& c )
472 {
473  const vector< 3, T > cb = b - c;
474  const vector< 3, T > ba = a - b;
475 
476  vector< 4, T > plane = cb.cross( ba );
477  plane.normalize();
478  plane.w() = -plane.x() * a.x() - plane.y() * a.y() - plane.z() * a.z();
479  return plane;
480 }
481 
482 template< size_t M, typename T >
483 vector< M, T >::vector( const T& _a )
484 {
485  for( iterator it = begin(), it_end = end(); it != it_end; ++it )
486  {
487  *it = _a;
488  }
489 }
490 
491 template< size_t M, typename T >
492 vector< M, T >::vector( const T& _x, const T& _y )
493 {
494  array[ 0 ] = _x;
495  array[ 1 ] = _y;
496 }
497 
498 
499 template< size_t M, typename T >
500 vector< M, T >::vector( const T& _x, const T& _y, const T& _z )
501 {
502  array[ 0 ] = _x;
503  array[ 1 ] = _y;
504  array[ 2 ] = _z;
505 }
506 
507 
508 
509 template< size_t M, typename T >
510 vector< M, T >::vector( const T& _x, const T& _y, const T& _z, const T& _w )
511 {
512  array[ 0 ] = _x;
513  array[ 1 ] = _y;
514  array[ 2 ] = _z;
515  array[ 3 ] = _w;
516 }
517 
518 
519 template< size_t M, typename T >
520 vector< M, T >::vector( const T* values )
521 {
522  memcpy( array, values, M * sizeof( T ));
523 }
524 
525 #ifdef __OSG_MATH
526 template< size_t M, typename T >
527 template< typename OSGVEC3 >
528 vector< M, T >::vector( const OSGVEC3& from,
529  typename enable_if< M == 3, OSGVEC3 >::type* )
530 {
531  array[ 0 ] = from.x();
532  array[ 1 ] = from.y();
533  array[ 2 ] = from.z();
534 }
535 #endif
536 
537 #ifndef SWIG
538 template< size_t M, typename T >
539 // initializes the first M-1 values from vector_, the last from last_
540 vector< M, T >::vector( const vector< M-1, T >& vector_, T last_ )
541 {
542  typename vector< M-1, T >::const_iterator
543  it = vector_.begin(), it_end = vector_.end();
544 
545  iterator my_it = begin();
546 
547  for( ; it != it_end; ++it, ++my_it )
548  {
549  (*my_it) = *it;
550  }
551  (*my_it) = last_;
552 }
553 #endif
554 
555 
556 
557 // to-homogenous-coordinates ctor
558 template< size_t M, typename T >
559 template< size_t N >
560 vector< M, T >::vector( const vector< N, T >& source_,
561  typename enable_if< N == M - 1 >::type* )
562 {
563  (*this) = source_;
564 }
565 
566 
567 
568 
569 // from-homogenous-coordinates ctor
570 template< size_t M, typename T >
571 template< size_t N >
572 vector< M, T >::vector( const vector< N, T >& source_,
573  typename enable_if< N == M + 1 >::type* )
574 {
575  (*this) = source_;
576 }
577 
578 
579 template< size_t M, typename T >
580 template< typename U >
581 vector< M, T >::vector( const vector< M, U >& source_ )
582 {
583  (*this) = source_;
584 }
585 
586 
587 
588 template< size_t M, typename T > void vector< M, T >::set( T _a )
589 {
590  for( iterator it = begin(), it_end = end(); it != it_end; ++it )
591  {
592  *it = _a;
593  }
594 }
595 
596 
597 #ifndef SWIG
598 template< size_t M, typename T >
599 void vector< M, T >::set( const vector< M-1, T >& v, T _a )
600 {
601  memcpy( array, v.array, sizeof( T ) * (M-1) );
602  at( M-1 ) = _a;
603 }
604 #endif
605 
606 template< size_t M, typename T > template< size_t N >
607 void vector< M, T >::set( const vector< N, T >& v )
608 {
609  size_t minimum = M;
610  if (N < M) minimum = N;
611  memcpy( array, v.array, sizeof( T ) * minimum );
612 }
613 
614 template< size_t M, typename T >
615 void vector< M, T >::set( T _x, T _y )
616 {
617  array[ 0 ] = _x;
618  array[ 1 ] = _y;
619 }
620 
621 
622 template< size_t M, typename T >
623 void vector< M, T >::set( T _x, T _y, T _z )
624 {
625  array[ 0 ] = _x;
626  array[ 1 ] = _y;
627  array[ 2 ] = _z;
628 }
629 
630 
631 
632 template< size_t M, typename T >
633 void vector< M, T >::set( T _x, T _y, T _z, T _w )
634 {
635  array[ 0 ] = _x;
636  array[ 1 ] = _y;
637  array[ 2 ] = _z;
638  array[ 3 ] = _w;
639 }
640 
641 
642 template< size_t M, typename T >
643 inline T&
644 vector< M, T >::operator()( size_t index )
645 {
646  return at( index );
647 }
648 
649 
650 
651 template< size_t M, typename T >
652 inline const T&
653 vector< M, T >::operator()( size_t index ) const
654 {
655  return at( index );
656 }
657 
658 
659 
660 template< size_t M, typename T >
661 inline T&
662 vector< M, T >::at( size_t index )
663 {
664  #ifdef VMMLIB_SAFE_ACCESSORS
665  if ( index >= M )
666  {
667  VMMLIB_ERROR( "at() - index out of bounds", VMMLIB_HERE );
668  }
669  #endif
670  return array[ index ];
671 }
672 
673 
674 
675 template< size_t M, typename T >
676 inline const T&
677 vector< M, T >::at( size_t index ) const
678 {
679  #ifdef VMMLIB_SAFE_ACCESSORS
680  if ( index >= M )
681  {
682  VMMLIB_ERROR( "at() - index out of bounds", VMMLIB_HERE );
683  }
684  #endif
685  return array[ index ];
686 }
687 
688 
689 #ifndef VMMLIB_NO_CONVERSION_OPERATORS
690 
691 template< size_t M, typename T >
692 vector< M, T >::operator T*()
693 {
694  return array;
695 }
696 
697 
698 
699 template< size_t M, typename T >
700 vector< M, T >::operator const T*() const
701 {
702  return array;
703 }
704 #else
705 
706 template< size_t M, typename T >
707 T&
708 vector< M, T >::operator[]( size_t index )
709 {
710  return at( index );
711 }
712 
713 template< size_t M, typename T >
714 const T&
715 vector< M, T >::operator[]( size_t index ) const
716 {
717  return at( index );
718 }
719 
720 
721 #endif
722 
723 
724 #if 0
725 template< size_t M, typename T >
726 inline T&
727 vector< M, T >::operator[]( size_t index )
728 {
729  return at( index );
730 }
731 
732 
733 
734 template< size_t M, typename T >
735 inline const T&
736 vector< M, T >::operator[]( size_t index ) const
737 {
738  return at( index );
739 }
740 #endif
741 
742 
743 template< size_t M, typename T >
744 vector< M, T >
745 vector< M, T >::operator*( const vector< M, T >& other ) const
746 {
747  vector< M, T > result;
748  for( size_t index = 0; index < M; ++index )
749  {
750  result.at( index ) = at( index ) * other.at( index );
751  }
752  return result;
753 }
754 
755 
756 
757 template< size_t M, typename T >
758 vector< M, T >
759 vector< M, T >::operator/( const vector< M, T >& other ) const
760 {
761  vector< M, T > result;
762  for( size_t index = 0; index < M; ++index )
763  {
764  result.at( index ) = at( index ) / other.at( index );
765  }
766  return result;
767 }
768 
769 
770 
771 template< size_t M, typename T >
772 vector< M, T >
773 vector< M, T >::operator+( const vector< M, T >& other ) const
774 {
775  vector< M, T > result;
776  for( size_t index = 0; index < M; ++index )
777  {
778  result.at( index ) = at( index ) + other.at( index );
779  }
780  return result;
781 }
782 
783 
784 
785 template< size_t M, typename T >
786 vector< M, T >
787 vector< M, T >::operator-( const vector< M, T >& other ) const
788 {
789  vector< M, T > result;
790  for( size_t index = 0; index < M; ++index )
791  {
792  result.at( index ) = at( index ) - other.at( index );
793  }
794  return result;
795 }
796 
797 
798 
799 
800 template< size_t M, typename T >
801 void
802 vector< M, T >::operator*=( const vector< M, T >& other )
803 {
804  for( size_t index = 0; index < M; ++index )
805  {
806  at( index ) *= other.at( index );
807  }
808 }
809 
810 
811 
812 template< size_t M, typename T >
813 void
814 vector< M, T >::operator/=( const vector< M, T >& other )
815 {
816  for( size_t index = 0; index < M; ++index )
817  {
818  at( index ) /= other.at( index );
819  }
820 }
821 
822 
823 
824 template< size_t M, typename T >
825 void
826 vector< M, T >::operator+=( const vector< M, T >& other )
827 {
828  for( size_t index = 0; index < M; ++index )
829  {
830  at( index ) += other.at( index );
831  }
832 }
833 
834 
835 
836 template< size_t M, typename T >
837 void
838 vector< M, T >::operator-=( const vector< M, T >& other )
839 {
840  for( size_t index = 0; index < M; ++index )
841  {
842  at( index ) -= other.at( index );
843  }
844 }
845 
846 
847 
848 template< size_t M, typename T >
849 vector< M, T >
850 vector< M, T >::operator*( const T other ) const
851 {
852  vector< M, T > result;
853  for( size_t index = 0; index < M; ++index )
854  {
855  result.at( index ) = at( index ) * other;
856  }
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  {
869  result.at( index ) = at( index ) / other;
870  }
871  return result;
872 }
873 
874 
875 
876 template< size_t M, typename T >
877 vector< M, T >
878 vector< M, T >::operator+( const T other ) const
879 {
880  vector< M, T > result;
881  for( size_t index = 0; index < M; ++index )
882  {
883  result.at( index ) = at( index ) + other;
884  }
885  return result;
886 }
887 
888 
889 
890 template< size_t M, typename T >
891 vector< M, T >
892 vector< M, T >::operator-( const T other ) const
893 {
894  vector< M, T > result;
895  for( size_t index = 0; index < M; ++index )
896  {
897  result.at( index ) = at( index ) - other;
898  }
899  return result;
900 }
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  {
911  at( index ) *= other;
912  }
913 }
914 
915 
916 
917 template< size_t M, typename T >
918 void
919 vector< M, T >::operator/=( const T other )
920 {
921  for( size_t index = 0; index < M; ++index )
922  {
923  at( index ) /= other;
924  }
925 }
926 
927 
928 
929 template< size_t M, typename T >
930 void
931 vector< M, T >::operator+=( const T other )
932 {
933  for( size_t index = 0; index < M; ++index )
934  {
935  at( index ) += other;
936  }
937 }
938 
939 
940 
941 template< size_t M, typename T >
942 void
943 vector< M, T >::operator-=( const T other )
944 {
945  for( size_t index = 0; index < M; ++index )
946  {
947  at( index ) -= other;
948  }
949 }
950 
951 
952 
953 template< size_t M, typename T >
954 vector< M, T >
955 vector< M, T >::operator-() const
956 {
957  vector< M, T > v( *this );
958  return v.negate();
959 }
960 
961 
962 
963 template< size_t M, typename T >
964 const vector< M, T >&
965 vector< M, T >::negate()
966 {
967  for( size_t index = 0; index < M; ++index )
968  {
969  array[ index ] = -array[ index ];
970  }
971  return *this;
972 }
973 
974 
975 
976 template< size_t M, typename T >
977 inline T&
978 vector< M, T >::x()
979 {
980  return array[ 0 ];
981 }
982 
983 
984 
985 template< size_t M, typename T >
986 inline T&
987 vector< M, T >::y()
988 {
989  return array[ 1 ];
990 }
991 
992 
993 
994 template< size_t M, typename T >
995 inline T&
996 vector< M, T >::z()
997 {
998  return array[ 2 ];
999 }
1000 
1001 
1002 
1003 template< size_t M, typename T >
1004 inline T&
1005 vector< M, T >::w()
1006 {
1007  return array[ 3 ];
1008 }
1009 
1010 
1011 
1012 template< size_t M, typename T >
1013 inline const T&
1014 vector< M, T >::x() const
1015 {
1016  return array[ 0 ];
1017 }
1018 
1019 
1020 
1021 template< size_t M, typename T >
1022 inline const T&
1023 vector< M, T >::y() const
1024 {
1025  return array[ 1 ];
1026 }
1027 
1028 
1029 
1030 template< size_t M, typename T >
1031 inline const T&
1032 vector< M, T >::z() const
1033 {
1034  return array[ 2 ];
1035 }
1036 
1037 
1038 
1039 template< size_t M, typename T >
1040 inline const T&
1041 vector< M, T >::w() const
1042 {
1043  return array[ 3 ];
1044 }
1045 
1046 
1047 template< size_t M, typename T >
1048 inline T&
1049 vector< M, T >::r()
1050 {
1051  return array[ 0 ];
1052 }
1053 
1054 
1055 
1056 template< size_t M, typename T >
1057 inline T&
1058 vector< M, T >::g()
1059 {
1060  return array[ 1 ];
1061 }
1062 
1063 
1064 
1065 template< size_t M, typename T >
1066 inline T&
1067 vector< M, T >::b()
1068 {
1069  return array[ 2 ];
1070 }
1071 
1072 
1073 
1074 template< size_t M, typename T >
1075 inline T&
1076 vector< M, T >::a()
1077 {
1078  return array[ 3 ];
1079 }
1080 
1081 
1082 
1083 template< size_t M, typename T >
1084 inline const T&
1085 vector< M, T >::r() const
1086 {
1087  return array[ 0 ];
1088 }
1089 
1090 
1091 
1092 template< size_t M, typename T >
1093 inline const T&
1094 vector< M, T >::g() const
1095 {
1096  return array[ 1 ];
1097 }
1098 
1099 
1100 
1101 template< size_t M, typename T >
1102 inline const T&
1103 vector< M, T >::b() const
1104 {
1105  return array[ 2 ];
1106 }
1107 
1108 
1109 
1110 template< size_t M, typename T >
1111 inline const T&
1112 vector< M, T >::a() const
1113 {
1114  return array[ 3 ];
1115 }
1116 
1117 // result = vec1.cross( vec2 ) => result = vec1 x vec2
1118 template< size_t M, typename T >
1119 template< typename TT >
1120 inline vector< M, T > vector< M, T >::cross( const vector< M, TT >& rhs,
1121  typename enable_if< M == 3, TT >::type* ) const
1122 {
1123  vector< M, T > result;
1124  result.cross( *this, rhs );
1125  return result;
1126 }
1127 
1128 
1129 
1130 // result.cross( vec1, vec2 ) => (this) = vec1 x vec2
1131 template< size_t M, typename T >
1132 template< typename TT >
1133 void vector< M, T >::cross( const vector< M, TT >& aa,
1134  const vector< M, TT >& bb,
1135  typename enable_if< M == 3, TT >::type* )
1136 {
1137  array[ 0 ] = aa.y() * bb.z() - aa.z() * bb.y();
1138  array[ 1 ] = aa.z() * bb.x() - aa.x() * bb.z();
1139  array[ 2 ] = aa.x() * bb.y() - aa.y() * bb.x();
1140 }
1141 
1142 
1143 
1144 template< size_t M, typename T >
1145 inline T vector< M, T >::dot( const vector< M, T >& other ) const
1146 {
1147  T tmp = 0.0;
1148  for( size_t index = 0; index < M; ++index )
1149  tmp += at( index ) * other.at( index );
1150 
1151  return tmp;
1152 }
1153 
1154 
1155 template< size_t M, typename T >
1156 inline T vector< M, T >::normalize()
1157 {
1158  const T len = length();
1159 
1160  if ( len == 0 )
1161  return 0;
1162 
1163  const T tmp = 1.0 / len;
1164  (*this) *= tmp;
1165  return len;
1166 }
1167 
1168 template< size_t M, typename T >
1169 inline T vector< M, T >::length() const
1170 {
1171  return std::sqrt( squared_length() );
1172 }
1173 
1174 template< size_t M, typename T >
1175 inline T vector< M, T >::squared_length() const
1176 {
1177  T _squared_length = 0.0;
1178  for( const_iterator it = begin(), it_end = end(); it != it_end; ++it )
1179  _squared_length += (*it) * (*it);
1180 
1181  return _squared_length;
1182 }
1183 
1184 
1185 
1186 template< size_t M, typename T >
1187 inline T
1188 vector< M, T >::distance( const vector< M, T >& other_vector_ ) const
1189 {
1190  return std::sqrt( squared_distance( other_vector_ ) );
1191 }
1192 
1193 
1194 
1195 template< size_t M, typename T >
1196 inline T
1197 vector< M, T >::squared_distance( const vector< M, T >& other_vector_ ) const
1198 {
1199  vector< M, T > tmp( *this );
1200  tmp -= other_vector_;
1201  return tmp.squared_length();
1202 }
1203 
1204 
1205 
1206 template< size_t M, typename T >
1207 void
1208 vector< M, T >::compute_normal(
1209  const vector< M, T >& aa,
1210  const vector< M, T >& bb,
1211  const vector< M, T >& cc
1212  )
1213 {
1214  vector< M, T > u,v;
1215  // right hand system, CCW triangle
1216  u = bb - aa;
1217  v = cc - aa;
1218  cross( u, v );
1219  normalize();
1220 }
1221 
1222 
1223 
1224 template< size_t M, typename T >
1225 vector< M, T >
1226 vector< M, T >::compute_normal( const vector< M, T >& bb,
1227  const vector< M, T >& cc ) const
1228 {
1229  vector< M, T > tmp;
1230  tmp.compute_normal( *this, bb, cc);
1231  return tmp;
1232 }
1233 
1234 template< size_t M, typename T >
1235 template< typename TT >
1236 vector< 3, T > vector< M, T >::rotate( const T theta, vector< M, TT > axis,
1237  typename enable_if< M == 3, TT >::type* ) const
1238 {
1239  axis.normalize();
1240 
1241  const T costheta = std::cos( theta );
1242  const T sintheta = std::sin( theta );
1243 
1244  return vector< 3, T >(
1245  (costheta + ( 1.0f - costheta ) * axis.x() * axis.x() ) * x() +
1246  (( 1 - costheta ) * axis.x() * axis.y() - axis.z() * sintheta ) * y() +
1247  (( 1 - costheta ) * axis.x() * axis.z() + axis.y() * sintheta ) * z(),
1248 
1249  (( 1 - costheta ) * axis.x() * axis.y() + axis.z() * sintheta ) * x() +
1250  ( costheta + ( 1 - costheta ) * axis.y() * axis.y() ) * y() +
1251  (( 1 - costheta ) * axis.y() * axis.z() - axis.x() * sintheta ) * z(),
1252 
1253  (( 1 - costheta ) * axis.x() * axis.z() - axis.y() * sintheta ) * x() +
1254  (( 1 - costheta ) * axis.y() * axis.z() + axis.x() * sintheta ) * y() +
1255  ( costheta + ( 1 - costheta ) * axis.z() * axis.z() ) * z() );
1256 }
1257 
1258 
1259 // sphere layout: center xyz, radius w
1260 template< size_t M, typename T >
1261 template< typename TT >
1262 inline vector< 3, T >
1263 vector< M, T >::
1264 project_point_onto_sphere( const vector< 3, TT >& point,
1265  typename enable_if< M == 4, TT >::type* ) const
1266 {
1267  const vector< 3, T >& _center = get_sub_vector< 3 >( 0 );
1268 
1269  vector< 3, T > projected_point( point );
1270  projected_point -= _center;
1271  projected_point.normalize();
1272  projected_point *= w();
1273  return _center + projected_point;
1274 }
1275 
1276 
1277 
1278 // sphere layout: center xyz, radius w
1279 template< size_t M, typename T >
1280 template< typename TT >
1281 inline T
1282 vector< M, T >::
1283 distance_to_sphere( const vector< 3, TT >& point,
1284  typename enable_if< M == 4, TT >::type* ) const
1285 {
1286  const vector< 3, T >& center_ = get_sub_vector< 3 >( 0 );
1287  return ( point - center_ ).length() - w();
1288 }
1289 
1290 template< size_t M, typename T >
1291 template< size_t N >
1292 void
1293 vector< M, T >::get_sub_vector( vector< N, T >& sub_v, size_t offset,
1294  typename enable_if< M >= N >::type* )
1295 {
1296  assert( offset <= M - N );
1297  sub_v = reinterpret_cast< vector< N, T >& >( *( begin() + offset ) );
1298 }
1299 
1300 
1301 
1302 template< size_t M, typename T >
1303 template< size_t N >
1304 inline vector< N, T >&
1305 vector< M, T >::get_sub_vector( size_t offset,
1306  typename enable_if< M >= N >::type* )
1307 {
1308  assert( offset <= M - N );
1309  return reinterpret_cast< vector< N, T >& >( *( begin() + offset ) );
1310 }
1311 
1312 
1313 
1314 template< size_t M, typename T >
1315 template< size_t N >
1316 inline const vector< N, T >&
1317 vector< M, T >::get_sub_vector( size_t offset,
1318  typename enable_if< M >= N >::type* ) const
1319 {
1320  assert( offset <= M - N );
1321  return reinterpret_cast< const vector< N, T >& >( *( begin() + offset ) );
1322 }
1323 
1324 
1325 // plane: normal xyz, distance w
1326 template< size_t M, typename T > template< typename TT >
1327 inline T vector< M, T >::distance_to_plane( const vector< 3, TT >& point,
1328  typename enable_if< M == 4, TT >::type* ) const
1329 {
1330  const vector< 3, T >& normal = get_sub_vector< 3 >( 0 );
1331  return normal.dot( point ) + w();
1332 }
1333 
1334 
1335 
1336 // plane: normal xyz, distance w
1337 template< size_t M, typename T >
1338 template< typename TT >
1339 vector< 3, T >
1340 vector< M, T >::project_point_onto_plane( const vector< 3, TT >& point,
1341  typename enable_if< M == 4, TT >::type* ) const
1342 {
1343  const vector< 3, T >& normal = get_sub_vector< 3 >( 0 );
1344  return point - ( normal * distance_to_plane( point ) );
1345 }
1346 
1347 
1348 
1349 template< size_t M, typename T >
1350 bool vector< M, T >::operator==( const vector< M, T >& other ) const
1351 {
1352  return memcmp( array, other.array, sizeof( array )) == 0;
1353 }
1354 
1355 
1356 template< size_t M, typename T >
1357 bool vector< M, T >::operator!=( const vector< M, T >& other ) const
1358 {
1359  return ! this->operator==( other );
1360 }
1361 
1362 
1363 template< size_t M, typename T >
1364 bool
1365 vector< M, T >::
1366 equals( const vector< M, T >& other, T tolerance ) const
1367 {
1368  bool is_ok = true;
1369  for( size_t index = 0; is_ok && index < M; ++index )
1370  {
1371  is_ok = fabs( at( index ) - other( index ) ) < tolerance;
1372  }
1373  return is_ok;
1374 
1375 }
1376 
1377 
1378 template< size_t M, typename T >
1379 bool
1380 vector< M, T >::operator<( const vector< M, T >& other ) const
1381 {
1382  for(size_t index = 0; index < M; ++index )
1383  {
1384  if (at( index ) < other.at( index )) return true;
1385  if (other.at( index ) < at( index )) return false;
1386  }
1387  return false;
1388 }
1389 
1390 
1391 // to-homogenous-coordinates assignment operator
1392 // non-chainable because of sfinae
1393 template< size_t M, typename T > template< size_t N >
1394 typename enable_if< N == M - 1 >::type*
1395 vector< M, T >::operator=( const vector< N, T >& source_ )
1396 {
1397  std::copy( source_.begin(), source_.end(), begin() );
1398  at( M - 1 ) = static_cast< T >( 1.0 );
1399  return 0;
1400 }
1401 
1402 
1403 // from-homogenous-coordinates assignment operator
1404 // non-chainable because of sfinae
1405 template< size_t M, typename T > template< size_t N >
1406 typename enable_if< N == M + 1 >::type*
1407 vector< M, T >::operator=( const vector< N, T >& source_ )
1408 {
1409  const T w_reci = static_cast< T >( 1.0 ) / source_( M );
1410  iterator it = begin(), it_end = end();
1411  for( size_t index = 0; it != it_end; ++it, ++index )
1412  {
1413  *it = source_( index ) * w_reci;
1414  }
1415  return 0;
1416 }
1417 
1418 
1419 template< size_t M, typename T >
1420 vector< M, T >& vector< M, T >::operator=( const T* c_array )
1421 {
1422  iter_set( c_array, c_array + M );
1423  return *this;
1424 }
1425 
1426 
1427 
1428 template< size_t M, typename T >
1429 T vector< M, T >::operator=( T filler_value )
1430 {
1431  for( size_t index = 0; index < M; ++index )
1432  {
1433  at( index ) = filler_value;
1434  }
1435  return filler_value;
1436 }
1437 
1438 
1439 
1440 
1441 template< size_t M, typename T >
1442 vector< M, T >& vector< M, T >::operator=( const vector< M, T >& other )
1443 {
1444  if( this != &other )
1445  memcpy( array, other.array, M * sizeof( T ) );
1446  return *this;
1447 }
1448 
1449 
1450 
1451 // returns void to avoid 'silent' loss of precision when chaining
1452 template< size_t M, typename T > template< typename U >
1453 void vector< M, T >::operator=( const vector< M, U >& source_ )
1454 {
1455  typedef typename vector< M, U >::const_iterator u_c_iter;
1456  u_c_iter it = source_.begin(), it_end = source_.end();
1457  for( iterator my_it = begin(); it != it_end; ++it, ++my_it )
1458  {
1459  *my_it = static_cast< T >( *it );
1460  }
1461 }
1462 
1463 
1464 
1465 template< size_t M, typename T >
1466 template< typename input_iterator_t >
1467 void
1468 vector< M, T >::iter_set( input_iterator_t begin_, input_iterator_t end_ )
1469 {
1470  input_iterator_t in_it = begin_;
1471  iterator it = begin(), it_end = end();
1472  for( ; it != it_end && in_it != end_; ++it, ++in_it )
1473  {
1474  (*it) = static_cast< T >( *in_it );
1475  }
1476 
1477 }
1478 
1479 template< size_t M, typename T >
1480 void vector< M, T >::clamp( const T& min, const T& max )
1481 {
1482  for( size_t i = 0; i < M; ++i )
1483  {
1484  if( array[i] < min )
1485  array[i] = min;
1486  if( array[i] > max )
1487  array[i] = max;
1488  }
1489 }
1490 
1491 
1492 
1493 template< size_t M, typename T >
1494 template< typename TT >
1495 void
1496 vector< M, T >::scale_to( vector< M, TT >& result_,
1497  T min_value, T max_value ) const
1498 {
1499  T range = max_value-min_value;
1500  T half_range = range * 0.5;
1501  T scale = ( 1.0 / range ) * static_cast< T >( std::numeric_limits< TT >::max() );
1502 
1503  for( size_t index = 0; index < M; ++index )
1504  {
1505  result_.at( index )
1506  = static_cast< TT >( ( at( index ) + half_range ) * scale );
1507  }
1508 
1509 }
1510 
1511 
1512 
1513 template< size_t M, typename T >
1514 inline size_t
1515 vector< M, T >::size()
1516 {
1517  return M;
1518 }
1519 
1520 
1521 
1522 template< size_t M, typename T >
1523 size_t
1524 vector< M, T >::find_min_index() const
1525 {
1526  return std::min_element( begin(), end() ) - begin();
1527 }
1528 
1529 
1530 
1531 template< size_t M, typename T >
1532 size_t
1533 vector< M, T >::find_max_index() const
1534 {
1535  return std::max_element( begin(), end() ) - begin();
1536 }
1537 
1538 
1539 
1540 template< size_t M, typename T >
1541 size_t
1542 vector< M, T >::find_abs_min_index() const
1543 {
1544  return std::min_element( begin(), end(), vmml::math::abs_less< T >() ) - begin();
1545 }
1546 
1547 
1548 
1549 template< size_t M, typename T >
1550 size_t
1551 vector< M, T >::find_abs_max_index() const
1552 {
1553  return std::max_element( begin(), end(), vmml::math::abs_greater< T >() ) - begin();
1554 }
1555 
1556 
1557 
1558 template< size_t M, typename T >
1559 T&
1560 vector< M, T >::find_min()
1561 {
1562  return *std::min_element( begin(), end() );
1563 }
1564 
1565 
1566 
1567 template< size_t M, typename T >
1568 const T&
1569 vector< M, T >::find_min() const
1570 {
1571  return *std::min_element( begin(), end() );
1572 }
1573 
1574 
1575 
1576 template< size_t M, typename T >
1577 T&
1578 vector< M, T >::find_max()
1579 {
1580  return *std::max_element( begin(), end() );
1581 }
1582 
1583 
1584 
1585 template< size_t M, typename T >
1586 const T&
1587 vector< M, T >::find_max() const
1588 {
1589  return *std::max_element( begin(), end() );
1590 }
1591 
1592 
1593 template< size_t M, typename T >
1594 inline typename vector< M, T >::iterator
1595 vector< M, T >::begin()
1596 {
1597  return array;
1598 }
1599 
1600 
1601 template< size_t M, typename T >
1602 inline typename vector< M, T >::iterator
1603 vector< M, T >::end()
1604 {
1605  return array + M; ;
1606 }
1607 
1608 
1609 template< size_t M, typename T >
1610 inline typename vector< M, T >::const_iterator
1611 vector< M, T >::begin() const
1612 {
1613  return array;
1614 }
1615 
1616 
1617 template< size_t M, typename T >
1618 inline typename vector< M, T >::const_iterator
1619 vector< M, T >::end() const
1620 {
1621  return array + M; ;
1622 }
1623 
1624 
1625 
1626 template< size_t M, typename T >
1627 inline typename vector< M, T >::reverse_iterator
1628 vector< M, T >::rbegin()
1629 {
1630  return array + M - 1;
1631 }
1632 
1633 
1634 template< size_t M, typename T >
1635 inline typename vector< M, T >::reverse_iterator
1636 vector< M, T >::rend()
1637 {
1638  return array - 1;
1639 }
1640 
1641 
1642 template< size_t M, typename T >
1643 inline typename vector< M, T >::const_reverse_iterator
1644 vector< M, T >::rbegin() const
1645 {
1646  return array + M - 1;
1647 }
1648 
1649 
1650 template< size_t M, typename T >
1651 inline typename vector< M, T >::const_reverse_iterator
1652 vector< M, T >::rend() const
1653 {
1654  return array - 1;
1655 }
1656 
1657 
1658 
1659 template< size_t M, typename T >
1660 bool
1661 vector< M, T >::is_unit_vector() const
1662 {
1663  const_iterator it = begin(), it_end = end();
1664  bool one = false;
1665  for( ; it != it_end; ++it )
1666  {
1667  if ( *it == 1.0 )
1668  {
1669  if ( one )
1670  return false;
1671  one = true;
1672  }
1673  else if ( *it != 0.0 )
1674  {
1675  return false;
1676  }
1677  }
1678  return one;
1679 }
1680 
1681 
1682 
1683 
1684 template< size_t M, typename T >
1685 void
1686 vector< M, T >::perturb( T perturbation )
1687 {
1688  for( iterator it = begin(), it_end = end(); it != it_end; ++it )
1689  {
1690  (*it) += ( rand() & 1u ) ? perturbation : -perturbation;
1691  }
1692 
1693 }
1694 
1695 template< size_t M, typename T >
1696 void
1697 vector< M, T >::sqrt_elementwise()
1698 {
1699  for( iterator it = begin(), it_end = end(); it != it_end; ++it )
1700  {
1701  (*it) = std::sqrt(*it);
1702  }
1703 }
1704 
1705 
1706 
1707 template< size_t M, typename T >
1708 void
1709 vector< M, T >::reciprocal()
1710 {
1711  for( iterator it = begin(), it_end = end(); it != it_end; ++it )
1712  {
1713  (*it) = static_cast< T >( 1.0 ) / (*it);
1714  }
1715 }
1716 
1717 
1718 
1719 template< size_t M, typename T >
1720 void
1721 vector< M, T >::reciprocal_safe()
1722 {
1723  for( iterator it = begin(), it_end = end(); it != it_end; ++it )
1724  {
1725  T& v = *it;
1726 
1727  if ( v == static_cast< T >( 0 ) )
1728  v = std::numeric_limits< T >::max();
1729  else
1730  v = static_cast< T >( 1.0 ) / v;
1731  }
1732 }
1733 
1734 
1735 
1736 template< size_t M, typename T >
1737 template< typename TT >
1738 void
1739 vector< M, T >::cast_from( const vector< M, TT >& other )
1740 {
1741  typedef vmml::vector< M, TT > vector_tt_type ;
1742  typedef typename vector_tt_type::const_iterator tt_const_iterator;
1743 
1744  iterator it = begin(), it_end = end();
1745  tt_const_iterator other_it = other.begin();
1746  for( ; it != it_end; ++it, ++other_it )
1747  {
1748  *it = static_cast< T >( *other_it );
1749  }
1750 }
1751 
1752 template< size_t M, typename T >
1753 size_t
1754 vector< M, T >::nnz() const
1755 {
1756  size_t counter = 0;
1757 
1758  const_iterator it = begin(),
1759  it_end = end();
1760  for( ; it != it_end; ++it)
1761  {
1762  if ( *it != 0 ) {
1763  ++counter;
1764  }
1765  }
1766 
1767  return counter;
1768 }
1769 
1770 
1771 template< size_t M, typename T >
1772 double
1773 vector< M, T >::norm( ) const
1774 {
1775  double norm_v = 0.0;
1776 
1777  const_iterator it = begin(), it_end = end();
1778  for( ; it != it_end; ++it )
1779  {
1780  norm_v += *it * *it;
1781  }
1782 
1783  return std::sqrt(norm_v);
1784 }
1785 
1786 template< size_t M, typename T >
1787 void
1788 vector< M, T >::set_random( int seed )
1789 {
1790  if ( seed >= 0 )
1791  srand( seed );
1792 
1793  for( size_t i = 0; i < M; ++i )
1794  {
1795  const double fillValue = double( rand( )) / double( RAND_MAX );
1796  at( i ) = -1.0 + 2.0 * fillValue;
1797  }
1798 }
1799 
1800 
1801 } // namespace vmml
1802 
1803 #endif
vmml::enable_if
Definition: enable_if.hpp:46
vmml::math::abs_greater
Definition: math.hpp:92
vmml::math::abs_less
Definition: math.hpp:82