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