Line data Source code
1 : /*
2 : * Copyright (c) 2006-2016, 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/enable_if.hpp>
37 :
38 : #include <algorithm>
39 : #include <cmath>
40 : #include <cstring>
41 : #include <iomanip>
42 : #include <iostream>
43 : #include <limits>
44 : #include <stdexcept>
45 : #include <stdint.h>
46 : #include <string>
47 : #include <vector>
48 :
49 : namespace vmml
50 : {
51 : template< size_t M, typename T > class vector
52 : {
53 : public:
54 : typedef T value_type;
55 : typedef T* iterator;
56 : typedef const T* const_iterator;
57 : typedef std::reverse_iterator< iterator > reverse_iterator;
58 : typedef std::reverse_iterator< const_iterator > const_reverse_iterator;
59 :
60 : static const size_t DIMENSION = M;
61 :
62 : // constructors
63 101 : vector() : array() {} // http://stackoverflow.com/questions/5602030
64 : explicit vector( const T& a ); // sets all components to a;
65 : vector( const T& x, const T& y );
66 : vector( const T& x, const T& y, const T& z );
67 : vector( const T& x, const T& y, const T& z, const T& w );
68 :
69 : #ifndef SWIG
70 : // initializes the first M-1 values from vector_, the last from last_
71 : vector( const vector< M-1, T >& vector_, T last_ );
72 : #endif
73 :
74 : explicit vector( const T* values );
75 :
76 : #ifdef __OSG_MATH
77 : template< typename OSGVEC3 >
78 : explicit vector( const OSGVEC3& from,
79 : typename enable_if< M == 3, OSGVEC3 >::type* = 0 );
80 : #endif
81 :
82 : // vec< M > with homogeneous coordinates <-> vec< M-1 > conversion ctor
83 : // to-homogenous-coordinates ctor
84 : template< size_t N >
85 : vector( const vector< N, T >& source_,
86 : typename enable_if< N == M - 1 >::type* = 0 );
87 :
88 : // from-homogenous-coordinates vector
89 : template< size_t N >
90 : vector( const vector< N, T >& source_,
91 : typename enable_if< N == M + 1 >::type* = 0 );
92 :
93 : template< typename U > vector( const vector< M, U >& source_ );
94 :
95 : // iterators
96 : inline iterator begin();
97 : inline iterator end();
98 : inline const_iterator begin() const;
99 : inline const_iterator end() const;
100 : inline reverse_iterator rbegin();
101 : inline reverse_iterator rend();
102 : inline const_reverse_iterator rbegin() const;
103 : inline const_reverse_iterator rend() const;
104 :
105 : # ifndef VMMLIB_NO_CONVERSION_OPERATORS
106 : // conversion operators
107 : inline operator T*();
108 : inline operator const T*() const;
109 : # else
110 : inline T& operator[]( size_t index );
111 : inline const T& operator[]( size_t index ) const;
112 : # endif
113 :
114 : // accessors
115 : inline T& operator()( size_t index );
116 : inline const T& operator()( size_t index ) const;
117 :
118 : inline T& at( size_t index );
119 : inline const T& at( size_t index ) const;
120 :
121 : // element accessors for M <= 4;
122 : inline T& x();
123 : inline T& y();
124 : inline T& z();
125 : inline T& w();
126 : inline const T& x() const;
127 : inline const T& y() const;
128 : inline const T& z() const;
129 : inline const T& w() const;
130 :
131 : // pixel color element accessors for M<= 4
132 : inline T& r();
133 : inline T& g();
134 : inline T& b();
135 : inline T& a();
136 : inline const T& r() const;
137 : inline const T& g() const;
138 : inline const T& b() const;
139 : inline const T& a() const;
140 :
141 : bool operator==( const vector& other ) const;
142 : bool operator!=( const vector& other ) const;
143 : bool equals( const vector& other,
144 : T tolerance = std::numeric_limits< T >::epsilon( )) const;
145 : bool operator<( const vector& other ) const;
146 :
147 : // remember kids: c_arrays are dangerous and evil!
148 : vector& operator=( const T* c_array );
149 : T operator=( T filler );
150 :
151 : vector& operator=( const vector& other );
152 :
153 : // returns void to avoid 'silent' loss of precision when chaining
154 : template< typename U > void operator=( const vector< M, U >& other );
155 :
156 : // to-homogenous-coordinates assignment operator
157 : // non-chainable because of sfinae
158 : template< size_t N >
159 : typename enable_if< N == M - 1 >::type*
160 : operator=( const vector< N, T >& source_ );
161 :
162 : // from-homogenous-coordinates assignment operator
163 : // non-chainable because of sfinae
164 : template< size_t N >
165 : typename enable_if< N == M + 1 >::type*
166 : operator=( const vector< N, T >& source_ );
167 :
168 : vector operator*( const vector& other ) const;
169 : vector operator/( const vector& other ) const;
170 : vector operator+( const vector& other ) const;
171 : vector operator-( const vector& other ) const;
172 :
173 : void operator*=( const vector& other );
174 : void operator/=( const vector& other );
175 : void operator+=( const vector& other );
176 : void operator-=( const vector& other );
177 :
178 : vector operator*( const T other ) const;
179 : vector operator/( const T other ) const;
180 : vector operator+( const T other ) const;
181 : vector operator-( const T other ) const;
182 :
183 : void operator*=( const T other );
184 : void operator/=( const T other );
185 : void operator+=( const T other );
186 : void operator-=( const T other );
187 :
188 : vector operator-() const;
189 :
190 : const vector& negate();
191 :
192 : void set( T a ); // sets all components to a;
193 : #ifndef SWIG
194 : void set( const vector< M-1, T >& v, T a );
195 : #endif
196 : template< size_t N >
197 : void set( const vector< N, T >& v );
198 :
199 : // sets the first few components to a certain value
200 : void set( T x, T y );
201 : void set( T x, T y, T z );
202 : void set( T x, T y, T z, T w );
203 :
204 : template< typename input_iterator_t >
205 : void iter_set( input_iterator_t begin_, input_iterator_t end_ );
206 :
207 : // compute the cross product of two vectors
208 : // note: there's also a free function:
209 : // vector<> cross( const vector<>, const vector<> )
210 : template< typename TT >
211 : vector< M, T >& cross( const vector< M, TT >& b,
212 : typename enable_if< M == 3, TT >::type* = 0 );
213 :
214 : // compute the dot product of two vectors
215 : // note: there's also a free function:
216 : // T dot( const vector<>, const vector<> );
217 : inline T dot( const vector& other ) const;
218 :
219 : // normalize the vector
220 : // note: there's also a free function:
221 : // vector<> normalize( const vector<> );
222 : inline T normalize();
223 :
224 : //sets all vector components to random values
225 : //remember to set srand( seed );
226 : //if seed is set to -1, srand( seed ) was set outside set_random
227 : //otherwise srand( seed ) will be called with the given seed
228 : void set_random( int seed = -1 );
229 :
230 : inline T length() const;
231 : inline T squared_length() const;
232 :
233 : inline T distance( const vector& other ) const;
234 : inline T squared_distance( const vector& other ) const;
235 :
236 : /** @return the product of all elements of this vector */
237 : T product() const;
238 :
239 : template< typename TT >
240 : vector< 3, T >& rotate( T theta, vector< M, TT > axis,
241 : typename enable_if< M == 3, TT >::type* = 0 );
242 :
243 : /** @return the sub vector of the given length at the given offset. */
244 : template< size_t N, size_t O >
245 : vector< N, T > get_sub_vector( typename enable_if< M >= N+O >::type* = 0 )
246 : const;
247 :
248 : /** Set the sub vector of the given length at the given offset. */
249 : template< size_t N, size_t O >
250 : void set_sub_vector( const vector< N, T >& sub,
251 : typename enable_if< M >= N+O >::type* = 0 );
252 :
253 : // sphere functions - sphere layout: center xyz, radius w
254 : template< typename TT >
255 : inline vector< 3, T > project_point_onto_sphere(
256 : const vector< 3, TT >& point,
257 : typename enable_if< M == 4, TT >::type* = 0 ) const;
258 :
259 : // returns a negative distance if the point lies in the sphere
260 : template< typename TT >
261 : inline T distance_to_sphere( const vector< 3, TT >& point,
262 : typename enable_if< M == 4, TT >::type* = 0 ) const;
263 :
264 : // plane functions - plane layout; normal xyz, distance w
265 : template< typename TT >
266 : inline T distance_to_plane( const vector< 3, TT >& point,
267 : typename enable_if< M == 4, TT >::type* = 0 ) const;
268 :
269 : template< typename TT >
270 : inline vector< 3, T > project_point_onto_plane(
271 : const vector< 3, TT >& point,
272 : typename enable_if< M == 4, TT >::type* = 0 ) const;
273 :
274 : // returns the index of the minimal resp. maximal value in the vector
275 : size_t find_min_index() const;
276 : size_t find_max_index() const;
277 :
278 : // returns minimal resp. maximal value in the vector
279 : T& find_min();
280 : T& find_max();
281 : const T& find_min() const;
282 : const T& find_max() const;
283 :
284 : void clamp( const T& min = 0.0, const T& max = 1.0 );
285 :
286 : inline static size_t size(); // returns M
287 :
288 : bool is_unit_vector() const;
289 :
290 : // perturbs each component by randomly + or - the perturbation parameter
291 : void perturb( T perturbation = 0.0001 );
292 :
293 : void sqrt_elementwise();
294 : double norm() const; //l2 norm
295 :
296 : // computes the reciprocal value for each component, x = 1/x;
297 : // WARNING: might result in nans if division by 0!
298 : void reciprocal();
299 : // computes the reciprocal value for each component, x = 1/x;
300 : // checks every component for 0, sets to max value if zero.
301 : void reciprocal_safe();
302 :
303 : template< typename TT >
304 : void cast_from( const vector< M, TT >& other );
305 :
306 : size_t nnz() const;
307 :
308 : // test each component of the vector for isnan and isinf
309 : // inline bool is_valid() const; -> moved to class validator
310 :
311 0 : friend std::ostream& operator<< ( std::ostream& os, const vector& vector_ )
312 : {
313 0 : const std::ios::fmtflags flags = os.flags();
314 0 : const int prec = os.precision();
315 :
316 0 : os.setf( std::ios::right, std::ios::adjustfield );
317 0 : os.precision( 5 );
318 0 : os << "[ ";
319 0 : for( size_t index = 0; index < M; ++index )
320 0 : os << std::setw(10) << vector_.at( index ) << " ";
321 0 : os << "]";
322 0 : os.precision( prec );
323 0 : os.setf( flags );
324 0 : return os;
325 : }
326 :
327 : T array[ M ]; //!< storage
328 :
329 : #ifndef SWIG
330 : // Vector3 defaults
331 : static const vector FORWARD;
332 : static const vector BACKWARD;
333 : static const vector UP;
334 : static const vector DOWN;
335 : static const vector LEFT;
336 : static const vector RIGHT;
337 :
338 : static const vector ONE;
339 : static const vector ZERO;
340 :
341 : // Unit vectors
342 : static const vector UNIT_X;
343 : static const vector UNIT_Y;
344 : static const vector UNIT_Z;
345 : #endif
346 :
347 : }; // class vector
348 :
349 : //
350 : // typedefs and statics
351 : //
352 : #ifndef SWIG
353 : template< size_t M, typename T >
354 2 : const vector< M, T > vector< M, T >::FORWARD( 0, 0, -1 );
355 : template< size_t M, typename T >
356 1 : const vector< M, T > vector< M, T >::BACKWARD( 0, 0, 1 );
357 : template< size_t M, typename T >
358 2 : const vector< M, T > vector< M, T >::UP( 0, 1, 0 );
359 : template< size_t M, typename T >
360 : const vector< M, T > vector< M, T >::DOWN( 0, -1, 0 );
361 : template< size_t M, typename T >
362 1 : const vector< M, T > vector< M, T >::LEFT( -1, 0, 0 );
363 : template< size_t M, typename T >
364 : const vector< M, T > vector< M, T >::RIGHT( 1, 0, 0 );
365 : template< size_t M, typename T >
366 :
367 : const vector< M, T > vector< M, T >::ONE( static_cast< T >( 1 ));
368 : template< size_t M, typename T >
369 1 : const vector< M, T > vector< M, T >::ZERO( static_cast< T >( 0 ));
370 : template< size_t M, typename T >
371 :
372 : const vector< M, T > vector< M, T >::UNIT_X( 1, 0, 0 );
373 : template< size_t M, typename T >
374 : const vector< M, T > vector< M, T >::UNIT_Y( 0, 1, 0 );
375 : template< size_t M, typename T >
376 : const vector< M, T > vector< M, T >::UNIT_Z( 0, 0, 1 );
377 : #endif
378 :
379 : //
380 : // some free functions for convenience
381 : //
382 :
383 : template< size_t M, typename T >
384 : bool equals( const vector< M, T >& a, const vector< M, T >& b )
385 : {
386 : return a.equals( b );
387 : }
388 :
389 : // allows float * vector, not only vector * float
390 : template< size_t M, typename T >
391 1 : static vector< M, T > operator* ( T factor, const vector< M, T >& vector_ )
392 : {
393 1 : return vector_ * factor;
394 : }
395 :
396 : template< size_t M, typename T >
397 3 : inline T dot( const vector< M, T >& first, const vector< M, T >& second )
398 : {
399 3 : return first.dot( second );
400 : }
401 :
402 : template< size_t M, typename T >
403 10 : inline vector< M, T > cross( vector< M, T > a, const vector< M, T >& b )
404 : {
405 10 : return a.cross( b );
406 : }
407 :
408 : template< size_t M, typename T >
409 1 : vector< M, T > compute_normal( const vector< M, T >& a, const vector< M, T >& b,
410 : const vector< M, T >& c )
411 : {
412 : // right hand system, CCW triangle
413 1 : const vector< M, T > u = b - a;
414 1 : const vector< M, T > v = c - a;
415 1 : vector< M, T > w = cross( u, v );
416 1 : w.normalize();
417 1 : return w;
418 : }
419 :
420 : template< typename T >
421 2 : vector< 3, T > rotate( vector< 3, T > vec, const T theta,
422 : const vector< 3, T >& axis )
423 : {
424 2 : return vec.rotate( theta, axis );
425 : }
426 :
427 :
428 : template< size_t M, typename T >
429 2 : inline vector< M, T > normalize( vector< M, T > vector_ )
430 : {
431 2 : vector_.normalize();
432 2 : return vector_;
433 : }
434 :
435 : template< typename T >
436 6 : inline vector< 4, T > compute_plane( const vector< 3, T >& a,
437 : const vector< 3, T >& b,
438 : const vector< 3, T >& c )
439 : {
440 6 : const vector< 3, T > cb = b - c;
441 6 : const vector< 3, T > ba = a - b;
442 :
443 6 : vector< 4, T > plane = vector< 4, T >( cross( cb, ba ));
444 6 : plane.normalize();
445 6 : plane.w() = -plane.x() * a.x() - plane.y() * a.y() - plane.z() * a.z();
446 6 : return plane;
447 : }
448 :
449 : template< size_t M, typename T >
450 4 : vector< M, T >::vector( const T& _a )
451 : {
452 17 : for( iterator it = begin(), it_end = end(); it != it_end; ++it )
453 : {
454 13 : *it = _a;
455 : }
456 4 : }
457 :
458 : template< size_t M, typename T >
459 1 : vector< M, T >::vector( const T& _x, const T& _y )
460 : {
461 1 : array[ 0 ] = _x;
462 1 : array[ 1 ] = _y;
463 1 : }
464 :
465 : template< size_t M, typename T >
466 68 : vector< M, T >::vector( const T& _x, const T& _y, const T& _z )
467 : {
468 68 : array[ 0 ] = _x;
469 68 : array[ 1 ] = _y;
470 68 : array[ 2 ] = _z;
471 68 : }
472 :
473 : template< size_t M, typename T >
474 19 : vector< M, T >::vector( const T& _x, const T& _y, const T& _z, const T& _w )
475 : {
476 19 : array[ 0 ] = _x;
477 19 : array[ 1 ] = _y;
478 19 : array[ 2 ] = _z;
479 19 : array[ 3 ] = _w;
480 19 : }
481 :
482 : template< size_t M, typename T >
483 7 : vector< M, T >::vector( const T* values )
484 : {
485 7 : memcpy( array, values, M * sizeof( T ));
486 7 : }
487 :
488 : #ifdef __OSG_MATH
489 : template< size_t M, typename T >
490 : template< typename OSGVEC3 >
491 : vector< M, T >::vector( const OSGVEC3& from,
492 : typename enable_if< M == 3, OSGVEC3 >::type* )
493 : {
494 : array[ 0 ] = from.x();
495 : array[ 1 ] = from.y();
496 : array[ 2 ] = from.z();
497 : }
498 : #endif
499 :
500 : #ifndef SWIG
501 : template< size_t M, typename T >
502 : // initializes the first M-1 values from vector_, the last from last_
503 1 : vector< M, T >::vector( const vector< M-1, T >& vector_, T last_ )
504 : {
505 : typename vector< M-1, T >::const_iterator
506 1 : it = vector_.begin(), it_end = vector_.end();
507 :
508 1 : iterator my_it = begin();
509 :
510 4 : for( ; it != it_end; ++it, ++my_it )
511 : {
512 3 : (*my_it) = *it;
513 : }
514 1 : (*my_it) = last_;
515 1 : }
516 : #endif
517 :
518 : // to-homogenous-coordinates ctor
519 : template< size_t M, typename T >
520 : template< size_t N >
521 42 : vector< M, T >::vector( const vector< N, T >& source_,
522 : typename enable_if< N == M - 1 >::type* )
523 : {
524 42 : (*this) = source_;
525 42 : }
526 :
527 : // from-homogenous-coordinates ctor
528 : template< size_t M, typename T >
529 : template< size_t N >
530 1 : vector< M, T >::vector( const vector< N, T >& source_,
531 : typename enable_if< N == M + 1 >::type* )
532 : {
533 1 : (*this) = source_;
534 1 : }
535 :
536 : template< size_t M, typename T >
537 : template< typename U >
538 1 : vector< M, T >::vector( const vector< M, U >& source_ )
539 : {
540 1 : (*this) = source_;
541 1 : }
542 :
543 1 : template< size_t M, typename T > void vector< M, T >::set( T _a )
544 : {
545 5 : for( iterator it = begin(), it_end = end(); it != it_end; ++it )
546 4 : *it = _a;
547 1 : }
548 :
549 : #ifndef SWIG
550 : template< size_t M, typename T >
551 : void vector< M, T >::set( const vector< M-1, T >& v, T _a )
552 : {
553 : memcpy( array, v.array, sizeof( T ) * (M-1) );
554 : at( M-1 ) = _a;
555 : }
556 : #endif
557 :
558 : template< size_t M, typename T > template< size_t N >
559 : void vector< M, T >::set( const vector< N, T >& v )
560 : {
561 : size_t minimum = M;
562 : if (N < M) minimum = N;
563 : memcpy( array, v.array, sizeof( T ) * minimum );
564 : }
565 :
566 : template< size_t M, typename T >
567 : void vector< M, T >::set( T _x, T _y )
568 : {
569 : array[ 0 ] = _x;
570 : array[ 1 ] = _y;
571 : }
572 :
573 : template< size_t M, typename T >
574 : void vector< M, T >::set( T _x, T _y, T _z )
575 : {
576 : array[ 0 ] = _x;
577 : array[ 1 ] = _y;
578 : array[ 2 ] = _z;
579 : }
580 :
581 : template< size_t M, typename T >
582 1 : void vector< M, T >::set( T _x, T _y, T _z, T _w )
583 : {
584 1 : array[ 0 ] = _x;
585 1 : array[ 1 ] = _y;
586 1 : array[ 2 ] = _z;
587 1 : array[ 3 ] = _w;
588 1 : }
589 :
590 : template< size_t M, typename T >
591 : inline T&
592 26 : vector< M, T >::operator()( size_t index )
593 : {
594 26 : return at( index );
595 : }
596 :
597 : template< size_t M, typename T >
598 : inline const T&
599 48 : vector< M, T >::operator()( size_t index ) const
600 : {
601 48 : return at( index );
602 : }
603 :
604 : template< size_t M, typename T >
605 : inline T&
606 441 : vector< M, T >::at( size_t index )
607 : {
608 441 : if( index >= M )
609 0 : throw std::runtime_error( "at() - index out of bounds" );
610 441 : return array[ index ];
611 : }
612 :
613 : template< size_t M, typename T >
614 : inline const T&
615 694 : vector< M, T >::at( size_t index ) const
616 : {
617 694 : if ( index >= M )
618 0 : throw std::runtime_error( "at() - index out of bounds" );
619 694 : return array[ index ];
620 : }
621 :
622 : #ifndef VMMLIB_NO_CONVERSION_OPERATORS
623 :
624 : template< size_t M, typename T >
625 0 : vector< M, T >::operator T*()
626 : {
627 0 : return array;
628 : }
629 :
630 : template< size_t M, typename T >
631 87 : vector< M, T >::operator const T*() const
632 : {
633 87 : return array;
634 : }
635 : #else
636 :
637 : template< size_t M, typename T >
638 : T&
639 : vector< M, T >::operator[]( size_t index )
640 : {
641 : return at( index );
642 : }
643 :
644 : template< size_t M, typename T >
645 : const T&
646 : vector< M, T >::operator[]( size_t index ) const
647 : {
648 : return at( index );
649 : }
650 :
651 : #endif
652 :
653 : #if 0
654 : template< size_t M, typename T >
655 : inline T&
656 : vector< M, T >::operator[]( size_t index )
657 : {
658 : return at( index );
659 : }
660 :
661 : template< size_t M, typename T >
662 : inline const T&
663 : vector< M, T >::operator[]( size_t index ) const
664 : {
665 : return at( index );
666 : }
667 : #endif
668 :
669 : template< size_t M, typename T >
670 : vector< M, T >
671 1 : vector< M, T >::operator*( const vector< M, T >& other ) const
672 : {
673 1 : vector< M, T > result;
674 5 : for( size_t index = 0; index < M; ++index )
675 4 : result.at( index ) = at( index ) * other.at( index );
676 1 : return result;
677 : }
678 :
679 : template< size_t M, typename T >
680 : vector< M, T >
681 1 : vector< M, T >::operator/( const vector< M, T >& other ) const
682 : {
683 1 : vector< M, T > result;
684 5 : for( size_t index = 0; index < M; ++index )
685 4 : result.at( index ) = at( index ) / other.at( index );
686 1 : return result;
687 : }
688 :
689 : template< size_t M, typename T >
690 : vector< M, T >
691 12 : vector< M, T >::operator+( const vector< M, T >& other ) const
692 : {
693 12 : vector< M, T > result;
694 52 : for( size_t index = 0; index < M; ++index )
695 40 : result.at( index ) = at( index ) + other.at( index );
696 12 : return result;
697 : }
698 :
699 : template< size_t M, typename T >
700 : vector< M, T >
701 27 : vector< M, T >::operator-( const vector< M, T >& other ) const
702 : {
703 27 : vector< M, T > result;
704 112 : for( size_t index = 0; index < M; ++index )
705 85 : result.at( index ) = at( index ) - other.at( index );
706 27 : return result;
707 : }
708 :
709 : template< size_t M, typename T >
710 : void
711 1 : vector< M, T >::operator*=( const vector< M, T >& other )
712 : {
713 5 : for( size_t index = 0; index < M; ++index )
714 4 : at( index ) *= other.at( index );
715 1 : }
716 :
717 : template< size_t M, typename T >
718 : void
719 1 : vector< M, T >::operator/=( const vector< M, T >& other )
720 : {
721 5 : for( size_t index = 0; index < M; ++index )
722 4 : at( index ) /= other.at( index );
723 1 : }
724 :
725 : template< size_t M, typename T >
726 : void
727 1 : vector< M, T >::operator+=( const vector< M, T >& other )
728 : {
729 5 : for( size_t index = 0; index < M; ++index )
730 4 : at( index ) += other.at( index );
731 1 : }
732 :
733 : template< size_t M, typename T >
734 : void
735 1 : vector< M, T >::operator-=( const vector< M, T >& other )
736 : {
737 5 : for( size_t index = 0; index < M; ++index )
738 4 : at( index ) -= other.at( index );
739 1 : }
740 :
741 : template< size_t M, typename T >
742 : vector< M, T >
743 14 : vector< M, T >::operator*( const T other ) const
744 : {
745 14 : vector< M, T > result;
746 58 : for( size_t index = 0; index < M; ++index )
747 44 : result.at( index ) = at( index ) * other;
748 14 : return result;
749 : }
750 :
751 : template< size_t M, typename T >
752 : vector< M, T >
753 1 : vector< M, T >::operator/( const T other ) const
754 : {
755 1 : vector< M, T > result;
756 5 : for( size_t index = 0; index < M; ++index )
757 4 : result.at( index ) = at( index ) / other;
758 1 : return result;
759 : }
760 :
761 : template< size_t M, typename T >
762 : vector< M, T >
763 1 : vector< M, T >::operator+( const T other ) const
764 : {
765 1 : vector< M, T > result;
766 5 : for( size_t index = 0; index < M; ++index )
767 4 : result.at( index ) = at( index ) + other;
768 1 : return result;
769 : }
770 :
771 : template< size_t M, typename T >
772 : vector< M, T >
773 1 : vector< M, T >::operator-( const T other ) const
774 : {
775 1 : vector< M, T > result;
776 5 : for( size_t index = 0; index < M; ++index )
777 4 : result.at( index ) = at( index ) - other;
778 1 : return result;
779 : }
780 :
781 : template< size_t M, typename T >
782 : void
783 20 : vector< M, T >::operator*=( const T other )
784 : {
785 88 : for( size_t index = 0; index < M; ++index )
786 68 : at( index ) *= other;
787 20 : }
788 :
789 : template< size_t M, typename T >
790 : void
791 1 : vector< M, T >::operator/=( const T other )
792 : {
793 5 : for( size_t index = 0; index < M; ++index )
794 4 : at( index ) /= other;
795 1 : }
796 :
797 : template< size_t M, typename T >
798 : void
799 1 : vector< M, T >::operator+=( const T other )
800 : {
801 5 : for( size_t index = 0; index < M; ++index )
802 4 : at( index ) += other;
803 1 : }
804 :
805 : template< size_t M, typename T >
806 : void
807 1 : vector< M, T >::operator-=( const T other )
808 : {
809 5 : for( size_t index = 0; index < M; ++index )
810 4 : at( index ) -= other;
811 1 : }
812 :
813 : template< size_t M, typename T >
814 : vector< M, T >
815 2 : vector< M, T >::operator-() const
816 : {
817 2 : vector< M, T > v( *this );
818 2 : return v.negate();
819 : }
820 :
821 : template< size_t M, typename T >
822 : const vector< M, T >&
823 2 : vector< M, T >::negate()
824 : {
825 8 : for( size_t index = 0; index < M; ++index )
826 6 : array[ index ] = -array[ index ];
827 2 : return *this;
828 : }
829 :
830 : template< size_t M, typename T >
831 : inline T&
832 98 : vector< M, T >::x()
833 : {
834 98 : return array[ 0 ];
835 : }
836 :
837 : template< size_t M, typename T >
838 : inline T&
839 99 : vector< M, T >::y()
840 : {
841 99 : return array[ 1 ];
842 : }
843 :
844 : template< size_t M, typename T >
845 : inline T&
846 94 : vector< M, T >::z()
847 : {
848 94 : return array[ 2 ];
849 : }
850 :
851 : template< size_t M, typename T >
852 : inline T&
853 13 : vector< M, T >::w()
854 : {
855 13 : return array[ 3 ];
856 : }
857 :
858 : template< size_t M, typename T >
859 : inline const T&
860 187 : vector< M, T >::x() const
861 : {
862 187 : return array[ 0 ];
863 : }
864 :
865 : template< size_t M, typename T >
866 : inline const T&
867 181 : vector< M, T >::y() const
868 : {
869 181 : return array[ 1 ];
870 : }
871 :
872 : template< size_t M, typename T >
873 : inline const T&
874 179 : vector< M, T >::z() const
875 : {
876 179 : return array[ 2 ];
877 : }
878 :
879 : template< size_t M, typename T >
880 : inline const T&
881 100 : vector< M, T >::w() const
882 : {
883 100 : return array[ 3 ];
884 : }
885 :
886 : template< size_t M, typename T >
887 : inline T&
888 : vector< M, T >::r()
889 : {
890 : return array[ 0 ];
891 : }
892 :
893 : template< size_t M, typename T >
894 : inline T&
895 : vector< M, T >::g()
896 : {
897 : return array[ 1 ];
898 : }
899 :
900 : template< size_t M, typename T >
901 : inline T&
902 : vector< M, T >::b()
903 : {
904 : return array[ 2 ];
905 : }
906 :
907 : template< size_t M, typename T >
908 : inline T&
909 : vector< M, T >::a()
910 : {
911 : return array[ 3 ];
912 : }
913 :
914 : template< size_t M, typename T >
915 : inline const T&
916 : vector< M, T >::r() const
917 : {
918 : return array[ 0 ];
919 : }
920 :
921 : template< size_t M, typename T >
922 : inline const T&
923 : vector< M, T >::g() const
924 : {
925 : return array[ 1 ];
926 : }
927 :
928 : template< size_t M, typename T >
929 : inline const T&
930 : vector< M, T >::b() const
931 : {
932 : return array[ 2 ];
933 : }
934 :
935 : template< size_t M, typename T >
936 : inline const T&
937 : vector< M, T >::a() const
938 : {
939 : return array[ 3 ];
940 : }
941 :
942 : template< size_t M, typename T > template< typename TT >
943 11 : vector< M, T >& vector< M, T >::cross( const vector< M, TT >& rhs,
944 : typename enable_if< M == 3, TT >::type* )
945 : {
946 11 : const T x_ = y() * rhs.z() - z() * rhs.y();
947 11 : const T y_ = z() * rhs.x() - x() * rhs.z();
948 11 : const T z_ = x() * rhs.y() - y() * rhs.x();
949 11 : x() = x_;
950 11 : y() = y_;
951 11 : z() = z_;
952 11 : return *this;
953 : }
954 :
955 : template< size_t M, typename T >
956 38 : inline T vector< M, T >::dot( const vector< M, T >& other ) const
957 : {
958 38 : T tmp = 0.0;
959 186 : for( size_t index = 0; index < M; ++index )
960 148 : tmp += at( index ) * other.at( index );
961 :
962 38 : return tmp;
963 : }
964 :
965 19 : template< size_t M, typename T > inline T vector< M, T >::normalize()
966 : {
967 19 : const T len = length();
968 19 : if ( len <= std::numeric_limits< T >::epsilon( ))
969 0 : return 0;
970 :
971 19 : const T tmp = 1.0 / len;
972 19 : (*this) *= tmp;
973 19 : return len;
974 : }
975 :
976 : template< size_t M, typename T >
977 27 : inline T vector< M, T >::length() const
978 : {
979 27 : return std::sqrt( squared_length() );
980 : }
981 :
982 : template< size_t M, typename T >
983 28 : inline T vector< M, T >::squared_length() const
984 : {
985 28 : T _squared_length = 0.0;
986 122 : for( const_iterator it = begin(), it_end = end(); it != it_end; ++it )
987 94 : _squared_length += (*it) * (*it);
988 :
989 28 : return _squared_length;
990 : }
991 :
992 : template< size_t M, typename T >
993 : inline T
994 : vector< M, T >::distance( const vector< M, T >& other ) const
995 : {
996 : return std::sqrt( squared_distance( other ) );
997 : }
998 :
999 : template< size_t M, typename T >
1000 : inline T vector< M, T >::squared_distance( const vector< M, T >& other ) const
1001 : {
1002 : vector< M, T > tmp( *this );
1003 : tmp -= other;
1004 : return tmp.squared_length();
1005 : }
1006 :
1007 2 : template< size_t M, typename T > inline T vector< M, T >::product() const
1008 : {
1009 2 : T result = at( 0 );
1010 6 : for( size_t i = 1; i < M; ++i )
1011 4 : result *= at( i );
1012 2 : return result;
1013 : }
1014 :
1015 : template< size_t M, typename T > template< typename TT >
1016 3 : vector< 3, T >& vector< M, T >::rotate( const T theta, vector< M, TT > axis,
1017 : typename enable_if< M==3, TT >::type* )
1018 : {
1019 3 : const T costheta = std::cos( theta );
1020 3 : const T sintheta = std::sin( theta );
1021 :
1022 3 : axis.normalize();
1023 105 : return *this = vector< 3, T >(
1024 9 : (costheta + ( 1 - costheta ) * axis.x() * axis.x() ) * x() +
1025 12 : (( 1 - costheta ) * axis.x() * axis.y() - axis.z() * sintheta ) * y() +
1026 12 : (( 1 - costheta ) * axis.x() * axis.z() + axis.y() * sintheta ) * z(),
1027 :
1028 12 : (( 1 - costheta ) * axis.x() * axis.y() + axis.z() * sintheta ) * x() +
1029 9 : ( costheta + ( 1 - costheta ) * axis.y() * axis.y() ) * y() +
1030 12 : (( 1 - costheta ) * axis.y() * axis.z() - axis.x() * sintheta ) * z(),
1031 :
1032 12 : (( 1 - costheta ) * axis.x() * axis.z() - axis.y() * sintheta ) * x() +
1033 12 : (( 1 - costheta ) * axis.y() * axis.z() + axis.x() * sintheta ) * y() +
1034 15 : ( costheta + ( 1 - costheta ) * axis.z() * axis.z() ) * z( ));
1035 : }
1036 :
1037 : // sphere layout: center xyz, radius w
1038 : template< size_t M, typename T > template< typename TT > inline vector< 3, T >
1039 : vector< M, T >::project_point_onto_sphere( const vector< 3, TT >& point,
1040 : typename enable_if< M == 4, TT >::type* ) const
1041 : {
1042 : const vector< 3, T >& center_ = get_sub_vector< 3 >();
1043 :
1044 : vector< 3, T > projected_point( point );
1045 : projected_point -= center_;
1046 : projected_point.normalize();
1047 : projected_point *= w();
1048 : return center_ + projected_point;
1049 : }
1050 :
1051 : // sphere layout: center xyz, radius w
1052 : template< size_t M, typename T > template< typename TT > inline T
1053 : vector< M, T >::distance_to_sphere( const vector< 3, TT >& point,
1054 : typename enable_if< M == 4, TT >::type* )
1055 : const
1056 : {
1057 : const vector< 3, T >& center_ = get_sub_vector< 3 >();
1058 : return ( point - center_ ).length() - w();
1059 : }
1060 :
1061 : template< size_t M, typename T > template< size_t N, size_t O >
1062 7 : vector< N, T > vector< M, T >::get_sub_vector(
1063 : typename enable_if< M >= N+O >::type* ) const
1064 : {
1065 7 : return vector< N, T >( array + O );
1066 : }
1067 :
1068 : template< size_t M, typename T > template< size_t N, size_t O >
1069 1 : void vector< M, T >::set_sub_vector( const vector< N, T >& sub,
1070 : typename enable_if< M >= N+O >::type* )
1071 : {
1072 1 : ::memcpy( array + O, sub.array, N * sizeof( T ));
1073 1 : }
1074 :
1075 : // plane: normal xyz, distance w
1076 : template< size_t M, typename T > template< typename TT >
1077 : inline T vector< M, T >::distance_to_plane( const vector< 3, TT >& point,
1078 : typename enable_if< M == 4, TT >::type* ) const
1079 : {
1080 : const vector< 3, T >& normal = get_sub_vector< 3 >();
1081 : return normal.dot( point ) + w();
1082 : }
1083 :
1084 : // plane: normal xyz, distance w
1085 : template< size_t M, typename T > template< typename TT > vector< 3, T >
1086 : vector< M, T >::project_point_onto_plane( const vector< 3, TT >& point,
1087 : typename enable_if< M == 4, TT >::type* ) const
1088 : {
1089 : const vector< 3, T >& normal = get_sub_vector< 3 >();
1090 : return point - ( normal * distance_to_plane( point ) );
1091 : }
1092 :
1093 : template< size_t M, typename T >
1094 28 : bool vector< M, T >::operator==( const vector< M, T >& other ) const
1095 : {
1096 28 : return memcmp( array, other.array, sizeof( array )) == 0;
1097 : }
1098 :
1099 : template< size_t M, typename T >
1100 : bool vector< M, T >::operator!=( const vector< M, T >& other ) const
1101 : {
1102 : return ! this->operator==( other );
1103 : }
1104 :
1105 : template< size_t M, typename T >
1106 2 : bool vector< M, T >::equals( const vector< M, T >& other, T tolerance ) const
1107 : {
1108 8 : for( size_t index = 0; index < M; ++index )
1109 6 : if( fabs( at( index ) - other( index ) ) >= tolerance )
1110 0 : return false;
1111 2 : return true;
1112 :
1113 : }
1114 :
1115 : template< size_t M, typename T >
1116 : bool
1117 : vector< M, T >::operator<( const vector< M, T >& other ) const
1118 : {
1119 : for(size_t index = 0; index < M; ++index )
1120 : {
1121 : if (at( index ) < other.at( index )) return true;
1122 : if (other.at( index ) < at( index )) return false;
1123 : }
1124 : return false;
1125 : }
1126 :
1127 : // to-homogenous-coordinates assignment operator
1128 : // non-chainable because of sfinae
1129 : template< size_t M, typename T > template< size_t N >
1130 : typename enable_if< N == M - 1 >::type*
1131 42 : vector< M, T >::operator=( const vector< N, T >& source_ )
1132 : {
1133 42 : std::copy( source_.begin(), source_.end(), begin() );
1134 42 : at( M - 1 ) = static_cast< T >( 1.0 );
1135 42 : return 0;
1136 : }
1137 :
1138 : // from-homogenous-coordinates assignment operator
1139 : // non-chainable because of sfinae
1140 : template< size_t M, typename T > template< size_t N >
1141 : typename enable_if< N == M + 1 >::type*
1142 2 : vector< M, T >::operator=( const vector< N, T >& source_ )
1143 : {
1144 2 : const T w_reci = static_cast< T >( 1.0 ) / source_( M );
1145 2 : iterator it = begin(), it_end = end();
1146 8 : for( size_t index = 0; it != it_end; ++it, ++index )
1147 6 : *it = source_( index ) * w_reci;
1148 2 : return 0;
1149 : }
1150 :
1151 : template< size_t M, typename T >
1152 19 : vector< M, T >& vector< M, T >::operator=( const T* c_array )
1153 : {
1154 19 : iter_set( c_array, c_array + M );
1155 19 : return *this;
1156 : }
1157 :
1158 : template< size_t M, typename T >
1159 : T vector< M, T >::operator=( T filler_value )
1160 : {
1161 : for( size_t index = 0; index < M; ++index )
1162 : at( index ) = filler_value;
1163 : return filler_value;
1164 : }
1165 :
1166 : template< size_t M, typename T >
1167 34 : vector< M, T >& vector< M, T >::operator=( const vector< M, T >& other )
1168 : {
1169 34 : if( this != &other )
1170 34 : memcpy( array, other.array, M * sizeof( T ) );
1171 34 : return *this;
1172 : }
1173 :
1174 : // returns void to avoid 'silent' loss of precision when chaining
1175 : template< size_t M, typename T > template< typename U >
1176 2 : void vector< M, T >::operator=( const vector< M, U >& source_ )
1177 : {
1178 : typedef typename vector< M, U >::const_iterator u_c_iter;
1179 2 : u_c_iter it = source_.begin(), it_end = source_.end();
1180 8 : for( iterator my_it = begin(); it != it_end; ++it, ++my_it )
1181 6 : *my_it = static_cast< T >( *it );
1182 2 : }
1183 :
1184 : template< size_t M, typename T >
1185 : template< typename input_iterator_t >
1186 : void
1187 23 : vector< M, T >::iter_set( input_iterator_t begin_, input_iterator_t end_ )
1188 : {
1189 23 : input_iterator_t in_it = begin_;
1190 23 : iterator it = begin(), it_end = end();
1191 110 : for( ; it != it_end && in_it != end_; ++it, ++in_it )
1192 87 : (*it) = static_cast< T >( *in_it );
1193 23 : }
1194 :
1195 : template< size_t M, typename T >
1196 : void vector< M, T >::clamp( const T& min, const T& max )
1197 : {
1198 : for( size_t i = 0; i < M; ++i )
1199 : {
1200 : if( array[i] < min )
1201 : array[i] = min;
1202 : if( array[i] > max )
1203 : array[i] = max;
1204 : }
1205 : }
1206 :
1207 : template< size_t M, typename T >
1208 : inline size_t
1209 : vector< M, T >::size()
1210 : {
1211 : return M;
1212 : }
1213 :
1214 : template< size_t M, typename T >
1215 : size_t
1216 2 : vector< M, T >::find_min_index() const
1217 : {
1218 2 : return std::min_element( begin(), end() ) - begin();
1219 : }
1220 :
1221 : template< size_t M, typename T >
1222 : size_t
1223 2 : vector< M, T >::find_max_index() const
1224 : {
1225 2 : return std::max_element( begin(), end() ) - begin();
1226 : }
1227 :
1228 : template< size_t M, typename T >
1229 : T&
1230 2 : vector< M, T >::find_min()
1231 : {
1232 2 : return *std::min_element( begin(), end() );
1233 : }
1234 :
1235 : template< size_t M, typename T >
1236 : const T&
1237 : vector< M, T >::find_min() const
1238 : {
1239 : return *std::min_element( begin(), end() );
1240 : }
1241 :
1242 : template< size_t M, typename T >
1243 : T&
1244 2 : vector< M, T >::find_max()
1245 : {
1246 2 : return *std::max_element( begin(), end() );
1247 : }
1248 :
1249 : template< size_t M, typename T >
1250 : const T&
1251 : vector< M, T >::find_max() const
1252 : {
1253 : return *std::max_element( begin(), end() );
1254 : }
1255 :
1256 : template< size_t M, typename T >
1257 : inline typename vector< M, T >::iterator
1258 83 : vector< M, T >::begin()
1259 : {
1260 83 : return array;
1261 : }
1262 :
1263 : template< size_t M, typename T >
1264 : inline typename vector< M, T >::iterator
1265 37 : vector< M, T >::end()
1266 : {
1267 37 : return array + M; ;
1268 : }
1269 :
1270 : template< size_t M, typename T >
1271 : inline typename vector< M, T >::const_iterator
1272 82 : vector< M, T >::begin() const
1273 : {
1274 82 : return array;
1275 : }
1276 :
1277 : template< size_t M, typename T >
1278 : inline typename vector< M, T >::const_iterator
1279 78 : vector< M, T >::end() const
1280 : {
1281 78 : return array + M; ;
1282 : }
1283 :
1284 : template< size_t M, typename T >
1285 : inline typename vector< M, T >::reverse_iterator
1286 : vector< M, T >::rbegin()
1287 : {
1288 : return array + M - 1;
1289 : }
1290 :
1291 : template< size_t M, typename T >
1292 : inline typename vector< M, T >::reverse_iterator
1293 : vector< M, T >::rend()
1294 : {
1295 : return array - 1;
1296 : }
1297 :
1298 : template< size_t M, typename T >
1299 : inline typename vector< M, T >::const_reverse_iterator
1300 : vector< M, T >::rbegin() const
1301 : {
1302 : return array + M - 1;
1303 : }
1304 :
1305 : template< size_t M, typename T >
1306 : inline typename vector< M, T >::const_reverse_iterator
1307 : vector< M, T >::rend() const
1308 : {
1309 : return array - 1;
1310 : }
1311 :
1312 : template< size_t M, typename T >
1313 : bool
1314 : vector< M, T >::is_unit_vector() const
1315 : {
1316 : const_iterator it = begin(), it_end = end();
1317 : bool one = false;
1318 : for( ; it != it_end; ++it )
1319 : {
1320 : if ( *it == 1.0 )
1321 : {
1322 : if ( one )
1323 : return false;
1324 : one = true;
1325 : }
1326 : else if ( *it != 0.0 )
1327 : {
1328 : return false;
1329 : }
1330 : }
1331 : return one;
1332 : }
1333 :
1334 : template< size_t M, typename T >
1335 : void
1336 : vector< M, T >::perturb( T perturbation )
1337 : {
1338 : for( iterator it = begin(), it_end = end(); it != it_end; ++it )
1339 : {
1340 : (*it) += ( rand() & 1u ) ? perturbation : -perturbation;
1341 : }
1342 :
1343 : }
1344 :
1345 : template< size_t M, typename T >
1346 : void
1347 1 : vector< M, T >::sqrt_elementwise()
1348 : {
1349 5 : for( iterator it = begin(), it_end = end(); it != it_end; ++it )
1350 : {
1351 4 : (*it) = std::sqrt(*it);
1352 : }
1353 1 : }
1354 :
1355 1 : template< size_t M, typename T > void vector< M, T >::reciprocal()
1356 : {
1357 5 : for( iterator it = begin(), it_end = end(); it != it_end; ++it )
1358 4 : (*it) = static_cast< T >( 1 ) / (*it);
1359 1 : }
1360 :
1361 : template< size_t M, typename T > void vector< M, T >::reciprocal_safe()
1362 : {
1363 : for( iterator it = begin(), it_end = end(); it != it_end; ++it )
1364 : {
1365 : T& v = *it;
1366 :
1367 : if( v == 0 )
1368 : v = std::numeric_limits< T >::max();
1369 : else
1370 : v = static_cast< T >( 1 ) / v;
1371 : }
1372 : }
1373 :
1374 : template< size_t M, typename T >
1375 : template< typename TT >
1376 : void
1377 : vector< M, T >::cast_from( const vector< M, TT >& other )
1378 : {
1379 : typedef vmml::vector< M, TT > vector_tt_type ;
1380 : typedef typename vector_tt_type::const_iterator tt_const_iterator;
1381 :
1382 : iterator it = begin(), it_end = end();
1383 : tt_const_iterator other_it = other.begin();
1384 : for( ; it != it_end; ++it, ++other_it )
1385 : {
1386 : *it = static_cast< T >( *other_it );
1387 : }
1388 : }
1389 :
1390 : template< size_t M, typename T >
1391 : size_t
1392 : vector< M, T >::nnz() const
1393 : {
1394 : size_t counter = 0;
1395 :
1396 : const_iterator it = begin(),
1397 : it_end = end();
1398 : for( ; it != it_end; ++it)
1399 : {
1400 : if ( *it != 0 ) {
1401 : ++counter;
1402 : }
1403 : }
1404 :
1405 : return counter;
1406 : }
1407 :
1408 : template< size_t M, typename T >
1409 : double
1410 1 : vector< M, T >::norm( ) const
1411 : {
1412 1 : double norm_v = 0.0;
1413 :
1414 1 : const_iterator it = begin(), it_end = end();
1415 5 : for( ; it != it_end; ++it )
1416 : {
1417 4 : norm_v += *it * *it;
1418 : }
1419 :
1420 1 : return std::sqrt(norm_v);
1421 : }
1422 :
1423 : template< size_t M, typename T >
1424 : void
1425 : vector< M, T >::set_random( int seed )
1426 : {
1427 : if ( seed >= 0 )
1428 : srand( seed );
1429 :
1430 : for( size_t i = 0; i < M; ++i )
1431 : {
1432 : const double fillValue = double( rand( )) / double( RAND_MAX );
1433 : at( i ) = -1.0 + 2.0 * fillValue;
1434 : }
1435 : }
1436 :
1437 : } // namespace vmml
1438 :
1439 : #endif
|
