1 /** 
2  * This module contains the `Mat` templated struct representing matrices and
3  * their operations.
4  */
5 module dvec.matrix;
6 
7 import std.traits : isNumeric, isFloatingPoint;
8 import dvec.vector;
9 import dvec.vector_types;
10 
11 /** 
12  * Generic struct that represents a matrix with `rowCount` rows and `colCount`
13  * columns, holding elements of type `T`. A matrix must be at least 1x1.
14  */
15 struct Mat(T, size_t rowCount, size_t colCount) if (isNumeric!T && rowCount > 0 && colCount > 0) {
16     /** 
17      * The width of this matrix.
18      */
19     public static const size_t WIDTH = colCount;
20 
21     /** 
22      * The height of this matrix.
23      */
24     public static const size_t HEIGHT = rowCount;
25 
26     /** 
27      * The internal static array storing the elements in row-major form.
28      */
29     public T[rowCount * colCount] data;
30 
31     /** 
32      * Constructs a matrix from the given elements.
33      * Params:
34      *   elements = The elements to place in the matrix. Should be in row-major
35      *              form.
36      */
37     public this(T[rowCount * colCount] elements) {
38         setData(elements);
39     }
40 
41     /** 
42      * Constructs a matrix from the given elements.
43      * Params:
44      *   elements = The elements to place in the matrix. Should be in row-major
45      *              form.
46      */
47     public this(T[] elements...) {
48         data[0 .. $] = elements[0 .. data.length];
49     }
50 
51     /** 
52      * Constructs a matrix from the given 2-dimensional array.
53      * Params:
54      *   elements = The elements to place in the matrix.
55      */
56     public this(T[colCount][rowCount] elements) {
57         static foreach (i; 0 .. rowCount) {
58             static foreach (j; 0 .. colCount) {
59                 this[i, j] = elements[i][j];
60             }
61         }
62     }
63 
64     /** 
65      * Constructs a matrix as a copy of another.
66      * Params:
67      *   other = The matrix to copy.
68      */
69     public this(Mat!(T, rowCount, colCount) other) {
70         setData(other.data);
71     }
72 
73     /** 
74      * Constructs a matrix with all elements initialized to the given value.
75      * Params:
76      *   value = The value to set all elements to.
77      */
78     public this(T value) {
79         static foreach (i; 0 .. data.length) data[i] = value;
80     }
81 
82     /** 
83      * Helper method to convert a 2-dimensional `[i, j]` index into a
84      * 1-dimensional `[i]` index.
85      * Params:
86      *   i = The row number.
87      *   j = The column number.
88      * Returns: A 1-dimensional index.
89      */
90     private static size_t convertToIndex(size_t i, size_t j) {
91         return colCount * i + j;
92     }
93 
94     /** 
95      * Gets the value at a specified location in the matrix. For example,
96      * with a `Mat!(float, 2, 2) m`, we can say `float v = m[0, 1];`
97      * Params:
98      *   i = The row number, starting from 0.
99      *   j = The column number, starting from 0.
100      * Returns: The value at the specified location.
101      */
102     public T opIndex(size_t i, size_t j) const {
103         return data[convertToIndex(i, j)];
104     }
105 
106     /** 
107      * Sets the value at a specified location in the matrix.
108      * Params:
109      *   value = The value to assign.
110      *   i = The row number, starting from 0.
111      *   j = The column number, starting from 0.
112      */
113     public void opIndexAssign(T value, size_t i, size_t j) {
114         data[convertToIndex(i, j)] = value;
115     }
116 
117     /** 
118      * Gets a specified row as a vector.
119      * Params:
120      *   row = The row number, starting from 0.
121      * Returns: The row.
122      */
123     public Vec!(T, colCount) getRow(size_t row) {
124         size_t idx = convertToIndex(row, 0);
125         return Vec!(T, colCount)(data[idx .. idx + colCount]);
126     }
127 
128     /** 
129      * Sets a specified row to the given vector of elements.
130      * Params:
131      *   row = The row number, starting from 0.
132      *   vector = The elements to set in the row.
133      * Returns: A reference to this matrix, for method chaining.
134      */
135     public ref Mat!(T, rowCount, colCount) setRow(size_t row, Vec!(T, colCount) vector) {
136         size_t idx = convertToIndex(row, 0);
137         data[idx .. idx + colCount] = vector.data;
138         return this;
139     }
140 
141     /** 
142      * Gets a specified column as a vector.
143      * Params:
144      *   col = The column number, starting from 0.
145      * Returns: The column.
146      */
147     public Vec!(T, rowCount) getCol(size_t col) const {
148         Vec!(T, rowCount) v;
149         static foreach (i; 0 .. rowCount) {
150             v[i] = this[i, col];
151         }
152         return v;
153     }
154 
155     /** 
156      * Sets a specified column to the given vector of elements.
157      * Params:
158      *   col = The column number, starting from 0.
159      *   vector = The elements to set in the column.
160      * Returns: A reference to this matrix, for method chaining.
161      */
162     public ref Mat!(T, rowCount, colCount) setCol(size_t col, Vec!(T, rowCount) vector) {
163         static foreach (i; 0 .. rowCount) {
164             this[i, col] = vector[i];
165         }
166         return this;
167     }
168 
169     /** 
170      * Sets all elements of this matrix using the given elements.
171      * Params:
172      *   elements = The elements to set.
173      * Returns: A reference to this matrix, for method chaining.
174      */
175     public ref Mat!(T, rowCount, colCount) setData(T[rowCount * colCount] elements) {
176         data[0 .. $] = elements[0 .. $];
177         return this;
178     }
179 
180     /** 
181      * Gets a copy of this matrix.
182      * Returns: A copy of this matrix.
183      */
184     public Mat!(T, rowCount, colCount) copy() const {
185         return Mat!(T, rowCount, colCount)(this);
186     }
187 
188     /** 
189      * Adds a given matrix to this one.
190      * Params:
191      *   other = The other matrix to add to this one.
192      * Returns: A reference to this matrix, for method chaining.
193      */
194     public ref Mat!(T, rowCount, colCount) add(Mat!(T, rowCount, colCount) other) {
195         static foreach (i; 0 .. data.length) data[i] += other.data[i];
196         return this;
197     }
198 
199     /** 
200      * Subtracts a given matrix from this one.
201      * Params:
202      *   other = The other matrix to subtract from this one.
203      * Returns: A reference to this matrix, for method chaining.
204      */
205     public ref Mat!(T, rowCount, colCount) sub(Mat!(T, rowCount, colCount) other) {
206         static foreach (i; 0 .. data.length) data[i] -= other.data[i];
207         return this;
208     }
209 
210     /** 
211      * Multiplies this matrix by a factor.
212      * Params:
213      *   factor = The factor to muliply by.
214      * Returns: A reference to this matrix, for method chaining.
215      */
216     public ref Mat!(T, rowCount, colCount) mul(T factor) {
217         static foreach (i; 0 .. data.length) data[i] *= factor;
218         return this;
219     }
220 
221     /** 
222      * Divides this matrix by a factor.
223      * Params:
224      *   factor = The factor to divide by.
225      * Returns: A reference to this matrix, for method chaining.
226      */
227     public ref Mat!(T, rowCount, colCount) div(T factor) {
228         static foreach (i; 0 .. data.length) data[i] /= factor;
229         return this;
230     }
231 
232     /** 
233      * Gets the [transpose](https://en.wikipedia.org/wiki/Transpose) of this
234      * matrix, and stores it in this matrix.
235      * Returns: A reference to this matrix, for method chaining.
236      */
237     public ref Mat!(T, rowCount, colCount) transpose() {
238         Mat!(T, colCount, rowCount) m;
239         static foreach (i; 0 .. rowCount) {
240             static foreach (j; 0 .. colCount) {
241                 m[j, i] = this[i, j];
242             }
243         }
244         return setData(m.data);
245     }
246 
247     /** 
248      * Computes the [matrix multiplication](https://en.wikipedia.org/wiki/Matrix_multiplication)
249      * of `this * other`.
250      * Params:
251      *   other = The matrix to multiply with this one.
252      * Returns: A new matrix of size `this.rowCount x other.colCount`.
253      */
254     public Mat!(T, rowCount, otherColCount) mul(T, size_t otherRowCount, size_t otherColCount)
255         (Mat!(T, otherRowCount, otherColCount) other) const {
256         Mat!(T, rowCount, otherColCount) m;
257         T sum;
258         static foreach (i; 0 .. rowCount) {
259             static foreach (j; 0 .. otherColCount) {
260                 sum = 0;
261                 static foreach (k; 0 .. colCount) {
262                     sum += this[i, k] * other[k, j];
263                 }
264                 m[i, j] = sum;
265             }
266         }
267         return m;
268     }
269 
270     /** 
271      * Multiplies a vector against this matrix.
272      * Params:
273      *   vector = The vector to multiply.
274      * Returns: The resultant transformed vector.
275      */
276     public Vec!(T, rowCount) mul(Vec!(T, colCount) vector) const {
277         Vec!(T, rowCount) result;
278         T sum;
279         static foreach (i; 0 .. rowCount) {
280             sum = 0;
281             static foreach (j; 0 .. colCount) {
282                 sum += this[i, j] * vector[j];
283             }
284             result[i] = sum;
285         }
286         return result;
287     }
288 
289     /** 
290      * Switches two rows in this matrix.
291      * Params:
292      *   rowI = A row to swap, starting from 0.
293      *   rowJ = A row to swap, starting from 0.
294      * Returns: A reference to this matrix, for method chaining.
295      */
296     public ref Mat!(T, rowCount, colCount) rowSwitch(size_t rowI, size_t rowJ) {
297         auto r = getRow(rowI);
298         setRow(rowI, getRow(rowJ));
299         setRow(rowJ, r);
300         return this;
301     }
302 
303     /** 
304      * Multiplies a row by a factor.
305      * Params:
306      *   row = The row number, starting from 0.
307      *   factor = The factor to multiply by.
308      * Returns: A reference to this matrix, for method chaining.
309      */
310     public ref Mat!(T, rowCount, colCount) rowMultiply(size_t row, T factor) {
311         size_t idx = convertToIndex(row, 0);
312         static foreach (i; 0 .. colCount) {
313             data[idx + i] *= factor;
314         }
315         return this;
316     }
317 
318     /** 
319      * Adds a row, multiplied by a factor, to another row.
320      * Params:
321      *   rowI = The row number to add to.
322      *   factor = The factor to multiply `rowJ` by.
323      *   rowJ = The row to add to `rowI`.
324      * Returns: A reference to this matrix, for method chaining.
325      */
326     public ref Mat!(T, rowCount, colCount) rowAdd(size_t rowI, T factor, size_t rowJ) {
327         auto row = getRow(rowJ);
328         row.mul(factor);
329         setRow(rowI, row);
330         return this;
331     }
332 
333     /** 
334      * Gets a submatrix of this matrix, with the given rows and columns removed.
335      * Params:
336      *   rows = The set of rows to remove.
337      *   cols = The set of columns to remove.
338      * Returns: A matrix with the given rows and columns removed.
339      */
340     public Mat!(T, rowCount - n, colCount - m) subMatrix(size_t n, size_t m)(size_t[n] rows, size_t[m] cols) const @nogc
341         if (rowCount - n > 0 && colCount - m > 0) {
342         Mat!(T, rowCount - n, colCount - m) sub;
343         size_t subIdx = 0;
344         size_t row;
345         size_t col;
346         bool skip;
347         static foreach (idx; 0 .. data.length) {
348             row = idx / colCount;
349             col = idx % colCount;
350             skip = false;
351             foreach (r; rows) {
352                 if (r == row) {
353                     skip = true;
354                     break;
355                 }
356             }
357             if (!skip) {
358                 foreach (c; cols) {
359                     if (c == col) {
360                         skip = true;
361                         break;
362                     }
363                 }
364             }
365             if (!skip) {
366                 sub.data[subIdx++] = this[row, col];
367             }
368         }
369         return sub;
370     }
371     unittest {
372         import dvec.matrix_types;
373         Mat3i m1 = Mat3i([
374             1, 2, 3,
375             4, 5, 6,
376             7, 8, 9
377         ]);
378         assert(m1.subMatrix([0], [0]) == Mat2i(5, 6, 8, 9));
379         assert(m1.subMatrix([0], [1]) == Mat2i(4, 6, 7, 9));
380         assert(m1.subMatrix([0], [2]) == Mat2i(4, 5, 7, 8));
381         assert(m1.subMatrix([1], [1]) == Mat2i(1, 3, 7, 9));
382     }
383 
384     /** 
385      * Converts this matrix to a string.
386      * Returns: A string representation of the matrix.
387      */
388     public string toString() const {
389         import std.conv : to;
390         import std.algorithm : max;
391         import std.array : appender, replicate;
392         string[data.length] values;
393         size_t[colCount] columnWidths;
394         foreach (i; 0 .. data.length) {
395             values[i] = data[i].to!string;
396             size_t colIdx = i % colCount;
397             columnWidths[colIdx] = max(columnWidths[colIdx], values[i].length + 1);
398         }
399         auto s = appender!string;
400         foreach (r; 0 .. rowCount) {
401             s ~= "| ";
402             foreach (c; 0 .. colCount) {
403                 size_t idx = convertToIndex(r, c);
404                 size_t padAmount = columnWidths[c] - values[idx].length;
405                 string padding = replicate(" ", padAmount);
406                 s ~= padding ~ values[idx];
407                 if (c < colCount - 1) s ~= ", ";
408             }
409             s ~= " |";
410             if (r < rowCount - 1) s ~= "\n";
411         }
412         return s[];
413     }
414 
415     // Special methods for square matrices.
416     static if (rowCount == colCount) {
417         alias N = rowCount;
418 
419         /** 
420          * Gets an [identity matrix](https://en.wikipedia.org/wiki/Identity_matrix).
421          * Returns: An identity matrix.
422          */
423         public static Mat!(T, rowCount, colCount) identity() {
424             Mat!(T, rowCount, colCount) m;
425             static foreach (i; 0 .. N) {
426                 static foreach (j; 0 .. N) {
427                     m[i, j] = i == j ? 1 : 0;
428                 }
429             }
430             return m;
431         }
432 
433         /** 
434          * Gets the [determinant](https://en.wikipedia.org/wiki/Determinant)
435          * of this matrix.
436          * Returns: The determinant of this matrix.
437          */
438         public T det() const {
439             static if (N == 1) {
440                 return data[0];
441             } else static if (N == 2) {
442                 return data[0] * data[3] - data[1] * data[2];
443             } else {
444                 // Laplace expansion, taking i = 0.
445                 T sum = 0;
446                 static foreach (j; 0 .. N) {
447                     sum += (j % 2 == 0 ? 1 : -1) * this[0, j] * this.subMatrix([0], [j]).det();
448                 }
449                 return sum;
450             }
451         }
452 
453         /** 
454          * Determines if this matrix is invertible, which simply means a
455          * non-zero determinant.
456          * Returns: True if this matrix is invertible.
457          */
458         public bool invertible() const {
459             return det() != 0;
460         }
461 
462         /** 
463          * Gets a [cofactor matrix](https://en.wikipedia.org/wiki/Minor_(linear_algebra)#Cofactor_expansion_of_the_determinant).
464          * Returns: A reference to this matrix, for method chaining.
465          */
466         public ref Mat!(T, rowCount, colCount) cofactor() {
467             static if (N > 1) {
468                 Mat!(T, rowCount, colCount) c;
469                 static foreach (i; 0 .. N) {
470                     static foreach (j; 0 .. N) {
471                         c[i, j] = ((i + j) % 2 == 0 ? 1 : -1) * this.subMatrix([i], [j]).det();
472                     }
473                 }
474                 setData(c.data);
475             }
476             return this;
477         }
478 
479         /** 
480          * Gets an [adjugate matrix](https://en.wikipedia.org/wiki/Adjugate_matrix).
481          * Returns: A reference to this matrix, for method chaining.
482          */
483         public ref Mat!(T, rowCount, colCount) adjugate() {
484             cofactor();
485             transpose();
486             return this;
487         }
488 
489         /** 
490          * Gets the inverse of this matrix.
491          * Returns: A reference to this matrix, for method chaining.
492          */
493         public ref Mat!(T, rowCount, colCount) inv() {
494             T d = det();
495             adjugate();
496             div(d);
497             return this;
498         }
499 
500         // Special case for 3x3 floating-point matrices: linear transformations
501         static if (N == 3 && isFloatingPoint!T) {
502             /** 
503              * Applies a 2D rotation to this matrix.
504              * Params:
505              *   theta = The angle in radians.
506              * Returns: A reference to this matrix, for method chaining.
507              */
508             public ref Mat!(T, rowCount, colCount) rotate(T theta) {
509                 import std.math : cos, sin;
510                 this = mul(Mat!(T, rowCount, colCount)(
511                     cos(theta), -sin(theta), 0,
512                     sin(theta), cos(theta),  0,
513                     0,          0,           1
514                 ));
515                 return this;
516             }
517 
518             /** 
519              * Applies a 2D translation to this matrix.
520              * Params:
521              *   dx = The translation on the x-axis.
522              *   dy = The translation on the y-axis.
523              * Returns: A reference to this matrix, for method chaining.
524              */
525             public ref Mat!(T, rowCount, colCount) translate(T dx, T dy) {
526                 this = mul(Mat!(T, rowCount, colCount)(
527                     1, 0, dx,
528                     0, 1, dy,
529                     0, 0, 1
530                 ));
531                 return this;
532             }
533 
534             /** 
535              * Applies a 2D scaling to this matrix.
536              * Params:
537              *   sx = The scale factor on the x-axis.
538              *   sy = The scale factor on the y-axis.
539              * Returns: A reference to this matrix, for method chaining.
540              */
541             public ref Mat!(T, rowCount, colCount) scale(T sx, T sy) {
542                 this = mul(Mat!(T, rowCount, colCount)(
543                     sx, 0,  0,
544                     0,  sy, 0,
545                     0,  0,  1
546                 ));
547                 return this;
548             }
549 
550             /** 
551              * Applies a uniform 2D scaling to this matrix on all axes.
552              * Params:
553              *   s = The scale factor to apply.
554              * Returns: A reference to this matrix, for method chaining.
555              */
556             public ref Mat!(T, rowCount, colCount) scale(T s) {
557                 return scale(s, s);
558             }
559 
560             /** 
561              * Applies a 2D shear to this matrix.
562              * Params:
563              *   sx = The shear factor on the x-axis.
564              *   sy = The shear factor on the y-axis.
565              * Returns: A reference to this matrix, for method chaining.
566              */
567             public ref Mat!(T, rowCount, colCount) shear(T sx, T sy) {
568                 this = mul(Mat!(T, rowCount, colCount)(
569                     1,  sx, 0,
570                     sy, 1,  0,
571                     0,  0,  1
572                 ));
573                 return this;
574             }
575 
576             /** 
577              * Maps the given 2D vector using this matrix.
578              * Params:
579              *   v = The vector to map.
580              * Returns: The resultant vector.
581              */
582             public Vec!(T, 2) map(Vec!(T, 2) v) {
583                 Vec!(T, 3) v1 = this.mul(Vec!(T, 3)(v[0], v[1], 1));
584                 return Vec!(T, 2)(v1[0], v1[1]);
585             }
586         }
587 
588         // Special case for 4x4 floating-point matrices: linear transformations
589         static if (N == 4 && isFloatingPoint!T) {
590             /** 
591              * Applies a 3D translation to this matrix.
592              * Params:
593              *   dx = The translation on the x-axis.
594              *   dy = The translation on the y-axis.
595              *   dz = The translation on the z-axis.
596              * Returns: A reference to this matrix, for method chaining.
597              */
598             public ref Mat!(T, rowCount, colCount) translate(T dx, T dy, T dz) {
599                 this = mul(Mat!(T, rowCount, colCount)(
600                     1, 0, 0, dx,
601                     0, 1, 0, dy,
602                     0, 0, 1, dz,
603                     0, 0, 0, 1
604                 ));
605                 return this;
606             }
607 
608             /** 
609              * Applies a 3D scaling to this matrix.
610              * Params:
611              *   sx = The scale factor on the x-axis.
612              *   sy = The scale factor on the y-axis.
613              *   sz = The scale factor on the z-axis.
614              * Returns: A reference to this matrix, for method chaining.
615              */
616             public ref Mat!(T, rowCount, colCount) scale(T sx, T sy, T sz) {
617                 this = mul(Mat!(T, rowCount, colCount)(
618                     sx, 0,  0,  0,
619                     0,  sy, 0,  0,
620                     0,  0,  sz, 0,
621                     0,  0,  0,  1
622                 ));
623                 return this;
624             }
625 
626             /** 
627              * Applies a uniform 3D scaling to this matrix on all axes.
628              * Params:
629              *   s = The scale factor to apply.
630              * Returns: A reference to this matrix, for method chaining.
631              */
632             public ref Mat!(T, rowCount, colCount) scale(T s) {
633                 return scale(s, s, s);
634             }
635 
636             /** 
637              * Applies a rotation to this matrix about the x-axis.
638              * Params:
639              *   theta = The angle in radians.
640              * Returns: A reference to this matrix, for method chaining.
641              */
642             public ref Mat!(T, rowCount, colCount) rotateX(T theta) {
643                 import std.math : cos, sin;
644                 this = mul(Mat!(T, rowCount, colCount)(
645                     1, 0,          0,           0,
646                     0, cos(theta), -sin(theta), 0,
647                     0, sin(theta), cos(theta),  0,
648                     0, 0,          0,           1
649                 ));
650                 return this;
651             }
652 
653             /** 
654              * Applies a rotation to this matrix about the y-axis.
655              * Params:
656              *   theta = The angle in radians.
657              * Returns: A reference to this matrix, for method chaining.
658              */
659             public ref Mat!(T, rowCount, colCount) rotateY(T theta) {
660                 import std.math : cos, sin;
661                 this = mul(Mat!(T, rowCount, colCount)(
662                     cos(theta),  0, sin(theta), 0,
663                     0,           1, 0,          0,
664                     -sin(theta), 0, cos(theta), 0,
665                     0,           0, 0,          1
666                 ));
667                 return this;
668             }
669 
670             /** 
671              * Applies a rotation to this matrix about the z-axis.
672              * Params:
673              *   theta = The angle in radians.
674              * Returns: A reference to this matrix, for method chaining.
675              */
676             public ref Mat!(T, rowCount, colCount) rotateZ(T theta) {
677                 import std.math : cos, sin;
678                 this = mul(Mat!(T, rowCount, colCount)(
679                     cos(theta), -sin(theta), 0, 0,
680                     sin(theta), cos(theta),  0, 0,
681                     0,          0,           1, 0,
682                     0,          0,           0, 1
683                 ));
684                 return this;
685             }
686 
687             /** 
688              * Applies a 3D rotation to this matrix about the x, y, and then z
689              * axes.
690              * Params:
691              *   x = The angle to rotate about the x-axis.
692              *   y = The angle to rotate about the y-axis.
693              *   z = The angle to rotate about the z-axis.
694              * Returns: A reference to this matrix, for method chaining.
695              */
696             public ref Mat!(T, rowCount, colCount) rotate(T x, T y, T z) {
697                 rotateX(x);
698                 rotateY(y);
699                 rotateZ(z);
700                 return this;
701             }
702         }
703     }
704 }
705 
706 unittest {
707     import std.stdio;
708     import dvec.vector;
709     import dvec.matrix_types;
710 
711     auto m1 = Mat3d();
712     assert(m1.data.length == 9);
713 
714     auto m2 = Mat!(double, 2, 3)([1, 2, 3, 0, -6, 7]);
715     assert(m2.transpose.data == [1, 0, 2, -6, 3, 7]);
716 
717     auto m4 = Mat2d([1, 2, 3, 4]);
718     assert(m4.getRow(0).data == [1, 2]);
719     assert(m4.getRow(1).data == [3, 4]);
720     assert(m4.getCol(0).data == [1, 3]);
721     assert(m4.getCol(1).data == [2, 4]);
722     auto m5 = m4.mul(Mat2d([0, 1, 0, 0]));
723     assert(m5.data == [0, 1, 0, 3]);
724 
725     assert(Mat2d([[1, 2], [3, 4]]).data == [1, 2, 3, 4]);
726     assert(Mat!(double, 2, 3)([[1, 2, 3], [-1, -2, -3]]).data == [1, 2, 3, -1, -2, -3]);
727 
728     assert(Mat2d.identity.data == [1, 0, 0, 1]);
729     assert(Mat3d.identity.data == [1, 0, 0, 0, 1, 0, 0, 0, 1]);
730 
731     auto m8 = Mat!(double, 2, 3)([1, -1, 2, 0, -3, 1]);
732     Vec3d v1 = Vec3d(2, 1, 0);
733     assert(m8.mul(v1).data == [1, -3]);
734 
735     Vec3f p = Vec3f(0, 0, 1);
736     Mat3f tx = Mat3f([
737         [1, 0, 42],
738         [0, 1, 64],
739         [0, 0, 1]
740     ]);
741     auto transformed = tx.mul(p);
742     assert(transformed.data == [42, 64, 1]);
743 
744     auto p2 = Vec2f(0, 0);
745     auto tx2 = Mat3f.identity();
746     tx2.translate(42, 64);
747     auto transformed2 = tx2.map(p2);
748     assert(transformed2.data == [42, 64]);
749 
750     auto m9 = Mat!(double, 3, 4)([
751         1, 2, 3, 4,
752         5, 6, 7, 8,
753         9, 10, 11, 12
754     ]);
755     auto m10 = m9.subMatrix([2], [1]);
756     assert(m10.data == [
757         1, 3, 4,
758         5, 7, 8
759     ]);
760     assert(m10.subMatrix(cast(size_t[])[], [0]).data == [
761         3, 4,
762         7, 8
763     ]);
764 
765     assert(Mat2d([3, 7, 1, -4]).det == -19);
766     assert(Mat2d([1, 2, 3, 4]).det == -2);
767     assert(Mat3d([1, 2, 3, 4, 5, 6, 7, 8, 9]).det == 0);
768 
769     assert(Mat2d.identity.inv() == Mat2d.identity());
770     assert(Mat3f.identity.inv() == Mat3f.identity());
771     assert(Mat2d(4, 7, 2, 6).inv() == Mat2d(0.6, -0.7, -0.2, 0.4));
772     assert(Mat2d(-3, 1, 5, -2).inv() == Mat2d(-2, -1, -5, -3));
773     assert(Mat3d(1, 3, 3, 1, 4, 3, 1, 3, 4).inv() == Mat3d(7, -3, -3, -1, 1, 0, -1, 0, 1));
774 }