Hosting by Solid Eight Studios ,PhotoTangler Collage Maker

This section of the archives stores flipcode's complete Developer Toolbox collection, featuring a variety of mini-articles and source code contributions from our readers.

 Matrix Template with Rare Features   Submitted by

This is a much improved version of a matrix library that I previously posted to COTD. The previous version supported potentially ambiguous expressions by having operator* represent both scalar and vector product, generating the longest recorded discussion for a COTD. That feature have been removed, the more conventional dot() and cross() are now used. In addition to the temporary-free evaluation that was present in the previous version, there have been performance improvements and some new features and SSE support have been added. One new feature that I would assume most of you have never seen before is "fancy indexing". New "fake" matrices may be composed of elements from an existing matrix.

 ```matrix<4 A; vector<4 col(1,0,0,0);A[c0] = col; // set the first column to 1,0,0,0 A[r1] << 0,1,0,0; // set second row to 0,1,0,0 A[r12,c12] = zero; // set the center four elements to zero ```

The same infrastructure is used to enable swizzling of vector elements, transpose and cofactor operations.
 ```matrix<1,4 transvec; vector<4 vec; vector<3 vec3;vec = transpose(transvec); vec = transvec[tr]; // same as above vec[tr] = transvec; // also the same vec3 = vec[xyw]; // select elements vec[xyw] = vec3; // works both ways vec = vec3[z|x|y|y]; // custom swizzling ```

Natural operator syntax is still used for the most common operations.
 ```matrix<4 A, B, C; vector<4 u, v, w; // matrix multiplication (vectors are column matrices) A *= B; C = A * B; u = A * v;// elementwise operations for vectors u *= v; w = v + u; u /= 5; ```

To clear out this before the discussion starts: The operator syntax does not require any extra (nontrivial) temporaries and is as efficient as the ugly but apparently popular

 `C.MatrixMult(A, B); // not supported, not needed `

The only quirk is that you need a decent compiler. SSE is now supported when using GCC 3.1, more than doubling the performance for simple expressions. Without SSE, performance has been increased by 20% since the last COTD version. I would appreciate if you people helped me check which compilers can handle the code, extensive support for templates is needed. All I know is that GCC 3 works. Don't even bother with Visual C++ 5 and 6, the last time I tried they choked within the first 50 lines of the declaration header. If somebody reports that Intel C++ works, I will make sure SSE will be supported there as well. The source code is public domain. The latest version will be available at the SourceForge project page. (Not up to date yet.)

Currently browsing [mathlib_cotd.zip] (40,608 bytes) - [mathlib/include/mathlib/box.h] - (603 bytes)

 ```#ifndef MATHLIB_BOX_H #define MATHLIB_BOX_H#include namespace mathlib {// axis aligned bounding box template class box { public: vector min; vector max; box() {} box(const vector &min_, const vector &max_) : min(min_), max(max_) {} box &operator|=(const vector &v) { for (int i=0; i!=D; ++i) if (v(i) < min(i)) min(i) = v(i); else if (v(i) > max(i)) max(i) = v(i); return *this; } };template inline box operator|(const box &b, const vector &v) { box r = b; return r |= v; }} // namespace mathlib #endif ```

Currently browsing [mathlib_cotd.zip] (40,608 bytes) - [mathlib/include/mathlib/complex.h] - (163 bytes)

 ```#ifndef MATHLIB_COMPLEX_H #define MATHLIB_COMPLEX_H#include namespace mathlib {typedef std::complex complex;} // namespace mathlib #endif ```

Currently browsing [mathlib_cotd.zip] (40,608 bytes) - [mathlib/include/mathlib/line.h] - (295 bytes)

 ```#ifndef MATHLIB_LINE_H #define MATHLIB_LINE_H#include namespace mathlib {template class line { public: vector base; vector direction; line() {} line(const vector &r, const vector &v) : base(r), direction(v) {} };} // namespace mathlib #endif ```

Currently browsing [mathlib_cotd.zip] (40,608 bytes) - [mathlib/include/mathlib/mathconf.h] - (708 bytes)

 ```#ifndef MATHLIB_MATHCONF_H #define MATHLIB_MATHCONF_H// detects defects and features of the compiler // g++ #if __GNUG__ // unaliased references are supported in g++, use them # define MATHLIB_RESTRICT __restrict # if __GNUG__ < 3 # ifndef MATHLIB_NO_USING # define MATHLIB_NO_USING 1 # endif # ifndef MATHLIB_LIBSTDCPP2 # define MATHLIB_LIBSTDCPP2 1 # endif # endif #endif// work around libstdc++-v2 defects #if MATHLIB_LIBSTDCPP2 # ifndef MATHLIB_NO_LIMITS # define MATHLIB_NO_LIMITS 1 # endif # ifndef MATHLIB_NO_OSTREAM # define MATHLIB_NO_OSTREAM 1 # endif #endif// unaliased references does not appear to be supported #ifndef MATHLIB_RESTRICT # define MATHLIB_RESTRICT #endif#endif ```

Currently browsing [mathlib_cotd.zip] (40,608 bytes) - [mathlib/include/mathlib/mathdefs.h] - (697 bytes)

 ```#ifndef MATHLIB_MATHDEFS_H #define MATHLIB_MATHDEFS_H// constants #include #if MATHLIB_NO_LIMITS # include #else # include #endifnamespace mathlib {const float pi = 3.14159265358979323846; const float ln_2 = 0.693147180559945309417; #if MATHLIB_NO_LIMITS const float epsilon = FLT_EPSILON; #else const float epsilon = std::numeric_limits::epsilon(); #endifstruct identity_t {}; struct zero_t { operator float() const { return 0.f; } };// absence of geometry. struct nothing {};// typeless multiplicative identity const identity_t identity = identity_t();// typeless zero const zero_t zero = zero_t();} // namespace mathlib #endif ```

Currently browsing [mathlib_cotd.zip] (40,608 bytes) - [mathlib/include/mathlib/matrix.h] - (10,908 bytes)

 ```#ifndef MATHLIB_MATRIX_H #define MATHLIB_MATRIX_H// class definitions for matrix and vector #include #include #include namespace mathlib {class float_keeper; template struct assign_dispatch; template struct matrix_indexer; template struct matrix_unop; template struct unop_index; template struct unop_cindex; template class tmatrix; template class matrix;// representation of matrices template class matrix_storage { float m_data[Cols][Rows]; public: float &operator()(int r, int c) { return m_data[c][r]; } float operator()(int r, int c) const { return m_data[c][r]; } float *data() { return &m_data[0][0]; } const float *data() const { return &m_data[0][0]; } }; template<> class matrix_storage<3, 3> { float m_data[3][3]; public: float &operator()(int r, int c) { return m_data[c][r]; } float operator()(int r, int c) const { return m_data[c][r]; } float *data() { return &m_data[0][0]; } const float *data() const { return &m_data[0][0]; } #if MATHLIB_SIMD4 private: float m_filler; #endif } MATHLIB_SIMD4_ALIGN; template<> class matrix_storage<4, 4> { float m_data[4][4]; public: float &operator()(int r, int c) { return m_data[c][r]; } float operator()(int r, int c) const { return m_data[c][r]; } float *data() { return &m_data[0][0]; } const float *data() const { return &m_data[0][0]; } } MATHLIB_SIMD4_ALIGN;// specializations for vectors template class matrix_storage { float m_data[Rows]; public: float &operator()(int r, int) { return m_data[r]; } float operator()(int r, int) const { return m_data[r]; } float &operator()(int r) { return m_data[r]; } float operator()(int r) const { return m_data[r]; } float *data() { return m_data; } const float *data() const { return m_data; } }; template<> class matrix_storage<2, 1> { public: union { float x; float s; }; union { float y; float t; }; float &operator()(int r, int) { return (&x)[r]; } float operator()(int r, int) const { return (&x)[r]; } float &operator()(int r) { return (&x)[r]; } float operator()(int r) const { return (&x)[r]; } float *data() { return &x; } const float *data() const { return &x; } void assign(float a, float b) { x = a; y = b; } }; template<> class matrix_storage<3, 1> { public: union { float x; float s; float r; }; union { float y; float t; float g; }; union { float z; float u; float b; }; float &operator()(int r, int) { return (&x)[r]; } float operator()(int r, int) const { return (&x)[r]; } float &operator()(int r) { return (&x)[r]; } float operator()(int r) const { return (&x)[r]; } float *data() { return &x; } const float *data() const { return &x; } void assign(float a, float b, float c) { x = a; y = b; z = c; }#if MATHLIB_SIMD4 matrix_storage() {} matrix_storage(const matrix_storage &); matrix_storage &operator=(const matrix_storage &); private: float m_filler; #endif } MATHLIB_SIMD4_ALIGN; template<> class matrix_storage<4, 1> { public: union { float x; float s; float r; }; union { float y; float t; float g; }; union { float z; float u; float b; }; union { float w; float v; float a; }; float &operator()(int r, int) { return (&x)[r]; } float operator()(int r, int) const { return (&x)[r]; } float &operator()(int r) { return (&x)[r]; } float operator()(int r) const { return (&x)[r]; } float *data() { return &x; } const float *data() const { return &x; } void assign(float a, float b, float c, float d) { x = a; y = b; z = c; w = d; }#if MATHLIB_SIMD4 matrix_storage() {} matrix_storage(const matrix_storage &); matrix_storage &operator=(const matrix_storage &); #endif } MATHLIB_SIMD4_ALIGN;// the type most matrix operations work with. subclass tmatrix // to make your own matrix type. do not use basic_matrix directly template class basic_matrix : public matrix_storage { public: static const int rows = Rows; static const int columns = Cols; static const bool is_storage = true; static const bool is_lvalue = true; static const bool error = false; static const bool element_op = true; typedef basic_matrix base_type; basic_matrix() {} // used for implicit evaluation by the select_* templates in // matrix_operators.tpp template basic_matrix(const T &m); template float element() const { enum { i = R+C*Rows }; return data()[i]; } template float &element() { enum { i = R+C*Rows }; return data()[i]; } mathlib::matrix &matrix() { return *static_cast *>(this); } const mathlib::matrix &matrix() const { return *static_cast *>(this); } };// should be wrapped around matrix expressions that are strictly for // initialization to avoid having to write two constructor/assignment // pairs per initializer template class tmatrix_ndim : public Base { public: tmatrix_ndim() {} template explicit tmatrix_ndim(const A &a) : Base(a) {} template tmatrix_ndim(const A &a, const B &b) : Base(a, b) {} };// wrapped around matrix expressions by tmatrix_type to make overloading based // on size possible. the type most operators accept template class tmatrix_dim : public Base { public: tmatrix_dim() {} template explicit tmatrix_dim(A &a) : Base(a) {} template explicit tmatrix_dim(const A &a) : Base(a) {} template tmatrix_dim(const A &a, const B &b) : Base(a, b) {} tmatrix_dim &operator*=(float); tmatrix_dim &operator/=(float); tmatrix_dim &operator=(zero_t); tmatrix_dim &operator=(identity_t); template tmatrix_dim &operator=(const tmatrix_dim &); template tmatrix_dim &operator=(const tmatrix_ndim &); template tmatrix_dim &operator*=(const tmatrix_dim &); template tmatrix_dim &operator/=(const tmatrix_dim &); template tmatrix_dim &operator*=(const tmatrix_dim &); template tmatrix_dim &operator+=(const tmatrix_dim &); template tmatrix_dim &operator-=(const tmatrix_dim &); template tmatrix, Base> > operator[](matrix_indexer); template const tmatrix, Base> > operator[](matrix_indexer) const; };// wrapped around expressions by tmatrix. used for early detection of // illegal operations template class tmatrix_type { // error - cannot be used }; template class tmatrix_type : public tmatrix_dim { typedef tmatrix_dim m_base; public: // not scalar, is l-value tmatrix_type() {} template explicit tmatrix_type(A &a) : m_base(a) {} template explicit tmatrix_type(const A &a) : m_base(a) {} template tmatrix_type(const A &a, const B &b) : m_base(a, b) {}#if MATHLIB_NO_USING template tmatrix_type &operator=(const A &a) { m_base::operator=(a); return *this; } #else using m_base::operator=; #endif }; template class tmatrix_type : public tmatrix_dim { typedef tmatrix_dim m_base; public: // not scalar, not l-value tmatrix_type() {} template explicit tmatrix_type(A &a) : m_base(a) {} template explicit tmatrix_type(const A &a) : m_base(a) {} template tmatrix_type(const A &a, const B &b) : m_base(a, b) {} }; template class tmatrix_type : public Base { typedef Base m_base; public: // scalar, l-value. disable matrix operations tmatrix_type() {} template explicit tmatrix_type(A &a) : m_base(a) {} template explicit tmatrix_type(const A &a) : m_base(a) {} template tmatrix_type(const A &a, const B &b) : m_base(a, b) {} float &operator=(float f) { return Base::template element<0,0>() = f; } #if MATHLIB_NO_USING template float &operator=(const A &f) { float t = const_cast(f); return *this = t; } #endif operator float &() { return Base::template element<0,0>(); } }; template class tmatrix_type : public Base { typedef Base m_base; public: // scalar, not l-value. disable matrix operations tmatrix_type() {} template explicit tmatrix_type(A &a) : m_base(a) {} template explicit tmatrix_type(const A &a) : m_base(a) {} template tmatrix_type(const A &a, const B &b) : m_base(a, b) {} operator float() const { basic_matrix<1,1> t(*this); return *t.data(); } };// selects superclasses to provide (only) the supported operations for // expressions template class tmatrix : public tmatrix_type { typedef tmatrix_type m_base; public: tmatrix() {} template explicit tmatrix(A &a) : m_base(a) {} template explicit tmatrix(const A &a) : m_base(a) {} template tmatrix(const A &a, const B &b) : m_base(a, b) {}#if MATHLIB_NO_USING template tmatrix &operator=(const A &a) { m_base::operator=(a); return *this; } #else using m_base::operator=; #endif };// friendly matrix template class matrix : public tmatrix > { typedef tmatrix > m_base; public: matrix() {} matrix(zero_t); matrix(identity_t); template matrix(const tmatrix_dim &); template matrix(const tmatrix_ndim &);#if MATHLIB_NO_USING template matrix &operator=(const A &a) { m_base::operator=(a); return *this; } #else using m_base::operator=; #endif };// friendly column matrix template class vector : public tmatrix > { typedef tmatrix > m_base; public: vector() {} explicit vector(float filler); vector(float x, float y) { assign(x, y); } vector(float x, float y, float z) { assign(x, y, z); } vector(float x, float y, float z, float w) { assign(x, y, z, w); } vector(zero_t); template vector(const tmatrix_dim &); template vector(const tmatrix_ndim &);#if MATHLIB_NO_USING template vector &operator=(const A &a) { m_base::operator=(a); return *this; } #else using m_base::operator=; #endif };} // namespace mathlib #endif ```

Currently browsing [mathlib_cotd.zip] (40,608 bytes) - [mathlib/include/mathlib/plane.h] - (593 bytes)

 ```#ifndef MATHLIB_PLANE_H #define MATHLIB_PLANE_H#include namespace mathlib {template class plane { public: vector n; float D; plane() {} plane(const vector &n_, float d) : n(n_), D(d) {} };template class parameter_plane { public: vector v, u; parameter_plane(const vector &v_, const vector &u_) : v(v_), u(u_) {} parameter_plane(const plane &); operator plane() const; };template inline float alignment(const plane &p, const vector &v) { return dot(p.n, v) + p.D; }} // namespace mathlib #endif ```

Currently browsing [mathlib_cotd.zip] (40,608 bytes) - [mathlib/include/mathlib/quaternion.h] - (1,751 bytes)

 ```#ifndef MATHLIB_QUATERNION_H #define MATHLIB_QUATERNION_H#include namespace mathlib {class quaternion { public: float a, b, c, d; quaternion() {} quaternion(float a_, float b_ = 0.f, float c_ = 0.f, float d_ = 0.f) : a(a_), b(b_), c(c_), d(d_) {} float real() const { return a; } quaternion unreal() const { return quaternion(0.f, b, c, d); } quaternion &operator+=(float f) { a += f; return *this; } quaternion &operator-=(float f) { a -= f; return *this; } quaternion &operator*=(float f) { a *= f; b *= f; c *= f; d *= f; return *this; } quaternion &operator/=(float f) { a /= f; b /= f; c /= f; d /= f; return *this; } quaternion &operator+=(complex c) { a += c.real(); b += c.imag(); return *this; } quaternion &operator-=(complex c) { a -= c.real(); b -= c.imag(); return *this; } quaternion &operator*=(complex); quaternion &operator/=(complex); quaternion &operator+=(const quaternion &q) { a += q.a; b += q.b; c += q.c; d += q.d; return *this; } quaternion &operator-=(const quaternion &q) { a -= q.a; b -= q.b; c -= q.c; d -= q.d; return *this; } quaternion &operator*=(const quaternion &q) { return *this = *this * q; } quaternion &operator/=(const quaternion &q) { return *this = *this / q; } friend quaternion operator*(const quaternion &, const quaternion &); friend quaternion operator/(const quaternion &, const quaternion &); friend float real(const quaternion &q) { return q.real(); } friend quaternion unreal(const quaternion &q) { return q.unreal(); } friend float abs(const quaternion &q) { return sqrt(q.a*q.a + q.b*q.b + q.c*q.c + q.d*q.d); } friend quaternion normalize(const quaternion &q) { return q / abs(q); } };} // namespace mathlib #endif ```

Currently browsing [mathlib_cotd.zip] (40,608 bytes) - [mathlib/include/mathlib/segment.h] - (242 bytes)

 ```#ifndef MATHLIB_SEGMENT_H #define MATHLIB_SEGMENT_Hnamespace mathlib {template class segment { public: vector v1, v2; segment() {} segment(const vector &a, const vector &b) : v1(a), v2(b) {} };} // namespace mathlib ```

Currently browsing [mathlib_cotd.zip] (40,608 bytes) - [mathlib/include/mathlib/simd/config.h] - (517 bytes)

 ```#ifndef MATHLIB_SIMD_CONFIG_H #define MATHLIB_SIMD_CONFIG_H// type definitions and defines for simd instructions namespace mathlib {#if defined(__GNUG__) && defined(__SSE__) || 1 // g++ with sse instructions enabled # define MATHLIB_SIMD 1 # define MATHLIB_SIMD4 1 # define MATHLIB_GCC_SSE 1# define MATHLIB_SIMD4_ALIGN __attribute__((aligned(__alignof__(v4sf))))// vector of 4 floats typedef int v4sf __attribute__ ((mode(V4SF)));#else # define MATHLIB_SIMD4_ALIGN #endif} // namespace mathlib #endif ```

Currently browsing [mathlib_cotd.zip] (40,608 bytes) - [mathlib/include/mathlib/sphere.h] - (294 bytes)

 ```#ifndef MATHLIB_SPHERE_H #define MATHLIB_SPHERE_H#include namespace mathlib {template class sphere { public: vector position; float radius; sphere() {} sphere(const vector &pos, float r) : position(pos), radius(r) {} };} // namespace mathlib #endif ```

Currently browsing [mathlib_cotd.zip] (40,608 bytes) - [mathlib/include/mathlib/triangle.h] - (317 bytes)

 ```#ifndef MATHLIB_TRIANGLE_H #define MATHLIB_TRIANGLE_H#include namespace mathlib {template class triangle { public: vector a, b, c; triangle() {} triangle(const vector &a_, const vector &b_, const vector &c_) : a(a_), b(b_), c(c_) {} };} // namespace mathlib #endif ```

Currently browsing [mathlib_cotd.zip] (40,608 bytes) - [mathlib/include/mathlib/utility.h] - (405 bytes)

 ```#ifndef MATHLIB_UTILITY_H #define MATHLIB_UTILITY_H// random mathematics functions #include #include namespace mathlib {inline float sqr(float a) { return a*a; }inline bool close(float a, float b) { return std::abs(a - b) < epsilon; }inline int sign(float a) { if (a > epsilon) return 1; if (a < -epsilon) return -1; return 0; }} // namespace mathlib #endif ```

Currently browsing [mathlib_cotd.zip] (40,608 bytes) - [mathlib/src/matrix2.cpp] - (407 bytes)

 ```#include namespace mathlib {typedef basic_matrix<2, 2> m2;void matrix_assign_handler::go(m2 &dest, matrix_unop src) { float det_src = det_bm(src.a); float d = 1.f / det_src; // dest = adjoint(src.a) / det(src.a) dest(0,0) = d * +src.a(1,1); dest(1,0) = d * -src.a(1,0); dest(0,1) = d * -src.a(0,1); dest(1,1) = d * +src.a(0,0); }} // namespace mathlib ```

Currently browsing [mathlib_cotd.zip] (40,608 bytes) - [mathlib/src/matrix3.cpp] - (1,380 bytes)

 ```#include namespace mathlib {typedef basic_matrix<3, 3> m3; typedef basic_matrix<3, 1> v3;float det_bm(const m3 &m) { return m(0,0) * det(cofactor_bm<0,0>(m)) - m(1,0) * det(cofactor_bm<1,0>(m)) + m(2,0) * det(cofactor_bm<2,0>(m)); }void matrix_assign_handler::go(m3 &dest, matrix_unop src) { dest(0,0) = +det(cofactor_bm<0,0>(src.a)); dest(1,0) = -det(cofactor_bm<0,1>(src.a)); dest(2,0) = +det(cofactor_bm<0,2>(src.a)); dest(0,1) = -det(cofactor_bm<1,0>(src.a)); dest(1,1) = +det(cofactor_bm<1,1>(src.a)); dest(2,1) = -det(cofactor_bm<1,2>(src.a)); dest(0,2) = +det(cofactor_bm<2,0>(src.a)); dest(1,2) = -det(cofactor_bm<2,1>(src.a)); dest(2,2) = +det(cofactor_bm<2,2>(src.a)); }void matrix_assign_handler::go(m3 &dest, matrix_unop src) { float det_src = det_bm(src.a); float d = 1.f / det_src; // dest = adjoint(src.a) / det(src.a) dest(0,0) = d * +det(cofactor_bm<0,0>(src.a)); dest(1,0) = d * -det(cofactor_bm<0,1>(src.a)); dest(2,0) = d * +det(cofactor_bm<0,2>(src.a)); dest(0,1) = d * -det(cofactor_bm<1,0>(src.a)); dest(1,1) = d * +det(cofactor_bm<1,1>(src.a)); dest(2,1) = d * -det(cofactor_bm<1,2>(src.a)); dest(0,2) = d * +det(cofactor_bm<2,0>(src.a)); dest(1,2) = d * -det(cofactor_bm<2,1>(src.a)); dest(2,2) = d * +det(cofactor_bm<2,2>(src.a)); }} // namespace mathlib ```

Currently browsing [mathlib_cotd.zip] (40,608 bytes) - [mathlib/src/matrix4.cpp] - (5,164 bytes)

 ```#include namespace mathlib {const matrix<4> frustum(float l, float r, float b, float t, float n, float f) { float n2 = 2 * n; float rl = r - l; float tb = t - b; float fn = f - n; matrix<4> dest = zero; dest(0,0) = n2/rl; dest(0,2) = (r+l)/rl; dest(1,1) = n2/tb; dest(1,2) = (t+b)/tb; dest(2,2) = -(f+n)/fn; dest(2,3) = -2*f*n/fn; dest(3,2) = -1; return dest; }const matrix<4> ortho(float l, float r, float b, float t, float n, float f) { float rl = r - l; float tb = t - b; float fn = f - n; matrix<4> dest = zero; dest(0,0) = 2/rl; dest(0,3) = (r+l)/rl; dest(1,1) = 2/tb; dest(1,3) = -(t+b)/tb; dest(2,2) = -2/fn; dest(2,3) = (f+n)/fn; dest(3,3) = 1; return dest; }void matrix_assign_handler::go(bm44 &dest, matrix_unop m) { dest(0,3) = m.a.x; dest(1,3) = m.a.y; dest(2,3) = m.a.z; dest(0,1) = dest(0,2) = 0; dest(1,0) = dest(1,2) = 0; dest(2,0) = dest(2,1) = 0; dest(3,0) = dest(3,1) = dest(3,2) = 0; dest(0,0) = dest(1,1) = dest(2,2) = dest(3,3) = 1; }void matrix_assign_handler::go(bm44 &dest, matrix_unop src) { dest(0,0) = src.a.x; dest(1,1) = src.a.y; dest(2,2) = src.a.z; dest(3,3) = 1; dest(0,1) = dest(0,2) = dest(0,3) = 0; dest(1,0) = dest(1,2) = dest(1,3) = 0; dest(2,0) = dest(2,1) = dest(2,3) = 0; dest(3,0) = dest(3,1) = dest(3,2) = 0; }void matrix_assign_handler::go(bm44 &dest, matrix_unop src) { dest(0,0) = dest(1,1) = dest(2,2) = src.a.f; dest(3,3) = 1; dest(0,1) = dest(0,2) = dest(0,3) = 0; dest(1,0) = dest(1,2) = dest(1,3) = 0; dest(2,0) = dest(2,1) = dest(2,3) = 0; dest(3,0) = dest(3,1) = dest(3,2) = 0; }void matrix_assign_handler::go(bm44 &dest, matrix_binop src) { vector<3> u = norm(src.a.matrix()); matrix<3> S; S << 0, -u.z, u.y, u.z, 0, -u.x, -u.y, u.x, 0; matrix<3> uut = u*transpose(u); matrix<3> I = identity; matrix<4> & MATHLIB_RESTRICT R = dest.matrix(); float a = src.b.f; R[r02,c02] = uut + cos(a)*(I - uut) + sin(a)*S; R[r3,c02] = zero; R[r02,c3] = zero; R(3,3) = 1; }void matrix_assign_handler::go(bm44 &dest, matrix_binop src) { float c = cos(src.b.f), s = sin(src.b.f); dest(1,1) = c; dest(2,2) = c; dest(1,2) = -s; dest(2,1) = s; dest(0,1) = dest(0,2) = 0; dest(1,0) = dest(2,0) = 0; dest(0,0) = 1; dest(0,3) = dest(1,3) = dest(2,3) = dest(3,0) = dest(3,1) = dest(3,2) = 0; dest(3,3) = 1; }void matrix_assign_handler::go(bm44 &dest, matrix_binop src) { float c = cos(src.b.f), s = sin(src.b.f); dest(0,0) = c; dest(2,2) = c; dest(2,0) = -s; dest(0,2) = s; dest(0,1) = dest(1,2) = 0; dest(1,0) = dest(2,1) = 0; dest(1,1) = 1; dest(0,3) = dest(1,3) = dest(2,3) = dest(3,0) = dest(3,1) = dest(3,2) = 0; dest(3,3) = 1; }void matrix_assign_handler::go(bm44 &dest, matrix_binop src) { float c = cos(src.b.f), s = sin(src.b.f); dest(0,0) = c; dest(1,1) = c; dest(0,1) = -s; dest(1,0) = s; dest(0,2) = dest(1,2) = 0; dest(2,0) = dest(2,1) = 0; dest(2,2) = 1; dest(0,3) = dest(1,3) = dest(2,3) = dest(3,0) = dest(3,1) = dest(3,2) = 0; dest(3,3) = 1; }float det_bm(const bm44 &m) { return m(0,0) * det(cofactor_bm<0,0>(m)) - m(1,0) * det(cofactor_bm<1,0>(m)) + m(2,0) * det(cofactor_bm<2,0>(m)) - m(3,0) * det(cofactor_bm<3,0>(m)); }void matrix_assign_handler::go(bm44 &dest, matrix_unop src) { matrix<4> & MATHLIB_RESTRICT R = dest.matrix(); const matrix<4> & MATHLIB_RESTRICT S = src.a.matrix(); R(0,0) = det(cofactor<0,0>(S)); R(1,0) = -det(cofactor<0,1>(S)); R(2,0) = det(cofactor<0,2>(S)); R(3,0) = -det(cofactor<0,3>(S)); R(0,1) = -det(cofactor<1,0>(S)); R(1,1) = det(cofactor<1,1>(S)); R(2,1) = -det(cofactor<1,2>(S)); R(3,1) = det(cofactor<1,3>(S)); R(0,2) = det(cofactor<2,0>(S)); R(1,2) = -det(cofactor<2,1>(S)); R(2,2) = det(cofactor<2,2>(S)); R(3,2) = -det(cofactor<2,3>(S)); R(0,3) = -det(cofactor<3,0>(S)); R(1,3) = det(cofactor<3,1>(S)); R(2,3) = -det(cofactor<3,2>(S)); R(3,3) = det(cofactor<3,3>(S)); }void matrix_assign_handler::go(bm44 &dest, matrix_unop src) { matrix<4> & MATHLIB_RESTRICT R = dest.matrix(); const matrix<4> & MATHLIB_RESTRICT S = src.a.matrix(); float det_s = det(S); float d = 1 / det_s; // dest = adjoint(src.a) / det(src.a) R(0,0) = d * det(cofactor<0,0>(S)); R(1,0) = d * -det(cofactor<0,1>(S)); R(2,0) = d * det(cofactor<0,2>(S)); R(3,0) = d * -det(cofactor<0,3>(S)); R(0,1) = d * -det(cofactor<1,0>(S)); R(1,1) = d * det(cofactor<1,1>(S)); R(2,1) = d * -det(cofactor<1,2>(S)); R(3,1) = d * det(cofactor<1,3>(S)); R(0,2) = d * det(cofactor<2,0>(S)); R(1,2) = d * -det(cofactor<2,1>(S)); R(2,2) = d * det(cofactor<2,2>(S)); R(3,2) = d * -det(cofactor<2,3>(S)); R(0,3) = d * -det(cofactor<3,0>(S)); R(1,3) = d * det(cofactor<3,1>(S)); R(2,3) = d * -det(cofactor<3,2>(S)); R(3,3) = d * det(cofactor<3,3>(S)); }} // namespace mathlib ```

Currently browsing [mathlib_cotd.zip] (40,608 bytes) - [mathlib/src/quaternion.cpp] - (680 bytes)

 ```#include namespace mathlib {quaternion operator*(const quaternion &a, const quaternion &b) { quaternion r; r.a = a.a*b.a - a.b*b.b - a.c*b.c - a.d*b.d; r.b = a.a*b.b + a.b*b.a + a.c*b.d - a.d*b.c; r.c = a.a*b.c - a.b*b.d + a.c*b.a + a.d*b.b; r.d = a.a*b.d + a.b*b.c - a.c*b.b + a.d*b.a; return r; }quaternion operator/(const quaternion &a, const quaternion &b) { quaternion r; float d = abs(b); r.a = (+a.a*b.a + a.b*b.b + a.c*b.c + a.d*b.d) / d; r.b = (-a.a*b.b + a.b*b.a - a.c*b.d + a.d*b.c) / d; r.c = (-a.a*b.c + a.b*b.d + a.c*b.a - a.d*b.b) / d; r.d = (-a.a*b.d - a.b*b.c + a.c*b.b + a.d*b.a) / d; return r; }} // namespace mathlib ```

Currently browsing [mathlib_cotd.zip] (40,608 bytes) - [mathlib/src/simd/gcc_sse.cpp] - (2,277 bytes)

 ```#include #if MATHLIB_GCC_SSE#include #include namespace mathlib {void simd_assign_handler::go(bm44 &dest, matrix_binop, bm44, bm44> src) { v4sf v, f1, f2, f3, f4; v4sf c1, c2, c3, c4; v4sf ac1, ac2, ac3, ac4; v4sf s1, s2; ac1 = ssev(src.a)[0]; ac2 = ssev(src.a)[1]; ac3 = ssev(src.a)[2]; ac4 = ssev(src.a)[3]; v = ssev(src.b)[0]; f1 = __builtin_ia32_shufps(v, v, 0x00); f2 = __builtin_ia32_shufps(v, v, 0x55); f3 = __builtin_ia32_shufps(v, v, 0xaa); f4 = __builtin_ia32_shufps(v, v, 0xff); c1 = __builtin_ia32_mulps(f1, ac1); c2 = __builtin_ia32_mulps(f2, ac2); c3 = __builtin_ia32_mulps(f3, ac3); c4 = __builtin_ia32_mulps(f4, ac4); s1 = __builtin_ia32_addps(c1, c2); s2 = __builtin_ia32_addps(c3, c4); ssev(dest)[0] = __builtin_ia32_addps(s1, s2); v = ssev(src.b)[1]; f1 = __builtin_ia32_shufps(v, v, 0x00); f2 = __builtin_ia32_shufps(v, v, 0x55); f3 = __builtin_ia32_shufps(v, v, 0xaa); f4 = __builtin_ia32_shufps(v, v, 0xff); c1 = __builtin_ia32_mulps(f1, ac1); c2 = __builtin_ia32_mulps(f2, ac2); c3 = __builtin_ia32_mulps(f3, ac3); c4 = __builtin_ia32_mulps(f4, ac4); s1 = __builtin_ia32_addps(c1, c2); s2 = __builtin_ia32_addps(c3, c4); ssev(dest)[1] = __builtin_ia32_addps(s1, s2); v = ssev(src.b)[2]; f1 = __builtin_ia32_shufps(v, v, 0x00); f2 = __builtin_ia32_shufps(v, v, 0x55); f3 = __builtin_ia32_shufps(v, v, 0xaa); f4 = __builtin_ia32_shufps(v, v, 0xff); c1 = __builtin_ia32_mulps(f1, ac1); c2 = __builtin_ia32_mulps(f2, ac2); c3 = __builtin_ia32_mulps(f3, ac3); c4 = __builtin_ia32_mulps(f4, ac4); s1 = __builtin_ia32_addps(c1, c2); s2 = __builtin_ia32_addps(c3, c4); ssev(dest)[2] = __builtin_ia32_addps(s1, s2); v = ssev(src.b)[3]; f1 = __builtin_ia32_shufps(v, v, 0x00); f2 = __builtin_ia32_shufps(v, v, 0x55); f3 = __builtin_ia32_shufps(v, v, 0xaa); f4 = __builtin_ia32_shufps(v, v, 0xff); c1 = __builtin_ia32_mulps(f1, ac1); c2 = __builtin_ia32_mulps(f2, ac2); c3 = __builtin_ia32_mulps(f3, ac3); c4 = __builtin_ia32_mulps(f4, ac4); s1 = __builtin_ia32_addps(c1, c2); s2 = __builtin_ia32_addps(c3, c4); ssev(dest)[3] = __builtin_ia32_addps(s1, s2); }} // namespace mathlib #endif // MATHLIB_GCC_SSE ```

The zip file viewer built into the Developer Toolbox made use of the zlib library, as well as the zlibdll source additions.