Theoretica
Mathematical Library
Loading...
Searching...
No Matches
vec.h
Go to the documentation of this file.
1
5
6#ifndef THEORETICA_VECTOR_H
7#define THEORETICA_VECTOR_H
8
9#ifndef THEORETICA_NO_PRINT
10#include <sstream>
11#include <ostream>
12#endif
13
14#include "../core/error.h"
15#include "../core/real_analysis.h"
16#include "./algebra.h"
17#include <vector>
18
19
20namespace theoretica {
21
22
28 template<typename Vector, typename ReturnType = vector_element_t<Vector>&>
30
31 private:
32
34 Vector& vector;
35
37 size_t i;
38
39 public:
40
41 using iterator_category = std::forward_iterator_tag;
42 using value_type = vector_element_t<Vector>;
43 using pointer = value_type*;
44 using reference = value_type&;
45
48 vec_iterator(Vector& vector, size_t index) : vector(vector), i(index) {}
49
52 ReturnType operator*() {
53 return vector[i];
54 }
55
58 ++i;
59 return *this;
60 }
61
63 // vec_iterator& operator--() {
64 // --i;
65 // return *this;
66 // }
67
68
70 size_t index() {
71 return i;
72 }
73
74
76 bool operator==(const vec_iterator& other) const {
77 return i == other.i;
78 }
79
80 bool operator!=(const vec_iterator& other) const {
81 return !(*this == other);
82 }
83 };
84
85
91 template<typename Type = real, unsigned int N = 0>
92 class vec {
93
94 private:
95 Type data[N];
96
97 public:
98
99 // Vector size template argument
100 static constexpr size_t size_argument = N;
101
104 vec() {
105 algebra::vec_zeroes(*this);
106 }
107
108
111 vec(Type val) {
112
113 for (unsigned int i = 0; i < N; ++i)
114 data[i] = val;
115 }
116
117
121 vec(unsigned int size, Type val) {
122
123 resize(size);
124
125 for (unsigned int i = 0; i < N; ++i)
126 data[i] = val;
127 }
128
129
131 template<unsigned int M>
132 vec(const vec<Type, M>& other) {
133 algebra::vec_copy(*this, other);
134 }
135
136
138 template<typename Vector>
139 vec<Type, N>& operator=(const Vector& other) {
140 return algebra::vec_copy(*this, other);
141 }
142
143
145 vec(std::initializer_list<Type> l) {
146
147 if(l.size() != N) {
148 TH_MATH_ERROR("vec::vec(initializer_list<Type>)", l.size(), INVALID_ARGUMENT);
149 algebra::vec_error(*this);
150 return;
151 }
152
153 std::copy(l.begin(), l.end(), &data[0]);
154 }
155
156 ~vec() = default;
157
158
160 inline vec<Type, N> operator+() const {
161 return *this;
162 }
163
164
166 inline vec<Type, N> operator+(const vec<Type, N>& other) const {
167
168 vec<Type, N> result;
169 algebra::vec_sum(result, *this, other);
170 return result;
171 }
172
173
175 inline vec<Type, N> operator-() const {
176 return *this * (Type) -1;
177 }
178
179
181 inline vec<Type, N> operator-(const vec<Type, N>& other) const {
182
183 vec<Type, N> result;
184 algebra::vec_diff(result, *this, other);
185 return result;
186 }
187
188
190 inline vec<Type, N> operator*(Type scalar) const {
191 vec<Type, N> result;
192
193 for (unsigned int i = 0; i < N; ++i)
194 result.data[i] = scalar * data[i];
195
196 return result;
197 }
198
199
201 inline vec<Type, N> operator/(Type scalar) const {
202 vec<Type, N> result;
203
204 for (unsigned int i = 0; i < N; ++i)
205 result.data[i] = data[i] / scalar;
206
207 return result;
208 }
209
210
212 template<typename Vector>
213 inline Type dot(const Vector& other) const {
214 return algebra::dot(*this, other);
215 }
216
217
219 template<typename Vector>
220 inline Type operator*(const Vector& other) const {
221 return dot(other);
222 }
223
224
226 inline vec<Type, N> cross(const vec<Type, N>& other) const {
227 static_assert(N == 3, "The vector must be three dimensional");
228 return algebra::cross(*this, other);
229 }
230
231
233 template<typename Vector>
234 inline vec<Type, N> cross(const Vector& other) const {
235 return algebra::cross(*this, other);
236 }
237
238
240 template<typename Vector>
241 inline vec<Type, N>& operator+=(const Vector& other) {
242
243 for (unsigned int i = 0; i < N; ++i)
244 data[i] += other.data[i];
245
246 return *this;
247 }
248
249
251 template<typename Vector>
252 inline vec<Type, N>& operator-=(const Vector& other) {
253
254 for (unsigned int i = 0; i < N; ++i)
255 data[i] -= other.data[i];
256
257 return *this;
258 }
259
260
262 inline vec<Type, N>& operator*=(Type scalar) {
263
264 for (unsigned int i = 0; i < N; ++i)
265 data[i] *= scalar;
266
267 return *this;
268 }
269
270
272 inline vec<Type, N>& operator/=(Type scalar) {
273
274 if(abs(scalar) < MACH_EPSILON) {
275 TH_MATH_ERROR("vec::operator/=", scalar, DIV_BY_ZERO);
276 algebra::vec_error(*this);
277 return *this;
278 }
279
280 for (unsigned int i = 0; i < N; ++i)
281 data[i] /= scalar;
282
283 return *this;
284 }
285
286
288 inline Type norm() const {
289 return algebra::norm(*this);
290 }
291
292
294 inline Type sqr_norm() const {
295 return algebra::sqr_norm(*this);
296 }
297
298
300 inline Type& operator[](unsigned int i) {
301 return data[i];
302 }
303
304
306 inline const Type& operator[](unsigned int i) const {
307 return data[i];
308 }
309
310
315 inline Type& at(unsigned int i) {
316
317 if (i >= N) {
318 throw std::out_of_range(
319 "The element index in vec::at() is out of bounds"
320 );
321 }
322
323 return data[i];
324 }
325
326
331 inline Type at(unsigned int i) const {
332
333 if (i >= N) {
334 throw std::out_of_range(
335 "The element index in vec::at() is out of bounds"
336 );
337 }
338
339 return data[i];
340 }
341
342
344 using iterator = vec_iterator<vec<Type, N>, Type&>;
345
347 using const_iterator = vec_iterator<const vec<Type, N>, const Type&>;
348
349
352 inline auto begin() {
353 return iterator(*this, 0);
354 }
355
356
359 inline auto end() {
360 return iterator(*this, size());
361 }
362
363
365 inline auto begin() const {
366 return const_iterator(*this, 0);
367 }
368
369
372 inline auto end() const {
373 return const_iterator(*this, size());
374 }
375
376
378 inline void normalize() {
380 }
381
382
384 inline vec<Type, N> normalized() const {
385 return algebra::normalize(*this);
386 }
387
388
390 template<typename Vector>
391 inline bool operator==(const Vector& other) const {
392
393 if(size() != other.size())
394 return false;
395
396 for (unsigned int i = 0; i < N; ++i)
397 if(data[i] != other[i])
398 return false;
399
400 return true;
401 }
402
403
405 template<typename Vector>
406 inline bool operator!=(const Vector& other) const {
407 return !(*this == other);
408 }
409
410
412 inline TH_CONSTEXPR unsigned int size() const {
413 return N;
414 }
415
416
422 inline void resize(size_t n) const {
423
424 if(N != n) {
425 TH_MATH_ERROR("vec::resize", N, INVALID_ARGUMENT);
426 }
427 }
428
429
432 inline static vec<Type, N> euclidean_base(
433 unsigned int i, unsigned int n = N) {
434
435 if(i >= n) {
436 TH_MATH_ERROR("vec::euclidean_base", i, INVALID_ARGUMENT);
437 return vec<Type, N>(Type(nan()));
438 }
439
440 vec<Type, N> e_i = vec<Type, N>(n, Type(0.0));
441 e_i.resize(n);
442 e_i[i] = 1;
443
444 return e_i;
445 }
446
447
450 inline friend vec<Type, N> operator*(Type a, const vec<Type, N>& v) {
451 return v * a;
452 }
453
454
455#ifndef THEORETICA_NO_PRINT
456
458 inline std::string to_string(
459 const std::string& separator = ", ",
460 bool parenthesis = true) const {
461
462 std::stringstream res;
463
464 if(parenthesis)
465 res << "(";
466
467 for (unsigned int i = 0; i < N; ++i) {
468 res << data[i];
469 if(i != N - 1)
470 res << separator;
471 }
472
473 if(parenthesis)
474 res << ")";
475
476 return res.str();
477 }
478
479
481 inline operator std::string() {
482 return to_string();
483 }
484
485
487 inline friend std::ostream& operator<<(
488 std::ostream& out, const vec<Type, N>& obj) {
489 return out << obj.to_string();
490 }
491
492#endif
493
494 };
495
496
502 template<typename Type>
503 class vec<Type, 0> {
504
505 // Container type for storage (alias for std::vector)
506 template<typename T>
507 using Container = std::vector<T>;
508
509 private:
510 Container<Type> data;
511
512 public:
513
514 // Vector size template argument
515 static constexpr size_t size_argument = 0;
516
518 vec() {}
519
520
523 vec(unsigned int n) {
524 resize(n);
525 algebra::vec_zeroes(*this);
526 }
527
528
531 vec(unsigned int n, Type a) {
532 data = std::vector<Type>(n, a);
533 }
534
535
537 template<unsigned int M>
538 vec(const vec<Type, M>& other) {
539 algebra::vec_copy(*this, other);
540 }
541
542
544 template<typename Vector>
545 vec<Type>& operator=(const Vector& other) {
546 return algebra::vec_copy(*this, other);
547 }
548
549
551 vec(std::initializer_list<Type> l) {
552
553 resize(l.size());
554 std::copy(l.begin(), l.end(), &data[0]);
555 }
556
557 ~vec() = default;
558
559
561 inline vec<Type> operator+() const {
562 return *this;
563 }
564
565
567 template<typename Vector>
568 inline vec<Type> operator+(const Vector& other) const {
569
570 vec<Type> result;
571 result.resize(size());
572 algebra::vec_sum(result, *this, other);
573 return result;
574 }
575
576
578 inline vec<Type> operator-() const {
579 return *this * (Type) -1;
580 }
581
582
584 template<typename Vector>
585 inline vec<Type> operator-(const Vector& other) const {
586
587 vec<Type> result;
588 result.resize(size());
589 algebra::vec_diff(result, *this, other);
590 return result;
591 }
592
593
595 inline vec<Type> operator*(Type scalar) const {
596
597 vec<Type> result;
598 result.resize(size());
599
600 for (unsigned int i = 0; i < size(); ++i)
601 result.data[i] = scalar * data[i];
602
603 return result;
604 }
605
606
608 inline vec<Type> operator/(Type scalar) const {
609
610 vec<Type> result;
611 result.resize(size());
612
613 for (unsigned int i = 0; i < size(); ++i)
614 result.data[i] = data[i] / scalar;
615
616 return result;
617 }
618
619
621 template<typename Vector>
622 inline Type dot(const Vector& other) const {
623 return algebra::dot(*this, other);
624 }
625
626
628 template<typename Vector>
629 inline Type operator*(const Vector& other) const {
630 return dot(other);
631 }
632
633
635 template<typename Vector>
636 inline vec<Type> cross(const Vector& other) const {
637 return algebra::cross(*this, other);
638 }
639
640
642 template<typename Vector>
643 inline vec<Type>& operator+=(const Vector& other) {
644
645 // If the vector is uninitialized,
646 // initialize it to be a zero vector
647 // with the target size
648 if(size() == 0)
649 resize(other.size());
650
651 if(size() != other.size()) {
652 TH_MATH_ERROR("vec::operator+=", size(), INVALID_ARGUMENT);
653 return (*this = vec<Type>(max(size(), 1), nan()));
654 }
655
656 for (unsigned int i = 0; i < size(); ++i)
657 data[i] += other.data[i];
658
659 return *this;
660 }
661
662
664 template<typename Vector>
665 inline vec<Type>& operator-=(const Vector& other) {
666
667 if(size() != other.size()) {
668 TH_MATH_ERROR("vec::operator-=", size(), INVALID_ARGUMENT);
669 return (*this = vec<Type>(max(size(), 1), nan()));
670 }
671
672 for (unsigned int i = 0; i < size(); ++i)
673 data[i] -= other.data[i];
674
675 return *this;
676 }
677
678
680 inline vec<Type>& operator*=(Type scalar) {
681
682 for (unsigned int i = 0; i < size(); ++i)
683 data[i] *= scalar;
684
685 return *this;
686 }
687
688
690 inline vec<Type>& operator/=(Type scalar) {
691
692 if(abs(scalar) < MACH_EPSILON) {
693 TH_MATH_ERROR("vec::operator/=", scalar, DIV_BY_ZERO);
694 *this = vec<Type>(max(size(), 1), nan());
695 return *this;
696 }
697
698 for (unsigned int i = 0; i < size(); ++i)
699 data[i] /= scalar;
700
701 return *this;
702 }
703
704
706 inline Type norm() const {
707 return algebra::norm(*this);
708 }
709
710
712 inline Type sqr_norm() const {
713 return algebra::sqr_norm(*this);
714 }
715
716
718 inline Type& operator[](unsigned int i) {
719 return data[i];
720 }
721
722
724 inline const Type& operator[](unsigned int i) const {
725 return data[i];
726 }
727
728
733 inline Type& at(unsigned int i) {
734 return data.at(i);
735 }
736
737
742 inline Type at(unsigned int i) const {
743 return data.at(i);
744 }
745
746
747 using iterator = typename Container<Type>::iterator;
748
749
752 inline auto begin() {
753 return data.begin();
754 }
755
756
759 inline auto begin() const {
760 return data.cbegin();
761 }
762
763
766 inline auto end() {
767 return data.end();
768 }
769
770
773 inline auto end() const {
774 return data.cend();
775 }
776
777
779 inline void normalize() {
781 }
782
783
785 inline vec<Type> normalized() const {
786 return algebra::normalize(*this);
787 }
788
789
791 template<typename Vector>
792 inline bool operator==(const Vector& other) const {
793
794 if(size() != other.size())
795 return false;
796
797 for (unsigned int i = 0; i < size(); ++i)
798 if(data[i] != other[i])
799 return false;
800
801 return true;
802 }
803
804
806 template<typename Vector>
807 inline bool operator!=(const Vector& other) const {
808 return !(*this == other);
809 }
810
811
813 inline TH_CONSTEXPR unsigned int size() const {
814 return data.size();
815 }
816
817
819 inline void resize(size_t n) {
820 data.resize(n);
821 }
822
823
826 inline void push(const Type& x) {
827 data.push_back(x);
828 }
829
830
833 inline void push(Type&& x) {
834 data.push_back(x);
835 }
836
837
840 inline static vec<Type> euclidean_base(
841 unsigned int i, unsigned int n) {
842
843 if(i >= n) {
844 TH_MATH_ERROR("vec::euclidean_base", i, INVALID_ARGUMENT);
845 return vec<Type>(n, nan());
846 }
847
848 vec<Type> e_i = vec<Type>(n, Type(0.0));
849 e_i.resize(n);
850 e_i[i] = 1;
851
852 return e_i;
853 }
854
855
858 inline friend vec<Type> operator*(Type a, const vec<Type>& v) {
859 return v * a;
860 }
861
862
863#ifndef THEORETICA_NO_PRINT
864
866 inline std::string to_string(
867 const std::string& separator = ", ",
868 bool parenthesis = true) const {
869
870 std::stringstream res;
871
872 if(parenthesis)
873 res << "(";
874
875 for (unsigned int i = 0; i < size(); ++i) {
876 res << data[i];
877 if(i != size() - 1)
878 res << separator;
879 }
880
881 if(parenthesis)
882 res << ")";
883
884 return res.str();
885 }
886
887
889 inline operator std::string() {
890 return to_string();
891 }
892
893
895 inline friend std::ostream& operator<<(std::ostream& out, const vec<Type>& obj) {
896 return out << obj.to_string();
897 }
898
899#endif
900
901 };
902
903
910 template<typename ElementType, typename Type, typename ...Args>
911 void make_vec(vec<ElementType>& v, size_t index, Type last) {
912 v[index] = last;
913 }
914
915
924 template<typename ElementType, typename Type, typename ...Args>
925 void make_vec(vec<ElementType>& v, size_t index, Type first, Args... elements) {
926
927 v[index] = first;
928 make_vec<ElementType>(v, index + 1, elements...);
929 }
930
931
934 template<typename Type, typename ...Args>
935 vec<Type> make_vec(Type first, Args... elements) {
936
937 vec<Type> v;
938 v.resize(sizeof...(elements) + 1);
939
940 v[0] = first;
941 make_vec<Type>(v, 1, elements...);
942
943 return v;
944 }
945
946}
947
948
949#endif
Linear algebra routines.
A sequential iterator for traversing vector-like containers.
Definition vec.h:29
ReturnType operator*()
Dereference the iterator to get the current element.
Definition vec.h:52
vec_iterator(Vector &vector, size_t index)
Construct the iterator from a pointer to the elements and a starting index.
Definition vec.h:48
size_t index()
Move to the previous element in the vector.
Definition vec.h:70
vec_iterator & operator++()
Move to the next element in the vector.
Definition vec.h:57
bool operator==(const vec_iterator &other) const
Comparison operators.
Definition vec.h:76
A statically allocated N-dimensional vector with elements of the given type.
Definition vec.h:92
const Type & operator[](unsigned int i) const
Get the i-th component by value.
Definition vec.h:306
vec< Type, N > cross(const vec< Type, N > &other) const
Cross product between vectors.
Definition vec.h:226
Type & operator[](unsigned int i)
Access i-th component by reference.
Definition vec.h:300
bool operator!=(const Vector &other) const
Check whether all elements of both vectors are unequal.
Definition vec.h:406
Type operator*(const Vector &other) const
Dot product between vectors (v * w = v.x * w.x + ...)
Definition vec.h:220
auto begin() const
Get a const iterator to the first element of the vector.
Definition vec.h:365
vec_iterator< vec< Type, N >, Type & > iterator
Sequential iterator for statically allocated vectors.
Definition vec.h:344
void resize(size_t n) const
Compatibility function to allow for allocation or resizing of dynamic vectors.
Definition vec.h:422
vec()
Construct a vector with all elements equal to zero.
Definition vec.h:104
auto end()
Get an iterator to one plus the last element of the vector.
Definition vec.h:359
static vec< Type, N > euclidean_base(unsigned int i, unsigned int n=N)
Returns an N-dimensional euclidean base unit vector with the i-th element set to 1.
Definition vec.h:432
vec(const vec< Type, M > &other)
Copy constructor.
Definition vec.h:132
vec_iterator< const vec< Type, N >, const Type & > const_iterator
Const sequential iterator for statically allocated vectors.
Definition vec.h:347
vec(unsigned int size, Type val)
Construct a vector with all elements equal to the given value, checking that the given size matches t...
Definition vec.h:121
vec< Type, N > & operator*=(Type scalar)
Multiply the vector itself by a scalar.
Definition vec.h:262
Type sqr_norm() const
Compute the square norm of the vector (v * v)
Definition vec.h:294
std::string to_string(const std::string &separator=", ", bool parenthesis=true) const
Convert the vector to string representation.
Definition vec.h:458
vec< Type, N > & operator+=(const Vector &other)
Sum a vector the the vector itself.
Definition vec.h:241
vec< Type, N > operator-(const vec< Type, N > &other) const
Vector subtraction.
Definition vec.h:181
TH_CONSTEXPR unsigned int size() const
Returns the size of the vector (N)
Definition vec.h:412
vec< Type, N > & operator-=(const Vector &other)
Subtract a vector the the vector itself.
Definition vec.h:252
vec< Type, N > cross(const Vector &other) const
Cross product between vectors.
Definition vec.h:234
vec< Type, N > & operator=(const Vector &other)
Copy from other.
Definition vec.h:139
Type norm() const
Compute the norm of the vector (sqrt(v * v))
Definition vec.h:288
friend vec< Type, N > operator*(Type a, const vec< Type, N > &v)
Friend operator to enable equations of the form (Type) * (vec)
Definition vec.h:450
void normalize()
Vector normalization (v / |v|)
Definition vec.h:378
vec(std::initializer_list< Type > l)
Initialize from a list, e.g. {1, 2, 3}.
Definition vec.h:145
Type dot(const Vector &other) const
Dot product between vectors (v * w = v.x * w.x + ...)
Definition vec.h:213
Type & at(unsigned int i)
Access i-th element by reference, with bound checking.
Definition vec.h:315
vec< Type, N > normalized() const
Return the normalized vector (v / |v|)
Definition vec.h:384
auto begin()
Get an iterator to the first element of the vector.
Definition vec.h:352
vec< Type, N > operator+(const vec< Type, N > &other) const
Vector sum (v + w = (v.x + w.x, ...))
Definition vec.h:166
friend std::ostream & operator<<(std::ostream &out, const vec< Type, N > &obj)
Stream the vector in string representation to an output stream (std::ostream)
Definition vec.h:487
bool operator==(const Vector &other) const
Check whether all elements of both vectors are equal.
Definition vec.h:391
auto end() const
Get a const iterator to one plus the last element of the vector.
Definition vec.h:372
vec< Type, N > operator/(Type scalar) const
Scalar division (v / a = (v.x / a, ...))
Definition vec.h:201
Type at(unsigned int i) const
Get the i-th element by value, with bound checking.
Definition vec.h:331
vec(Type val)
Construct a vector with all elements equal to the given value.
Definition vec.h:111
vec< Type, N > & operator/=(Type scalar)
Divide the vector itself by a scalar.
Definition vec.h:272
vec< Type, N > operator+() const
Identity.
Definition vec.h:160
vec< Type, N > operator*(Type scalar) const
Scalar multiplication (av = (v.x * a, ...))
Definition vec.h:190
vec< Type, N > operator-() const
Opposite vector.
Definition vec.h:175
#define TH_CONSTEXPR
Enable constexpr in function declarations if C++14 is supported.
Definition constants.h:161
#define TH_MATH_ERROR(F_NAME, VALUE, EXCEPTION)
TH_MATH_ERROR is a macro which throws exceptions or modifies errno (depending on which compiling opti...
Definition error.h:225
Vector1 cross(const Vector1 &v1, const Vector2 &v2)
Compute the cross product between two tridimensional vectors.
Definition algebra.h:373
Vector & make_normalized(Vector &v)
Normalize a given vector overwriting it.
Definition algebra.h:327
Vector1 & vec_copy(Vector1 &dest, const Vector2 &src)
Copy a vector by overwriting another.
Definition algebra.h:143
auto dot(const Vector1 &v, const Vector2 &w)
Computes the dot product between two vectors, using the Hermitian form if needed.
Definition algebra.h:351
auto sqr_norm(const Vector &v)
Returns the square of the Euclidean/Hermitian norm of the given vector.
Definition algebra.h:275
Vector2 & vec_sum(Vector1 &v1, const Vector2 &v2)
Sum two vectors and store the result in the first vector.
Definition algebra.h:1135
Vector2 & vec_diff(Vector1 &v1, const Vector2 &v2)
Subtract two vectors and store the result in the first vector.
Definition algebra.h:1183
Vector & vec_zeroes(Vector &v)
Overwrite a vector with all zeroes.
Definition algebra.h:109
Vector & vec_error(Vector &v)
Overwrite the given vector with the error vector with NaN values.
Definition algebra.h:62
auto norm(const Vector &v)
Returns the Euclidean/Hermitian norm of the given vector.
Definition algebra.h:292
Vector normalize(const Vector &v)
Returns the normalized vector.
Definition algebra.h:302
Main namespace of the library which contains all functions and objects.
Definition algebra.h:27
dual2 abs(dual2 x)
Compute the absolute value of a second order dual number.
Definition dual2_functions.h:198
auto max(const Vector &X)
Finds the maximum value inside a dataset.
Definition dataset.h:330
TH_CONSTEXPR real nan()
Return a quiet NaN number in floating point representation.
Definition error.h:54
void make_vec(vec< ElementType > &v, size_t index, Type last)
Populates a vector with a single element at the specified index.
Definition vec.h:911
constexpr real MACH_EPSILON
Machine epsilon for the real type.
Definition constants.h:207