1 /**
2 dlsl.matrix
3 
4 
5 Authors: Peter Particle ( based on gl3n by David Herberth )
6 License: MIT
7 
8 Note: All methods marked with pure are weakly pure since, they all access an instance member.
9 All static methods are strongly pure.
10 */
11 
12 module dlsl.matrix;
13 
14 import dlsl.vector;
15 
16 import std.conv	: to;
17 import std.math	: sqrt, sin, cos, tan, PI;
18 import std.array : join;
19 import std.traits : isNumeric, isFloatingPoint, isArray;
20 import std..string : format, rightJustify;
21 import std.algorithm : max, min, reduce;
22 
23 
24 /// predefined matrix types, the first number represents the number of cols
25 /// and the second the number of columns, if there's just one it's a nxn matrix.
26 /// All of these matrices are floating - point matrices.
27 alias Matrix!( float, 2, 2 ) mat2;
28 alias Matrix!( float, 2, 3 ) mat23;
29 alias Matrix!( float, 2, 4 ) mat24;
30 
31 alias Matrix!( float, 3, 2 ) mat32;
32 alias Matrix!( float, 3, 3 ) mat3;
33 alias Matrix!( float, 3, 4 ) mat34;
34 
35 alias Matrix!( float, 4, 2 ) mat42;
36 alias Matrix!( float, 4, 3 ) mat43;
37 alias Matrix!( float, 4, 4 ) mat4;
38 
39 
40 enum float	deg2rad	= 0.0174532925199432957692369076849f;
41 enum float	rad2deg	= 57.295779513082320876798154814105f;
42 
43 version( NoReciprocalMul ) {
44 	private enum rmul = false;
45 }	else {
46 	private enum rmul = true;
47 }
48 
49 /// If T is a matrix, this evaluates to true, otherwise false
50 template isMatrix( T ) {  enum isMatrix = is( typeof( isMatrixImpl( T.init )));  }
51 private void isMatrixImpl( T, int cols, int rows )( Matrix!( T, cols, rows ) mat ) {}
52 
53 
54 /// Base template for all matrix types. See predefined matrix types
55 /// Params:
56 /// type = the value type of each matrix element
57 /// cols = count of columns of the matrix
58 /// rows = count of rows of the matrix
59 /// - - -
60 struct Matrix( type, int numCols, int numRows ) if(( numCols > 1 ) && ( numRows > 1 )) {
61 	alias type valueType;						/// Holds the internal primitive type of the matrix;
62 	alias Vector!( type, rows ) vectorType;		/// Holds the internal vectorType of the matrix;
63 
64 	/// Holds the matrix $( RED column - major ) in memory.
65 	vectorType[ cols ] data;					/// Each Column is Vector of length rows
66 	alias data this;							/// Treat this matrix as an nxm array
67 	alias cols = numCols;						/// Holds the number of cols;
68 	alias rows = numRows;						/// Holds the number of columns;
69 
70 
71 	unittest {		// Construction through aliased this static array of Vectors
72 
73 		mat2 m2 = mat2( 0.0f, 1.0f, 2.0f, 3.0f );
74 		m2 = mat2( [ 0.0f, 1.0f ], [ 2.0f, 3.0f ] );
75 		m2 = mat2( vec2( 0.0f, 1.0f ), vec2( 2.0f, 3.0f ));
76 		assert( m2[ 0 ][ 0 ] == 0.0f );
77 		assert( m2[ 0 ][ 1 ] == 1.0f );
78 		assert( m2[ 1 ][ 0 ] == 2.0f );
79 		assert( m2[ 1 ][ 1 ] == 3.0f );
80 
81 		m2[ 0 ] = [ 2.0, 2.0 ];
82 		m2[ 0 ] = vec2( 2.0, 2.0 );
83 		//m2[ 0 .. 1 ] = [ [ 2.0f, 2.0f ] ];
84 		m2[ 0 .. 1 ] = [ vec2( 2.0, 2.0 ) ];
85 
86 		assert( m2 == [ [ 2.0f, 2.0f ], [ 2.0f, 3.0f ] ] );
87 
88 		mat3 m3 = mat3( 0.0f, 0.1f, 0.2f, 1.0f, 1.1f, 1.2f, 2.0f, 2.1f, 2.2f );
89 		assert( m3[ 0 ][ 1 ] == 0.1f );
90 		assert( m3[ 2 ][ 0 ] == 2.0f );
91 		assert( m3[ 1 ][ 2 ] == 1.2f );
92 
93 		m3[ 0 ][ 0 .. $ ] = 0.0f;
94 		assert(  m3 == [[ 0.0f, 0.0f, 0.0f ],
95 						[ 1.0f, 1.1f, 1.2f ],
96 						[ 2.0f, 2.1f, 2.2f ] ] );
97 		//m3[ 1 .. 3 ]  = [ [ 1, 1, 1 ], [ 2, 2, 2 ] ];
98 		m3[ 1 .. 3 ] =  [ vec3( 1, 1, 1 ), vec3( 2, 2, 2 ) ];
99 		assert(  m3 == [[ 0.0f, 0.0f, 0.0f ],
100 						[ 1.0f, 1.0f, 1.0f ],
101 						[ 2.0f, 2.0f, 2.0f ]] );
102 
103 		mat4 m4 = mat4( 0.0f, 0.1f, 0.2f, 0.3f,
104 							 1.0f, 1.1f, 1.2f, 1.3f,
105 							 2.0f, 2.1f, 2.2f, 2.3f,
106 							 3.0f, 3.1f, 3.2f, 3.3f );
107 		assert( m4[ 0 ][ 3 ] == 0.3f );
108 		assert( m4[ 1 ][ 1 ] == 1.1f );
109 		assert( m4[ 2 ][ 0 ] == 2.0f );
110 		assert( m4[ 3 ][ 2 ] == 3.2f );
111 
112 		m4[ 2 ][ 1 .. 3 ] = [ 1.0f, 2.0f ];
113 		assert(  m4 == [[ 0.0f, 0.1f, 0.2f, 0.3f ],
114 						[ 1.0f, 1.1f, 1.2f, 1.3f ],
115 						[ 2.0f, 1.0f, 2.0f, 2.3f ],
116 						[ 3.0f, 3.1f, 3.2f, 3.3f ]] );
117 
118 	}
119 
120 
121 	/// Returns the pointer to the stored values as OpenGL requires it.
122 	/// Note this will return a pointer to a $( RED column - major ) matrix,
123 	/// $( RED this is the OpneGL convention and expected in programs via Uniforms or UBOs ).
124 	@property auto ptr() { return data[ 0 ].ptr; }
125 
126 	/// Returns the current matrix formatted as flat string.
127 	@property string asString() {  return format( "%s", data );  }
128 	alias asString toString;
129 
130 	/// Returns the current matrix as pretty formatted string.
131 	/// TODO : Check This
132 	@property string asPrettyString() {
133 		string fmtr = "%s";
134 
135 		size_t rjust = max( format( fmtr, reduce!( max )( data[] )).length,
136 								  format( fmtr, reduce!( min )( data[] )).length ) - 1;
137 
138 		string[] outer_parts;
139 		foreach( valueType[] col; data ) {
140 			string[] inner_parts;
141 			foreach( valueType row; col ) {
142 				inner_parts ~= rightJustify( format( fmtr, row ), rjust );
143 			}
144 			outer_parts ~= " [ " ~ join( inner_parts, ", " ) ~ " ]";
145 		}
146 
147 		return "[ " ~ join( outer_parts, "\n" )[ 1 .. $ ] ~ " ]";
148 	} alias asPrettyString toPrettyString;
149 
150 	static @property Matrix identity() if( cols == rows ) {
151 		Matrix result;
152 		foreach( ref r; rows )
153 			foreach( ref c; cols )
154 				result[r][c] = r == c ? 1 : 0;
155 		return result;
156 	}
157 
158 	@safe pure nothrow :
159 	template isCompatibleMatrix( T ) {  enum isCompatibleMatrix = is( typeof( isCompatibleMatrixImpl( T.init )));  }
160 	static void isCompatibleMatrixImpl( int col, int row )( Matrix!( valueType, col, row ) mat ) {}
161 
162 	template isCompatibleVector( T ) {  enum isCompatibleVector = is( typeof( isCompatibleVectorImpl( T.init )));  }
163 	static void isCompatibleVectorImpl( int dim )( Vector!( valueType, dim ) vec ) {}
164 
165 
166 	/// TODO: Fix Construction with combination of numeric, static and dynamic array, vector and matrix
167 	private void construct( int i, T, Tail... )( T head, Tail tail ) {
168 		static if( i >= cols * rows ) {
169 			static assert( false, "constructor has too many arguments" );
170 		}	else static if( is( T : valueType )) {
171 			data[ i / rows ][ i % rows ] = head;
172 			construct!( i + 1 )( tail );
173 		}	else static if( is( T == Vector!( valueType, rows ))) {
174 			static if( i % rows == 0 ) {
175 				data[ i / rows ] = head;
176 				construct!( i + T.dimension )( tail );
177 			}	else {
178 				static assert( false, "Can't convert Vector into the matrix. Maybe it doesn't align to the columns correctly or dimension doesn't fit" );
179 			}
180 		}	/*else {
181 			static assert( false, "Matrix constructor argument must be of type " ~ valueType.stringof ~ " or Vector, not " ~ T.stringof );
182 		}*/
183 	}
184 
185 	private void construct( int i )() {}	// terminate
186 
187 
188 	/// Constructs the matrix:
189 	/// If a single value is passed, the matrix will be cleared with this value ( each column in each col will contain this value ).
190 	/// If a matrix with more cols and columns is passed, the matrix will be the upper left nxm matrix.
191 	/// If a matrix with less cols and columns is passed, the passed matrix will be stored in the upper left of an identity matrix.
192 	/// It's also allowed to pass vectors and scalars at a time, but the vectors dimension must match the number of columns and align correctly.
193 	/// Examples:
194 	/// - - -
195 	/// mat2 m2 = mat2( 0.0f ); // mat2 m2 = mat2( 0.0f, 0.0f, 0.0f, 0.0f );
196 	/// mat3 m3 = mat3( m2 ); // mat3 m3 = mat3( 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f );
197 	/// mat3 m3_2 = mat3( vec3( 1.0f, 2.0f, 3.0f ), 4.0f, 5.0f, 6.0f, vec3( 7.0f, 8.0f, 9.0f ));
198 	/// mat4 m4 = mat4.identity; // just an identity matrix
199 	/// mat3 m3_3 = mat3( m4 ); // mat3 m3_3 = mat3.identity
200 	/// - - -
201 	this( Args... )( Args args ) {  construct!( 0 )( args );  }
202 
203 	/// Construct a Matrix from another Matrix, equal sized or bigger
204 //	this( T )( T mat ) if( isMatrix!T && ( T.rows >= rows ) && ( T.cols >= cols )) {
205 //		for( int c = 0; c < cols; ++c ) {
206 //			for( int r = 0; r < rows; ++r ) {
207 //				data[ c ][ r ] = mat.data[ c ][ r ];
208 //			}
209 //		}
210 //	}
211 
212 
213 	/// Construct a Matrix from another Matrix
214 	this( T )( T mat ) if( isMatrix!T ) {
215 		this( 1 );	// Could be optimized with not initializing and setting missing values within static if range checks
216 		//alias Cols = cols < T.cols ? cols : T.cols;
217 		//alias Rows = rows < T.rows ? rows : T.rows;
218 		foreach( c; 0 .. cols < T.cols ? cols : T.cols ) {
219 			foreach( r; 0 .. rows < T.rows ? rows : T.rows ) {
220 				data[ c ][ r ] = mat.data[ c ][ r ];
221 			}
222 		}
223 	}
224 
225 	/// Construct a Matrix from a single value, for non- and square matrices
226 	/// see GLSL 4.5 Spec, section 5.4.2 Vector and Matrix Constructors
227 	this()( valueType value ) {
228 		clear( 0 );
229 		foreach( rc; 0 .. cols < rows ? cols : rows ) {
230 			data[ rc ][ rc ] = value;
231 		}
232 	}
233 
234 	/// Sets all values of the matrix to value ( each column in each col will contain this value ).
235 	void clear( valueType value ) {
236 		foreach( ref vec; data )
237 			vec[] = value;
238 	}
239 
240 	unittest {		/// Matrix.isvalid, Constructing
241 
242 		mat2 m2 = mat2( 1.0f, 1.0f, vec2( 2.0f, 2.0f ));
243 		assert( m2.data == [[ 1.0f, 1.0f ], [ 2.0f, 2.0f ]] );
244 		m2.clear( 3.0f );
245 		assert( m2.data == [[ 3.0f, 3.0f ], [ 3.0f, 3.0f ]] );
246 		assert( m2.isvalid );
247 		m2.clear( float.nan );
248 		assert( !m2.isvalid );
249 		m2.clear( float.infinity );
250 		assert( !m2.isvalid );
251 		m2.clear( 0.0f );
252 		assert( m2.isvalid );
253 
254 		mat3 m3 = mat3( 1.0f );
255 		assert( m3.data == [[ 1.0f, 0.0f, 0.0f ],
256 							[ 0.0f, 1.0f, 0.0f ],
257 							[ 0.0f, 0.0f, 1.0f ]] );
258 
259 		mat4 m4 = mat4( vec4( 	1.0f, 1.0f, 1.0f, 1.0f ),
260 								2.0f, 2.0f, 2.0f, 2.0f,
261 								3.0f, 3.0f, 3.0f, 3.0f,
262 							 vec4( 4.0f, 4.0f, 4.0f, 4.0f ));
263 		assert( m4.data == [[ 1.0f, 1.0f, 1.0f, 1.0f ],
264 							[ 2.0f, 2.0f, 2.0f, 2.0f ],
265 							[ 3.0f, 3.0f, 3.0f, 3.0f ],
266 							[ 4.0f, 4.0f, 4.0f, 4.0f ]] );
267 		assert( mat3( m4 ).data == [[ 1.0f, 1.0f, 1.0f ],
268 									[ 2.0f, 2.0f, 2.0f ],
269 									[ 3.0f, 3.0f, 3.0f ]] );
270 		assert( mat2( mat3( m4 )).data == [[ 1.0f, 1.0f ],
271 														[ 2.0f, 2.0f ] ] );
272 		assert( mat2( m4 ).data == mat2( mat3( m4 )).data );
273 		assert( mat4( mat3( m4 )).data ==  [[ 1.0f, 1.0f, 1.0f, 0.0f ],
274 											[ 2.0f, 2.0f, 2.0f, 0.0f ],
275 											[ 3.0f, 3.0f, 3.0f, 0.0f ],
276 											[ 0.0f, 0.0f, 0.0f, 1.0f ]] );
277 
278 		Matrix!( float, 2, 3 ) mt1 = Matrix!( float, 2, 3 )( 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f );
279 		Matrix!( float, 3, 2 ) mt2 = Matrix!( float, 3, 2 )( 6.0f, - 1.0f, 3.0f, 2.0f, 0.0f, - 3.0f );
280 
281 		assert( mt1.data == [[ 1.0f, 2.0f, 3.0f ], [ 4.0f, 5.0f, 6.0f ]] );
282 		assert( mt2.data == [[ 6.0f, - 1.0f ], [ 3.0f, 2.0f ], [ 0.0f, - 3.0f ]] );
283 	}
284 
285 
286 	//void opAssign( T )( T array )  if( is( T array == U[][], U )) { //if( isArray( A )) {
287 	//	pragma( msg, "Instantiat opAssign" );
288 	//	data[] = array[];
289 	//}
290 
291 	static if( cols == rows ) {
292 		/// static create identity matrix
293 		static @property Matrix identity() {
294 			return Matrix( 1 );
295 		}
296 
297 		/// Returns a identity matrix.
298 		@property bool isIdentity() {
299 			return data == Matrix( 1 ).data;
300 		}
301 
302 		/// Unittest identity, isIdentity
303 		unittest {
304 			mat2 m2 = mat2( 1.0f );
305 			assert( m2.isIdentity );
306 			assert( m2.data == mat2( 1 ).data );
307 
308 			m2 = mat2.identity;
309 			assert( m2.data == [[ 1.0f, 0.0f ],
310 								[ 0.0f, 1.0f ]] );
311 
312 			mat3 m3 = mat3.identity;
313 			assert( m3.isIdentity );
314 			assert( m3.data == mat3( 1.0f ).data );
315 
316 			mat4 m4 = mat4( 2.0f );
317 			assert( m4.data == mat4( 2.0f ).data );
318 
319 			m4 = mat4.identity;
320 			assert( m4.isIdentity );
321 			assert( m4.data == [[ 1.0f, 0.0f, 0.0f, 0.0f ],
322 								[ 0.0f, 1.0f, 0.0f, 0.0f ],
323 								[ 0.0f, 0.0f, 1.0f, 0.0f ],
324 								[ 0.0f, 0.0f, 0.0f, 1.0f ]] );
325 			assert( m4.data == m4.identity.data );
326 		}
327 
328 
329 
330 
331 		/////////////////
332 		// Translation //
333 		/////////////////
334 		static if( cols == 2 ) {			// homogeneous 1D translations
335 			/// static construction of a translation matrix
336 			static Matrix translation( valueType x ) {
337 				Matrix result = Matrix.identity;
338 				result.data[ 1 ][ 0 ] = x;
339 				return result;
340 			}
341 
342 			/// translate an existing matrix
343 			Matrix translate( valueType x ) {
344 				this = Matrix.translation( x ) * this;
345 				return this;
346 			}
347 		}
348 
349 		else static if( cols == 3 ) {		// homogeneous 2D translations
350 			/// static construction of a translation matrix from an array/vector
351 			static Matrix translation( valueType[ 2 ] vec )  {
352 				Matrix result = Matrix.identity;
353 				result.data[ 2 ][ 0..2 ] = vec;
354 				return result;
355 			}
356 
357 			/// static construction of a translation matrix from two coordinates
358 			static Matrix translation( valueType x, valueType y )  {
359 				return translation( [ x, y ] );
360 			}
361 
362 			/// translate an existing matrix with an array/vector
363 			Matrix translate( valueType[ 2 ] vec ) {
364 				this = Matrix.translation( vec ) * this;
365 				return this;
366 			}
367 
368 			/// translate an existing matrix with two scalars
369 			Matrix translate( valueType x, valueType y ) {
370 				this = Matrix.translation( [ x, y ] ) * this;
371 				return this;
372 			}
373 		}
374 
375 		else static if( cols == 4 ) {		// homogeneous 3D translations
376 			/// static construction of a translation matrix from an array/vector
377 			static Matrix translation( valueType[ 3 ] vec )  {
378 				Matrix result = Matrix.identity;
379 				result.data[ 3 ][ 0..3 ] = vec;
380 				return result;
381 			}
382 
383 			/// static construction of a translation matrix from three coordinates
384 			static Matrix translation( valueType x, valueType y, valueType z )  {
385 				return translation( [ x, y, z ] );
386 			}
387 
388 			/// translate an existing matrix with a vector
389 			Matrix translate( valueType[ 3 ] vec ) {
390 				this = Matrix.translation( vec ) * this;
391 				return this;
392 			}
393 
394 			/// translate an existing matrix with three scalars
395 			Matrix translate( valueType x, valueType y, valueType z ) {
396 				this = Matrix.translation( [ x, y, z ] ) * this;
397 				return this;
398 			}
399 
400 			unittest {
401 				mat4 m4 = mat4( 1.0f );
402 				assert( m4.translation( 1.0f, 2.0f, 3.0f ).data == mat4.translation( 1.0f, 2.0f, 3.0f ).data );
403 				assert( mat4.translation( 1.0f, 2.0f, 3.0f ).data  ==  [[ 1.0f, 0.0f, 0.0f, 0.0f ],
404 																		[ 0.0f, 1.0f, 0.0f, 0.0f ],
405 																		[ 0.0f, 0.0f, 1.0f, 0.0f ],
406 																		[ 1.0f, 2.0f, 3.0f, 1.0f ] ] );
407 				assert( mat4.identity.translate( 0.0f, 1.0f, 2.0f ).data == mat4.translation( 0.0f, 1.0f, 2.0f ).data );
408 			}
409 		}
410 
411 		/// extract translation part of the matrix in a copy
412 		Matrix translation() {
413 			Matrix result = Matrix.identity;
414 			result[ $-1 ][ 0 .. $-1 ] = data[ $-1 ][ 0 .. $-1 ];
415 			return result;
416 		}
417 
418 		unittest  {
419 			mat3 m3 = mat3( 0.0f, 1.0f, 2.0f,
420 							3.0f, 4.0f, 5.0f,
421 							6.0f, 7.0f, 1.0f );
422 			assert( m3.translation.data == [ [ 1.0f, 0.0f, 2.0f ], [ 0.0f, 1.0f, 5.0f ], [ 0.0f, 0.0f, 1.0f ] ] );
423 
424 
425 			mat4 m4 = mat4(  0.0f,  1.0f,  2.0f,  3.0f,
426 							 4.0f,  5.0f,  6.0f,  7.0f,
427 							 8.0f,  9.0f, 10.0f, 11.0f,
428 							12.0f, 13.0f, 14.0f,  1.0f );
429 			assert( m4.translation.data == [ [ 1.0f, 0.0f, 0.0f,  3.0f ],
430 											 [ 0.0f, 1.0f, 0.0f,  7.0f ],
431 											 [ 0.0f, 0.0f, 1.0f, 11.0f ],
432 											 [ 0.0f, 0.0f, 0.0f,  1.0f ] ] );
433 		}
434 
435 
436 
437 		//////////////
438 		// Rotation // - integer rotations might not be useful, but should they be disabled?
439 		//////////////
440 		static if( cols == 2 || cols == 3 ) {		// non- and homogeneous 2D rotations
441 			/// static construction of a rotation matrix
442 			static Matrix rotation( real angle ) {
443 				Matrix result = Matrix.identity;
444 				valueType cosAngle = to!valueType( cos( angle ));
445 				valueType sinAngle = to!valueType( sin( angle ));
446 				result.data[ 0 ][ 0 ] = cosAngle;
447 				result.data[ 1 ][ 0 ] = - sinAngle;
448 				result.data[ 0 ][ 1 ] = sinAngle;
449 				result.data[ 1 ][ 1 ] = cosAngle;
450 				return result;
451 			}
452 
453 			/// rotate an existing matrix
454 			Matrix rotate( real angle ) {
455 				this = rotation( angle ) * this;
456 				return this;
457 			}
458 		}
459 
460 		static if( cols >= 3 ) {						// non- and homogeneous 3D rotations
461 			/// static construction of a rotation matrix angle and axis
462 			static Matrix rotation( real angle, Vector!( valueType, 3 ) axis ) {
463 				Matrix result = Matrix.identity;
464 
465 				auto length = axis.length;
466 				if( length != 1 ) {
467 					version( NoReciprocalMul ) {
468 						axis /= length;
469 					}	else	{
470 						auto invLength = 1.0 / length;
471 						axis *= invLength;
472 					}
473 				}
474 
475 				const real cosAngle = cos( angle );
476 				const real sinAngle = sin( angle );
477 
478 				Vector!( valueType, 3 ) temp = ( 1 - cosAngle ) * axis;
479 
480 				result.data[ 0 ][ 0 ] = to!valueType( cosAngle	+	temp.x * axis.x );
481 				result.data[ 0 ][ 1 ] = to!valueType(				temp.x * axis.y + sinAngle * axis.z );
482 				result.data[ 0 ][ 2 ] = to!valueType(				temp.x * axis.z - sinAngle * axis.y );
483 				result.data[ 1 ][ 0 ] = to!valueType(				temp.y * axis.x - sinAngle * axis.z );
484 				result.data[ 1 ][ 1 ] = to!valueType( cosAngle	+	temp.y * axis.y );
485 				result.data[ 1 ][ 2 ] = to!valueType(				temp.y * axis.z + sinAngle * axis.x );
486 				result.data[ 2 ][ 0 ] = to!valueType(				temp.z * axis.x + sinAngle * axis.y );
487 				result.data[ 2 ][ 1 ] = to!valueType(				temp.z * axis.y - sinAngle * axis.x );
488 				result.data[ 2 ][ 2 ] = to!valueType( cosAngle	+	temp.z * axis.z );
489 
490 				return result;
491 			}
492 
493 			/// static construction of a rotation matrix angle and axis coordinates
494 			static Matrix rotation( real angle, valueType x, valueType y, valueType z ) {
495 				return Matrix.rotation( angle, Vector!( valueType, 3 )( x, y, z ));
496 			}
497 
498 			/// rotate an existing matrix with an angle around an axis
499 			Matrix rotate( real angle, Vector!( valueType, 3 ) axis ) {
500 				this = rotation( angle, axis ) * this;
501 				return this;
502 			}
503 
504 			/// rotate an existing matrix with an angle around axis coordinates
505 			Matrix rotate( real angle, valueType x, valueType y, valueType z  ) {
506 				this = rotation( angle, x, y, z ) * this;
507 				return this;
508 			}
509 
510 			/// static construction of a rotation matrix with an angle around a canonical axis
511 			private static Matrix rotationA( ubyte a, ubyte b )( real angle ) {
512 				Matrix result = Matrix.identity;
513 
514 				valueType cosAngle = to!valueType( cos( angle ));
515 				valueType sinAngle = to!valueType( sin( angle ));
516 
517 				result.data[ a ][ a ] = cosAngle;
518 				result.data[ b ][ a ] = - sinAngle;
519 				result.data[ a ][ b ] = sinAngle;
520 				result.data[ b ][ b ] = cosAngle;
521 
522 				return result;
523 			}
524 
525 			/// static construction of a rotation matrix with an angle around X axis
526 			static Matrix rotationX( real angle ) { return rotationA!( 1, 2 )( angle );	} /// A-canonical = X
527 
528 			/// static construction of a rotation matrix with an angle around Y axis
529 			static Matrix rotationY( real angle ) { return rotationA!( 2, 0 )( angle );	} /// A-canonical = Y
530 
531 			/// static construction of a rotation matrix with an angle around Z axis
532 			static Matrix rotationZ( real angle ) { return rotationA!( 0, 1 )( angle );	} /// A-canonical = Z
533 
534 			/// rotate an existing matrix with an angle around X axis
535 			Matrix rotateX( real angle ) {
536 				this = rotationX( angle ) * this;
537 				return this;
538 			}
539 
540 			/// rotate an existing matrix with an angle around Y axis
541 			Matrix rotateY( real angle ) {
542 				this = rotationY( angle ) * this;
543 				return this;
544 			}
545 
546 			/// rotate an existing matrix with an angle around Z axis
547 			Matrix rotateZ( real angle ) {
548 				this = rotationZ( angle ) * this;
549 				return this;
550 			}
551 
552 			unittest {
553 				assert( mat4.rotationX( 0 ).data == [[ 1.0f,   0.0f,   0.0f, 0.0f ],
554 													 [ 0.0f,   1.0f, - 0.0f, 0.0f ],
555 													 [ 0.0f,   0.0f,   1.0f, 0.0f ],
556 													 [ 0.0f,   0.0f,   0.0f, 1.0f ]] );
557 				assert( mat4.rotationY( 0 ).data == [[ 1.0f,   0.0f,   0.0f, 0.0f ],
558 													 [ 0.0f,   1.0f,   0.0f, 0.0f ],
559 													 [ 0.0f,   0.0f,   1.0f, 0.0f ],
560 													 [ 0.0f,   0.0f,   0.0f, 1.0f ]] );
561 				assert( mat4.rotationZ( 0 ).data == [[ 1.0f, - 0.0f,   0.0f, 0.0f ],
562 													 [ 0.0f,   1.0f,   0.0f, 0.0f ],
563 													 [ 0.0f,   0.0f,   1.0f, 0.0f ],
564 													 [ 0.0f,   0.0f,   0.0f, 1.0f ]] );
565 
566 				//mat4 rotX = mat4.identity;
567 				//rotX.rotateX( 1 );
568 				//assert( mat4.rotationX( 1 ).data == rotX.data );
569 				//assert( rotX.data == mat4.identity.rotateX( 1 ).data );
570 				//assert( rotX.data == mat4.rotation( 1, vec3( 1.0f, 0.0f, 0.0f )).data );
571 
572 				//mat4 rotY = mat4.identity;
573 				//rotY.rotateY( 2 );
574 				//assert( mat4.rotationY( 2 ).data == rotY.data );
575 				//assert( rotY.data == mat4.identity.rotateY( 2 ).data );
576 				//assert( rotY.data == mat4.rotation( 2, vec3( 0.0f, 1.0f, 0.0f )).data );
577 
578 				//mat4 rotZ = mat4.identity;
579 				//rotZ.rotateZ( 3 );
580 				//assert( mat4.rotationZ( 3 ).data == rotZ.data );
581 				//assert( rotZ.data == mat4.identity.rotateZ( 3 ).data );
582 				//assert( rotZ.data == mat4.rotation( 3, vec3( 0.0f, 0.0f, 1.0f )).data );
583 			}
584 		}
585 
586 		/// extract rotation part of the matrix in a copy
587 		// TODO(pp): fix proper extraction, currently works only right if determinant = 1
588 		Matrix rotation() {
589 			Matrix result = data;
590 			result[ $-1 ][ 0..$-1 ] = 0;
591 			result[ 0..$-1 ][ $-1 ] = 0;
592 			result[ $-1 ][ $-1 ] = 1;
593 			return result;
594 		}
595 
596 		unittest {
597 			mat3 m3 = mat3( 0.0f, 1.0f, 2.0f,
598 							3.0f, 4.0f, 5.0f,
599 							6.0f, 7.0f, 1.0f );
600 			assert( m3.rotation.data == [[ 0.0f, 1.0f, 0.0f ],
601 										 [ 3.0f, 4.0f, 0.0f ],
602 										 [ 0.0f, 0.0f, 1.0f ]] );
603 
604 			mat4 m4 = mat4(  0.0f,  1.0f,  2.0f,  3.0f,
605 							 4.0f,  5.0f,  6.0f,  7.0f,
606 							 8.0f,  9.0f, 10.0f, 11.0f,
607 							12.0f, 13.0f, 14.0f,  1.0f );
608 
609 			assert( m4.rotation.data == [[  0.0f,  1.0f,  2.0f, 0.0f ],
610 										 [  4.0f,  5.0f,  6.0f, 0.0f ],
611 										 [  8.0f,  9.0f, 10.0f, 0.0f ],
612 										 [  0.0f,  0.0f,  0.0f, 1.0f ]] );
613 		}
614 
615 
616 
617 		/////////////
618 		// Scaling //
619 		/////////////
620 		static if( cols == 2 ) {						// homogeneous 1D scaling
621 			/// static construction of a scaling matrix
622 			static Matrix scaling( valueType x ) {
623 				Matrix result = Matrix.identity;
624 				result.data[ 0 ][ 0 ] = x;
625 				return result;
626 			}
627 
628 			/// scale an existing matrix
629 			Matrix scale( valueType x ) {
630 				this = Matrix.scaling( x ) * this;
631 				return this;
632 			}
633 		}
634 
635 		static if( cols == 2 || cols == 3 ) {		// non- and homogeneous 2D scaling
636 			/// static construction of a scaling matrix from an array/vector
637 			static Matrix scaling( valueType[ 2 ] vec ) {
638 				Matrix result = Matrix.identity;
639 				result.data[ 0 ][ 0 ] = vec[ 0 ];
640 				result.data[ 1 ][ 1 ] = vec[ 1 ];
641 				return result;
642 			}
643 
644 			/// static construction of a scaling matrix from two scalars
645 			static Matrix scaling( valueType x, valueType y ) {
646 				return scaling( [ x, y ] );
647 			}
648 
649 			/// scale an existing matrix with an array/vector
650 			Matrix scale( valueType[ 2 ] vec ) {
651 				this = Matrix.scaling( vec ) * this;
652 				return this;
653 			}
654 
655 			/// scale an existing matrix with two scalars
656 			Matrix scale( valueType x, valueType y ) {
657 				this = Matrix.scaling( [ x, y ] ) * this;
658 				return this;
659 			}
660 		}
661 
662 		else static if( cols >= 3 ) {
663 			/// static construction of a scaling matrix from an array/vector
664 			static Matrix scaling( valueType[ 3 ] vec ) {
665 				Matrix result = Matrix.identity;
666 				result.data[ 0 ][ 0 ] = vec[ 0 ];
667 				result.data[ 1 ][ 1 ] = vec[ 1 ];
668 				result.data[ 2 ][ 2 ] = vec[ 2 ];
669 				return result;
670 			}
671 
672 			/// static construction of a scaling matrix from three scalars
673 			static Matrix scaling( valueType x, valueType y, valueType z ) {
674 				return scaling( [ x, y, z ] );
675 			}
676 
677 			/// scale an existing matrix with an array/vector
678 			Matrix scale( valueType[ 3 ] vec ) {
679 				this = Matrix.scaling( vec ) * this;
680 				return this;
681 			}
682 
683 			/// scale an existing matrix with three scalars
684 			Matrix scale( valueType x, valueType y, valueType z ) {
685 				this = Matrix.scaling( [ x, y, z ] ) * this;
686 				return this;
687 			}
688 
689 			unittest {
690 				mat4 m4 = mat4( 1.0f );
691 				assert( m4.scaling( 0.0f, 1.0f, 2.0f ).data == mat4.scaling( 0.0f, 1.0f, 2.0f ).data );
692 				assert( mat4.scaling( 0.0f, 1.0f, 2.0f ).data  ==  [[ 0.0f, 0.0f, 0.0f, 0.0f ],
693 																	[ 0.0f, 1.0f, 0.0f, 0.0f ],
694 																	[ 0.0f, 0.0f, 2.0f, 0.0f ],
695 																	[ 0.0f, 0.0f, 0.0f, 1.0f ] ] );
696 				assert( mat4.identity.scale( 0.0f, 1.0f, 2.0f ).data == mat4.scaling( 0.0f, 1.0f, 2.0f ).data );
697 			}
698 		}
699 
700 		/// extract scaling part of the matrix in a copy
701 		Matrix scale() {
702 			Matrix result = Matrix.identity;
703 			foreach( rc; 0 .. cols - 1 )
704 				result[ rc ][ rc ] = data[ rc ][ rc ];
705 			return result;
706 		}
707 
708 		unittest {
709 			mat3 m3 = mat3( 0.0f, 1.0f, 2.0f,
710 							3.0f, 4.0f, 5.0f,
711 							6.0f, 7.0f, 1.0f );
712 
713 			assert( m3.scale.data  ==  [[ 0.0f, 0.0f, 0.0f ],
714 										[ 0.0f, 4.0f, 0.0f ],
715 										[ 0.0f, 0.0f, 1.0f ] ] );
716 
717 
718 			mat4 m4 = mat4(  0.0f,  1.0f,  2.0f,  3.0f,
719 							 4.0f,  5.0f,  6.0f,  7.0f,
720 							 8.0f,  9.0f, 10.0f, 11.0f,
721 							12.0f, 13.0f, 14.0f,  1.0f );
722 
723 			assert( m4.scale.data  ==  [[ 0.0f, 0.0f,  0.0f, 0.0f ],
724 										[ 0.0f, 5.0f,  0.0f, 0.0f ],
725 										[ 0.0f, 0.0f, 10.0f, 0.0f ],
726 										[ 0.0f, 0.0f,  0.0f, 1.0f ]] );
727 		}
728 	}	// end static if( cols == rows )
729 
730 
731 	/// Component-wise binary matrix-scalar operation: addition, subtraction, multiplication, division
732 	auto opBinary( string op, T )( T s ) const if( isNumeric!T && (( op == "+" ) || ( op == "-" ) || ( op == "*" ) || ( op == "/" ))) {
733 		Matrix result;
734 		mixin( "result[0] = data[0] " ~ op ~ " s;" );
735 		mixin( "result[1] = data[1] " ~ op ~ " s;" );
736 		static if( cols >= 3 )  mixin( "result[2] = data[2] " ~ op ~ " s;" );
737 		static if( cols == 4 )  mixin( "result[3] = data[3] " ~ op ~ " s;" );
738 		return result;
739 	}
740 
741 	/// Component-wise binary scalar-matrix operation: addition, subtraction, multiplication, division
742 	auto opBinaryRight( string op, T )( T s ) const if( isNumeric!T && (( op == "+" ) || ( op == "-" ) || ( op == "*" ) || ( op == "/" ))) {
743 		Matrix result;
744 		mixin( "result[0] = s " ~ op ~ "data[0];" );
745 		mixin( "result[1] = s " ~ op ~ "data[1];" );
746 		static if( cols >= 3 )  mixin( "result[2] = s " ~ op ~ "data[2];" );
747 		static if( cols == 4 )  mixin( "result[3] = s " ~ op ~ "data[3];" );
748 		return result;
749 	}
750 
751 	unittest {		/// Matrix-Scalar Operations
752 		mat2 m2 = mat2( 1.0f, 2.0f, 3.0f, 4.0f );
753 		assert(( m2 * 2.0f ).data == [ [ 2.0f, 4.0f ], [ 6.0f, 8.0f ] ] );
754 		assert(( 2 * m2 ).data == ( m2 * 2 ).data );
755 
756 		mat3 m3 = mat3( 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f );
757 		assert(( m3 * 2 ).data == [ [ 2.0f, 4.0f, 6.0f ], [ 8.0f, 10.0f, 12.0f ], [ 14.0f, 16.0f, 18.0f ] ] );
758 		assert(( 2 * m3 ).data == ( m3 * 2 ).data );
759 
760 		//TODO: tests for mat4, mat34
761 	}
762 
763 
764 	/// matrix-vector multiplication
765 	auto opBinary( string op : "*", V )( V vec ) const if( isVector!V && V.dimension == cols ) {
766 		vectorType result;
767 		result.clear( 0 );
768 		foreach( r; 0 .. rows )
769 			foreach( c; 0 .. cols )
770 				result[ r ] += data[ c ][ r ] * vec[ c ];
771 		return result;
772 	}
773 
774 	/// vector-Matrix multiplication, optimized instead of transposing matrix and multiplying
775 	auto opBinaryRight( string op : "*", V )( V vec ) const if( isVector!V && V.dimension == rows ) {
776 		Vector!( valueType, cols ) result;
777 		result.clear( 0 );
778 		foreach( c; 0 .. cols )
779 			foreach( r; 0 .. rows )
780 				result[ c ] += data[ c ][ r ] * vec[ r ];
781 		return result;
782 	}
783 
784 	unittest {		/// matrix-vector, vector-matrix multiplication
785 		mat2 m2 = mat2( 2.0f, 4.0f, 6.0f, 8.0f );
786 		vec2 v2 = vec2( 2.0f, 2.0f );
787 		assert(( m2 * v2 ) == [ 16.0f, 24.0f ] );
788 		assert(( v2 * m2 ) == [ 12.0f, 28.0f ] );
789 
790 		mat3 m3 = mat3( 2.0f, 4.0f, 6.0f, 8.0f, 10.0f, 12.0f, 14.0f, 16.0f, 18.0f );
791 		vec3 v3 = vec3( 2.0f, 2.0f, 2.0f );
792 		assert(( m3 * v3 ) == [ 48.0f, 60.0f, 72.0f ] );
793 		assert(( v3 * m3 ) == [ 24.0f, 60.0f, 96.0f ] );
794 
795 		mat4 m4 = mat4( 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32 );
796 		vec4 v4 = vec4( 2, 2, 2, 2 );
797 		assert(( m4 * v4 ) == [ 112, 128, 144, 160 ] );
798 		assert(( v4 * m4 ) == [  40, 104, 168, 232 ] );
799 
800 		//TODO: tests for mat4, mat34, mat43, mat24, mat 42
801 	}
802 
803 
804 	/// matrix-matrix component-wise operations addition, subtraction and division, using vector-vector operations of all colums
805 	Matrix opBinary( string op )( Matrix mat ) const if( op == "+" || op == "-" || op == "/" ) {
806 		Matrix result;
807 		mixin( "result[0] = data[0] " ~ op ~ "ma t[0];" );
808 		mixin( "result[1] = data[1] " ~ op ~ "ma t[1];" );
809 		static if( cols >= 3 )  mixin( "result[2] = data[2] " ~ op ~ "ma t[2];" );
810 		static if( cols == 4 )  mixin( "result[3] = data[3] " ~ op ~ "ma t[3];" );
811 		return result;
812 	}
813 
814 	/// matrix-matrix multiplication, using matrix-vector multiplication for each column of mat
815 	Matrix opBinary( string op : "*", M )( M mat ) const if( isCompatibleMatrix!M && cols == M.rows ) {
816 		Matrix!( valueType, cols, M.cols ) result;
817 		result[0] = this * mat[0];
818 		result[1] = this * mat[1];
819 		static if( M.cols >= 3 )  result[2] = this * mat[2];
820 		static if( M.cols == 4 )  result[3] = this * mat[3];
821 		return result;
822 	}
823 
824 
825 	unittest {		/// Matrix-Matrix Operations
826 		mat2 m2 = mat2( 2.0f, 4.0f, 6.0f, 8.0f );
827 		assert(( m2 * m2 ) == [ [ 28.0f, 40.0f ], [ 60.0f, 88.0f ] ] );
828 		assert(( m2 - m2 ) == [ [ 0.0f, 0.0f ], [ 0.0f, 0.0f ] ] );
829 		mat2 m2_2 = m2 + m2;
830 		assert(  m2_2  == [ [ 4.0f, 8.0f ], [ 12.0f, 16.0f ] ] );
831 		assert(( m2_2 / m2 ) == [ vec2( 2.0f, 2.0f ), vec2( 2.0f, 2.0f ) ] );
832 
833 		mat3 m3 = mat3( 2.0f, 4.0f, 6.0f, 8.0f, 10.0f, 12.0f, 14.0f, 16.0f, 18.0f );
834 		assert(( m3 * m3 ) == [ [ 120.0f, 144.0f, 168.0f ], [ 264.0f, 324.0f, 384.0f ], [ 408.0f, 504.0f, 600.0f ] ] );
835 		assert(( m3 - m3 ) == [ [ 0.0f, 0.0f, 0.0f ], [ 0.0f, 0.0f, 0.0f ], [ 0.0f, 0.0f, 0.0f ] ] );
836 		assert(( m3 + m3 ) == [ [ 4.0f, 8.0f, 12.0f ], [ 16.0f, 20.0f, 24.0f ], [ 28.0f, 32.0f, 36.0f ] ] );
837 		assert(( m3 / m3 ) == [ [ 1.0f, 1.0f, 1.0f ], [ 1.0f, 1.0f, 1.0f ], vec3( 1.0f, 1.0f, 1.0f ) ] );
838 
839 		//TODO: tests for mat4, mat34
840 	}
841 
842 
843 	void opOpAssign( string op, T )( T val ) {
844 		mixin( "this = this " ~ op ~  "val;" );
845 	}
846 
847 	unittest {		/// Matrix Unary Operations
848 		mat2  m2 = mat2( 2.0f, 4.0f, 6.0f, 8.0f );
849 		m2 += m2; assert( m2 == [ [ 4.0f, 8.0f ], [ 12.0f, 16.0f ] ] );
850 		m2 /= m2; assert( m2 == [ [ 1.0f, 1.0f ], [ 1.0f, 1.0f ] ] );
851 		m2 -= m2; assert( m2 == [ [ 0.0f, 0.0f ], [ 0.0f, 0.0f ] ] );
852 
853 		mat3  m3 = mat3( 2.0f, 4.0f, 6.0f, 8.0f, 10.0f, 12.0f, 14.0f, 16.0f, 18.0f );
854 		m3 += m3; assert( m3 == [ [ 4.0f, 8.0f, 12.0f ], [ 16.0f, 20.0f, 24.0f ], [ 28.0f, 32.0f, 36.0f ] ] );
855 		m3 /= m3; assert( m3 == [ [ 1.0f, 1.0f, 1.0f ], [ 1.0f, 1.0f, 1.0f ], [ 1.0f, 1.0f, 1.0f ] ] );
856 		m3 -= m3; assert( m3 == [ [ 0.0f, 0.0f, 0.0f ], [ 0.0f, 0.0f, 0.0f ], [ 0.0f, 0.0f, 0.0f ] ] );
857 
858 		//TODO: tests for mat4, mat34
859 	}
860 
861 
862 	bool opCast( T : bool )() const {
863 		return isvalid;
864 	}
865 
866 	unittest {
867 		assert( mat2( 1.0f, 2.0f, 1.0f, 1.0f ) == mat2( 1.0f, 2.0f, 1.0f, 1.0f ));
868 		assert( mat2( 1.0f, 2.0f, 1.0f, 1.0f ) != mat2( 1.0f, 1.0f, 1.0f, 1.0f ));
869 
870 		assert( mat3( 1.0f ) == mat3( 1.0f ));
871 		assert( mat3( 1.0f ) != mat3( 2.0f ));
872 
873 		assert( mat4( 1.0f ) == mat4( 1.0f ));
874 		assert( mat4( 1.0f ) != mat4( 2.0f ));
875 
876 		assert( !( mat4( float.nan )));
877 		if( mat4( 1.0f )) { }
878 		else { assert( false ); }
879 	}
880 
881 }
882 
883 
884 /////////////////
885 /// Transpose ///
886 /////////////////
887 auto transpose( Matrix )( Matrix mat ) if( isMatrix!Matrix ) {
888 	Matrix!( Matrix.valueType, Matrix.rows, Matrix.cols ) result;
889 
890 	foreach( c; 0 .. Matrix.cols ) {
891 		foreach( r; 0 .. Matrix.rows ) {
892 			result[ r ][ c ] = mat[ c ][ r ];
893 		}
894 	}
895 
896 	return result;
897 }
898 
899 /// Unittest transpose
900 unittest {
901 	mat2 m2 = mat2( 1.0f );
902 	assert( m2.transpose.data == mat2( 1.0f ).data );
903 
904 	assert( m2.transpose.data  ==  [[ 1.0f, 0.0f ],
905 									[ 0.0f, 1.0f ]] );
906 	assert( m2.data == mat2.identity.data );
907 
908 	mat3 m3 = mat3( 1.1f, 1.2f, 1.3f,
909 					2.1f, 2.2f, 2.3f,
910 					3.1f, 3.2f, 3.3f );
911 	m3 = m3.transpose;
912 	assert( m3.data == [[ 1.1f, 2.1f, 3.1f ],
913 						[ 1.2f, 2.2f, 3.2f ],
914 						[ 1.3f, 2.3f, 3.3f ]] );
915 
916 	mat4 m4 = mat4( 2.0f );
917 	assert( m4.transpose.data == mat4( 2.0f ).data );
918 }
919 
920 
921 
922 
923 /////////////////
924 // Determinant //
925 /////////////////
926 auto determinant( Matrix )( Matrix mat ) if( isMatrix!Matrix && Matrix.cols == Matrix.rows ) {
927 	static if( mat.cols == 2 ) {
928 		return ( mat[ 0 ][ 0 ] * mat[ 1 ][ 1 ] - mat[ 0 ][ 1 ] * mat[ 1 ][ 0 ] );
929 	}
930 
931 	else static if( mat.cols == 3 ) {
932 		return (  mat[ 0 ][ 0 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 2 ]
933 				+ mat[ 0 ][ 1 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 0 ]
934 				+ mat[ 0 ][ 2 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 1 ]
935 				- mat[ 0 ][ 2 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 0 ]
936 				- mat[ 0 ][ 1 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 2 ]
937 				- mat[ 0 ][ 0 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 1 ] );
938 	}
939 
940 	else {
941 		return (  mat[ 0 ][ 3 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 0 ] - mat[ 0 ][ 2 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 0 ]
942 				- mat[ 0 ][ 3 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 0 ] + mat[ 0 ][ 1 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 0 ]
943 				+ mat[ 0 ][ 2 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 0 ] - mat[ 0 ][ 1 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 0 ]
944 				- mat[ 0 ][ 3 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 1 ] + mat[ 0 ][ 2 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 1 ]
945 				+ mat[ 0 ][ 3 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 1 ] - mat[ 0 ][ 0 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 1 ]
946 				- mat[ 0 ][ 2 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 1 ] + mat[ 0 ][ 0 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 1 ]
947 				+ mat[ 0 ][ 3 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 2 ] - mat[ 0 ][ 1 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 2 ]
948 				- mat[ 0 ][ 3 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 2 ] + mat[ 0 ][ 0 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 2 ]
949 				+ mat[ 0 ][ 1 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 2 ] - mat[ 0 ][ 0 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 2 ]
950 				- mat[ 0 ][ 2 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 3 ] + mat[ 0 ][ 1 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 3 ]
951 				+ mat[ 0 ][ 2 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 3 ] - mat[ 0 ][ 0 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 3 ]
952 				- mat[ 0 ][ 1 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 3 ] + mat[ 0 ][ 0 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 3 ] );
953 	}
954 }
955 
956 ////////////
957 // Invert //
958 ////////////
959 auto invert( Matrix )( Matrix mat ) if( isMatrix!Matrix && Matrix.cols == Matrix.rows ) {
960 
961 	static if( isFloatingPoint!( Matrix.valueType ) && rmul ) {
962 		Matrix.valueType d = 1 / mat.determinant;
963 		enum op = "*";
964 	}	else {
965 		Matrix.valueType d = mat.determinant;
966 		enum op = "/";
967 	}
968 
969 	Matrix result;
970 	alias vectorType = Matrix.vectorType;
971 
972 	static if( Matrix.cols == 2 )
973 		mixin( "
974 			result = [
975 				vectorType(   mat[ 1 ][ 1 ], - mat[ 0 ][ 1 ] ) " ~ op ~ " d,
976 				vectorType( - mat[ 1 ][ 0 ],   mat[ 0 ][ 0 ] ) " ~ op ~ " d ];" );
977 
978 	else static if( Matrix.cols == 3 )
979 		mixin( "
980 			result = [
981 				vectorType((  mat[ 1 ][ 1 ] * mat[ 2 ][ 2 ] - mat[ 1 ][ 2 ] * mat[ 2 ][ 1 ] ) " ~ op ~ " d,
982 							( mat[ 0 ][ 2 ] * mat[ 2 ][ 1 ] - mat[ 0 ][ 1 ] * mat[ 2 ][ 2 ] ) " ~ op ~ " d,
983 							( mat[ 0 ][ 1 ] * mat[ 1 ][ 2 ] - mat[ 0 ][ 2 ] * mat[ 1 ][ 1 ] ) " ~ op ~ " d ),
984 				vectorType((  mat[ 1 ][ 2 ] * mat[ 2 ][ 0 ] - mat[ 1 ][ 0 ] * mat[ 2 ][ 2 ] ) " ~ op ~ " d,
985 							( mat[ 0 ][ 0 ] * mat[ 2 ][ 2 ] - mat[ 0 ][ 2 ] * mat[ 2 ][ 0 ] ) " ~ op ~ " d,
986 							( mat[ 0 ][ 2 ] * mat[ 1 ][ 0 ] - mat[ 0 ][ 0 ] * mat[ 1 ][ 2 ] ) " ~ op ~ " d ),
987 				vectorType((  mat[ 1 ][ 0 ] * mat[ 2 ][ 1 ] - mat[ 1 ][ 1 ] * mat[ 2 ][ 0 ] ) " ~ op ~ " d,
988 							( mat[ 0 ][ 1 ] * mat[ 2 ][ 0 ] - mat[ 0 ][ 0 ] * mat[ 2 ][ 1 ] ) " ~ op ~ " d,
989 							( mat[ 0 ][ 0 ] * mat[ 1 ][ 1 ] - mat[ 0 ][ 1 ] * mat[ 1 ][ 0 ] ) " ~ op ~ " d ) ];" );
990 
991 	else
992 		mixin( "
993 			result = [
994 				vectorType((  mat[ 1 ][ 1 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 3 ] + mat[ 1 ][ 2 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 1 ] + mat[ 1 ][ 3 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 2 ]
995 							- mat[ 1 ][ 1 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 2 ] - mat[ 1 ][ 2 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 3 ] - mat[ 1 ][ 3 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 1 ] ) " ~ op ~ " d,
996 							( mat[ 0 ][ 1 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 2 ] + mat[ 0 ][ 2 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 3 ] + mat[ 0 ][ 3 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 1 ]
997 							- mat[ 0 ][ 1 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 3 ] - mat[ 0 ][ 2 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 1 ] - mat[ 0 ][ 3 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 2 ] ) " ~ op ~ " d,
998 							( mat[ 0 ][ 1 ] * mat[ 1 ][ 2 ] * mat[ 3 ][ 3 ] + mat[ 0 ][ 2 ] * mat[ 1 ][ 3 ] * mat[ 3 ][ 1 ] + mat[ 0 ][ 3 ] * mat[ 1 ][ 1 ] * mat[ 3 ][ 2 ]
999 							- mat[ 0 ][ 1 ] * mat[ 1 ][ 3 ] * mat[ 3 ][ 2 ] - mat[ 0 ][ 2 ] * mat[ 1 ][ 1 ] * mat[ 3 ][ 3 ] - mat[ 0 ][ 3 ] * mat[ 1 ][ 2 ] * mat[ 3 ][ 1 ] ) " ~ op ~ " d,
1000 							( mat[ 0 ][ 1 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 2 ] + mat[ 0 ][ 2 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 3 ] + mat[ 0 ][ 3 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 1 ]
1001 							- mat[ 0 ][ 1 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 3 ] - mat[ 0 ][ 2 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 1 ] - mat[ 0 ][ 3 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 2 ] ) " ~ op ~ " d ),
1002 				vectorType((  mat[ 1 ][ 0 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 2 ] + mat[ 1 ][ 2 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 3 ] + mat[ 1 ][ 3 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 0 ]
1003 							- mat[ 1 ][ 0 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 3 ] - mat[ 1 ][ 2 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 0 ] - mat[ 1 ][ 3 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 2 ] ) " ~ op ~ " d,
1004 							( mat[ 0 ][ 0 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 3 ] + mat[ 0 ][ 2 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 0 ] + mat[ 0 ][ 3 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 2 ]
1005 							- mat[ 0 ][ 0 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 2 ] - mat[ 0 ][ 2 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 3 ] - mat[ 0 ][ 3 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 0 ] ) " ~ op ~ " d,
1006 							( mat[ 0 ][ 0 ] * mat[ 1 ][ 3 ] * mat[ 3 ][ 2 ] + mat[ 0 ][ 2 ] * mat[ 1 ][ 0 ] * mat[ 3 ][ 3 ] + mat[ 0 ][ 3 ] * mat[ 1 ][ 2 ] * mat[ 3 ][ 0 ]
1007 							- mat[ 0 ][ 0 ] * mat[ 1 ][ 2 ] * mat[ 3 ][ 3 ] - mat[ 0 ][ 2 ] * mat[ 1 ][ 3 ] * mat[ 3 ][ 0 ] - mat[ 0 ][ 3 ] * mat[ 1 ][ 0 ] * mat[ 3 ][ 2 ] ) " ~ op ~ " d,
1008 							( mat[ 0 ][ 0 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 3 ] + mat[ 0 ][ 2 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 0 ] + mat[ 0 ][ 3 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 2 ]
1009 							- mat[ 0 ][ 0 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 2 ] - mat[ 0 ][ 2 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 3 ] - mat[ 0 ][ 3 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 0 ] ) " ~ op ~ " d ),
1010 				vectorType((  mat[ 1 ][ 0 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 3 ] + mat[ 1 ][ 1 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 0 ] + mat[ 1 ][ 3 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 1 ]
1011 							- mat[ 1 ][ 0 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 1 ] - mat[ 1 ][ 1 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 3 ] - mat[ 1 ][ 3 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 0 ] ) " ~ op ~ " d,
1012 							( mat[ 0 ][ 0 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 1 ] + mat[ 0 ][ 1 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 3 ] + mat[ 0 ][ 3 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 0 ]
1013 							- mat[ 0 ][ 0 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 3 ] - mat[ 0 ][ 1 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 0 ] - mat[ 0 ][ 3 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 1 ] ) " ~ op ~ " d,
1014 							( mat[ 0 ][ 0 ] * mat[ 1 ][ 1 ] * mat[ 3 ][ 3 ] + mat[ 0 ][ 1 ] * mat[ 1 ][ 3 ] * mat[ 3 ][ 0 ] + mat[ 0 ][ 3 ] * mat[ 1 ][ 0 ] * mat[ 3 ][ 1 ]
1015 							- mat[ 0 ][ 0 ] * mat[ 1 ][ 3 ] * mat[ 3 ][ 1 ] - mat[ 0 ][ 1 ] * mat[ 1 ][ 0 ] * mat[ 3 ][ 3 ] - mat[ 0 ][ 3 ] * mat[ 1 ][ 1 ] * mat[ 3 ][ 0 ] ) " ~ op ~ " d,
1016 							( mat[ 0 ][ 0 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 1 ] + mat[ 0 ][ 1 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 3 ] + mat[ 0 ][ 3 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 0 ]
1017 							- mat[ 0 ][ 0 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 3 ] - mat[ 0 ][ 1 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 0 ] - mat[ 0 ][ 3 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 1 ] ) " ~ op ~ " d ),
1018 				vectorType((  mat[ 1 ][ 0 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 1 ] + mat[ 1 ][ 1 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 2 ] + mat[ 1 ][ 2 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 0 ]
1019 							- mat[ 1 ][ 0 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 2 ] - mat[ 1 ][ 1 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 0 ] - mat[ 1 ][ 2 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 1 ] ) " ~ op ~ " d,
1020 							( mat[ 0 ][ 0 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 2 ] + mat[ 0 ][ 1 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 0 ] + mat[ 0 ][ 2 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 1 ]
1021 							- mat[ 0 ][ 0 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 1 ] - mat[ 0 ][ 1 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 2 ] - mat[ 0 ][ 2 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 0 ] ) " ~ op ~ " d,
1022 							( mat[ 0 ][ 0 ] * mat[ 1 ][ 2 ] * mat[ 3 ][ 1 ] + mat[ 0 ][ 1 ] * mat[ 1 ][ 0 ] * mat[ 3 ][ 2 ] + mat[ 0 ][ 2 ] * mat[ 1 ][ 1 ] * mat[ 3 ][ 0 ]
1023 							- mat[ 0 ][ 0 ] * mat[ 1 ][ 1 ] * mat[ 3 ][ 2 ] - mat[ 0 ][ 1 ] * mat[ 1 ][ 2 ] * mat[ 3 ][ 0 ] - mat[ 0 ][ 2 ] * mat[ 1 ][ 0 ] * mat[ 3 ][ 1 ] ) " ~ op ~ " d,
1024 							( mat[ 0 ][ 0 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 2 ] + mat[ 0 ][ 1 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 0 ] + mat[ 0 ][ 2 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 1 ]
1025 							- mat[ 0 ][ 0 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 1 ] - mat[ 0 ][ 1 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 2 ] - mat[ 0 ][ 2 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 0 ] ) " ~ op ~ " d ) ];" );
1026 
1027 	return result;
1028 }
1029 
1030 
1031 ///////////////////////////////////////////////
1032 // Invert purely Translate-Rotate 4*4 matrix //
1033 ///////////////////////////////////////////////
1034 auto invertTR( Matrix )( Matrix mat ) if( isMatrix!Matrix && Matrix.cols == 4 && Matrix.rows == 4 ) {
1035 	Matrix result = mat.rotation.transpose;
1036 	result[3] = - result.opBinary!( "*", vec4 )( mat[3] );	// why does result * mat3[3] not work ???
1037 	return result;
1038 }
1039 
1040 
1041 unittest {
1042 	mat2 m2 = mat2( 1.0f, 2.0f, vec2( 3.0f, 4.0f ));
1043 	assert( m2.determinant == - 2.0f );
1044 	assert( m2.invert.data ==  [[ - 2.0f,   1.0f ],
1045 								[   1.5f, - 0.5f ]  );
1046 
1047 	mat3 m3 = mat3( 1.0f, - 2.0f,   3.0f,
1048 					7.0f, - 1.0f,   0.0f,
1049 					3.0f,   2.0f, - 4.0f );
1050 	assert( m3.determinant == - 1.0f );
1051 	assert( m3.invert.data ==  [[ -  4.0f,  2.0f, -  3.0f ],
1052 								[ - 28.0f, 13.0f, - 21.0f ],
1053 								[ - 17.0f,  8.0f, - 13.0f ]] );
1054 
1055 	mat4 m4 = mat4(	  1.0f,   2.0f,   3.0f,   4.0f,
1056 					- 2.0f,   1.0f,   5.0f, - 2.0f,
1057 					  2.0f, - 1.0f,   7.0f,   1.0f,
1058 					  3.0f, - 3.0f,   2.0f,   0.0f );
1059 	assert( m4.determinant == - 8.0f );
1060 	assert( m4.invert.data ==  [[   6.875f,   7.875f, - 11.75f,  11.125f ],
1061 								[   6.625f,   7.625f, - 11.25f,  10.375f ],
1062 								[ - 0.375f, - 0.375f,    0.75f, - 0.625f ],
1063 								[ - 4.5f,   - 5.5f,      8.0f,  - 7.5f ]] );
1064 }
1065 
1066 
1067 /// query if any entry is nan
1068 bool isnan( genType )( const ref genType mat ) if( isMatrix!genType && isFloatingPoint!( genType.valueType )) {
1069 	foreach( const ref vec; mat )
1070 		if( vec.isnan )
1071 			return true;
1072 	return false;
1073 }
1074 
1075 
1076 /// query if any entry is inf
1077 bool isinf( genType )( const ref genType mat ) if( isMatrix!genType && isFloatingPoint!( genType.valueType )) {
1078 	foreach( const ref vec; mat )
1079 		if( vec.isinf )
1080 			return true;
1081 	return false;
1082 }
1083 
1084 
1085 /// query if all entries are not nan and not inf
1086 bool isvalid( genType )( const ref genType mat ) if( isMatrix!genType && isFloatingPoint!( genType.valueType )) {
1087 	return !( mat.isinf || mat.isnan );
1088 }
1089 
1090 
1091 unittest {		/// Matrix.isvalid
1092 
1093 	mat2 m2 = mat2( 1.0f, 1.0f, vec2( 2.0f, 2.0f ));
1094 	assert( m2.data == [[ 1.0f, 1.0f ], [ 2.0f, 2.0f ]] );
1095 	m2.clear( 3.0f );
1096 	assert( m2.data == [[ 3.0f, 3.0f ], [ 3.0f, 3.0f ]] );
1097 	assert( m2.isvalid );
1098 	m2.clear( float.nan );
1099 	assert( !m2.isvalid );
1100 	m2.clear( float.infinity );
1101 	assert( !m2.isvalid );
1102 	m2.clear( 0.0f );
1103 	assert( m2.isvalid );
1104 
1105 	mat3 m3 = mat3( 1.0f );
1106 	assert( m3.data == [[ 1.0f, 0.0f, 0.0f ],
1107 						[ 0.0f, 1.0f, 0.0f ],
1108 						[ 0.0f, 0.0f, 1.0f ]] );
1109 
1110 	mat4 m4 = mat4( vec4( 	1.0f, 1.0f, 1.0f, 1.0f ),
1111 							2.0f, 2.0f, 2.0f, 2.0f,
1112 							3.0f, 3.0f, 3.0f, 3.0f,
1113 						 vec4( 4.0f, 4.0f, 4.0f, 4.0f ));
1114 	assert( m4.data == [[ 1.0f, 1.0f, 1.0f, 1.0f ],
1115 						[ 2.0f, 2.0f, 2.0f, 2.0f ],
1116 						[ 3.0f, 3.0f, 3.0f, 3.0f ],
1117 						[ 4.0f, 4.0f, 4.0f, 4.0f ]] );
1118 	assert( mat3( m4 ).data == [[ 1.0f, 1.0f, 1.0f ],
1119 								[ 2.0f, 2.0f, 2.0f ],
1120 								[ 3.0f, 3.0f, 3.0f ]] );
1121 	assert( mat2( mat3( m4 )).data == [[ 1.0f, 1.0f ],
1122 													[ 2.0f, 2.0f ] ] );
1123 	assert( mat2( m4 ).data == mat2( mat3( m4 )).data );
1124 	assert( mat4( mat3( m4 )).data ==  [[ 1.0f, 1.0f, 1.0f, 0.0f ],
1125 										[ 2.0f, 2.0f, 2.0f, 0.0f ],
1126 										[ 3.0f, 3.0f, 3.0f, 0.0f ],
1127 										[ 0.0f, 0.0f, 0.0f, 1.0f ]] );
1128 
1129 	Matrix!( float, 2, 3 ) mt1 = Matrix!( float, 2, 3 )( 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f );
1130 	Matrix!( float, 3, 2 ) mt2 = Matrix!( float, 3, 2 )( 6.0f, - 1.0f, 3.0f, 2.0f, 0.0f, - 3.0f );
1131 
1132 	assert( mt1.data == [[ 1.0f, 2.0f, 3.0f ], [ 4.0f, 5.0f, 6.0f ]] );
1133 	assert( mt2.data == [[ 6.0f, - 1.0f ], [ 3.0f, 2.0f ], [ 0.0f, - 3.0f ]] );
1134 }
1135 
1136 
1137 
1138