1 module gfm.math.matrix;
2 
3 import std.math,
4        std.typetuple,
5        std.traits,
6        std.string,
7        std.typecons,
8        std.conv;
9 
10 import gfm.math.vector,
11        gfm.math.shapes,
12        gfm.math.quaternion;
13 
14 static if( __VERSION__ < 2066 ) private enum nogc = 1;
15 
16 /// Generic non-resizeable matrix with R rows and C columns.
17 /// Intended for 3D use (size 3x3 and 4x4).
18 /// Important: <b>Matrices here are in row-major order whereas OpenGL is column-major.</b>
19 /// Params:
20 ///   T = type of elements
21 ///   R = number of rows
22 ///   C = number of columns
23 align(1) struct Matrix(T, int R, int C)
24 {
25     align(1):
26     public
27     {
28         static assert(R >= 1 && C >= 1);
29 
30         alias Vector!(T, C) row_t;
31         alias Vector!(T, R) column_t;
32 
33         enum bool isSquare = (R == C);
34 
35         // fields definition
36         union
37         {
38             T[C*R] v;        // all elements
39             row_t[R] rows;   // all rows
40             T[C][R] c;       // components
41         }
42 
43         @nogc this(U...)(U values) pure nothrow
44         {
45             static if ((U.length == C*R) && allSatisfy!(isTConvertible, U))
46             {
47                 // construct with components
48                 foreach(int i, x; values)
49                     v[i] = x;
50             }
51             else static if ((U.length == 1) && (isAssignable!(U[0])) && (!is(U[0] : Matrix)))
52             {
53                 // construct with assignment
54                 opAssign!(U[0])(values[0]);
55             }
56             else static assert(false, "cannot create a matrix from given arguments");
57         }
58 
59         /// Construct a matrix from columns.
60         @nogc static Matrix fromColumns(column_t[] columns) pure nothrow
61         {
62             assert(columns.length == C);
63             Matrix res;
64             for (int i = 0; i < R; ++i)
65                 for (int j = 0; j < C; ++j)
66                 {
67                    res.c[i][j] = columns[j][i];
68                 }
69             return res;
70         }
71 
72         /// Construct a matrix from rows.
73         @nogc static Matrix fromRows(row_t[] rows) pure nothrow
74         {
75             assert(rows.length == R);
76             Matrix res;
77             res.rows[] = rows[];
78             return res;
79         }
80 
81         /// Construct matrix with a scalar.
82         @nogc this(U)(T x) pure nothrow
83         {
84             for (int i = 0; i < _N; ++i)
85                 v[i] = x;
86         }
87 
88         /// Assign with a samey matrice.
89         @nogc ref Matrix opAssign(U : Matrix)(U x) pure nothrow
90         {
91             for (int i = 0; i < R * C; ++i)
92                 v[i] = x.v[i];
93             return this;
94         }
95 
96         /// Assign from other small matrices (same size, compatible type).
97         @nogc ref Matrix opAssign(U)(U x) pure nothrow
98             if (is(typeof(U._isMatrix))
99                 && is(U._T : _T)
100                 && (!is(U: Matrix))
101                 && (U._R == R) && (U._C == C))
102         {
103             for (int i = 0; i < R * C; ++i)
104                 v[i] = x.v[i];
105             return this;
106         }
107 
108         /// Assign with a static array of size R * C.
109         @nogc ref Matrix opAssign(U)(U x) pure nothrow
110             if ((isStaticArray!U)
111                 && is(typeof(x[0]) : T)
112                 && (U.length == R * C))
113         {
114             for (int i = 0; i < R * C; ++i)
115                 v[i] = x[i];
116             return this;
117         }
118 
119         /// Assign with a dynamic array of size R * C.
120         @nogc ref Matrix opAssign(U)(U x) pure nothrow
121             if ((isDynamicArray!U)
122                 && is(typeof(x[0]) : T))
123         {
124             assert(x.length == R * C);
125             for (int i = 0; i < R * C; ++i)
126                 v[i] = x[i];
127             return this;
128         }
129 
130         /// Return a pointer to content.
131         @nogc inout(T)* ptr() pure inout nothrow @property
132         {
133             return v.ptr;
134         }
135 
136         /// Returns: column j as a vector.
137         @nogc column_t column(int j) pure const nothrow
138         {
139             column_t res = void;
140             for (int i = 0; i < R; ++i)
141                 res.v[i] = c[i][j];
142             return res;
143         }
144 
145         /// Returns: row i as a vector.
146         @nogc row_t row(int i) pure const nothrow
147         {
148             return rows[i];
149         }
150 
151         /// Covnerts to pretty string.
152         string toString() const nothrow
153         {
154             try
155                 return format("%s", v);
156             catch (Exception e)
157                 assert(false); // should not happen since format is right
158         }
159 
160         /// Matrix * vector multiplication.
161         @nogc column_t opBinary(string op)(row_t x) pure const nothrow if (op == "*")
162         {
163             column_t res = void;
164             for (int i = 0; i < R; ++i)
165             {
166                 T sum = 0;
167                 for (int j = 0; j < C; ++j)
168                 {
169                     sum += c[i][j] * x.v[j];
170                 }
171                 res.v[i] = sum;
172             }
173             return res;
174         }
175 
176         /// Matrix * matrix multiplication.
177         @nogc auto opBinary(string op, U)(U x) pure const nothrow
178             if (is(typeof(U._isMatrix)) && (U._R == C) && (op == "*"))
179         {
180             Matrix!(T, R, U._C) result = void;
181 
182             for (int i = 0; i < R; ++i)
183             {
184                 for (int j = 0; j < U._C; ++j)
185                 {
186                     T sum = 0;
187                     for (int k = 0; k < C; ++k)
188                         sum += c[i][k] * x.c[k][j];
189                     result.c[i][j] = sum;
190                 }
191             }
192             return result;
193         }
194 
195         /// Matrix add and substraction.
196         @nogc Matrix opBinary(string op, U)(U other) pure const nothrow
197             if (is(U : Matrix) && (op == "+" || op == "-"))
198         {
199             Matrix result = void;
200 
201             for (int i = 0; i < R; ++i)
202             {
203                 for (int j = 0; j < C; ++j)
204                 {
205                     mixin("result.c[i][j] = c[i][j] " ~ op ~ " other.c[i][j];");
206                 }
207             }
208             return result;
209         }
210 
211         /// Assignment operator with another samey matrix.
212         @nogc ref Matrix opOpAssign(string op, U)(U operand) pure nothrow if (is(U : Matrix))
213         {
214             mixin("Matrix result = this " ~ op ~ " operand;");
215             return opAssign!Matrix(result);
216         }
217 
218         /// Assignment operator with another samey matrix.
219         @nogc ref Matrix opOpAssign(string op, U)(U operand) pure nothrow if (isConvertible!U)
220         {
221             Matrix conv = operand;
222             return opOpAssign!op(conv);
223         }
224 
225         /// Cast to other matrix types.
226         /// If the size are different, the resulting matrix is truncated
227         /// and/or filled with identity coefficients.
228         @nogc U opCast(U)() pure const nothrow if (is(typeof(U._isMatrix)))
229         {
230             U res = U.identity();
231             enum minR = R < U._R ? R : U._R;
232             enum minC = C < U._C ? C : U._C;
233             for (int i = 0; i < minR; ++i)
234                 for (int j = 0; j < minC; ++j)
235                 {
236                     res.c[i][j] = cast(U._T)(c[i][j]);
237                 }
238             return res;
239         }
240 
241         @nogc bool opEquals(U)(U other) pure const nothrow if (is(U : Matrix))
242         {
243             for (int i = 0; i < R * C; ++i)
244                 if (v[i] != other.v[i])
245                     return false;
246             return true;
247         }
248 
249         @nogc bool opEquals(U)(U other) pure const nothrow
250             if ((isAssignable!U) && (!is(U: Matrix)))
251         {
252             Matrix conv = other;
253             return opEquals(conv);
254         }
255 
256         // +matrix, -matrix, ~matrix, !matrix
257         @nogc Matrix opUnary(string op)() pure const nothrow if (op == "+" || op == "-" || op == "~" || op == "!")
258         {
259             Matrix res = void;
260             for (int i = 0; i < N; ++i)
261                 mixin("res.v[i] = " ~ op ~ "v[i];");
262             return res;
263         }
264 
265         /// Convert 3x3 rotation matrix to quaternion.
266         /// See_also: 3D Math Primer for Graphics and Game Development.
267         @nogc U opCast(U)() pure const nothrow if (is(typeof(U._isQuaternion))
268                                                    && is(U._T : _T)
269                                                    && (_R == 3) && (_C == 3))
270         {
271             T fourXSquaredMinus1 = c[0][0] - c[1][1] - c[2][2];
272             T fourYSquaredMinus1 = c[1][1] - c[0][0] - c[2][2];
273             T fourZSquaredMinus1 = c[2][2] - c[0][0] - c[1][1];
274             T fourWSquaredMinus1 = c[0][0] + c[1][1] + c[2][2];
275 
276             int biggestIndex = 0;
277             T fourBiggestSquaredMinus1 = fourWSquaredMinus1;
278 
279             if(fourXSquaredMinus1 > fourBiggestSquaredMinus1)
280             {
281                 fourBiggestSquaredMinus1 = fourXSquaredMinus1;
282                 biggestIndex = 1;
283             }
284 
285             if(fourYSquaredMinus1 > fourBiggestSquaredMinus1)
286             {
287                 fourBiggestSquaredMinus1 = fourYSquaredMinus1;
288                 biggestIndex = 2;
289             }
290 
291             if(fourZSquaredMinus1 > fourBiggestSquaredMinus1)
292             {
293                 fourBiggestSquaredMinus1 = fourZSquaredMinus1;
294                 biggestIndex = 3;
295             }
296 
297             T biggestVal = sqrt(fourBiggestSquaredMinus1 + 1) / 2;
298             T mult = 1 / (biggestVal * 4);
299 
300             U quat;
301             switch(biggestIndex)
302             {
303                 case 1:
304                     quat.w = (c[1][2] - c[2][1]) * mult;
305                     quat.x = biggestVal;
306                     quat.y = (c[0][1] + c[1][0]) * mult;
307                     quat.z = (c[2][0] + c[0][2]) * mult;
308                     break;
309 
310                 case 2:
311                     quat.w = (c[2][0] - c[0][2]) * mult;
312                     quat.x = (c[0][1] + c[1][0]) * mult;
313                     quat.y = biggestVal;
314                     quat.z = (c[1][2] + c[2][1]) * mult;
315                     break;
316 
317                 case 3:
318                     quat.w = (c[0][1] - c[1][0]) * mult;
319                     quat.x = (c[2][0] + c[0][2]) * mult;
320                     quat.y = (c[1][2] + c[2][1]) * mult;
321                     quat.z = biggestVal;
322                     break;
323 
324                 default: // biggestIndex == 0
325                     quat.w = biggestVal;
326                     quat.x = (c[1][2] - c[2][1]) * mult;
327                     quat.y = (c[2][0] - c[0][2]) * mult;
328                     quat.z = (c[0][1] - c[1][0]) * mult;
329                     break;
330             }
331 
332             return quat;
333         }
334 
335         /// Converts a 4x4 rotation matrix to quaternion.
336         @nogc U opCast(U)() pure const nothrow if (is(typeof(U._isQuaternion))
337                                                    && is(U._T : _T)
338                                                    && (_R == 4) && (_C == 4))
339         {
340             auto m3 = cast(mat3!T)(this);
341             return cast(U)(m3);
342         }
343 
344         /// Matrix inversion is provided for 2x2, 3x3 and 4x4 floating point matrices.
345 
346         static if (isSquare && isFloatingPoint!T && R == 2)
347         {
348             /// Returns: inverse of matrix.
349             @nogc Matrix inverse() pure const nothrow
350             {
351                 T invDet = 1 / (c[0][0] * c[1][1] - c[0][1] * c[1][0]);
352                 return Matrix( c[1][1] * invDet, -c[0][1] * invDet,
353                                    -c[1][0] * invDet,  c[0][0] * invDet);
354             }
355         }
356 
357         static if (isSquare && isFloatingPoint!T && R == 3)
358         {
359             /// Returns: inverse of matrix.
360             @nogc Matrix inverse() pure const nothrow
361             {
362                 T det = c[0][0] * (c[1][1] * c[2][2] - c[2][1] * c[1][2])
363                       - c[0][1] * (c[1][0] * c[2][2] - c[1][2] * c[2][0])
364                       + c[0][2] * (c[1][0] * c[2][1] - c[1][1] * c[2][0]);
365                 T invDet = 1 / det;
366 
367                 Matrix res = void;
368                 res.c[0][0] =  (c[1][1] * c[2][2] - c[2][1] * c[1][2]) * invDet;
369                 res.c[0][1] = -(c[0][1] * c[2][2] - c[0][2] * c[2][1]) * invDet;
370                 res.c[0][2] =  (c[0][1] * c[1][2] - c[0][2] * c[1][1]) * invDet;
371                 res.c[1][0] = -(c[1][0] * c[2][2] - c[1][2] * c[2][0]) * invDet;
372                 res.c[1][1] =  (c[0][0] * c[2][2] - c[0][2] * c[2][0]) * invDet;
373                 res.c[1][2] = -(c[0][0] * c[1][2] - c[1][0] * c[0][2]) * invDet;
374                 res.c[2][0] =  (c[1][0] * c[2][1] - c[2][0] * c[1][1]) * invDet;
375                 res.c[2][1] = -(c[0][0] * c[2][1] - c[2][0] * c[0][1]) * invDet;
376                 res.c[2][2] =  (c[0][0] * c[1][1] - c[1][0] * c[0][1]) * invDet;
377                 return res;
378             }
379         }
380 
381         static if (isSquare && isFloatingPoint!T && R == 4)
382         {
383             /// Returns: inverse of matrix.
384             @nogc Matrix inverse() pure const nothrow
385             {
386                 T det2_01_01 = c[0][0] * c[1][1] - c[0][1] * c[1][0];
387                 T det2_01_02 = c[0][0] * c[1][2] - c[0][2] * c[1][0];
388                 T det2_01_03 = c[0][0] * c[1][3] - c[0][3] * c[1][0];
389                 T det2_01_12 = c[0][1] * c[1][2] - c[0][2] * c[1][1];
390                 T det2_01_13 = c[0][1] * c[1][3] - c[0][3] * c[1][1];
391                 T det2_01_23 = c[0][2] * c[1][3] - c[0][3] * c[1][2];
392 
393                 T det3_201_012 = c[2][0] * det2_01_12 - c[2][1] * det2_01_02 + c[2][2] * det2_01_01;
394                 T det3_201_013 = c[2][0] * det2_01_13 - c[2][1] * det2_01_03 + c[2][3] * det2_01_01;
395                 T det3_201_023 = c[2][0] * det2_01_23 - c[2][2] * det2_01_03 + c[2][3] * det2_01_02;
396                 T det3_201_123 = c[2][1] * det2_01_23 - c[2][2] * det2_01_13 + c[2][3] * det2_01_12;
397 
398                 T det = - det3_201_123 * c[3][0] + det3_201_023 * c[3][1] - det3_201_013 * c[3][2] + det3_201_012 * c[3][3];
399                 T invDet = 1 / det;
400 
401                 T det2_03_01 = c[0][0] * c[3][1] - c[0][1] * c[3][0];
402                 T det2_03_02 = c[0][0] * c[3][2] - c[0][2] * c[3][0];
403                 T det2_03_03 = c[0][0] * c[3][3] - c[0][3] * c[3][0];
404                 T det2_03_12 = c[0][1] * c[3][2] - c[0][2] * c[3][1];
405                 T det2_03_13 = c[0][1] * c[3][3] - c[0][3] * c[3][1];
406                 T det2_03_23 = c[0][2] * c[3][3] - c[0][3] * c[3][2];
407                 T det2_13_01 = c[1][0] * c[3][1] - c[1][1] * c[3][0];
408                 T det2_13_02 = c[1][0] * c[3][2] - c[1][2] * c[3][0];
409                 T det2_13_03 = c[1][0] * c[3][3] - c[1][3] * c[3][0];
410                 T det2_13_12 = c[1][1] * c[3][2] - c[1][2] * c[3][1];
411                 T det2_13_13 = c[1][1] * c[3][3] - c[1][3] * c[3][1];
412                 T det2_13_23 = c[1][2] * c[3][3] - c[1][3] * c[3][2];
413 
414                 T det3_203_012 = c[2][0] * det2_03_12 - c[2][1] * det2_03_02 + c[2][2] * det2_03_01;
415                 T det3_203_013 = c[2][0] * det2_03_13 - c[2][1] * det2_03_03 + c[2][3] * det2_03_01;
416                 T det3_203_023 = c[2][0] * det2_03_23 - c[2][2] * det2_03_03 + c[2][3] * det2_03_02;
417                 T det3_203_123 = c[2][1] * det2_03_23 - c[2][2] * det2_03_13 + c[2][3] * det2_03_12;
418 
419                 T det3_213_012 = c[2][0] * det2_13_12 - c[2][1] * det2_13_02 + c[2][2] * det2_13_01;
420                 T det3_213_013 = c[2][0] * det2_13_13 - c[2][1] * det2_13_03 + c[2][3] * det2_13_01;
421                 T det3_213_023 = c[2][0] * det2_13_23 - c[2][2] * det2_13_03 + c[2][3] * det2_13_02;
422                 T det3_213_123 = c[2][1] * det2_13_23 - c[2][2] * det2_13_13 + c[2][3] * det2_13_12;
423 
424                 T det3_301_012 = c[3][0] * det2_01_12 - c[3][1] * det2_01_02 + c[3][2] * det2_01_01;
425                 T det3_301_013 = c[3][0] * det2_01_13 - c[3][1] * det2_01_03 + c[3][3] * det2_01_01;
426                 T det3_301_023 = c[3][0] * det2_01_23 - c[3][2] * det2_01_03 + c[3][3] * det2_01_02;
427                 T det3_301_123 = c[3][1] * det2_01_23 - c[3][2] * det2_01_13 + c[3][3] * det2_01_12;
428 
429                 Matrix res = void;
430                 res.c[0][0] = - det3_213_123 * invDet;
431                 res.c[1][0] = + det3_213_023 * invDet;
432                 res.c[2][0] = - det3_213_013 * invDet;
433                 res.c[3][0] = + det3_213_012 * invDet;
434 
435                 res.c[0][1] = + det3_203_123 * invDet;
436                 res.c[1][1] = - det3_203_023 * invDet;
437                 res.c[2][1] = + det3_203_013 * invDet;
438                 res.c[3][1] = - det3_203_012 * invDet;
439 
440                 res.c[0][2] = + det3_301_123 * invDet;
441                 res.c[1][2] = - det3_301_023 * invDet;
442                 res.c[2][2] = + det3_301_013 * invDet;
443                 res.c[3][2] = - det3_301_012 * invDet;
444 
445                 res.c[0][3] = - det3_201_123 * invDet;
446                 res.c[1][3] = + det3_201_023 * invDet;
447                 res.c[2][3] = - det3_201_013 * invDet;
448                 res.c[3][3] = + det3_201_012 * invDet;
449                 return res;
450             }
451         }
452 
453         /// Returns: transposed matrice.
454         @nogc Matrix!(T, C, R) transposed() pure const nothrow
455         {
456             Matrix!(T, C, R) res;
457             for (int i = 0; i < C; ++i)
458                 for (int j = 0; j < R; ++j)
459                     res.c[i][j] = c[j][i];
460             return res;
461         }
462 
463         static if (isSquare && R > 1)
464         {
465             /// In-place translate by (v, 1)
466             @nogc void translate(Vector!(T, R-1) v) pure nothrow
467             {
468                 for (int i = 0; i < R; ++i)
469                 {
470                     T dot = 0;
471                     for (int j = 0; j + 1 < C; ++j)
472                         dot += v.v[j] * c[i][j];
473 
474                     c[i][C-1] += dot;
475                 }
476             }
477 
478             /// Make a translation matrix.
479             @nogc static Matrix translation(Vector!(T, R-1) v) pure nothrow
480             {
481                 Matrix res = identity();
482                 for (int i = 0; i + 1 < R; ++i)
483                     res.c[i][C-1] += v.v[i];
484                 return res;
485             }
486 
487             /// In-place matrix scaling.
488             void scale(Vector!(T, R-1) v) pure nothrow
489             {
490                 for (int i = 0; i < R; ++i)
491                     for (int j = 0; j + 1 < C; ++j)
492                         c[i][j] *= v.v[j];
493             }
494 
495             /// Make a scaling matrix.
496             @nogc static Matrix scaling(Vector!(T, R-1) v) pure nothrow
497             {
498                 Matrix res = identity();
499                 for (int i = 0; i + 1 < R; ++i)
500                     res.c[i][i] = v.v[i];
501                 return res;
502             }
503         }
504 
505         // rotations are implemented for 3x3 and 4x4 matrices.
506         static if (isSquare && (R == 3 || R == 4) && isFloatingPoint!T)
507         {
508             @nogc public static Matrix rotateAxis(int i, int j)(T angle) pure nothrow
509             {
510                 Matrix res = identity();
511                 const T cosa = cos(angle);
512                 const T sina = sin(angle);
513                 res.c[i][i] = cosa;
514                 res.c[i][j] = -sina;
515                 res.c[j][i] = sina;
516                 res.c[j][j] = cosa;
517                 return res;
518             }
519 
520             /// Returns: rotation matrix along axis X
521             alias rotateAxis!(1, 2) rotateX;
522 
523             /// Returns: rotation matrix along axis Y
524             alias rotateAxis!(2, 0) rotateY;
525 
526             /// Returns: rotation matrix along axis Z
527             alias rotateAxis!(0, 1) rotateZ;
528 
529             /// Similar to the glRotate matrix, however the angle is expressed in radians
530             /// See_also: $(LINK http://www.cs.rutgers.edu/~decarlo/428/gl_man/rotate.html)
531             @nogc static Matrix rotation(T angle, vec3!T axis) pure nothrow
532             {
533                 Matrix res = identity();
534                 const T c = cos(angle);
535                 const oneMinusC = 1 - c;
536                 const T s = sin(angle);
537                 axis = axis.normalized();
538                 T x = axis.x,
539                   y = axis.y,
540                   z = axis.z;
541                 T xy = x * y,
542                   yz = y * z,
543                   xz = x * z;
544 
545                 res.c[0][0] = x * x * oneMinusC + c;
546                 res.c[0][1] = x * y * oneMinusC - z * s;
547                 res.c[0][2] = x * z * oneMinusC + y * s;
548                 res.c[1][0] = y * x * oneMinusC + z * s;
549                 res.c[1][1] = y * y * oneMinusC + c;
550                 res.c[1][2] = y * z * oneMinusC - x * s;
551                 res.c[2][0] = z * x * oneMinusC - y * s;
552                 res.c[2][1] = z * y * oneMinusC + x * s;
553                 res.c[2][2] = z * z * oneMinusC + c;
554                 return res;
555             }
556         }
557 
558         // 4x4 specific transformations for 3D usage
559         static if (isSquare && R == 4 && isFloatingPoint!T)
560         {
561             /// Returns: orthographic projection.
562             @nogc static Matrix orthographic(T left, T right, T bottom, T top, T near, T far) pure nothrow
563             {
564                 T dx = right - left,
565                   dy = top - bottom,
566                   dz = far - near;
567 
568                 T tx = -(right + left) / dx;
569                 T ty = -(top + bottom) / dy;
570                 T tz = -(far + near)   / dz;
571 
572                 return Matrix(2 / dx,   0,      0,    tx,
573                                 0,    2 / dy,   0,    ty,
574                                 0,      0,    2 / dz, tz,
575                                 0,      0,      0,     1);
576             }
577 
578             /// Returns: perspective projection.
579             @nogc static Matrix perspective(T FOVInRadians, T aspect, T zNear, T zFar) pure nothrow
580             {
581                 T f = 1 / tan(FOVInRadians / 2);
582                 T d = 1 / (zNear - zFar);
583 
584                 return Matrix(f / aspect, 0,                  0,                    0,
585                                        0, f,                  0,                    0,
586                                        0, 0, (zFar + zNear) * d, 2 * d * zFar * zNear,
587                                        0, 0,                 -1,                    0);
588             }
589 
590             /// Returns: "lookAt" projection.
591             /// Thanks to vuaru for corrections.
592             @nogc static Matrix lookAt(vec3!T eye, vec3!T target, vec3!T up) pure nothrow
593             {
594                 vec3!T Z = (eye - target).normalized();
595                 vec3!T X = cross(-up, Z).normalized();
596                 vec3!T Y = cross(Z, -X);
597 
598                 return Matrix(-X.x,        -X.y,        -X.z,      dot(X, eye),
599                                Y.x,         Y.y,         Y.z,     -dot(Y, eye),
600                                Z.x,         Z.y,         Z.z,     -dot(Z, eye),
601                                0,           0,           0,        1);
602             }
603 
604             /// Extract frustum from a 4x4 matrice.
605             @nogc Frustum!T frustum() pure const nothrow
606             {
607                 auto left   = Plane!T(row(3) + row(0));
608                 auto right  = Plane!T(row(3) - row(0));
609                 auto top    = Plane!T(row(3) - row(1));
610                 auto bottom = Plane!T(row(3) + row(1));
611                 auto near   = Plane!T(row(3) + row(2));
612                 auto far    = Plane!T(row(3) - row(2));
613                 return Frustum!T(left, right, top, bottom, near, far);
614             }
615 
616         }
617     }
618 
619     private
620     {
621         alias T _T;
622         enum _R = R;
623         enum _C = C;
624         enum bool _isMatrix = true;
625 
626         template isAssignable(T)
627         {
628             enum bool isAssignable =
629                 is(typeof(
630                 {
631                     T x;
632                     Matrix m = x;
633                 }()));
634         }
635 
636         template isConvertible(T)
637         {
638             enum bool isConvertible = (!is(T : Matrix)) && isAssignable!T;
639         }
640 
641         template isTConvertible(U)
642         {
643             enum bool isTConvertible = is(U : T);
644         }
645 
646         template isRowConvertible(U)
647         {
648             enum bool isRowConvertible = is(U : row_t);
649         }
650 
651         template isColumnConvertible(U)
652         {
653             enum bool isColumnConvertible = is(U : column_t);
654         }
655     }
656 
657     public
658     {
659         /// Returns: an identity matrice.
660         /// Note: the identity matrix, while only meaningful for square matrices,
661         /// is also defined for non-square ones.
662         @nogc static Matrix identity() pure nothrow
663         {
664             Matrix res = void;
665             for (int i = 0; i < R; ++i)
666                 for (int j = 0; j < C; ++j)
667                     res.c[i][j] = (i == j) ? 1 : 0;
668             return res;
669         }
670 
671         /// Returns: a constant matrice.
672         @nogc static Matrix constant(U)(U x) pure nothrow
673         {
674             Matrix res = void;
675 
676             for (int i = 0; i < R * C; ++i)
677                 res.v[i] = cast(T)x;
678             return res;
679         }
680     }
681 }
682 
683 // GLSL is a big inspiration here
684 // we defines types with more or less the same names
685 template mat2x2(T) { alias Matrix!(T, 2, 2) mat2x2; }
686 template mat3x3(T) { alias Matrix!(T, 3, 3) mat3x3; }
687 template mat4x4(T) { alias Matrix!(T, 4, 4) mat4x4; }
688 
689 // WARNING: in GLSL, first number is _columns_, second is rows
690 // It is the opposite here: first number is rows, second is columns
691 // With this convention mat2x3 * mat3x4 -> mat2x4.
692 template mat2x3(T) { alias Matrix!(T, 2, 3) mat2x3; }
693 template mat2x4(T) { alias Matrix!(T, 2, 4) mat2x4; }
694 template mat3x2(T) { alias Matrix!(T, 3, 2) mat3x2; }
695 template mat3x4(T) { alias Matrix!(T, 3, 4) mat3x4; }
696 template mat4x2(T) { alias Matrix!(T, 4, 2) mat4x2; }
697 template mat4x3(T) { alias Matrix!(T, 4, 3) mat4x3; }
698 
699 alias mat2x2 mat2;
700 alias mat3x3 mat3;  // shorter names for most common matrices
701 alias mat4x4 mat4;
702 
703 private string definePostfixAliases(string type)
704 {
705     return "alias " ~ type ~ "!byte "   ~ type ~ "b;\n"
706          ~ "alias " ~ type ~ "!ubyte "  ~ type ~ "ub;\n"
707          ~ "alias " ~ type ~ "!short "  ~ type ~ "s;\n"
708          ~ "alias " ~ type ~ "!ushort " ~ type ~ "us;\n"
709          ~ "alias " ~ type ~ "!int "    ~ type ~ "i;\n"
710          ~ "alias " ~ type ~ "!uint "   ~ type ~ "ui;\n"
711          ~ "alias " ~ type ~ "!long "   ~ type ~ "l;\n"
712          ~ "alias " ~ type ~ "!ulong "  ~ type ~ "ul;\n"
713          ~ "alias " ~ type ~ "!float "  ~ type ~ "f;\n"
714          ~ "alias " ~ type ~ "!double " ~ type ~ "d;\n";
715 }
716 
717 // define a lot of type names
718 mixin(definePostfixAliases("mat2"));
719 mixin(definePostfixAliases("mat3"));
720 mixin(definePostfixAliases("mat4"));
721 mixin(definePostfixAliases("mat2x2"));
722 mixin(definePostfixAliases("mat2x3"));
723 mixin(definePostfixAliases("mat2x4"));
724 mixin(definePostfixAliases("mat3x2"));
725 mixin(definePostfixAliases("mat3x3"));
726 mixin(definePostfixAliases("mat3x4"));
727 mixin(definePostfixAliases("mat4x2"));
728 mixin(definePostfixAliases("mat4x3"));
729 mixin(definePostfixAliases("mat4x4"));
730 
731 unittest
732 {
733     mat2i x = mat2i(0, 1,
734                     2, 3);
735     assert(x.c[0][0] == 0 && x.c[0][1] == 1 && x.c[1][0] == 2 && x.c[1][1] == 3);
736 
737     vec2i[2] cols = [vec2i(0, 2), vec2i(1, 3)];
738     mat2i y = mat2i.fromColumns(cols[]);
739     assert(y.c[0][0] == 0 && y.c[0][1] == 1 && y.c[1][0] == 2 && y.c[1][1] == 3);
740     y = mat2i.fromRows(cols[]);
741     assert(y.c[0][0] == 0 && y.c[1][0] == 1 && y.c[0][1] == 2 && y.c[1][1] == 3);
742     y = y.transposed();
743 
744     assert(x == y);
745     x = [0, 1, 2, 3];
746     assert(x == y);
747 
748 
749     mat2i z = x * y;
750     assert(z == mat2i([2, 3, 6, 11]));
751     vec2i vz = z * vec2i(2, -1);
752     assert(vz == vec2i(1, 1));
753 
754     mat2f a = z;
755     mat2d ad = a;
756     ad += a;
757     mat2f w = [4, 5, 6, 7];
758     z = cast(mat2i)w;
759     assert(w == z);
760 
761     {
762         mat2x3f A;
763         mat3x4f B;
764         mat2x4f C = A * B;
765     }
766 }