1 /**
2  * This module defines a matrix data structure and various operations
3  * on this data structure. All operations are @safe pure nothrow,
4  * so they can be used in any such function, and the results of any
5  * operation can be implicitly converted to an immutable type.
6  */
7 module dstruct.matrix;
8 
9 import std.traits;
10 
11 // A private implementation of matrix multiplication for use in types.
12 private auto matrixMultiply(ResultType, T, U)
13 (ref ResultType result, ref const(T) left, ref const(U) right) {
14     foreach(row; 0 .. result.rowCount) {
15         foreach(column; 0 .. result.columnCount) {
16             Unqual!(typeof(result[0, 0])) value = left[row, 0] * right[0, column];
17 
18             foreach(pivot; 1 .. left.columnCount) {
19                 value += left[row, pivot] * right[pivot, column];
20             }
21 
22             result[row, column] = value;
23         }
24     }
25 
26     return result;
27 }
28 
29 /**
30  * A matrix type. This is a 2D array of a guaranteed uniform size.
31  */
32 struct Matrix(Number) if(isNumeric!Number) {
33 private:
34     Number[] _data;
35     size_t _rowCount;
36     size_t _columnCount;
37 public:
38     /**
39      * Create an matrix from an array of data.
40      *
41      * Params:
42      *     rowCount = The number of rows for the matrix.
43      *     columnCount = The number of columns for the matrix.
44      *     data = The array of data for the matrix.
45      */
46     @safe pure nothrow
47     this(size_t rowCount, size_t columnCount, immutable(Number[]) data) immutable {
48         _data = data;
49         _rowCount = rowCount;
50         _columnCount = columnCount;
51     }
52 
53     // Copy-paste the constructors, because inout doesn't work with literals.
54 
55     /// ditto
56     @safe pure nothrow
57     this(size_t rowCount, size_t columnCount, const(Number[]) data) const {
58         _data = data;
59         _rowCount = rowCount;
60         _columnCount = columnCount;
61     }
62 
63     /// ditto
64     @safe pure nothrow
65     this(size_t rowCount, size_t columnCount, Number[] data) {
66         _data = data;
67         _rowCount = rowCount;
68         _columnCount = columnCount;
69     }
70 
71     /**
72      * Create a matrix of a given size.
73      *
74      * Params:
75      *     rowCount = The number of rows for the matrix.
76      *     columnCount = The number of columns for the matrix.
77      */
78     @safe pure nothrow
79     this(size_t rowCount, size_t columnCount) {
80         if (rowCount == 0 || columnCount == 0) {
81             return;
82         }
83 
84         _data = new Number[](rowCount * columnCount);
85 
86         _rowCount = rowCount;
87         _columnCount = columnCount;
88     }
89 
90     /**
91      * Returns: A new duplicate of this matrix.
92      */
93     @safe pure nothrow
94     Matrix!Number dup() const {
95         Matrix!Number mat;
96 
97         // We can't .dup in a nothrow function, but we can do this...
98         mat._data = new Number[](_rowCount * _columnCount);
99         mat._data[] = _data[];
100 
101         mat._rowCount = _rowCount;
102         mat._columnCount = _columnCount;
103 
104         return mat;
105     }
106 
107     /**
108      * Returns: A new immutable duplicate of this matrix.
109      */
110     @safe pure nothrow
111     immutable(Matrix!Number) idup() const {
112         return dup();
113     }
114 
115     /**
116      * When calling .idup on an already immutable matrix, the reference
117      * to the same immutable matrix is returned. It should be safe to
118      * share the immutable memory in this manner.
119      *
120      * Returns: A reference to this immutable matrix.
121      */
122     @safe pure nothrow
123     immutable(Matrix!Number) idup() immutable {
124         // There's no need to copy immutable to immutable, share it!
125         return this;
126     }
127 
128     unittest {
129         immutable m = immutable Matrix!int(1, 1);
130 
131         // Make sure this doesn't actually duplicate.
132         assert(m._data is m._data);
133 
134         auto o = Matrix!int(1, 1);
135 
136         // Make sure this still does.
137         assert(o.idup._data !is o._data);
138     }
139 
140     /// Returns: True if the matrix is empty.
141     @safe pure nothrow
142     @property bool empty() const {
143         return _data.length == 0;
144     }
145 
146     /// Returns: The number of rows in this matrix.
147     @safe pure nothrow
148     @property size_t rowCount() const {
149         return _rowCount;
150     }
151 
152     /// Returns: The number of columns in this matrix.
153     @safe pure nothrow
154     @property size_t columnCount() const {
155         return _columnCount;
156     }
157 
158     /// Returns: true if the matrix is a square matrix.
159     @safe pure nothrow
160     @property bool isSquare() const {
161         return _rowCount == _columnCount;
162     }
163 
164     /**
165      * Slice out a row from the matrix. Modifying this
166      * slice will modify the matrix, unless it is copied.
167      *
168      * Params:
169      *    row = A row index.
170      *
171      * Returns: A slice of the row of the matrix.
172      */
173     @trusted pure nothrow
174     inout(Number[]) opIndex(size_t row) inout
175     in {
176         assert(row <= rowCount, "row out of bounds!");
177     } body {
178         size_t offset = row * _columnCount;
179 
180         return _data[offset .. offset + _columnCount];
181     }
182 
183     /**
184      * Params:
185      *    row = A row index.
186      *    column = A column index.
187      *
188      * Returns: A value from the matrix
189      */
190     @safe pure nothrow
191     ref inout(Number) opIndex(size_t row, size_t column) inout
192     in {
193         assert(column <= columnCount, "column out of bounds!");
194     } body {
195         return _data[row * _columnCount + column];
196     }
197 
198     /**
199      * Overload for foreach(rowIndex, columnIndex, value; matrix) {}
200      */
201     @trusted
202     int opApply(int delegate(ref size_t, ref size_t, ref Number) dg) {
203         int result = 0;
204 
205         matrixLoop: foreach(row; 0 .. rowCount) {
206             size_t offset = row * columnCount;
207 
208             foreach(column; 0 .. columnCount) {
209                 result = dg(row, column, _data[offset + column]);
210 
211                 if (result) {
212                     break matrixLoop;
213                 }
214             }
215         }
216 
217         return result;
218     }
219 
220 
221     /**
222      * Modify this matrix, adding/subtracting values from another matrix.
223      *
224      * Example:
225      * ---
226      *     matrix += other_matrix;
227      *     matrix -= yet_another_matrix;
228      * ---
229      *
230      * Params:
231      *     other = Another matrix with an implicitly convertible numeric type.
232      */
233     @safe pure nothrow
234     void opOpAssign(string op, OtherNumber)
235     (const Matrix!OtherNumber other)
236     if((op == "+" || op == "-") && is(OtherNumber : Number))
237     in {
238         assert(this.rowCount == other.rowCount);
239         assert(this.columnCount == other.columnCount);
240     } body {
241         foreach(i; 0.._data.length) {
242             mixin(`_data[i]` ~ op ~ `= other._data[i];`);
243         }
244     }
245 
246     /**
247      * Add or subtract two matrices, yielding a new matrix.
248      *
249      * Params:
250      *     other = The other matrix.
251      *
252      * Returns: A new matrix.
253      */
254     @safe pure nothrow
255     Matrix!Number opBinary(string op, OtherNumber)
256     (ref const Matrix!OtherNumber other) const
257     if((op == "+" || op == "-") && is(OtherNumber : Number)) in {
258         assert(this.rowCount == other.rowCount);
259         assert(this.columnCount == other.columnCount);
260     } out(val) {
261         assert(this.rowCount == val.rowCount);
262         assert(this.rowCount == val.rowCount);
263     } body {
264         // Copy this matrix.
265         auto result = this.dup;
266 
267         mixin(`result ` ~ op ~ `= other;`);
268 
269         return result;
270     }
271     
272     /// ditto
273     @safe pure nothrow
274     Matrix!Number opBinary(string op, OtherNumber)
275     (const Matrix!OtherNumber other) const
276     if((op == "+" || op == "-") && is(OtherNumber : Number)) {
277         opBinary!(op, OtherNumber)(other);
278     }
279 
280     /**
281      * Multiply two matrices.
282      *
283      * Given a matrix of size (m, n) and a matrix of size (o, p).
284      * This operation can only work if n == o.
285      * The resulting matrix will be size (m, p).
286      *
287      * Params:
288      *     other = Another matrix
289      *
290      * Returns: The product of two matrices.
291      */
292     @safe pure nothrow
293     Matrix!Number opBinary(string op, OtherNumber)
294     (ref const Matrix!OtherNumber other) const
295     if((op == "*") && is(OtherNumber : Number)) in {
296         assert(this.columnCount == other.rowCount);
297     } out(val) {
298         assert(val.rowCount == this.rowCount);
299         assert(val.columnCount == other.columnCount);
300     } body {
301         auto result = Matrix!Number(this.rowCount, other.columnCount);
302 
303         matrixMultiply(result, this, other);
304 
305         return result;
306     }
307     
308     /// ditto
309     @safe pure nothrow
310     Matrix!Number opBinary(string op, OtherNumber)
311     (const Matrix!OtherNumber other) const
312     if((op == "*") && is(OtherNumber : Number)) in {
313         opBinary!(op, OtherNumber)(other);
314     }
315 
316     /// ditto
317     @safe pure nothrow
318     Matrix!OtherNumber opBinary(string op, OtherNumber)
319     (ref const Matrix!OtherNumber other) const
320     if((op == "*") && !is(Number == OtherNumber) && is(Number : OtherNumber)) in {
321         assert(this.columnCount == other.rowCount);
322     } out(val) {
323         assert(val.rowCount == this.rowCount);
324         assert(val.columnCount == other.columnCount);
325     } body {
326         auto result = Matrix!OtherNumber(this.rowCount, other.columnCount);
327 
328         matrixMultiply(result, this, other);
329 
330         return result;
331     }
332 
333     /// ditto
334     @safe pure nothrow
335     Matrix!OtherNumber opBinary(string op, OtherNumber)
336     (const Matrix!OtherNumber other) const 
337     if((op == "*") && !is(Number == OtherNumber) && is(Number : OtherNumber)) {
338         opBinary!(op, OtherNumber)(other);
339     }
340 
341     /**
342      * Modify this matrix with a scalar value.
343      */
344     @safe pure nothrow
345     void opOpAssign(string op, OtherNumber)(OtherNumber other)
346     if(op != "in" && op != "~" && is(OtherNumber : Number)) {
347         foreach(i; 0.._data.length) {
348             mixin(`_data[i]` ~ op ~ `= other;`);
349         }
350     }
351 
352     /**
353      * Returns: A new matrix produce by combining a matrix and a scalar value.
354      */
355     @safe pure nothrow
356     Matrix!Number opBinary(string op, OtherNumber)(OtherNumber other) const
357     if(op != "in" && op != "~" && is(OtherNumber : Number)) {
358         // Copy this matrix.
359         auto result = this.dup;
360 
361         mixin(`result ` ~ op ~ `= other;`);
362 
363         return result;
364     }
365 
366     /**
367      * Returns: true if two matrices are equal and have the same type.
368      */
369     @safe pure nothrow
370     bool opEquals(ref const Matrix!Number other) const {
371          return _rowCount == other._rowCount
372             && _columnCount == other._columnCount
373             && _data == other._data;
374     }
375     
376     /// ditto
377     @safe pure nothrow
378     bool opEquals(const Matrix!Number other) const {
379         return opEquals(other);
380     }
381 }
382 
383 // Test basic matrix initialisation and foreach.
384 unittest {
385     size_t rowCount = 4;
386     size_t columnCount = 3;
387     int expectedValue = 42;
388 
389     auto mat = Matrix!int(rowCount, columnCount, [
390         42, 42, 42,
391         42, 42, 42,
392         42, 42, 42,
393         42, 42, 42
394     ]);
395 
396     size_t expectedRow = 0;
397     size_t expectedColumn = 0;
398 
399     foreach(row, column, value; mat) {
400         assert(row == expectedRow);
401         assert(column == expectedColumn);
402         assert(value == expectedValue);
403 
404         if (++expectedColumn == columnCount) {
405             expectedColumn = 0;
406             ++expectedRow;
407         }
408     }
409 }
410 
411 // Test matrix referencing and copying
412 unittest {
413     auto mat = Matrix!int(3, 3);
414 
415     mat[0, 0] = 42;
416 
417     auto normalCopy = mat.dup;
418 
419     normalCopy[0, 0] = 27;
420 
421     assert(mat[0, 0] == 42, "Matrix .dup created a data reference!");
422 
423     immutable immutCopy = mat.idup;
424 }
425 
426 // Test modifying a matrix row externally.
427 unittest {
428     auto mat = Matrix!int(3, 3);
429 
430     auto row = mat[0];
431 
432     row[0] = 3;
433     row[1] = 4;
434     row[2] = 7;
435 
436     assert(mat[0] == [3, 4, 7]);
437 }
438 
439 // Test immutable initialisation for a matrix
440 unittest {
441     immutable mat = immutable Matrix!int(3, 3, [
442         1, 2, 3,
443         4, 5, 6,
444         7, 8, 9
445     ]);
446 }
447 
448 // Test matrix addition/subtraction.
449 unittest {
450     void runtest(string op)() {
451         import std.stdio;
452 
453         size_t rowCount = 4;
454         size_t columnCount = 4;
455 
456         int leftValue = 8;
457         byte rightValue = 10;
458 
459         auto left = Matrix!int(rowCount, columnCount, [
460             8, 8, 8, 8,
461             8, 8, 8, 8,
462             8, 8, 8, 8,
463             8, 8, 8, 8,
464         ]);
465 
466         auto right = Matrix!byte(rowCount, columnCount, [
467             10, 10, 10, 10,
468             10, 10, 10, 10,
469             10, 10, 10, 10,
470             10, 10, 10, 10,
471         ]);
472 
473         auto result = mixin(`left` ~ op ~ `right`);
474         auto expectedScalar = mixin(`leftValue` ~ op ~ `rightValue`);
475 
476         foreach(row, column, value; result) {
477             assert(value == expectedScalar, `Matrix op failed: ` ~ op);
478         }
479     }
480 
481     runtest!"+";
482     runtest!"-";
483 }
484 
485 // Text matrix-scalar operations.
486 unittest {
487     import std.stdio;
488 
489     void runtest(string op)() {
490         size_t rowCount = 2;
491         size_t columnCount = 3;
492 
493         // The results for these two values are always nonzero.
494         long matrixValue = 1_234_567;
495         int scalar = 11;
496 
497         auto matrix = Matrix!long(rowCount, columnCount, [
498             matrixValue, matrixValue, matrixValue,
499             matrixValue, matrixValue, matrixValue,
500         ]);
501 
502         auto result = mixin(`matrix` ~ op ~ `scalar`);
503 
504         auto expectedScalar = mixin(`matrixValue` ~ op ~ `scalar`);
505 
506         foreach(row, column, value; result) {
507             assert(value == expectedScalar, `Matirix scalar op failed: ` ~ op);
508         }
509     }
510 
511     runtest!"+";
512     runtest!"-";
513     runtest!"*";
514     runtest!"/";
515     runtest!"%";
516     runtest!"^^";
517     runtest!"&";
518     runtest!"|";
519     runtest!"^";
520     runtest!"<<";
521     runtest!">>";
522     runtest!">>>";
523 }
524 
525 unittest {
526     // Test matrix equality.
527 
528     auto left = Matrix!int(3, 3, [
529         1, 2, 3,
530         4, 5, 6,
531         7, 8, 9
532     ]);
533 
534     auto right = immutable Matrix!int(3, 3, [
535         1, 2, 3,
536         4, 5, 6,
537         7, 8, 9
538     ]);
539 
540     assert(left == right);
541 }
542 
543 // Test matrix multiplication
544 unittest {
545     // Let's test the Wikipedia example, why not?
546     auto left = Matrix!int(2, 3, [
547         2, 3, 4,
548         1, 0, 0
549     ]);
550 
551     auto right = Matrix!int(3, 2, [
552         0, 1000,
553         1, 100,
554         0, 10
555     ]);
556 
557     auto result = left * right;
558 
559     int[][] expected = [
560         [3, 2340],
561         [0, 1000]
562     ];
563 
564 
565     foreach(i; 0..2) {
566         foreach(j; 0..2) {
567             assert(result[i, j] == expected[i][j]);
568         }
569     }
570 }
571 
572 // Test matrix multiplication, with a numeric type on the 
573 // left which implicitly converts to the right.
574 unittest {
575     auto left = Matrix!int(1, 1);
576     auto right = Matrix!long(1, 1);
577 
578     Matrix!long result = left * right;
579 }
580 
581 /**
582  * This class defines a range of rows over a matrix.
583  */
584 struct Rows(Number) if(isNumeric!Number) {
585 private:
586     const(Number)[] _data;
587     size_t _columnCount;
588 public:
589     /**
590      * Create a new rows range for a given matrix.
591      */
592     @safe pure nothrow
593     this(ref const Matrix!Number matrix) {
594         _data = matrix._data;
595         _columnCount = matrix._columnCount;
596     }
597     
598     /// ditto
599     @safe pure nothrow
600     this(const Matrix!Number matrix) {
601         this(matrix);
602     }
603 
604     /// Returns: true if the range is empty.
605     @safe pure nothrow
606     @property bool empty() const {
607         return _data.length == 0;
608     }
609 
610     /// Advance to the next row.
611     @safe pure nothrow
612     void popFront() {
613         assert(!empty, "Attempted popFront on an empty Rows range!");
614         
615         _data = _data[_columnCount .. $];
616     }
617 
618     /// Returns: The current row.
619     @safe pure nothrow
620     @property const(Number[]) front() const {
621         assert(!empty, "Cannot get the front of an empty Rows range!");
622         
623         return this[0];
624     }
625 
626     /// Save a copy of this range.
627     @safe pure nothrow
628     Rows!Number save() const {
629         return this;
630     }
631 
632     /// Retreat a row backwards.
633     @safe pure nothrow
634     void popBack() {
635         assert(!empty, "Attempted popBack on an empty Rows range!");
636         
637         _data = _data[0 .. $ - _columnCount];
638     }
639 
640     /// Returns: The row at the end of the range.
641     @safe pure nothrow
642     @property const(Number[]) back() const {
643         assert(!empty, "Cannot get the back of an empty Rows range!");
644         
645         return this[$ - 1];
646     }
647 
648     /**
649      * Params:
650      *     index = An index for a row in the range.
651      *
652      * Returns: A row at an index in the range.
653      */
654     @safe pure nothrow
655     @property const(Number[]) opIndex(size_t index) const in {
656         assert(index >= 0, "Negative index given to Rows opIndex!");
657         assert(index < length, "Out of bounds index given to Rows opIndex!");
658     } body {
659         size_t offset = index * _columnCount;
660         
661         return _data[offset .. offset + _columnCount];
662     }
663 
664     /// Returns: The current length of the range.
665     @safe pure nothrow
666     @property size_t length() const {
667         if (_data.length == 0) {
668             return 0;
669         }
670 
671         return _data.length / _columnCount;
672     }
673     
674     /// ditto
675     @safe pure nothrow
676     @property size_t opDollar() const {
677         return length;
678     }
679 }
680 
681 /**
682  * Returns: A range through a matrix's rows.
683  */
684 @safe pure nothrow
685 Rows!Number rows(Number)(Matrix!Number matrix) {
686     return typeof(return)(matrix);
687 }
688 
689 unittest {
690     auto mat = Matrix!int(3, 3, [
691         1, 2, 3,
692         4, 5, 6,
693         7, 8, 9
694     ]);
695 
696     auto expected = [
697         [1, 2, 3],
698         [4, 5, 6],
699         [7, 8, 9]
700     ];
701 
702     size_t rowIndex = 0;
703 
704     // Test InputRange stuff
705     for (auto range = mat.rows; !range.empty; range.popFront) {
706         assert(range.front == expected[rowIndex++]);
707     }
708 
709     // Test ForwardRange
710 
711     auto range1 = mat.rows;
712 
713     range1.popFront;
714     range1.popBack;
715 
716     auto range2 = range1.save;
717 
718     range1.popFront;
719 
720     assert(range2.front == [4, 5, 6]);
721 
722     rowIndex = 3;
723 
724     // Test BidirectionalRange
725     for (auto range = mat.rows; !range.empty; range.popBack) {
726         assert(range.back == expected[--rowIndex]);
727     }
728 
729     // Test RandomAccessRange
730 
731     auto range3 = mat.rows;
732 
733     range3.popFront;
734     range3.popBack;
735 
736     assert(range3.length == 1);
737     assert(range3[0] == [4, 5, 6]);
738 }
739 
740 // Test 0 size Matrix rows
741 unittest {
742     Matrix!int mat;
743 
744     assert(mat.rows.length == 0); 
745 }
746 
747 /**
748  * A static matrix type. This is a value matrix value type created directly
749  * on the stack.
750  */
751 struct Matrix(Number, size_t _rowCount, size_t _columnCount)
752 if(isNumeric!Number && _rowCount > 0 && _columnCount > 0) {
753     /// The number of rows in this matrix.
754     enum rowCount = _rowCount;
755     /// The number of columns in this matrix.
756     enum columnCount = _columnCount;
757     /// true if this matrix is a zero-sized matrix.
758     enum empty = false;
759     /// true if this matrix is a square matrix.
760     enum isSquare = rowCount == columnCount;
761 
762     /// The data backing this matrix.
763     Number[columnCount][rowCount] array2D;
764 
765     alias array2D this;
766 
767     /**
768      * Construct this matrix from a 2 dimensional static array.
769      *
770      * Params:
771      *     array2D = A 2 dimension array of the same size.
772      */
773     @safe pure nothrow
774     this(ref const(Number[columnCount][rowCount]) data) inout {
775         array2D = data;
776     }
777 
778     /// ditto
779     @safe pure nothrow
780     this(const(Number[columnCount][rowCount]) data) inout {
781         array2D = data;
782     }
783 
784     /**
785      * Construct this matrix directly from a series of numbers.
786      * This constructor is designed to be executed at compile time.
787      *
788      * Params:
789      *     numbers... = A series of numbers to initialise the matrix with.
790      */
791     @safe pure nothrow
792     this(Number[rowCount * columnCount] numbers...) {
793         foreach(row; 0 .. rowCount) {
794             foreach(column; 0 .. columnCount) {
795                 array2D[row][column] = numbers[row * columnCount + column];
796             }
797         }
798     }
799 
800     /**
801      * Returns: A reference to this matrix's data as a 1D array.
802      */
803     @trusted pure nothrow
804     @property
805     ref inout(Number[rowCount * columnCount]) array1D() inout {
806         return (cast(Number*)array2D.ptr)[0 .. rowCount * columnCount];
807     }
808 
809     // Even with alias this, we still need this overload.
810     /**
811      * Params:
812      *    row = A row index.
813      *
814      * Returns: A row from the matrix.
815      */
816     @safe pure nothrow
817     ref inout(Number[columnCount]) opIndex(size_t row) inout {
818         return array2D[row];
819     }
820 
821     /**
822      * Params:
823      *    row = A row index.
824      *    column = A column index.
825      *
826      * Returns: A value from the matrix
827      */
828     @safe pure nothrow
829     ref inout(Number) opIndex(size_t row, size_t column) inout {
830         return array2D[row][column];
831     }
832 
833     /**
834      * Overload for foreach(rowIndex, columnIndex, value; matrix) {}
835      */
836     @trusted
837     int opApply(int delegate(ref size_t, ref size_t, ref Number) dg) {
838         int result = 0;
839 
840         matrixLoop: foreach(row, rowArray; array2D) {
841             foreach(column, value; rowArray) {
842                 result = dg(row, column, value);
843 
844                 if (result) {
845                     break matrixLoop;
846                 }
847             }
848         }
849 
850         return result;
851     }
852 
853     /**
854      * Modify this matrix, adding/subtracting values from another matrix.
855      *
856      * Example:
857      * ---
858      *     matrix += other_matrix;
859      *     matrix -= yet_another_matrix;
860      * ---
861      *
862      * Params:
863      *     other = Another matrix with an implicitly convertible numeric type.
864      */
865     @safe pure nothrow
866     void opOpAssign(string op, OtherNumber)
867     (ref const(Matrix!(OtherNumber, rowCount, columnCount)) other)
868     if((op == "+" || op == "-") && is(OtherNumber : Number)) {
869         foreach(i; 0 .. rowCount) {
870             foreach(j; 0 .. columnCount) {
871                 mixin(`array2D[i][j]` ~ op ~ `= other.array2D[i][j];`);
872             }
873         }
874     }
875 
876     /// ditto
877     @safe pure nothrow
878     void opOpAssign(string op, OtherNumber)
879     (const(Matrix!(OtherNumber, rowCount, columnCount)) other)
880     if((op == "+" || op == "-") && is(OtherNumber : Number)) {
881         opOpAssign(other);
882     }
883 
884     /**
885      * Add or subtract two matrices, yielding a new matrix.
886      *
887      * Params:
888      *     other = The other matrix.
889      *
890      * Returns: A new matrix.
891      */
892     @safe pure nothrow
893     Matrix!(Number, rowCount, columnCount) opBinary(string op, OtherNumber)
894     (ref const(Matrix!(OtherNumber, rowCount, columnCount)) other)
895     if((op == "+" || op == "-") && is(OtherNumber : Number)) {
896         // Copy this matrix.
897         typeof(return) result = this;
898 
899         mixin(`result ` ~ op ~ `= other;`);
900 
901         return result;
902     }
903 
904     /// ditto
905     @safe pure nothrow
906     Matrix!(Number, rowCount, columnCount) opBinary(string op, OtherNumber)
907     (const(Matrix!(OtherNumber, rowCount, columnCount)) other)
908     if((op == "+" || op == "-") && is(OtherNumber : Number)) {
909         return opBinary(other);
910     }
911 
912     /**
913      * Multiply two matrices.
914      *
915      * Given a matrix of size (m, n) and a matrix of size (o, p).
916      * This operation can only work if n == o.
917      * The resulting matrix will be size (m, p).
918      *
919      * Params:
920      *     other = Another matrix
921      *
922      * Returns: The product of two matrices.
923      */
924     @safe pure nothrow
925     Matrix!(Number, rowCount, otherColumnCount) 
926     opBinary(string op, OtherNumber, size_t otherRowCount, size_t otherColumnCount)
927     (ref const(Matrix!(OtherNumber, otherRowCount, otherColumnCount)) other) const
928     if((op == "*") && (columnCount == otherRowCount) && is(OtherNumber : Number)) {
929         typeof(return) result;
930 
931         matrixMultiply(result, this, other);
932 
933         return result;
934     }
935 
936     /// ditto
937     @safe pure nothrow
938     Matrix!(Number, rowCount, otherColumnCount) 
939     opBinary(string op, OtherNumber, size_t otherRowCount, size_t otherColumnCount)
940     (const(Matrix!(OtherNumber, otherRowCount, otherColumnCount)) other) const
941     if((op == "*") && (columnCount == otherRowCount) && is(OtherNumber : Number)) in {
942         return opBinary!(op, OtherNumber, otherRowCount, otherColumnCount)(other);
943     }
944 
945     @safe pure nothrow
946     Matrix!(OtherNumber, rowCount, otherColumnCount) 
947     opBinary(string op, OtherNumber, size_t otherRowCount, size_t otherColumnCount)
948     (ref const(Matrix!(OtherNumber, otherRowCount, otherColumnCount)) other) const
949     if((op == "*") && (columnCount == otherRowCount) && !is(Number == OtherNumber) && is(Number : OtherNumber)) {
950         typeof(return) result;
951 
952         matrixMultiply(result, this, other);
953 
954         return result;
955     }
956 
957     @safe pure nothrow
958     Matrix!(OtherNumber, rowCount, otherColumnCount) 
959     opBinary(string op, OtherNumber, size_t otherRowCount, size_t otherColumnCount)
960     (const(Matrix!(OtherNumber, otherRowCount, otherColumnCount)) other) const
961     if((op == "*") && (columnCount == otherRowCount) && !is(Number == OtherNumber) && is(Number : OtherNumber)) {
962         return opBinary!(op, OtherNumber, otherRowCount, otherColumnCount)(other);
963     }
964 }
965 
966 // 0 size matrices are special case.
967 
968 /// ditto
969 struct Matrix(Number, size_t _rowCount, size_t _columnCount)
970 if(isNumeric!Number && _rowCount == 0 && _columnCount == 0) {
971     /// The number of rows in this matrix.
972     enum rowCount = _rowCount;
973     /// The number of columns in this matrix.
974     enum columnCount = _columnCount;
975     /// True if this matrix is a zero-sized matrix.
976     enum empty = true;
977     /// true if this matrix is a square matrix.
978     enum isSquare = true;
979 }
980 
981 /**
982  * Alias all (M, 0), (0, N) size static matrices into one single zero-sized
983  * type of size (0, 0).
984  */
985 template Matrix(Number, size_t rowCount, size_t columnCount)
986 if ((rowCount > 0 && columnCount == 0) || (rowCount == 0 && columnCount > 0)) {
987     alias Matrix = Matrix!(Number, 0, 0);
988 }
989 
990 // Test copy constructor for 2D arrays.
991 unittest {
992     int[3][2] data = [
993         [1, 2, 3],
994         [4, 5, 6],
995     ];
996 
997     Matrix!(int, 2, 3) matrix = data;
998 
999     assert(data == matrix);
1000 }
1001 
1002 // Test move constructor for 2D arrays.
1003 unittest {
1004     auto matrix = Matrix!(int, 3, 3)([
1005         [1, 2, 3],
1006         [4, 5, 6],
1007         [7, 8, 9]
1008     ]);
1009 }
1010 
1011 // Test immutable too.
1012 unittest {
1013     immutable(int[3][3]) data = [
1014         [1, 2, 3],
1015         [4, 5, 6],
1016         [7, 8, 9]
1017     ];
1018 
1019     Matrix!(int, 3, 3) matrix = data;
1020 
1021     assert(data == matrix);
1022 }
1023 
1024 unittest {
1025     int[3][3] data = [
1026         [1, 2, 3],
1027         [4, 5, 6],
1028         [7, 8, 9]
1029     ];
1030 
1031     immutable(Matrix!(int, 3, 3,)) matrix = data;
1032 
1033     assert(data == matrix);
1034 }
1035 
1036 // Test 1D array matrix slicing.
1037 unittest {
1038     auto matrix = Matrix!(int, 3, 3)([
1039         [1, 2, 3],
1040         [4, 5, 6],
1041         [7, 8, 9]
1042     ]);
1043 
1044     assert(matrix.array1D == [1, 2, 3, 4, 5, 6, 7, 8, 9]);
1045 
1046     matrix.array1D[0] = 347;
1047 
1048     assert(matrix[0][0] == 347);
1049 }
1050 
1051 // Test that copy semantics work property for 1D arrays.
1052 unittest {
1053     auto matrix = Matrix!(int, 3, 3)([
1054         [1, 2, 3],
1055         [4, 5, 6],
1056         [7, 8, 9]
1057     ]);
1058 
1059     auto arr = matrix.array1D;
1060 
1061     arr[0] = 347;
1062 
1063     // The value should not have changed.
1064     assert(matrix[0][0] == 1);
1065 }
1066 
1067 // Test that reference semantics work property for 1D arrays.
1068 unittest {
1069     auto matrix = Matrix!(int, 3, 3)([
1070         [1, 2, 3],
1071         [4, 5, 6],
1072         [7, 8, 9]
1073     ]);
1074 
1075     void foo(ref int[9] arr) {
1076         arr[0] = 347;
1077     }
1078 
1079     foo(matrix.array1D);
1080 
1081     // The should have changed.
1082     assert(matrix[0][0] == 347);
1083 }
1084 
1085 // Test const and immutable 1D array, just in case.
1086 unittest {
1087     const matrix = Matrix!(int, 3, 3)([
1088         [1, 2, 3],
1089         [4, 5, 6],
1090         [7, 8, 9]
1091     ]);
1092 
1093     assert(matrix.array1D == [1, 2, 3, 4, 5, 6, 7, 8, 9]);
1094     assert(is(typeof(matrix.array1D) == const int[9]));
1095 }
1096 
1097 unittest {
1098     immutable matrix = Matrix!(int, 3, 3)([
1099         [1, 2, 3],
1100         [4, 5, 6],
1101         [7, 8, 9]
1102     ]);
1103 
1104     assert(matrix.array1D == [1, 2, 3, 4, 5, 6, 7, 8, 9]);
1105     assert(is(typeof(matrix.array1D) == immutable int[9]));
1106 }
1107 
1108 // Test compile time init with numbers.
1109 unittest {
1110     enum matrix = Matrix!(int, 3, 3)(
1111         1, 2, 3,
1112         4, 5, 6,
1113         7, 8, 9
1114     );
1115 
1116     assert(matrix.array2D[0][0] == 1);
1117     assert(matrix.array2D[0][1] == 2);
1118     assert(matrix.array2D[0][2] == 3);
1119     assert(matrix.array2D[1][0] == 4);
1120     assert(matrix.array2D[1][1] == 5);
1121     assert(matrix.array2D[1][2] == 6);
1122     assert(matrix.array2D[2][0] == 7);
1123     assert(matrix.array2D[2][1] == 8);
1124     assert(matrix.array2D[2][2] == 9);
1125 }
1126 
1127 // Test zero sized matrices
1128 unittest {
1129     Matrix!(int, 0, 0) nothing;
1130     Matrix!(int, 1, 1) scalar;
1131 
1132     assert(nothing.empty);
1133     assert(!scalar.empty);
1134 
1135     Matrix!(int, 0, 1) noRows;
1136     Matrix!(int, 1, 0) noColumns;
1137 
1138     assert(is(typeof(nothing) == typeof(noRows)));
1139     assert(is(typeof(noRows) == typeof(noColumns)));
1140 }
1141 
1142 // Test index and assignment for cells.
1143 unittest {
1144     Matrix!(int, 3, 3) matrix;
1145 
1146     assert(matrix[0, 0] == 0);
1147 
1148     matrix[0, 0] = 3;
1149 
1150     assert(matrix[0, 0] == 3);
1151 
1152     matrix[0, 0] *= 3;
1153 
1154     assert(matrix[0, 0] == 9);
1155 }
1156 
1157 /**
1158  * Transpose (flip) a matrix.
1159  *
1160  * Params:
1161  *     matrix = The matrix to produce a transpose for.
1162  *
1163  * Returns: A new matrix which is the transpose of the given matrix.
1164  */
1165 @safe pure nothrow
1166 Matrix!Number transpose(Number)(const Matrix!Number matrix)
1167 out(val) {
1168     assert(matrix.columnCount == val.rowCount);
1169     assert(matrix.rowCount == val.columnCount);
1170 } body {
1171     auto result = typeof(return)(matrix.columnCount, matrix.rowCount);
1172 
1173     foreach(row; 0 .. matrix.rowCount) {
1174         foreach(col; 0 .. matrix.columnCount) {
1175             result[col, row] = matrix[row, col];
1176         }
1177     }
1178 
1179     return result;
1180 }
1181 
1182 /// ditto
1183 @safe pure nothrow
1184 Matrix!(Number, columnCount, rowCount)
1185 transpose(Number, size_t rowCount, size_t columnCount)
1186 (ref const Matrix!(Number, rowCount, columnCount) matrix) {
1187     typeof(return) result;
1188 
1189     foreach(row; 0 .. matrix.rowCount) {
1190         foreach(col; 0 .. matrix.columnCount) {
1191             result[col, row] = matrix[row, col];
1192         }
1193     }
1194 
1195     return result;
1196 }
1197 
1198 /// ditto
1199 @safe pure nothrow
1200 Matrix!(Number, columnCount, rowCount)
1201 transpose(Number, size_t rowCount, size_t columnCount)
1202 (const Matrix!(Number, rowCount, columnCount) matrix) {
1203     return transpose(matrix);
1204 }
1205 
1206 unittest {
1207     auto matrix = Matrix!int(2, 3, [
1208         1, 2, 3,
1209         0, -6, 7
1210     ]);
1211 
1212     int[] expected = [1, 0, 2, -6, 3, 7];
1213 
1214     auto result = matrix.transpose;
1215 
1216     assert(result.rowCount == matrix.columnCount);
1217     assert(result.columnCount == matrix.rowCount);
1218     assert(result._data == expected);
1219 }
1220 
1221 unittest {
1222     // When transposed twice, we should get the same matrix.
1223     auto matrix = Matrix!int(2, 3, [
1224         1, 2, 3,
1225         0, -6, 7
1226     ]);
1227 
1228     assert(matrix == matrix.transpose.transpose);
1229 }
1230 
1231 unittest {
1232     auto matrix = Matrix!(int, 2, 3)(
1233         1, 2, 3,
1234         0, -6, 7
1235     );
1236 
1237     assert(matrix == matrix.transpose.transpose);
1238 }
1239 
1240 // Test foreach on static matrices.
1241 unittest {
1242     auto matrix = Matrix!(int, 2, 3)(
1243         1, 2, 3,
1244         0, -6, 7
1245     );
1246 
1247     foreach(row, col, value; matrix) {
1248         if (row == 0) {
1249             if (col == 0) {
1250                 assert(value == 1);
1251             } else if (col == 1) {
1252                 assert(value == 2);
1253             } else {
1254                 assert(value == 3);
1255             }
1256         } else {
1257             if (col == 0) {
1258                 assert(value == 0);
1259             } else if (col == 1) {
1260                 assert(value == -6);
1261             } else {
1262                 assert(value == 7);
1263             }
1264         }
1265     }
1266 }
1267 
1268 // Test binary modifying operations on static matrices.
1269 unittest {
1270     auto left = Matrix!(int, 2, 3)(
1271         1, 2, 3,
1272         4, 5, 6
1273     );
1274 
1275     auto right = Matrix!(int, 2, 3)(
1276         1, 2, 3,
1277         4, 5, 6
1278     );
1279 
1280     left -= right;
1281 
1282     foreach(row, col, value; left) {
1283         assert(value == 0);
1284     }
1285 
1286     left -= right;
1287     left += right;
1288 
1289     foreach(row, col, value; left) {
1290         assert(value == 0);
1291     }
1292 }
1293 
1294 // Test binary copying operations on static matrices.
1295 unittest {
1296     auto left = Matrix!(int, 2, 3)(
1297         1, 2, 3,
1298         4, 5, 6
1299     );
1300 
1301     auto right = Matrix!(int, 2, 3)(
1302         1, 2, 3,
1303         4, 5, 6
1304     );
1305 
1306     auto newMatrix = left - right;
1307 
1308     foreach(row, col, value; newMatrix) {
1309         assert(value == 0);
1310     }
1311 
1312     auto finalMatrix = newMatrix + left;
1313 
1314     finalMatrix -= left;
1315 
1316     foreach(row, col, value; finalMatrix) {
1317         assert(value == 0);
1318     }
1319 }
1320 
1321 // Test matrix multiplication on static matrices
1322 unittest {
1323     auto left = Matrix!(int, 2, 3)(
1324         2, 3, 4,
1325         1, 0, 0
1326     );
1327 
1328     auto right = Matrix!(int, 3, 2)(
1329         0, 1000,
1330         1, 100,
1331         0, 10
1332     );
1333 
1334     auto result = left * right;
1335 
1336     int[][] expected = [
1337         [3, 2340],
1338         [0, 1000]
1339     ];
1340 
1341     foreach(i; 0..2) {
1342         foreach(j; 0..2) {
1343             assert(result[i, j] == expected[i][j]);
1344         }
1345     }
1346 }
1347 
1348 // Test matrix multiplication, with a numeric type on the 
1349 // left which implicitly converts to the right.
1350 unittest {
1351     Matrix!(int, 1, 1) left;
1352     Matrix!(long, 1, 1) right;
1353 
1354     Matrix!(long, 1, 1) result = left * right;
1355 }