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