1 /**
2 gl3n.linalg
3 
4 Special thanks to:
5 $(UL
6   $(LI Tomasz Stachowiak (h3r3tic): allowed me to use parts of $(LINK2 https://bitbucket.org/h3r3tic/boxen/src/default/src/xf/omg, omg).)
7   $(LI Jakob Øvrum (jA_cOp): improved the code a lot!)
8   $(LI Florian Boesch (___doc__): helps me to understand opengl/complex maths better, see: $(LINK http://codeflow.org/).)
9   $(LI #D on freenode: answered general questions about D.)
10 )
11 
12 Authors: David Herberth
13 License: MIT
14 
15 Note: All methods marked with pure are weakly pure since, they all access an instance member.
16 All static methods are strongly pure.
17 */
18 
19 
20 module gl3n.linalg;
21 
22 private {
23     import std.math : isNaN, isInfinity;
24     import std.conv : to;
25     import std.traits : isIntegral, isFloatingPoint, isStaticArray, isDynamicArray, isImplicitlyConvertible, isArray;
26     import std.string : format, rightJustify;
27     import std.array : join;
28     import std.algorithm : max, min, reduce;
29     import gl3n.math : clamp, PI, sqrt, sin, cos, acos, tan, asin, atan2, almost_equal;
30     import gl3n.util : is_vector, is_matrix, is_quaternion, TupleRange;
31 }
32 
33 version(NoReciprocalMul) {
34     private enum rmul = false;
35 } else {
36     private enum rmul = true;
37 }
38 
39 /// Base template for all vector-types.
40 /// Params:
41 /// type = all values get stored as this type
42 /// dimension = specifies the dimension of the vector, can be 1, 2, 3 or 4
43 /// Examples:
44 /// ---
45 /// alias Vector!(int, 3) vec3i;
46 /// alias Vector!(float, 4) vec4;
47 /// alias Vector!(real, 2) vec2r;
48 /// ---
49 struct Vector(type, int dimension_) {
50     static assert(dimension > 0, "0 dimensional vectors don't exist.");
51 
52     alias type vt; /// Holds the internal type of the vector.
53     static const int dimension = dimension_; ///Holds the dimension of the vector.
54 
55     vt[dimension] vector; /// Holds all coordinates, length conforms dimension.
56 
57     /// Returns a pointer to the coordinates.
58     @property auto value_ptr() const { return vector.ptr; }
59 
60     /// Returns the current vector formatted as string, useful for printing the vector.
61     @property string as_string() {
62         return format("%s", vector);
63     }
64     alias as_string toString; /// ditto
65 
66     @safe pure nothrow:
67     ///
68     private @property ref inout(vt) get_(char coord)() inout {
69         return vector[coord_to_index!coord];
70     }
71 
72     alias get_!'x' x; /// static properties to access the values.
73     alias x u; /// ditto
74     alias x s; /// ditto
75     alias x r; /// ditto
76     static if(dimension >= 2) {
77         alias get_!'y' y; /// ditto
78         alias y v; /// ditto
79         alias y t; /// ditto
80         alias y g; /// ditto
81     }
82     static if(dimension >= 3) {
83         alias get_!'z' z; /// ditto
84         alias z b; /// ditto
85         alias z p; /// ditto
86     }
87     static if(dimension >= 4) {
88         alias get_!'w' w; /// ditto
89         alias w a; /// ditto
90         alias w q; /// ditto
91     }
92 
93     static if(dimension == 2) {
94         enum Vector e1 = Vector(1.to!vt, 0.to!vt); /// canonical basis for Euclidian space
95         enum Vector e2 = Vector(0.to!vt, 1.to!vt); /// ditto
96     } else static if(dimension == 3) {
97         enum Vector e1 = Vector(1.to!vt, 0.to!vt, 0.to!vt); /// canonical basis for Euclidian space
98         enum Vector e2 = Vector(0.to!vt, 1.to!vt, 0.to!vt); /// ditto
99         enum Vector e3 = Vector(0.to!vt, 0.to!vt, 1.to!vt); /// ditto
100     } else static if(dimension == 4) {
101         enum Vector e1 = Vector(1.to!vt, 0.to!vt, 0.to!vt, 0.to!vt); /// canonical basis for Euclidian space
102         enum Vector e2 = Vector(0.to!vt, 1.to!vt, 0.to!vt, 0.to!vt); /// ditto
103         enum Vector e3 = Vector(0.to!vt, 0.to!vt, 1.to!vt, 0.to!vt); /// ditto
104         enum Vector e4 = Vector(0.to!vt, 0.to!vt, 0.to!vt, 1.to!vt); /// ditto
105     }
106 
107     unittest {
108         assert(vec2.e1.vector == [1.0, 0.0]);
109         assert(vec2.e2.vector == [0.0, 1.0]);
110 
111         assert(vec3.e1.vector == [1.0, 0.0, 0.0]);
112         assert(vec3.e2.vector == [0.0, 1.0, 0.0]);
113         assert(vec3.e3.vector == [0.0, 0.0, 1.0]);
114 
115         assert(vec4.e1.vector == [1.0, 0.0, 0.0, 0.0]);
116         assert(vec4.e2.vector == [0.0, 1.0, 0.0, 0.0]);
117         assert(vec4.e3.vector == [0.0, 0.0, 1.0, 0.0]);
118         assert(vec4.e4.vector == [0.0, 0.0, 0.0, 1.0]);
119     }
120 
121     static void isCompatibleVectorImpl(int d)(Vector!(vt, d) vec) if(d <= dimension) {
122     }
123 
124     template isCompatibleVector(T) {
125         enum isCompatibleVector = is(typeof(isCompatibleVectorImpl(T.init)));
126     }
127 
128     static void isCompatibleMatrixImpl(int r, int c)(Matrix!(vt, r, c) m) {
129     }
130 
131     template isCompatibleMatrix(T) {
132         enum isCompatibleMatrix = is(typeof(isCompatibleMatrixImpl(T.init)));
133     }
134 
135     private void construct(int i, T, Tail...)(T head, Tail tail) {
136         static if(i >= dimension) {
137             static assert(false, "Too many arguments passed to constructor");
138         } else static if(is(T : vt)) {
139             vector[i] = head;
140             construct!(i + 1)(tail);
141         } else static if(isDynamicArray!T) {
142             static assert((Tail.length == 0) && (i == 0), "dynamic array can not be passed together with other arguments");
143             vector[] = head[];
144         } else static if(isStaticArray!T) {
145             vector[i .. i + T.length] = head[];
146             construct!(i + T.length)(tail);
147         } else static if(isCompatibleVector!T) {
148             vector[i .. i + T.dimension] = head.vector[];
149             construct!(i + T.dimension)(tail);
150         } else {
151             static assert(false, "Vector constructor argument must be of type " ~ vt.stringof ~ " or Vector, not " ~ T.stringof);
152         }
153     }
154 
155     private void construct(int i)() { // terminate
156         static assert(i == dimension, "Not enough arguments passed to constructor");
157     }
158 
159     /// Constructs the vector.
160     /// If a single value is passed the vector, the vector will be cleared with this value.
161     /// If a vector with a higher dimension is passed the vector will hold the first values up to its dimension.
162     /// If mixed types are passed they will be joined together (allowed types: vector, static array, $(I vt)).
163     /// Examples:
164     /// ---
165     /// vec4 v4 = vec4(1.0f, vec2(2.0f, 3.0f), 4.0f);
166     /// vec3 v3 = vec3(v4); // v3 = vec3(1.0f, 2.0f, 3.0f);
167     /// vec2 v2 = v3.xy; // swizzling returns a static array.
168     /// vec3 v3_2 = vec3(1.0f); // vec3 v3_2 = vec3(1.0f, 1.0f, 1.0f);
169     /// ---
170     this(Args...)(Args args) {
171         construct!(0)(args);
172     }
173 
174     /// ditto
175     this(T)(T vec) if(is_vector!T && is(T.vt : vt) && (T.dimension >= dimension)) {
176         foreach(i; TupleRange!(0, dimension)) {
177             vector[i] = vec.vector[i];
178         }
179     }
180 
181     /// ditto
182     this()(vt value) {
183         clear(value);
184     }
185 
186     /// Returns true if all values are not nan and finite, otherwise false.
187     @property bool isFinite() const {
188         static if(isIntegral!type) {
189             return true;
190         }
191         else {
192             foreach(v; vector) {
193                 if(isNaN(v) || isInfinity(v)) {
194                     return false;
195                 }
196             }
197             return true;
198         }
199     }
200     deprecated("Use isFinite instead of ok") alias ok = isFinite;
201 
202     /// Sets all values of the vector to value.
203     void clear(vt value) {
204         foreach(i; TupleRange!(0, dimension)) {
205             vector[i] = value;
206         }
207     }
208 
209     unittest {
210         vec3 vec_clear;
211         assert(!vec_clear.isFinite);
212         vec_clear.clear(1.0f);
213         assert(vec_clear.isFinite);
214         assert(vec_clear.vector == [1.0f, 1.0f, 1.0f]);
215         assert(vec_clear.vector == vec3(1.0f).vector);
216         vec_clear.clear(float.infinity);
217         assert(!vec_clear.isFinite);
218         vec_clear.clear(float.nan);
219         assert(!vec_clear.isFinite);
220         vec_clear.clear(1.0f);
221         assert(vec_clear.isFinite);
222 
223         vec4 b = vec4(1.0f, vec_clear);
224         assert(b.isFinite);
225         assert(b.vector == [1.0f, 1.0f, 1.0f, 1.0f]);
226         assert(b.vector == vec4(1.0f).vector);
227 
228         vec2 v2_1 = vec2(vec2(0.0f, 1.0f));
229         assert(v2_1.vector == [0.0f, 1.0f]);
230 
231         vec2 v2_2 = vec2(1.0f, 1.0f);
232         assert(v2_2.vector == [1.0f, 1.0f]);
233 
234         vec3 v3 = vec3(v2_1, 2.0f);
235         assert(v3.vector == [0.0f, 1.0f, 2.0f]);
236 
237         vec4 v4_1 = vec4(1.0f, vec2(2.0f, 3.0f), 4.0f);
238         assert(v4_1.vector == [1.0f, 2.0f, 3.0f, 4.0f]);
239         assert(vec3(v4_1).vector == [1.0f, 2.0f, 3.0f]);
240         assert(vec2(vec3(v4_1)).vector == [1.0f, 2.0f]);
241         assert(vec2(vec3(v4_1)).vector == vec2(v4_1).vector);
242         assert(v4_1.vector == vec4([1.0f, 2.0f, 3.0f, 4.0f]).vector);
243 
244         vec4 v4_2 = vec4(vec2(1.0f, 2.0f), vec2(3.0f, 4.0f));
245         assert(v4_2.vector == [1.0f, 2.0f, 3.0f, 4.0f]);
246         assert(vec3(v4_2).vector == [1.0f, 2.0f, 3.0f]);
247         assert(vec2(vec3(v4_2)).vector == [1.0f, 2.0f]);
248         assert(vec2(vec3(v4_2)).vector == vec2(v4_2).vector);
249         assert(v4_2.vector == vec4([1.0f, 2.0f, 3.0f, 4.0f]).vector);
250 
251         float[2] f2 = [1.0f, 2.0f];
252         float[3] f3 = [1.0f, 2.0f, 3.0f];
253         float[4] f4 = [1.0f, 2.0f, 3.0f, 4.0f];
254         assert(vec2(1.0f, 2.0f).vector == vec2(f2).vector);
255         assert(vec3(1.0f, 2.0f, 3.0f).vector == vec3(f3).vector);
256         assert(vec3(1.0f, 2.0f, 3.0f).vector == vec3(f2, 3.0f).vector);
257         assert(vec4(1.0f, 2.0f, 3.0f, 4.0f).vector == vec4(f4).vector);
258         assert(vec4(1.0f, 2.0f, 3.0f, 4.0f).vector == vec4(f3, 4.0f).vector);
259         assert(vec4(1.0f, 2.0f, 3.0f, 4.0f).vector == vec4(f2, 3.0f, 4.0f).vector);
260         // useful for: "vec4 v4 = […]" or "vec4 v4 = other_vector.rgba"
261 
262         assert(vec3(vec3i(1, 2, 3)) == vec3(1.0, 2.0, 3.0));
263         assert(vec3d(vec3(1.0, 2.0, 3.0)) == vec3d(1.0, 2.0, 3.0));
264 
265         static assert(!__traits(compiles, vec3(0.0f, 0.0f)));
266         static assert(!__traits(compiles, vec4(0.0f, 0.0f, 0.0f)));
267         static assert(!__traits(compiles, vec4(0.0f, vec2(0.0f, 0.0f))));
268         static assert(!__traits(compiles, vec4(vec3(0.0f, 0.0f, 0.0f))));
269     }
270 
271     template coord_to_index(char c) {
272         static if((c == 'x') || (c == 'r') || (c == 'u') || (c == 's')) {
273             enum coord_to_index = 0;
274         } else static if((c == 'y') || (c == 'g') || (c == 'v') || (c == 't')) {
275             enum coord_to_index = 1;
276         } else static if((c == 'z') || (c == 'b') || (c == 'p')) {
277             static assert(dimension >= 3, "the " ~ c ~ " property is only available on vectors with a third dimension.");
278             enum coord_to_index = 2;
279         } else static if((c == 'w') || (c == 'a') || (c == 'q')) {
280             static assert(dimension >= 4, "the " ~ c ~ " property is only available on vectors with a fourth dimension.");
281             enum coord_to_index = 3;
282         } else {
283             static assert(false, "accepted coordinates are x, s, r, u, y, g, t, v, z, p, b, w, q and a not " ~ c ~ ".");
284         }
285     }
286 
287     static if(dimension == 2) { void set(vt x, vt y) { vector[0] = x; vector[1] = y; } }
288     static if(dimension == 3) { void set(vt x, vt y, vt z) { vector[0] = x; vector[1] = y; vector[2] = z; } }
289     static if(dimension == 4) { void set(vt x, vt y, vt z, vt w) { vector[0] = x; vector[1] = y; vector[2] = z; vector[3] = w; } }
290 
291     /// Updates the vector with the values from other.
292     void update(Vector!(vt, dimension) other) {
293         vector = other.vector;
294     }
295 
296     unittest {
297         vec2 v2 = vec2(1.0f, 2.0f);
298         assert(v2.x == 1.0f);
299         assert(v2.y == 2.0f);
300         v2.x = 3.0f;
301         v2.x += 1;
302         v2.x -= 1;
303         assert(v2.vector == [3.0f, 2.0f]);
304         v2.y = 4.0f;
305         v2.y += 1;
306         v2.y -= 1;
307         assert(v2.vector == [3.0f, 4.0f]);
308         assert((v2.x == 3.0f) && (v2.x == v2.u) && (v2.x == v2.s) && (v2.x == v2.r));
309         assert(v2.y == 4.0f);
310         assert((v2.y == 4.0f) && (v2.y == v2.v) && (v2.y == v2.t) && (v2.y == v2.g));
311         v2.set(0.0f, 1.0f);
312         assert(v2.vector == [0.0f, 1.0f]);
313         v2.update(vec2(3.0f, 4.0f));
314         assert(v2.vector == [3.0f, 4.0f]);
315 
316         vec3 v3 = vec3(1.0f, 2.0f, 3.0f);
317         assert(v3.x == 1.0f);
318         assert(v3.y == 2.0f);
319         assert(v3.z == 3.0f);
320         v3.x = 3.0f;
321         v3.x += 1;
322         v3.x -= 1;
323         assert(v3.vector == [3.0f, 2.0f, 3.0f]);
324         v3.y = 4.0f;
325         v3.y += 1;
326         v3.y -= 1;
327         assert(v3.vector == [3.0f, 4.0f, 3.0f]);
328         v3.z = 5.0f;
329         v3.z += 1;
330         v3.z -= 1;
331         assert(v3.vector == [3.0f, 4.0f, 5.0f]);
332         assert((v3.x == 3.0f) && (v3.x == v3.s) && (v3.x == v3.r));
333         assert((v3.y == 4.0f) && (v3.y == v3.t) && (v3.y == v3.g));
334         assert((v3.z == 5.0f) && (v3.z == v3.p) && (v3.z == v3.b));
335         v3.set(0.0f, 1.0f, 2.0f);
336         assert(v3.vector == [0.0f, 1.0f, 2.0f]);
337         v3.update(vec3(3.0f, 4.0f, 5.0f));
338         assert(v3.vector == [3.0f, 4.0f, 5.0f]);
339 
340         vec4 v4 = vec4(1.0f, 2.0f, vec2(3.0f, 4.0f));
341         assert(v4.x == 1.0f);
342         assert(v4.y == 2.0f);
343         assert(v4.z == 3.0f);
344         assert(v4.w == 4.0f);
345         v4.x = 3.0f;
346         v4.x += 1;
347         v4.x -= 1;
348         assert(v4.vector == [3.0f, 2.0f, 3.0f, 4.0f]);
349         v4.y = 4.0f;
350         v4.y += 1;
351         v4.y -= 1;
352         assert(v4.vector == [3.0f, 4.0f, 3.0f, 4.0f]);
353         v4.z = 5.0f;
354         v4.z += 1;
355         v4.z -= 1;
356         assert(v4.vector == [3.0f, 4.0f, 5.0f, 4.0f]);
357         v4.w = 6.0f;
358         v4.w += 1;
359         v4.w -= 1;
360         assert(v4.vector == [3.0f, 4.0f, 5.0f, 6.0f]);
361         assert((v4.x == 3.0f) && (v4.x == v4.s) && (v4.x == v4.r));
362         assert((v4.y == 4.0f) && (v4.y == v4.t) && (v4.y == v4.g));
363         assert((v4.z == 5.0f) && (v4.z == v4.p) && (v4.z == v4.b));
364         assert((v4.w == 6.0f) && (v4.w == v4.q) && (v4.w == v4.a));
365         v4.set(0.0f, 1.0f, 2.0f, 3.0f);
366         assert(v4.vector == [0.0f, 1.0f, 2.0f, 3.0f]);
367         v4.update(vec4(3.0f, 4.0f, 5.0f, 6.0f));
368         assert(v4.vector == [3.0f, 4.0f, 5.0f, 6.0f]);
369     }
370 
371     private void dispatchImpl(int i, string s, int size)(ref vt[size] result) const {
372         static if(s.length > 0) {
373             result[i] = vector[coord_to_index!(s[0])];
374             dispatchImpl!(i + 1, s[1..$])(result);
375         }
376     }
377 
378     /// Implements dynamic swizzling.
379     /// Returns: a Vector
380     @property Vector!(vt, s.length) opDispatch(string s)() const {
381         vt[s.length] ret;
382         dispatchImpl!(0, s)(ret);
383         Vector!(vt, s.length) ret_vec;
384         ret_vec.vector = ret;
385         return ret_vec;
386     }
387 
388     unittest {
389         vec2 v2 = vec2(1.0f, 2.0f);
390         assert(v2.xytsy == [1.0f, 2.0f, 2.0f, 1.0f, 2.0f]);
391 
392         assert(vec3(1.0f, 2.0f, 3.0f).xybzyr == [1.0f, 2.0f, 3.0f, 3.0f, 2.0f, 1.0f]);
393         assert(vec4(v2, 3.0f, 4.0f).xyzwrgbastpq == [1.0f, 2.0f, 3.0f, 4.0f,
394                                                      1.0f, 2.0f, 3.0f, 4.0f,
395                                                      1.0f, 2.0f, 3.0f, 4.0f]);
396         assert(vec4(v2, 3.0f, 4.0f).wgyzax == [4.0f, 2.0f, 2.0f, 3.0f, 4.0f, 1.0f]);
397         assert(vec4(v2.xyst).vector == [1.0f, 2.0f, 1.0f, 2.0f]);
398     }
399 
400     /// Returns the squared magnitude of the vector.
401     @property real magnitude_squared() const {
402         real temp = 0;
403 
404         foreach(index; TupleRange!(0, dimension)) {
405             temp += vector[index]^^2;
406         }
407 
408         return temp;
409     }
410 
411     /// Returns the magnitude of the vector.
412     @property real magnitude() const {
413         return sqrt(magnitude_squared);
414     }
415 
416     alias magnitude_squared length_squared; /// ditto
417     alias magnitude length; /// ditto
418 
419     /// Normalizes the vector.
420     void normalize() {
421         real len = length;
422 
423         if(len != 0) {
424             foreach(index; TupleRange!(0, dimension)) {
425                 vector[index] = cast(type)(vector[index]/len);
426             }
427         }
428     }
429 
430     /// Returns a normalized copy of the current vector.
431     @property Vector normalized() const {
432         Vector ret;
433         ret.update(this);
434         ret.normalize();
435         return ret;
436     }
437 
438     Vector opUnary(string op : "-")() const {
439         Vector ret;
440 
441         foreach(index; TupleRange!(0, dimension)) {
442             ret.vector[index] = -vector[index];
443         }
444 
445         return ret;
446     }
447 
448     unittest {
449         assert(vec2(1.0f, 1.0f) == -vec2(-1.0f, -1.0f));
450         assert(vec2(-1.0f, 1.0f) == -vec2(1.0f, -1.0f));
451 
452         assert(-vec3(1.0f, 1.0f, 1.0f) == vec3(-1.0f, -1.0f, -1.0f));
453         assert(-vec3(-1.0f, 1.0f, -1.0f) == vec3(1.0f, -1.0f, 1.0f));
454 
455         assert(vec4(1.0f, 1.0f, 1.0f, 1.0f) == -vec4(-1.0f, -1.0f, -1.0f, -1.0f));
456         assert(vec4(-1.0f, 1.0f, -1.0f, 1.0f) == -vec4(1.0f, -1.0f, 1.0f, -1.0f));
457     }
458 
459     // let the math begin!
460     Vector opBinary(string op : "*")(vt r) const {
461         Vector ret;
462 
463         foreach(index; TupleRange!(0, dimension)) {
464             ret.vector[index] = vector[index] * r;
465         }
466 
467         return ret;
468     }
469 
470     Vector opBinary(string op : "/")(vt r) const {
471         Vector ret;
472 
473         foreach(index; TupleRange!(0, dimension)) {
474             ret.vector[index] = vector[index] / r;
475         }
476 
477         return ret;
478     }
479 
480     Vector opBinary(string op)(Vector r) const if((op == "+") || (op == "-")) {
481         Vector ret;
482 
483         foreach(index; TupleRange!(0, dimension)) {
484             ret.vector[index] = mixin("vector[index]" ~ op ~ "r.vector[index]");
485         }
486 
487         return ret;
488     }
489 
490     vt opBinary(string op : "*")(Vector r) const {
491         return dot(this, r);
492     }
493 
494     // vector * matrix (for matrix * vector -> struct Matrix)
495     Vector!(vt, T.cols) opBinary(string op : "*", T)(T inp) const if(isCompatibleMatrix!T && (T.rows == dimension)) {
496         Vector!(vt, T.cols) ret;
497         ret.clear(0);
498 
499         foreach(c; TupleRange!(0, T.cols)) {
500             foreach(r; TupleRange!(0, T.rows)) {
501                 ret.vector[c] += vector[r] * inp.matrix[r][c];
502             }
503         }
504 
505         return ret;
506     }
507 
508     auto opBinaryRight(string op, T)(T inp) const if(!is_vector!T && !is_matrix!T && !is_quaternion!T) {
509         return this.opBinary!(op)(inp);
510     }
511 
512     unittest {
513         vec2 v2 = vec2(1.0f, 3.0f);
514         auto v2times2 = 2 * v2;
515         assert((v2*2.5f).vector == [2.5f, 7.5f]);
516         assert((v2+vec2(3.0f, 1.0f)).vector == [4.0f, 4.0f]);
517         assert((v2-vec2(1.0f, 3.0f)).vector == [0.0f, 0.0f]);
518         assert((v2*vec2(2.0f, 2.0f)) == 8.0f);
519 
520         vec3 v3 = vec3(1.0f, 3.0f, 5.0f);
521         assert((v3*2.5f).vector == [2.5f, 7.5f, 12.5f]);
522         assert((v3+vec3(3.0f, 1.0f, -1.0f)).vector == [4.0f, 4.0f, 4.0f]);
523         assert((v3-vec3(1.0f, 3.0f, 5.0f)).vector == [0.0f, 0.0f, 0.0f]);
524         assert((v3*vec3(2.0f, 2.0f, 2.0f)) == 18.0f);
525 
526         vec4 v4 = vec4(1.0f, 3.0f, 5.0f, 7.0f);
527         assert((v4*2.5f).vector == [2.5f, 7.5f, 12.5f, 17.5]);
528         assert((v4+vec4(3.0f, 1.0f, -1.0f, -3.0f)).vector == [4.0f, 4.0f, 4.0f, 4.0f]);
529         assert((v4-vec4(1.0f, 3.0f, 5.0f, 7.0f)).vector == [0.0f, 0.0f, 0.0f, 0.0f]);
530         assert((v4*vec4(2.0f, 2.0f, 2.0f, 2.0f)) == 32.0f);
531 
532         mat2 m2 = mat2(1.0f, 2.0f, 3.0f, 4.0f);
533         vec2 v2_2 = vec2(2.0f, 2.0f);
534         assert((v2_2*m2).vector == [8.0f, 12.0f]);
535 
536         mat3 m3 = mat3(1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f);
537         vec3 v3_2 = vec3(2.0f, 2.0f, 2.0f);
538         assert((v3_2*m3).vector == [24.0f, 30.0f, 36.0f]);
539     }
540 
541     void opOpAssign(string op : "*")(vt r) {
542         foreach(index; TupleRange!(0, dimension)) {
543             vector[index] *= r;
544         }
545     }
546 
547     void opOpAssign(string op : "/")(vt r) {
548         foreach(index; TupleRange!(0, dimension)) {
549             vector[index] /= r;
550         }
551     }
552 
553     void opOpAssign(string op)(Vector r) if((op == "+") || (op == "-")) {
554         foreach(index; TupleRange!(0, dimension)) {
555             mixin("vector[index]" ~ op ~ "= r.vector[index];");
556         }
557     }
558 
559     unittest {
560         vec2 v2 = vec2(1.0f, 3.0f);
561         v2 *= 2.5f;
562         assert(v2.vector == [2.5f, 7.5f]);
563         v2 -= vec2(2.5f, 7.5f);
564         assert(v2.vector == [0.0f, 0.0f]);
565         v2 += vec2(1.0f, 3.0f);
566         assert(v2.vector == [1.0f, 3.0f]);
567         assert(almost_equal(v2.length, sqrt(10.0f)));
568         assert(v2.length_squared == 10.0f);
569         assert((v2.magnitude == v2.length) && (v2.magnitude_squared == v2.length_squared));
570         v2 /= 2.0f;
571         assert(v2.vector == [0.5f, 1.5f]);
572         assert(almost_equal(v2.normalized, vec2(1.0f/sqrt(10.0f), 3.0f/sqrt(10.0f))));
573 
574         vec3 v3 = vec3(1.0f, 3.0f, 5.0f);
575         v3 *= 2.5f;
576         assert(v3.vector == [2.5f, 7.5f, 12.5f]);
577         v3 -= vec3(2.5f, 7.5f, 12.5f);
578         assert(v3.vector == [0.0f, 0.0f, 0.0f]);
579         v3 += vec3(1.0f, 3.0f, 5.0f);
580         assert(v3.vector == [1.0f, 3.0f, 5.0f]);
581         assert(almost_equal(v3.length, sqrt(35.0f)));
582         assert(v3.length_squared == 35.0f);
583         assert((v3.magnitude == v3.length) && (v3.magnitude_squared == v3.length_squared));
584         v3 /= 2.0f;
585         assert(v3.vector == [0.5f, 1.5f, 2.5f]);
586         assert(almost_equal(v3.normalized, vec3(1.0f/sqrt(35.0f), 3.0f/sqrt(35.0f), 5.0f/sqrt(35.0f))));
587 
588         vec4 v4 = vec4(1.0f, 3.0f, 5.0f, 7.0f);
589         v4 *= 2.5f;
590         assert(v4.vector == [2.5f, 7.5f, 12.5f, 17.5]);
591         v4 -= vec4(2.5f, 7.5f, 12.5f, 17.5f);
592         assert(v4.vector == [0.0f, 0.0f, 0.0f, 0.0f]);
593         v4 += vec4(1.0f, 3.0f, 5.0f, 7.0f);
594         assert(v4.vector == [1.0f, 3.0f, 5.0f, 7.0f]);
595         assert(almost_equal(v4.length, sqrt(84.0f)));
596         assert(v4.length_squared == 84.0f);
597         assert((v4.magnitude == v4.length) && (v4.magnitude_squared == v4.length_squared));
598         v4 /= 2.0f;
599         assert(v4.vector == [0.5f, 1.5f, 2.5f, 3.5f]);
600         assert(almost_equal(v4.normalized, vec4(1.0f/sqrt(84.0f), 3.0f/sqrt(84.0f), 5.0f/sqrt(84.0f), 7.0f/sqrt(84.0f))));
601     }
602 
603     int opCmp(ref const Vector vec) const {
604         foreach(i, a; vector) {
605             if(a < vec.vector[i]) {
606                 return -1;
607             } else if(a > vec.vector[i]) {
608                 return 1;
609             }
610         }
611 
612         // Vectors are the same
613         return 0;
614     }
615 
616     bool opEquals(T)(const T vec) const if(!isArray!T && T.dimension == dimension) {
617         return vector == vec.vector;
618     }
619 
620     bool opEquals(T)(const(T)[] array) const if(!isArray!T && !is_vector!T) {
621         if(array.length != dimension) {
622             return false;
623         }
624 
625         foreach(index; TupleRange!(0, dimension)) {
626             if(vector[index] != array[index]) {
627                 return false;
628             }
629         }
630 
631         return true;
632     }
633 
634     bool opCast(T : bool)() const {
635         return isFinite;
636     }
637 
638     unittest {
639         assert(vec2(1.0f, 2.0f) == vec2(1.0f, 2.0f));
640         assert(vec2(1.0f, 2.0f) != vec2(1.0f, 1.0f));
641         assert(vec2(1.0f, 2.0f) == vec2d(1.0, 2.0));
642         assert(vec2(1.0f, 2.0f) != vec2d(1.0, 1.0));
643         assert(vec2(1.0f, 2.0f) == vec2(1.0f, 2.0f).vector);
644         assert(vec2(1.0f, 2.0f) != vec2(1.0f, 1.0f).vector);
645         assert(vec2(1.0f, 2.0f) == vec2d(1.0, 2.0).vector);
646         assert(vec2(1.0f, 2.0f) != vec2d(1.0, 1.0).vector);
647 
648         assert(vec3(1.0f, 2.0f, 3.0f) == vec3(1.0f, 2.0f, 3.0f));
649         assert(vec3(1.0f, 2.0f, 3.0f) != vec3(1.0f, 2.0f, 2.0f));
650         assert(vec3(1.0f, 2.0f, 3.0f) == vec3d(1.0, 2.0, 3.0));
651         assert(vec3(1.0f, 2.0f, 3.0f) != vec3d(1.0, 2.0, 2.0));
652         assert(vec3(1.0f, 2.0f, 3.0f) == vec3(1.0f, 2.0f, 3.0f).vector);
653         assert(vec3(1.0f, 2.0f, 3.0f) != vec3(1.0f, 2.0f, 2.0f).vector);
654         assert(vec3(1.0f, 2.0f, 3.0f) == vec3d(1.0, 2.0, 3.0).vector);
655         assert(vec3(1.0f, 2.0f, 3.0f) != vec3d(1.0, 2.0, 2.0).vector);
656 
657         assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) == vec4(1.0f, 2.0f, 3.0f, 4.0f));
658         assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) != vec4(1.0f, 2.0f, 3.0f, 3.0f));
659         assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) == vec4d(1.0, 2.0, 3.0, 4.0));
660         assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) != vec4d(1.0, 2.0, 3.0, 3.0));
661         assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) == vec4(1.0f, 2.0f, 3.0f, 4.0f).vector);
662         assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) != vec4(1.0f, 2.0f, 3.0f, 3.0f).vector);
663         assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) == vec4d(1.0, 2.0, 3.0, 4.0).vector);
664         assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) != vec4d(1.0, 2.0, 3.0, 3.0).vector);
665 
666         assert(!(vec4(float.nan)));
667         if(vec4(1.0f)) { }
668         else { assert(false); }
669     }
670 
671 }
672 
673 /// Calculates the product between two vectors.
674 T.vt dot(T)(const T veca, const T vecb) @safe pure nothrow if(is_vector!T) {
675     T.vt temp = 0;
676 
677     foreach(index; TupleRange!(0, T.dimension)) {
678         temp += veca.vector[index] * vecb.vector[index];
679     }
680 
681     return temp;
682 }
683 
684 /// Calculates the cross product of two 3-dimensional vectors.
685 T cross(T)(const T veca, const T vecb) @safe pure nothrow if(is_vector!T && (T.dimension == 3)) {
686    return T(veca.y * vecb.z - vecb.y * veca.z,
687             veca.z * vecb.x - vecb.z * veca.x,
688             veca.x * vecb.y - vecb.x * veca.y);
689 }
690 
691 /// Calculates the distance between two vectors.
692 T.vt distance(T)(const T veca, const T vecb) @safe pure nothrow if(is_vector!T) {
693     return (veca - vecb).length;
694 }
695 
696 unittest {
697     // dot is already tested in Vector.opBinary, so no need for testing with more vectors
698     vec3 v1 = vec3(1.0f, 2.0f, -3.0f);
699     vec3 v2 = vec3(1.0f, 3.0f, 2.0f);
700 
701     assert(dot(v1, v2) == 1.0f);
702     assert(dot(v1, v2) == (v1 * v2));
703     assert(dot(v1, v2) == dot(v2, v1));
704     assert((v1 * v2) == (v1 * v2));
705 
706     assert(cross(v1, v2).vector == [13.0f, -5.0f, 1.0f]);
707     assert(cross(v2, v1).vector == [-13.0f, 5.0f, -1.0f]);
708 
709     assert(distance(vec2(0.0f, 0.0f), vec2(0.0f, 10.0f)) == 10.0);
710 }
711 
712 /// reflect a vector using a surface normal
713 T reflect(T)(const T vec, const T norm) @safe pure nothrow if(is_vector!T) {
714     return (2 * (vec * norm) * norm) - vec;
715 }
716 
717 unittest
718 {
719     assert(vec2(1,1).reflect(vec2(0,1)) == vec2(-1,1));
720     assert(vec2(-1,1).reflect(vec2(0,1)) == vec2(1,1));
721     assert(vec2(2,1).reflect(vec2(0,1)) == vec2(-2,1));
722 
723     assert(vec3(1,1,1).reflect(vec3(0,1,0)) == vec3(-1,1,-1));
724 }
725 
726 /// Pre-defined vector types, the number represents the dimension and the last letter the type (none = float, d = double, i = int).
727 alias Vector!(float, 2) vec2;
728 alias Vector!(float, 3) vec3; /// ditto
729 alias Vector!(float, 4) vec4; /// ditto
730 
731 alias Vector!(double, 2) vec2d; /// ditto
732 alias Vector!(double, 3) vec3d; /// ditto
733 alias Vector!(double, 4) vec4d; /// ditto
734 
735 alias Vector!(int, 2) vec2i; /// ditto
736 alias Vector!(int, 3) vec3i; /// ditto
737 alias Vector!(int, 4) vec4i; /// ditto
738 
739 /*alias Vector!(ubyte, 2) vec2ub;
740 alias Vector!(ubyte, 3) vec3ub;
741 alias Vector!(ubyte, 4) vec4ub;*/
742 
743 
744 /// Base template for all matrix-types.
745 /// Params:
746 ///  type = all values get stored as this type
747 ///  rows_ = rows of the matrix
748 ///  cols_ = columns of the matrix
749 /// Examples:
750 /// ---
751 /// alias Matrix!(float, 4, 4) mat4;
752 /// alias Matrix!(double, 3, 4) mat34d;
753 /// alias Matrix!(real, 2, 2) mat2r;
754 /// ---
755 struct Matrix(type, int rows_, int cols_) if((rows_ > 0) && (cols_ > 0)) {
756     alias type mt; /// Holds the internal type of the matrix;
757     static const int rows = rows_; /// Holds the number of rows;
758     static const int cols = cols_; /// Holds the number of columns;
759 
760     /// Holds the matrix $(RED row-major) in memory.
761     mt[cols][rows] matrix; // In C it would be mt[rows][cols], D does it like this: (mt[foo])[bar]
762     alias matrix this;
763 
764     unittest {
765         mat2 m2 = mat2(0.0f, 1.0f, 2.0f, 3.0f);
766         assert(m2[0][0] == 0.0f);
767         assert(m2[0][1] == 1.0f);
768         assert(m2[1][0] == 2.0f);
769         assert(m2[1][1] == 3.0f);
770         m2[0..1] = [2.0f, 2.0f];
771         assert(m2 == [[2.0f, 2.0f], [2.0f, 3.0f]]);
772 
773         mat3 m3 = mat3(0.0f, 0.1f, 0.2f, 1.0f, 1.1f, 1.2f, 2.0f, 2.1f, 2.2f);
774         assert(m3[0][1] == 0.1f);
775         assert(m3[2][0] == 2.0f);
776         assert(m3[1][2] == 1.2f);
777         m3[0][0..$] = 0.0f;
778         assert(m3 == [[0.0f, 0.0f, 0.0f],
779                       [1.0f, 1.1f, 1.2f],
780                       [2.0f, 2.1f, 2.2f]]);
781 
782         mat4 m4 = mat4(0.0f, 0.1f, 0.2f, 0.3f,
783                        1.0f, 1.1f, 1.2f, 1.3f,
784                        2.0f, 2.1f, 2.2f, 2.3f,
785                        3.0f, 3.1f, 3.2f, 3.3f);
786        assert(m4[0][3] == 0.3f);
787        assert(m4[1][1] == 1.1f);
788        assert(m4[2][0] == 2.0f);
789        assert(m4[3][2] == 3.2f);
790        m4[2][1..3] = [1.0f, 2.0f];
791        assert(m4 == [[0.0f, 0.1f, 0.2f, 0.3f],
792                      [1.0f, 1.1f, 1.2f, 1.3f],
793                      [2.0f, 1.0f, 2.0f, 2.3f],
794                      [3.0f, 3.1f, 3.2f, 3.3f]]);
795 
796     }
797 
798     /// Returns the pointer to the stored values as OpenGL requires it.
799     /// Note this will return a pointer to a $(RED row-major) matrix,
800     /// $(RED this means you've to set the transpose argument to GL_TRUE when passing it to OpenGL).
801     /// Examples:
802     /// ---
803     /// // 3rd argument = GL_TRUE
804     /// glUniformMatrix4fv(programs.main.model, 1, GL_TRUE, mat4.translation(-0.5f, -0.5f, 1.0f).value_ptr);
805     /// ---
806     @property auto value_ptr() const { return matrix[0].ptr; }
807 
808     /// Returns the current matrix formatted as flat string.
809     @property string as_string() {
810         return format("%s", matrix);
811     }
812     alias as_string toString; /// ditto
813 
814     /// Returns the current matrix as pretty formatted string.
815     @property string as_pretty_string() {
816         string fmtr = "%s";
817 
818         size_t rjust = max(format(fmtr, reduce!(max)(matrix[])).length,
819                            format(fmtr, reduce!(min)(matrix[])).length) - 1;
820 
821         string[] outer_parts;
822         foreach(mt[] row; matrix) {
823             string[] inner_parts;
824             foreach(mt col; row) {
825                 inner_parts ~= rightJustify(format(fmtr, col), rjust);
826             }
827             outer_parts ~= " [" ~ join(inner_parts, ", ") ~ "]";
828         }
829 
830         return "[" ~ join(outer_parts, "\n")[1..$] ~ "]";
831     }
832     alias as_pretty_string toPrettyString; /// ditto
833 
834     @safe pure nothrow:
835     static void isCompatibleMatrixImpl(int r, int c)(Matrix!(mt, r, c) m) {
836     }
837 
838     template isCompatibleMatrix(T) {
839         enum isCompatibleMatrix = is(typeof(isCompatibleMatrixImpl(T.init)));
840     }
841 
842     static void isCompatibleVectorImpl(int d)(Vector!(mt, d) vec) {
843     }
844 
845     template isCompatibleVector(T) {
846         enum isCompatibleVector = is(typeof(isCompatibleVectorImpl(T.init)));
847     }
848 
849     private void construct(int i, T, Tail...)(T head, Tail tail) {
850         static if(i >= rows*cols) {
851             static assert(false, "Too many arguments passed to constructor");
852         } else static if(is(T : mt)) {
853             matrix[i / cols][i % cols] = head;
854             construct!(i + 1)(tail);
855         } else static if(is(T == Vector!(mt, cols))) {
856             static if(i % cols == 0) {
857                 matrix[i / cols] = head.vector;
858                 construct!(i + T.dimension)(tail);
859             } else {
860                 static assert(false, "Can't convert Vector into the matrix. Maybe it doesn't align to the columns correctly or dimension doesn't fit");
861             }
862         } else static if(isDynamicArray!T) {
863             foreach(j; 0..cols*rows)
864                 matrix[j / cols][j % cols] = head[j];
865         } else {
866             static assert(false, "Matrix constructor argument must be of type " ~ mt.stringof ~ " or Vector, not " ~ T.stringof);
867         }
868     }
869 
870     private void construct(int i)() { // terminate
871         static assert(i == rows*cols, "Not enough arguments passed to constructor");
872     }
873 
874     /// Constructs the matrix:
875     /// If a single value is passed, the matrix will be cleared with this value (each column in each row will contain this value).
876     /// If a matrix with more rows and columns is passed, the matrix will be the upper left nxm matrix.
877     /// If a matrix with less rows and columns is passed, the passed matrix will be stored in the upper left of an identity matrix.
878     /// It's also allowed to pass vectors and scalars at a time, but the vectors dimension must match the number of columns and align correctly.
879     /// Examples:
880     /// ---
881     /// mat2 m2 = mat2(0.0f); // mat2 m2 = mat2(0.0f, 0.0f, 0.0f, 0.0f);
882     /// mat3 m3 = mat3(m2); // mat3 m3 = mat3(0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f);
883     /// mat3 m3_2 = mat3(vec3(1.0f, 2.0f, 3.0f), 4.0f, 5.0f, 6.0f, vec3(7.0f, 8.0f, 9.0f));
884     /// mat4 m4 = mat4.identity; // just an identity matrix
885     /// mat3 m3_3 = mat3(m4); // mat3 m3_3 = mat3.identity
886     /// ---
887     this(Args...)(Args args) {
888         construct!(0)(args);
889     }
890 
891     /// ditto
892     this(T)(T mat) if(is_matrix!T && (T.cols >= cols) && (T.rows >= rows)) {
893         foreach(r; TupleRange!(0, rows)) {
894             foreach(c; TupleRange!(0, cols)) {
895                 matrix[r][c] = mat.matrix[r][c];
896             }
897         }
898     }
899 
900     /// ditto
901     this(T)(T mat) if(is_matrix!T && (T.cols < cols) && (T.rows < rows)) {
902         make_identity();
903 
904         foreach(r; TupleRange!(0, T.rows)) {
905             foreach(c; TupleRange!(0, T.cols)) {
906                 matrix[r][c] = mat.matrix[r][c];
907             }
908         }
909     }
910 
911     /// ditto
912     this()(mt value) {
913         clear(value);
914     }
915 
916     /// Returns true if all values are not nan and finite, otherwise false.
917     @property bool isFinite() const {
918         static if(isIntegral!type) {
919             return true;
920         }
921         else {
922             foreach(row; matrix) {
923                 foreach(col; row) {
924                     if(isNaN(col) || isInfinity(col)) {
925                         return false;
926                     }
927                 }
928             }
929             return true;
930         }
931 
932     }
933     deprecated("Use isFinite instead of ok") alias ok = isFinite;
934 
935     /// Sets all values of the matrix to value (each column in each row will contain this value).
936     void clear(mt value) {
937         foreach(r; TupleRange!(0, rows)) {
938             foreach(c; TupleRange!(0, cols)) {
939                 matrix[r][c] = value;
940             }
941         }
942     }
943 
944     unittest {
945         mat2 m2 = mat2(1.0f, 1.0f, vec2(2.0f, 2.0f));
946         assert(m2.matrix == [[1.0f, 1.0f], [2.0f, 2.0f]]);
947         m2.clear(3.0f);
948         assert(m2.matrix == [[3.0f, 3.0f], [3.0f, 3.0f]]);
949         assert(m2.isFinite);
950         m2.clear(float.nan);
951         assert(!m2.isFinite);
952         m2.clear(float.infinity);
953         assert(!m2.isFinite);
954         m2.clear(0.0f);
955         assert(m2.isFinite);
956 
957         mat3 m3 = mat3(1.0f);
958         assert(m3.matrix == [[1.0f, 1.0f, 1.0f],
959                              [1.0f, 1.0f, 1.0f],
960                              [1.0f, 1.0f, 1.0f]]);
961 
962         mat4 m4 = mat4(vec4(1.0f, 1.0f, 1.0f, 1.0f),
963                             2.0f, 2.0f, 2.0f, 2.0f,
964                             3.0f, 3.0f, 3.0f, 3.0f,
965                        vec4(4.0f, 4.0f, 4.0f, 4.0f));
966         assert(m4.matrix == [[1.0f, 1.0f, 1.0f, 1.0f],
967                              [2.0f, 2.0f, 2.0f, 2.0f],
968                              [3.0f, 3.0f, 3.0f, 3.0f],
969                              [4.0f, 4.0f, 4.0f, 4.0f]]);
970         assert(mat3(m4).matrix == [[1.0f, 1.0f, 1.0f],
971                                    [2.0f, 2.0f, 2.0f],
972                                    [3.0f, 3.0f, 3.0f]]);
973         assert(mat2(mat3(m4)).matrix == [[1.0f, 1.0f], [2.0f, 2.0f]]);
974         assert(mat2(m4).matrix == mat2(mat3(m4)).matrix);
975         assert(mat4(mat3(m4)).matrix == [[1.0f, 1.0f, 1.0f, 0.0f],
976                                          [2.0f, 2.0f, 2.0f, 0.0f],
977                                          [3.0f, 3.0f, 3.0f, 0.0f],
978                                          [0.0f, 0.0f, 0.0f, 1.0f]]);
979 
980         Matrix!(float, 2, 3) mt1 = Matrix!(float, 2, 3)(1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f);
981         Matrix!(float, 3, 2) mt2 = Matrix!(float, 3, 2)(6.0f, -1.0f, 3.0f, 2.0f, 0.0f, -3.0f);
982 
983         assert(mt1.matrix == [[1.0f, 2.0f, 3.0f], [4.0f, 5.0f, 6.0f]]);
984         assert(mt2.matrix == [[6.0f, -1.0f], [3.0f, 2.0f], [0.0f, -3.0f]]);
985 
986         static assert(!__traits(compiles, mat2(1, 2, 1)));
987         static assert(!__traits(compiles, mat3(1, 2, 3, 1, 2, 3, 1, 2)));
988         static assert(!__traits(compiles, mat4(1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3)));
989 
990         auto m5 = mat2([0.0f,1,2,3]);
991         assert(m5.matrix == [[0.0f, 1.0f], [2.0f, 3.0f]]);
992 
993         auto m6 = Matrix!(int, 2, 3)([0,1,2,3,4,5]);
994         assert(m6.matrix == [[0, 1, 2], [3, 4, 5]]);
995     }
996 
997     static if(rows == cols) {
998         /// Makes the current matrix an identity matrix.
999         void make_identity() {
1000             clear(0);
1001             foreach(r; TupleRange!(0, rows)) {
1002                 matrix[r][r] = 1;
1003             }
1004         }
1005 
1006         /// Returns a identity matrix.
1007         static @property Matrix identity() {
1008             Matrix ret;
1009             ret.clear(0);
1010 
1011             foreach(r; TupleRange!(0, rows)) {
1012                 ret.matrix[r][r] = 1;
1013             }
1014 
1015             return ret;
1016         }
1017 
1018         /// Transposes the current matrix;
1019         void transpose() {
1020             matrix = transposed().matrix;
1021         }
1022 
1023         unittest {
1024             mat2 m2 = mat2(1.0f);
1025             m2.transpose();
1026             assert(m2.matrix == mat2(1.0f).matrix);
1027             m2.make_identity();
1028             assert(m2.matrix == [[1.0f, 0.0f],
1029                                  [0.0f, 1.0f]]);
1030             m2.transpose();
1031             assert(m2.matrix == [[1.0f, 0.0f],
1032                                  [0.0f, 1.0f]]);
1033             assert(m2.matrix == m2.identity.matrix);
1034 
1035             mat3 m3 = mat3(1.1f, 1.2f, 1.3f,
1036                            2.1f, 2.2f, 2.3f,
1037                            3.1f, 3.2f, 3.3f);
1038             m3.transpose();
1039             assert(m3.matrix == [[1.1f, 2.1f, 3.1f],
1040                                  [1.2f, 2.2f, 3.2f],
1041                                  [1.3f, 2.3f, 3.3f]]);
1042 
1043             mat4 m4 = mat4(2.0f);
1044             m4.transpose();
1045             assert(m4.matrix == mat4(2.0f).matrix);
1046             m4.make_identity();
1047             assert(m4.matrix == [[1.0f, 0.0f, 0.0f, 0.0f],
1048                                  [0.0f, 1.0f, 0.0f, 0.0f],
1049                                  [0.0f, 0.0f, 1.0f, 0.0f],
1050                                  [0.0f, 0.0f, 0.0f, 1.0f]]);
1051             assert(m4.matrix == m4.identity.matrix);
1052         }
1053 
1054     }
1055 
1056     /// Returns a transposed copy of the matrix.
1057     @property Matrix!(mt, cols, rows) transposed() const {
1058         typeof(return) ret;
1059 
1060         foreach(r; TupleRange!(0, rows)) {
1061             foreach(c; TupleRange!(0, cols)) {
1062                 ret.matrix[c][r] = matrix[r][c];
1063             }
1064         }
1065 
1066         return ret;
1067     }
1068 
1069     // transposed already tested in last unittest
1070 
1071 
1072     static if((rows == 2) && (cols == 2)) {
1073         @property mt det() const {
1074             return (matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0]);
1075         }
1076 
1077         private Matrix invert(ref Matrix mat) const {
1078             static if(isFloatingPoint!mt && rmul) {
1079                 mt d = 1 / det;
1080 
1081                 mat.matrix = [[matrix[1][1]*d, -matrix[0][1]*d],
1082                               [-matrix[1][0]*d, matrix[0][0]*d]];
1083             } else {
1084                 mt d = det;
1085 
1086                 mat.matrix = [[matrix[1][1]/d, -matrix[0][1]/d],
1087                               [-matrix[1][0]/d, matrix[0][0]/d]];
1088             }
1089 
1090             return mat;
1091         }
1092 
1093         static Matrix scaling(mt x, mt y) {
1094             Matrix ret = Matrix.identity;
1095 
1096             ret.matrix[0][0] = x;
1097             ret.matrix[1][1] = y;
1098 
1099             return ret;
1100         }
1101 
1102         Matrix scale(mt x, mt y) {
1103             this = Matrix.scaling(x, y) * this;
1104             return this;
1105         }
1106 
1107         unittest {
1108             assert(mat2.scaling(3, 3).matrix == mat2.identity.scale(3, 3).matrix);
1109             assert(mat2.scaling(3, 3).matrix == [[3.0f, 0.0f], [0.0f, 3.0f]]);
1110         }
1111 
1112     } else static if((rows == 3) && (cols == 3)) {
1113         @property mt det() const {
1114             return (matrix[0][0] * matrix[1][1] * matrix[2][2]
1115                   + matrix[0][1] * matrix[1][2] * matrix[2][0]
1116                   + matrix[0][2] * matrix[1][0] * matrix[2][1]
1117                   - matrix[0][2] * matrix[1][1] * matrix[2][0]
1118                   - matrix[0][1] * matrix[1][0] * matrix[2][2]
1119                   - matrix[0][0] * matrix[1][2] * matrix[2][1]);
1120         }
1121 
1122         private Matrix invert(ref Matrix mat) const {
1123             static if(isFloatingPoint!mt && rmul) {
1124                 mt d = 1 / det;
1125                 enum op = "*";
1126             } else {
1127                 mt d = det;
1128                 enum op = "/";
1129             }
1130 
1131             mixin(`
1132             mat.matrix = [[(matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])`~op~`d,
1133                            (matrix[0][2] * matrix[2][1] - matrix[0][1] * matrix[2][2])`~op~`d,
1134                            (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])`~op~`d],
1135                           [(matrix[1][2] * matrix[2][0] - matrix[1][0] * matrix[2][2])`~op~`d,
1136                            (matrix[0][0] * matrix[2][2] - matrix[0][2] * matrix[2][0])`~op~`d,
1137                            (matrix[0][2] * matrix[1][0] - matrix[0][0] * matrix[1][2])`~op~`d],
1138                           [(matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0])`~op~`d,
1139                            (matrix[0][1] * matrix[2][0] - matrix[0][0] * matrix[2][1])`~op~`d,
1140                            (matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0])`~op~`d]];
1141             `);
1142 
1143             return mat;
1144         }
1145     } else static if((rows == 4) && (cols == 4)) {
1146         /// Returns the determinant of the current matrix (2x2, 3x3 and 4x4 matrices).
1147         @property mt det() const {
1148             return (matrix[0][3] * matrix[1][2] * matrix[2][1] * matrix[3][0] - matrix[0][2] * matrix[1][3] * matrix[2][1] * matrix[3][0]
1149                   - matrix[0][3] * matrix[1][1] * matrix[2][2] * matrix[3][0] + matrix[0][1] * matrix[1][3] * matrix[2][2] * matrix[3][0]
1150                   + matrix[0][2] * matrix[1][1] * matrix[2][3] * matrix[3][0] - matrix[0][1] * matrix[1][2] * matrix[2][3] * matrix[3][0]
1151                   - matrix[0][3] * matrix[1][2] * matrix[2][0] * matrix[3][1] + matrix[0][2] * matrix[1][3] * matrix[2][0] * matrix[3][1]
1152                   + matrix[0][3] * matrix[1][0] * matrix[2][2] * matrix[3][1] - matrix[0][0] * matrix[1][3] * matrix[2][2] * matrix[3][1]
1153                   - matrix[0][2] * matrix[1][0] * matrix[2][3] * matrix[3][1] + matrix[0][0] * matrix[1][2] * matrix[2][3] * matrix[3][1]
1154                   + matrix[0][3] * matrix[1][1] * matrix[2][0] * matrix[3][2] - matrix[0][1] * matrix[1][3] * matrix[2][0] * matrix[3][2]
1155                   - matrix[0][3] * matrix[1][0] * matrix[2][1] * matrix[3][2] + matrix[0][0] * matrix[1][3] * matrix[2][1] * matrix[3][2]
1156                   + matrix[0][1] * matrix[1][0] * matrix[2][3] * matrix[3][2] - matrix[0][0] * matrix[1][1] * matrix[2][3] * matrix[3][2]
1157                   - matrix[0][2] * matrix[1][1] * matrix[2][0] * matrix[3][3] + matrix[0][1] * matrix[1][2] * matrix[2][0] * matrix[3][3]
1158                   + matrix[0][2] * matrix[1][0] * matrix[2][1] * matrix[3][3] - matrix[0][0] * matrix[1][2] * matrix[2][1] * matrix[3][3]
1159                   - matrix[0][1] * matrix[1][0] * matrix[2][2] * matrix[3][3] + matrix[0][0] * matrix[1][1] * matrix[2][2] * matrix[3][3]);
1160         }
1161 
1162         private Matrix invert(ref Matrix mat) const {
1163             static if(isFloatingPoint!mt && rmul) {
1164                 mt d = 1 / det;
1165                 enum op = "*";
1166             } else {
1167                 mt d = det;
1168                 enum op = "/";
1169             }
1170 
1171             mixin(`
1172             mat.matrix = [[(matrix[1][1] * matrix[2][2] * matrix[3][3] + matrix[1][2] * matrix[2][3] * matrix[3][1] + matrix[1][3] * matrix[2][1] * matrix[3][2]
1173                           - matrix[1][1] * matrix[2][3] * matrix[3][2] - matrix[1][2] * matrix[2][1] * matrix[3][3] - matrix[1][3] * matrix[2][2] * matrix[3][1])`~op~`d,
1174                            (matrix[0][1] * matrix[2][3] * matrix[3][2] + matrix[0][2] * matrix[2][1] * matrix[3][3] + matrix[0][3] * matrix[2][2] * matrix[3][1]
1175                           - matrix[0][1] * matrix[2][2] * matrix[3][3] - matrix[0][2] * matrix[2][3] * matrix[3][1] - matrix[0][3] * matrix[2][1] * matrix[3][2])`~op~`d,
1176                            (matrix[0][1] * matrix[1][2] * matrix[3][3] + matrix[0][2] * matrix[1][3] * matrix[3][1] + matrix[0][3] * matrix[1][1] * matrix[3][2]
1177                           - matrix[0][1] * matrix[1][3] * matrix[3][2] - matrix[0][2] * matrix[1][1] * matrix[3][3] - matrix[0][3] * matrix[1][2] * matrix[3][1])`~op~`d,
1178                            (matrix[0][1] * matrix[1][3] * matrix[2][2] + matrix[0][2] * matrix[1][1] * matrix[2][3] + matrix[0][3] * matrix[1][2] * matrix[2][1]
1179                           - matrix[0][1] * matrix[1][2] * matrix[2][3] - matrix[0][2] * matrix[1][3] * matrix[2][1] - matrix[0][3] * matrix[1][1] * matrix[2][2])`~op~`d],
1180                           [(matrix[1][0] * matrix[2][3] * matrix[3][2] + matrix[1][2] * matrix[2][0] * matrix[3][3] + matrix[1][3] * matrix[2][2] * matrix[3][0]
1181                           - matrix[1][0] * matrix[2][2] * matrix[3][3] - matrix[1][2] * matrix[2][3] * matrix[3][0] - matrix[1][3] * matrix[2][0] * matrix[3][2])`~op~`d,
1182                            (matrix[0][0] * matrix[2][2] * matrix[3][3] + matrix[0][2] * matrix[2][3] * matrix[3][0] + matrix[0][3] * matrix[2][0] * matrix[3][2]
1183                           - matrix[0][0] * matrix[2][3] * matrix[3][2] - matrix[0][2] * matrix[2][0] * matrix[3][3] - matrix[0][3] * matrix[2][2] * matrix[3][0])`~op~`d,
1184                            (matrix[0][0] * matrix[1][3] * matrix[3][2] + matrix[0][2] * matrix[1][0] * matrix[3][3] + matrix[0][3] * matrix[1][2] * matrix[3][0]
1185                           - matrix[0][0] * matrix[1][2] * matrix[3][3] - matrix[0][2] * matrix[1][3] * matrix[3][0] - matrix[0][3] * matrix[1][0] * matrix[3][2])`~op~`d,
1186                            (matrix[0][0] * matrix[1][2] * matrix[2][3] + matrix[0][2] * matrix[1][3] * matrix[2][0] + matrix[0][3] * matrix[1][0] * matrix[2][2]
1187                           - matrix[0][0] * matrix[1][3] * matrix[2][2] - matrix[0][2] * matrix[1][0] * matrix[2][3] - matrix[0][3] * matrix[1][2] * matrix[2][0])`~op~`d],
1188                           [(matrix[1][0] * matrix[2][1] * matrix[3][3] + matrix[1][1] * matrix[2][3] * matrix[3][0] + matrix[1][3] * matrix[2][0] * matrix[3][1]
1189                           - matrix[1][0] * matrix[2][3] * matrix[3][1] - matrix[1][1] * matrix[2][0] * matrix[3][3] - matrix[1][3] * matrix[2][1] * matrix[3][0])`~op~`d,
1190                            (matrix[0][0] * matrix[2][3] * matrix[3][1] + matrix[0][1] * matrix[2][0] * matrix[3][3] + matrix[0][3] * matrix[2][1] * matrix[3][0]
1191                           - matrix[0][0] * matrix[2][1] * matrix[3][3] - matrix[0][1] * matrix[2][3] * matrix[3][0] - matrix[0][3] * matrix[2][0] * matrix[3][1])`~op~`d,
1192                            (matrix[0][0] * matrix[1][1] * matrix[3][3] + matrix[0][1] * matrix[1][3] * matrix[3][0] + matrix[0][3] * matrix[1][0] * matrix[3][1]
1193                           - matrix[0][0] * matrix[1][3] * matrix[3][1] - matrix[0][1] * matrix[1][0] * matrix[3][3] - matrix[0][3] * matrix[1][1] * matrix[3][0])`~op~`d,
1194                            (matrix[0][0] * matrix[1][3] * matrix[2][1] + matrix[0][1] * matrix[1][0] * matrix[2][3] + matrix[0][3] * matrix[1][1] * matrix[2][0]
1195                           - matrix[0][0] * matrix[1][1] * matrix[2][3] - matrix[0][1] * matrix[1][3] * matrix[2][0] - matrix[0][3] * matrix[1][0] * matrix[2][1])`~op~`d],
1196                           [(matrix[1][0] * matrix[2][2] * matrix[3][1] + matrix[1][1] * matrix[2][0] * matrix[3][2] + matrix[1][2] * matrix[2][1] * matrix[3][0]
1197                           - matrix[1][0] * matrix[2][1] * matrix[3][2] - matrix[1][1] * matrix[2][2] * matrix[3][0] - matrix[1][2] * matrix[2][0] * matrix[3][1])`~op~`d,
1198                            (matrix[0][0] * matrix[2][1] * matrix[3][2] + matrix[0][1] * matrix[2][2] * matrix[3][0] + matrix[0][2] * matrix[2][0] * matrix[3][1]
1199                           - matrix[0][0] * matrix[2][2] * matrix[3][1] - matrix[0][1] * matrix[2][0] * matrix[3][2] - matrix[0][2] * matrix[2][1] * matrix[3][0])`~op~`d,
1200                            (matrix[0][0] * matrix[1][2] * matrix[3][1] + matrix[0][1] * matrix[1][0] * matrix[3][2] + matrix[0][2] * matrix[1][1] * matrix[3][0]
1201                           - matrix[0][0] * matrix[1][1] * matrix[3][2] - matrix[0][1] * matrix[1][2] * matrix[3][0] - matrix[0][2] * matrix[1][0] * matrix[3][1])`~op~`d,
1202                            (matrix[0][0] * matrix[1][1] * matrix[2][2] + matrix[0][1] * matrix[1][2] * matrix[2][0] + matrix[0][2] * matrix[1][0] * matrix[2][1]
1203                           - matrix[0][0] * matrix[1][2] * matrix[2][1] - matrix[0][1] * matrix[1][0] * matrix[2][2] - matrix[0][2] * matrix[1][1] * matrix[2][0])`~op~`d]];
1204             `);
1205 
1206             return mat;
1207         }
1208 
1209         // some static fun ...
1210         // (1) glprogramming.com/red/appendixf.html - ortographic is broken!
1211         // (2) http://fly.cc.fer.hr/~unreal/theredbook/appendixg.html
1212         // (3) http://en.wikipedia.org/wiki/Orthographic_projection_(geometry)
1213 
1214         static if(isFloatingPoint!mt) {
1215             static private mt[6] cperspective(mt width, mt height, mt fov, mt near, mt far)
1216                 in { assert(height != 0); }
1217                 body {
1218                     mt aspect = width/height;
1219                     mt top = near * tan(fov*(PI/360.0));
1220                     mt bottom = -top;
1221                     mt right = top * aspect;
1222                     mt left = -right;
1223 
1224                     return [left, right, bottom, top, near, far];
1225                 }
1226 
1227             /// Returns a perspective matrix (4x4 and floating-point matrices only).
1228             static Matrix perspective(mt width, mt height, mt fov, mt near, mt far) {
1229                 mt[6] cdata = cperspective(width, height, fov, near, far);
1230                 return perspective(cdata[0], cdata[1], cdata[2], cdata[3], cdata[4], cdata[5]);
1231             }
1232 
1233             /// ditto
1234             static Matrix perspective(mt left, mt right, mt bottom, mt top, mt near, mt far)
1235                 in {
1236                     assert(right-left != 0);
1237                     assert(top-bottom != 0);
1238                     assert(far-near != 0);
1239                 }
1240                 body {
1241                     Matrix ret;
1242                     ret.clear(0);
1243 
1244                     ret.matrix[0][0] = (2*near)/(right-left);
1245                     ret.matrix[0][2] = (right+left)/(right-left);
1246                     ret.matrix[1][1] = (2*near)/(top-bottom);
1247                     ret.matrix[1][2] = (top+bottom)/(top-bottom);
1248                     ret.matrix[2][2] = -(far+near)/(far-near);
1249                     ret.matrix[2][3] = -(2*far*near)/(far-near);
1250                     ret.matrix[3][2] = -1;
1251 
1252                     return ret;
1253                 }
1254 
1255             /// Returns an inverse perspective matrix (4x4 and floating-point matrices only).
1256             static Matrix perspective_inverse(mt width, mt height, mt fov, mt near, mt far) {
1257                 mt[6] cdata = cperspective(width, height, fov, near, far);
1258                 return perspective_inverse(cdata[0], cdata[1], cdata[2], cdata[3], cdata[4], cdata[5]);
1259             }
1260 
1261             /// ditto
1262             static Matrix perspective_inverse(mt left, mt right, mt bottom, mt top, mt near, mt far)
1263                 in {
1264                     assert(near != 0);
1265                     assert(far != 0);
1266                 }
1267                 body {
1268                     Matrix ret;
1269                     ret.clear(0);
1270 
1271                     ret.matrix[0][0] = (right-left)/(2*near);
1272                     ret.matrix[0][3] = (right+left)/(2*near);
1273                     ret.matrix[1][1] = (top-bottom)/(2*near);
1274                     ret.matrix[1][3] = (top+bottom)/(2*near);
1275                     ret.matrix[2][3] = -1;
1276                     ret.matrix[3][2] = -(far-near)/(2*far*near);
1277                     ret.matrix[3][3] = (far+near)/(2*far*near);
1278 
1279                     return ret;
1280                 }
1281 
1282             // (2) and (3) say this one is correct
1283             /// Returns an orthographic matrix (4x4 and floating-point matrices only).
1284             static Matrix orthographic(mt left, mt right, mt bottom, mt top, mt near, mt far)
1285                 in {
1286                     assert(right-left != 0);
1287                     assert(top-bottom != 0);
1288                     assert(far-near != 0);
1289                 }
1290                 body {
1291                     Matrix ret;
1292                     ret.clear(0);
1293 
1294                     ret.matrix[0][0] = 2/(right-left);
1295                     ret.matrix[0][3] = -(right+left)/(right-left);
1296                     ret.matrix[1][1] = 2/(top-bottom);
1297                     ret.matrix[1][3] = -(top+bottom)/(top-bottom);
1298                     ret.matrix[2][2] = -2/(far-near);
1299                     ret.matrix[2][3] = -(far+near)/(far-near);
1300                     ret.matrix[3][3] = 1;
1301 
1302                     return ret;
1303                 }
1304 
1305             // (1) and (2) say this one is correct
1306             /// Returns an inverse ortographic matrix (4x4 and floating-point matrices only).
1307             static Matrix orthographic_inverse(mt left, mt right, mt bottom, mt top, mt near, mt far) {
1308                 Matrix ret;
1309                 ret.clear(0);
1310 
1311                 ret.matrix[0][0] = (right-left)/2;
1312                 ret.matrix[0][3] = (right+left)/2;
1313                 ret.matrix[1][1] = (top-bottom)/2;
1314                 ret.matrix[1][3] = (top+bottom)/2;
1315                 ret.matrix[2][2] = (far-near)/-2;
1316                 ret.matrix[2][3] = (far+near)/2;
1317                 ret.matrix[3][3] = 1;
1318 
1319                 return ret;
1320             }
1321 
1322             /// Returns a look at matrix (4x4 and floating-point matrices only).
1323             static Matrix look_at(Vector!(mt, 3) eye, Vector!(mt, 3) target, Vector!(mt, 3) up) {
1324                 alias Vector!(mt, 3) vec3mt;
1325                 vec3mt look_dir = (target - eye).normalized;
1326                 vec3mt up_dir = up.normalized;
1327 
1328                 vec3mt right_dir = cross(look_dir, up_dir).normalized;
1329                 vec3mt perp_up_dir = cross(right_dir, look_dir);
1330 
1331                 Matrix ret = Matrix.identity;
1332                 ret.matrix[0][0..3] = right_dir.vector[];
1333                 ret.matrix[1][0..3] = perp_up_dir.vector[];
1334                 ret.matrix[2][0..3] = (-look_dir).vector[];
1335 
1336                 ret.matrix[0][3] = -dot(eye, right_dir);
1337                 ret.matrix[1][3] = -dot(eye, perp_up_dir);
1338                 ret.matrix[2][3] = dot(eye, look_dir);
1339 
1340                 return ret;
1341             }
1342 
1343             unittest {
1344                 mt[6] cp = cperspective(600f, 900f, 60f, 1f, 100f);
1345                 assert(cp[4] == 1.0f);
1346                 assert(cp[5] == 100.0f);
1347                 assert(cp[0] == -cp[1]);
1348                 assert((cp[0] < -0.38489f) && (cp[0] > -0.38491f));
1349                 assert(cp[2] == -cp[3]);
1350                 assert((cp[2] < -0.577349f) && (cp[2] > -0.577351f));
1351 
1352                 assert(mat4.perspective(600f, 900f, 60.0, 1.0, 100.0) == mat4.perspective(cp[0], cp[1], cp[2], cp[3], cp[4], cp[5]));
1353                 float[4][4] m4p = mat4.perspective(600f, 900f, 60.0, 1.0, 100.0).matrix;
1354                 assert((m4p[0][0] < 2.598077f) && (m4p[0][0] > 2.598075f));
1355                 assert(m4p[0][2] == 0.0f);
1356                 assert((m4p[1][1] < 1.732052) && (m4p[1][1] > 1.732050));
1357                 assert(m4p[1][2] == 0.0f);
1358                 assert((m4p[2][2] < -1.020201) && (m4p[2][2] > -1.020203));
1359                 assert((m4p[2][3] < -2.020201) && (m4p[2][3] > -2.020203));
1360                 assert((m4p[3][2] < -0.9f) && (m4p[3][2] > -1.1f));
1361 
1362                 float[4][4] m4pi = mat4.perspective_inverse(600f, 900f, 60.0, 1.0, 100.0).matrix;
1363                 assert((m4pi[0][0] < 0.384901) && (m4pi[0][0] > 0.384899));
1364                 assert(m4pi[0][3] == 0.0f);
1365                 assert((m4pi[1][1] < 0.577351) && (m4pi[1][1] > 0.577349));
1366                 assert(m4pi[1][3] == 0.0f);
1367                 assert(m4pi[2][3] == -1.0f);
1368                 assert((m4pi[3][2] < -0.494999) && (m4pi[3][2] > -0.495001));
1369                 assert((m4pi[3][3] < 0.505001) && (m4pi[3][3] > 0.504999));
1370 
1371                 // maybe the next tests should be improved
1372                 float[4][4] m4o = mat4.orthographic(-1.0f, 1.0f, -1.0f, 1.0f, -1.0f, 1.0f).matrix;
1373                 assert(m4o == [[1.0f, 0.0f, 0.0f, 0.0f],
1374                                [0.0f, 1.0f, 0.0f, 0.0f],
1375                                [0.0f, 0.0f, -1.0f, 0.0f],
1376                                [0.0f, 0.0f, 0.0f, 1.0f]]);
1377 
1378                 float[4][4] m4oi = mat4.orthographic_inverse(-1.0f, 1.0f, -1.0f, 1.0f, -1.0f, 1.0f).matrix;
1379                 assert(m4oi == [[1.0f, 0.0f, 0.0f, 0.0f],
1380                                 [0.0f, 1.0f, 0.0f, 0.0f],
1381                                 [0.0f, 0.0f, -1.0f, 0.0f],
1382                                 [0.0f, 0.0f, 0.0f, 1.0f]]);
1383 
1384                 //TODO: look_at tests
1385             }
1386         }
1387     }
1388 
1389     static if((rows == cols) && (rows >= 3) && (rows <= 4)) {
1390         /// Returns a translation matrix (3x3 and 4x4 matrices).
1391         static Matrix translation(mt x, mt y, mt z) {
1392             Matrix ret = Matrix.identity;
1393 
1394             ret.matrix[0][cols-1] = x;
1395             ret.matrix[1][cols-1] = y;
1396             ret.matrix[2][cols-1] = z;
1397 
1398             return ret;
1399         }
1400 
1401         /// ditto
1402         static Matrix translation(Vector!(mt, 3) v) {
1403             Matrix ret = Matrix.identity;
1404             
1405             ret.matrix[0][cols-1] = v.x;
1406             ret.matrix[1][cols-1] = v.y;
1407             ret.matrix[2][cols-1] = v.z;
1408             
1409             return ret;
1410         }
1411 
1412         /// Applys a translation on the current matrix and returns $(I this) (3x3 and 4x4 matrices).
1413         Matrix translate(mt x, mt y, mt z) {
1414             this = Matrix.translation(x, y, z) * this;
1415             return this;
1416         }
1417 
1418         /// ditto
1419         Matrix translate(Vector!(mt, 3) v) {
1420             this = Matrix.translation(v) * this;
1421             return this;
1422         }
1423 
1424         /// Returns a scaling matrix (3x3 and 4x4 matrices);
1425         static Matrix scaling(mt x, mt y, mt z) {
1426             Matrix ret = Matrix.identity;
1427 
1428             ret.matrix[0][0] = x;
1429             ret.matrix[1][1] = y;
1430             ret.matrix[2][2] = z;
1431 
1432             return ret;
1433         }
1434 
1435         /// Applys a scale to the current matrix and returns $(I this) (3x3 and 4x4 matrices).
1436         Matrix scale(mt x, mt y, mt z) {
1437             this = Matrix.scaling(x, y, z) * this;
1438             return this;
1439         }
1440 
1441         unittest {
1442             mat3 m3 = mat3.identity;
1443             assert(m3.translate(1.0f, 2.0f, 3.0f).matrix == mat3.translation(1.0f, 2.0f, 3.0f).matrix);
1444             assert(mat3.translation(1.0f, 2.0f, 3.0f).matrix == [[1.0f, 0.0f, 1.0f],
1445                                                                  [0.0f, 1.0f, 2.0f],
1446                                                                  [0.0f, 0.0f, 3.0f]]);
1447             assert(mat3.identity.translate(0.0f, 1.0f, 2.0f).matrix == mat3.translation(0.0f, 1.0f, 2.0f).matrix);
1448 
1449             mat3 m31 = mat3.identity;
1450             assert(m31.translate(vec3(1.0f, 2.0f, 3.0f)).matrix == mat3.translation(vec3(1.0f, 2.0f, 3.0f)).matrix);
1451             assert(mat3.translation(vec3(1.0f, 2.0f, 3.0f)).matrix == [[1.0f, 0.0f, 1.0f],
1452                     [0.0f, 1.0f, 2.0f],
1453                     [0.0f, 0.0f, 3.0f]]);
1454             assert(mat3.identity.translate(vec3(0.0f, 1.0f, 2.0f)).matrix == mat3.translation(vec3(0.0f, 1.0f, 2.0f)).matrix);
1455 
1456             assert(m3.scaling(0.0f, 1.0f, 2.0f).matrix == mat3.scaling(0.0f, 1.0f, 2.0f).matrix);
1457             assert(mat3.scaling(0.0f, 1.0f, 2.0f).matrix == [[0.0f, 0.0f, 0.0f],
1458                                                              [0.0f, 1.0f, 0.0f],
1459                                                              [0.0f, 0.0f, 2.0f]]);
1460             assert(mat3.identity.scale(0.0f, 1.0f, 2.0f).matrix == mat3.scaling(0.0f, 1.0f, 2.0f).matrix);
1461 
1462             // same tests for 4x4
1463 
1464             mat4 m4 = mat4(1.0f);
1465             assert(m4.translation(1.0f, 2.0f, 3.0f).matrix == mat4.translation(1.0f, 2.0f, 3.0f).matrix);
1466             assert(mat4.translation(1.0f, 2.0f, 3.0f).matrix == [[1.0f, 0.0f, 0.0f, 1.0f],
1467                                                                  [0.0f, 1.0f, 0.0f, 2.0f],
1468                                                                  [0.0f, 0.0f, 1.0f, 3.0f],
1469                                                                  [0.0f, 0.0f, 0.0f, 1.0f]]);
1470             assert(mat4.identity.translate(0.0f, 1.0f, 2.0f).matrix == mat4.translation(0.0f, 1.0f, 2.0f).matrix);
1471 
1472             assert(m4.scaling(0.0f, 1.0f, 2.0f).matrix == mat4.scaling(0.0f, 1.0f, 2.0f).matrix);
1473             assert(mat4.scaling(0.0f, 1.0f, 2.0f).matrix == [[0.0f, 0.0f, 0.0f, 0.0f],
1474                                                              [0.0f, 1.0f, 0.0f, 0.0f],
1475                                                              [0.0f, 0.0f, 2.0f, 0.0f],
1476                                                              [0.0f, 0.0f, 0.0f, 1.0f]]);
1477             assert(mat4.identity.scale(0.0f, 1.0f, 2.0f).matrix == mat4.scaling(0.0f, 1.0f, 2.0f).matrix);
1478         }
1479     }
1480 
1481 
1482     static if((rows == cols) && (rows >= 3)) {
1483         static if(isFloatingPoint!mt) {
1484             /// Returns an identity matrix with an applied rotate_axis around an arbitrary axis (nxn matrices, n >= 3).
1485             static Matrix rotation(real alpha, Vector!(mt, 3) axis) {
1486                 Matrix mult = Matrix.identity;
1487 
1488                 if(axis.length != 1) {
1489                     axis.normalize();
1490                 }
1491 
1492                 real cosa = cos(alpha);
1493                 real sina = sin(alpha);
1494 
1495                 Vector!(mt, 3) temp = (1 - cosa)*axis;
1496 
1497                 mult.matrix[0][0] = to!mt(cosa + temp.x * axis.x);
1498                 mult.matrix[0][1] = to!mt(       temp.x * axis.y + sina * axis.z);
1499                 mult.matrix[0][2] = to!mt(       temp.x * axis.z - sina * axis.y);
1500                 mult.matrix[1][0] = to!mt(       temp.y * axis.x - sina * axis.z);
1501                 mult.matrix[1][1] = to!mt(cosa + temp.y * axis.y);
1502                 mult.matrix[1][2] = to!mt(       temp.y * axis.z + sina * axis.x);
1503                 mult.matrix[2][0] = to!mt(       temp.z * axis.x + sina * axis.y);
1504                 mult.matrix[2][1] = to!mt(       temp.z * axis.y - sina * axis.x);
1505                 mult.matrix[2][2] = to!mt(cosa + temp.z * axis.z);
1506 
1507                 return mult;
1508             }
1509 
1510             /// ditto
1511             static Matrix rotation(real alpha, mt x, mt y, mt z) {
1512                 return Matrix.rotation(alpha, Vector!(mt, 3)(x, y, z));
1513             }
1514 
1515             /// Returns an identity matrix with an applied rotation around the x-axis (nxn matrices, n >= 3).
1516             static Matrix xrotation(real alpha) {
1517                 Matrix mult = Matrix.identity;
1518 
1519                 mt cosamt = to!mt(cos(alpha));
1520                 mt sinamt = to!mt(sin(alpha));
1521 
1522                 mult.matrix[1][1] = cosamt;
1523                 mult.matrix[1][2] = -sinamt;
1524                 mult.matrix[2][1] = sinamt;
1525                 mult.matrix[2][2] = cosamt;
1526 
1527                 return mult;
1528             }
1529 
1530             /// Returns an identity matrix with an applied rotation around the y-axis (nxn matrices, n >= 3).
1531             static Matrix yrotation(real alpha) {
1532                 Matrix mult = Matrix.identity;
1533 
1534                 mt cosamt = to!mt(cos(alpha));
1535                 mt sinamt = to!mt(sin(alpha));
1536 
1537                 mult.matrix[0][0] = cosamt;
1538                 mult.matrix[0][2] = sinamt;
1539                 mult.matrix[2][0] = -sinamt;
1540                 mult.matrix[2][2] = cosamt;
1541 
1542                 return mult;
1543             }
1544 
1545             /// Returns an identity matrix with an applied rotation around the z-axis (nxn matrices, n >= 3).
1546             static Matrix zrotation(real alpha) {
1547                 Matrix mult = Matrix.identity;
1548 
1549                 mt cosamt = to!mt(cos(alpha));
1550                 mt sinamt = to!mt(sin(alpha));
1551 
1552                 mult.matrix[0][0] = cosamt;
1553                 mult.matrix[0][1] = -sinamt;
1554                 mult.matrix[1][0] = sinamt;
1555                 mult.matrix[1][1] = cosamt;
1556 
1557                 return mult;
1558             }
1559 
1560             Matrix rotate(real alpha, Vector!(mt, 3) axis) {
1561                 this = rotation(alpha, axis) * this;
1562                 return this;
1563             }
1564 
1565             /// Rotates the current matrix around the x-axis and returns $(I this) (nxn matrices, n >= 3).
1566             Matrix rotatex(real alpha) {
1567                 this = xrotation(alpha) * this;
1568                 return this;
1569             }
1570 
1571             /// Rotates the current matrix around the y-axis and returns $(I this) (nxn matrices, n >= 3).
1572             Matrix rotatey(real alpha) {
1573                 this = yrotation(alpha) * this;
1574                 return this;
1575             }
1576 
1577             /// Rotates the current matrix around the z-axis and returns $(I this) (nxn matrices, n >= 3).
1578             Matrix rotatez(real alpha) {
1579                 this = zrotation(alpha) * this;
1580                 return this;
1581             }
1582 
1583             unittest {
1584                 assert(mat4.xrotation(0).matrix == [[1.0f, 0.0f, 0.0f, 0.0f],
1585                                                     [0.0f, 1.0f, -0.0f, 0.0f],
1586                                                     [0.0f, 0.0f, 1.0f, 0.0f],
1587                                                     [0.0f, 0.0f, 0.0f, 1.0f]]);
1588                 assert(mat4.yrotation(0).matrix == [[1.0f, 0.0f, 0.0f, 0.0f],
1589                                                     [0.0f, 1.0f, 0.0f, 0.0f],
1590                                                     [0.0f, 0.0f, 1.0f, 0.0f],
1591                                                     [0.0f, 0.0f, 0.0f, 1.0f]]);
1592                 assert(mat4.zrotation(0).matrix == [[1.0f, -0.0f, 0.0f, 0.0f],
1593                                                     [0.0f, 1.0f, 0.0f, 0.0f],
1594                                                     [0.0f, 0.0f, 1.0f, 0.0f],
1595                                                     [0.0f, 0.0f, 0.0f, 1.0f]]);
1596                 mat4 xro = mat4.identity;
1597                 xro.rotatex(0);
1598                 assert(mat4.xrotation(0).matrix == xro.matrix);
1599                 assert(xro.matrix == mat4.identity.rotatex(0).matrix);
1600                 assert(xro.matrix == mat4.rotation(0, vec3(1.0f, 0.0f, 0.0f)).matrix);
1601                 mat4 yro = mat4.identity;
1602                 yro.rotatey(0);
1603                 assert(mat4.yrotation(0).matrix == yro.matrix);
1604                 assert(yro.matrix == mat4.identity.rotatey(0).matrix);
1605                 assert(yro.matrix == mat4.rotation(0, vec3(0.0f, 1.0f, 0.0f)).matrix);
1606                 mat4 zro = mat4.identity;
1607                 xro.rotatez(0);
1608                 assert(mat4.zrotation(0).matrix == zro.matrix);
1609                 assert(zro.matrix == mat4.identity.rotatez(0).matrix);
1610                 assert(zro.matrix == mat4.rotation(0, vec3(0.0f, 0.0f, 1.0f)).matrix);
1611             }
1612         } // isFloatingPoint
1613 
1614 
1615         /// Sets the translation of the matrix (nxn matrices, n >= 3).
1616         void set_translation(mt[] values...) // intended to be a property
1617             in { assert(values.length >= (rows-1)); }
1618             body {
1619                 foreach(r; TupleRange!(0, rows-1)) {
1620                     matrix[r][rows-1] = values[r];
1621                 }
1622             }
1623 
1624         /// Copyies the translation from mat to the current matrix (nxn matrices, n >= 3).
1625         void set_translation(Matrix mat) {
1626             foreach(r; TupleRange!(0, rows-1)) {
1627                 matrix[r][rows-1] = mat.matrix[r][rows-1];
1628             }
1629         }
1630 
1631         /// Returns an identity matrix with the current translation applied (nxn matrices, n >= 3)..
1632         Matrix get_translation() {
1633             Matrix ret = Matrix.identity;
1634 
1635             foreach(r; TupleRange!(0, rows-1)) {
1636                 ret.matrix[r][rows-1] = matrix[r][rows-1];
1637             }
1638 
1639             return ret;
1640         }
1641 
1642         unittest {
1643             mat3 m3 = mat3(0.0f, 1.0f, 2.0f,
1644                            3.0f, 4.0f, 5.0f,
1645                            6.0f, 7.0f, 1.0f);
1646             assert(m3.get_translation().matrix == [[1.0f, 0.0f, 2.0f], [0.0f, 1.0f, 5.0f], [0.0f, 0.0f, 1.0f]]);
1647             m3.set_translation(mat3.identity);
1648             assert(mat3.identity.matrix == m3.get_translation().matrix);
1649             m3.set_translation([2.0f, 5.0f]);
1650             assert(m3.get_translation().matrix == [[1.0f, 0.0f, 2.0f], [0.0f, 1.0f, 5.0f], [0.0f, 0.0f, 1.0f]]);
1651             assert(mat3.identity.matrix == mat3.identity.get_translation().matrix);
1652 
1653             mat4 m4 = mat4(0.0f, 1.0f, 2.0f, 3.0f,
1654                            4.0f, 5.0f, 6.0f, 7.0f,
1655                            8.0f, 9.0f, 10.0f, 11.0f,
1656                            12.0f, 13.0f, 14.0f, 1.0f);
1657             assert(m4.get_translation().matrix == [[1.0f, 0.0f, 0.0f, 3.0f],
1658                                        [0.0f, 1.0f, 0.0f, 7.0f],
1659                                        [0.0f, 0.0f, 1.0f, 11.0f],
1660                                        [0.0f, 0.0f, 0.0f, 1.0f]]);
1661             m4.set_translation(mat4.identity);
1662             assert(mat4.identity.matrix == m4.get_translation().matrix);
1663             m4.set_translation([3.0f, 7.0f, 11.0f]);
1664             assert(m4.get_translation().matrix == [[1.0f, 0.0f, 0.0f, 3.0f],
1665                                        [0.0f, 1.0f, 0.0f, 7.0f],
1666                                        [0.0f, 0.0f, 1.0f, 11.0f],
1667                                        [0.0f, 0.0f, 0.0f, 1.0f]]);
1668             assert(mat4.identity.matrix == mat4.identity.get_translation().matrix);
1669         }
1670 
1671         /// Sets the scale of the matrix (nxn matrices, n >= 3).
1672         void set_scale(mt[] values...)
1673             in { assert(values.length >= (rows-1)); }
1674             body {
1675                 foreach(r; TupleRange!(0, rows-1)) {
1676                     matrix[r][r] = values[r];
1677                 }
1678             }
1679 
1680         /// Copyies the scale from mat to the current matrix (nxn matrices, n >= 3).
1681         void set_scale(Matrix mat) {
1682             foreach(r; TupleRange!(0, rows-1)) {
1683                 matrix[r][r] = mat.matrix[r][r];
1684             }
1685         }
1686 
1687         /// Returns an identity matrix with the current scale applied (nxn matrices, n >= 3).
1688         Matrix get_scale() {
1689             Matrix ret = Matrix.identity;
1690 
1691             foreach(r; TupleRange!(0, rows-1)) {
1692                 ret.matrix[r][r] = matrix[r][r];
1693             }
1694 
1695             return ret;
1696         }
1697 
1698         unittest {
1699             mat3 m3 = mat3(0.0f, 1.0f, 2.0f,
1700                            3.0f, 4.0f, 5.0f,
1701                            6.0f, 7.0f, 1.0f);
1702             assert(m3.get_scale().matrix == [[0.0f, 0.0f, 0.0f], [0.0f, 4.0f, 0.0f], [0.0f, 0.0f, 1.0f]]);
1703             m3.set_scale(mat3.identity);
1704             assert(mat3.identity.matrix == m3.get_scale().matrix);
1705             m3.set_scale([0.0f, 4.0f]);
1706             assert(m3.get_scale().matrix == [[0.0f, 0.0f, 0.0f], [0.0f, 4.0f, 0.0f], [0.0f, 0.0f, 1.0f]]);
1707             assert(mat3.identity.matrix == mat3.identity.get_scale().matrix);
1708 
1709             mat4 m4 = mat4(0.0f, 1.0f, 2.0f, 3.0f,
1710                            4.0f, 5.0f, 6.0f, 7.0f,
1711                            8.0f, 9.0f, 10.0f, 11.0f,
1712                            12.0f, 13.0f, 14.0f, 1.0f);
1713             assert(m4.get_scale().matrix == [[0.0f, 0.0f, 0.0f, 0.0f],
1714                                        [0.0f, 5.0f, 0.0f, 0.0f],
1715                                        [0.0f, 0.0f, 10.0f, 0.0f],
1716                                        [0.0f, 0.0f, 0.0f, 1.0f]]);
1717             m4.set_scale(mat4.identity);
1718             assert(mat4.identity.matrix == m4.get_scale().matrix);
1719             m4.set_scale([0.0f, 5.0f, 10.0f]);
1720             assert(m4.get_scale().matrix == [[0.0f, 0.0f, 0.0f, 0.0f],
1721                                        [0.0f, 5.0f, 0.0f, 0.0f],
1722                                        [0.0f, 0.0f, 10.0f, 0.0f],
1723                                        [0.0f, 0.0f, 0.0f, 1.0f]]);
1724             assert(mat4.identity.matrix == mat4.identity.get_scale().matrix);
1725         }
1726 
1727         /// Copies rot into the upper left corner, the translation (nxn matrices, n >= 3).
1728         void set_rotation(Matrix!(mt, 3, 3) rot) {
1729             foreach(r; TupleRange!(0, 3)) {
1730                 foreach(c; TupleRange!(0, 3)) {
1731                     matrix[r][c] = rot[r][c];
1732                 }
1733             }
1734         }
1735 
1736         /// Returns an identity matrix with the current rotation applied (nxn matrices, n >= 3).
1737         Matrix!(mt, 3, 3) get_rotation() {
1738             Matrix!(mt, 3, 3) ret = Matrix!(mt, 3, 3).identity;
1739 
1740             foreach(r; TupleRange!(0, 3)) {
1741                 foreach(c; TupleRange!(0, 3)) {
1742                     ret.matrix[r][c] = matrix[r][c];
1743                 }
1744             }
1745 
1746             return ret;
1747         }
1748 
1749         unittest {
1750             mat3 m3 = mat3(0.0f, 1.0f, 2.0f,
1751                            3.0f, 4.0f, 5.0f,
1752                            6.0f, 7.0f, 1.0f);
1753             assert(m3.get_rotation().matrix == [[0.0f, 1.0f, 2.0f], [3.0f, 4.0f, 5.0f], [6.0f, 7.0f, 1.0f]]);
1754             m3.set_rotation(mat3.identity);
1755             assert(mat3.identity.matrix == m3.get_rotation().matrix);
1756             m3.set_rotation(mat3(0.0f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 1.0f));
1757             assert(m3.get_rotation().matrix == [[0.0f, 1.0f, 2.0f], [3.0f, 4.0f, 5.0f], [6.0f, 7.0f, 1.0f]]);
1758             assert(mat3.identity.matrix == mat3.identity.get_rotation().matrix);
1759 
1760             mat4 m4 = mat4(0.0f, 1.0f, 2.0f, 3.0f,
1761                            4.0f, 5.0f, 6.0f, 7.0f,
1762                            8.0f, 9.0f, 10.0f, 11.0f,
1763                            12.0f, 13.0f, 14.0f, 1.0f);
1764             assert(m4.get_rotation().matrix == [[0.0f, 1.0f, 2.0f], [4.0f, 5.0f, 6.0f], [8.0f, 9.0f, 10.0f]]);
1765             m4.set_rotation(mat3.identity);
1766             assert(mat3.identity.matrix == m4.get_rotation().matrix);
1767             m4.set_rotation(mat3(0.0f, 1.0f, 2.0f, 4.0f, 5.0f, 6.0f, 8.0f, 9.0f, 10.0f));
1768             assert(m4.get_rotation().matrix == [[0.0f, 1.0f, 2.0f], [4.0f, 5.0f, 6.0f], [8.0f, 9.0f, 10.0f]]);
1769             assert(mat3.identity.matrix == mat4.identity.get_rotation().matrix);
1770         }
1771 
1772     }
1773 
1774     static if((rows == cols) && (rows >= 2) && (rows <= 4)) {
1775         /// Returns an inverted copy of the current matrix (nxn matrices, 2 >= n <= 4).
1776         @property Matrix inverse() const {
1777             Matrix mat;
1778             invert(mat);
1779             return mat;
1780         }
1781 
1782         /// Inverts the current matrix (nxn matrices, 2 >= n <= 4).
1783         void invert() {
1784             // workaround Issue #11238
1785             // uses a temporary instead of invert(this)
1786             Matrix temp;
1787             invert(temp);
1788             this.matrix = temp.matrix;
1789         }
1790     }
1791 
1792     unittest {
1793         mat2 m2 = mat2(1.0f, 2.0f, vec2(3.0f, 4.0f));
1794         assert(m2.det == -2.0f);
1795         assert(m2.inverse.matrix == [[-2.0f, 1.0f], [1.5f, -0.5f]]);
1796 
1797         mat3 m3 = mat3(1.0f, -2.0f, 3.0f,
1798                        7.0f, -1.0f, 0.0f,
1799                        3.0f, 2.0f, -4.0f);
1800         assert(m3.det == -1.0f);
1801         assert(m3.inverse.matrix == [[-4.0f, 2.0f, -3.0f],
1802                                      [-28.0f, 13.0f, -21.0f],
1803                                      [-17.0f, 8.0f, -13.0f]]);
1804 
1805         mat4 m4 = mat4(1.0f, 2.0f, 3.0f, 4.0f,
1806                        -2.0f, 1.0f, 5.0f, -2.0f,
1807                        2.0f, -1.0f, 7.0f, 1.0f,
1808                        3.0f, -3.0f, 2.0f, 0.0f);
1809         assert(m4.det == -8.0f);
1810         assert(m4.inverse.matrix == [[6.875f, 7.875f, -11.75f, 11.125f],
1811                                      [6.625f, 7.625f, -11.25f, 10.375f],
1812                                      [-0.375f, -0.375f, 0.75f, -0.625f],
1813                                      [-4.5f, -5.5f, 8.0f, -7.5f]]);
1814     }
1815 
1816     private void mms(mt inp, ref Matrix mat) const { // mat * scalar
1817         for(int r = 0; r < rows; r++) {
1818             for(int c = 0; c < cols; c++) {
1819                 mat.matrix[r][c] = matrix[r][c] * inp;
1820             }
1821         }
1822     }
1823 
1824     private void masm(string op)(Matrix inp, ref Matrix mat) const { // mat + or - mat
1825         foreach(r; TupleRange!(0, rows)) {
1826             foreach(c; TupleRange!(0, cols)) {
1827                 mat.matrix[r][c] = mixin("inp.matrix[r][c]" ~ op ~ "matrix[r][c]");
1828             }
1829         }
1830     }
1831 
1832     Matrix!(mt, rows, T.cols) opBinary(string op : "*", T)(T inp) const if(isCompatibleMatrix!T && (T.rows == cols)) {
1833         Matrix!(mt, rows, T.cols) ret;
1834 
1835         foreach(r; TupleRange!(0, rows)) {
1836             foreach(c; TupleRange!(0, T.cols)) {
1837                 ret.matrix[r][c] = 0;
1838 
1839                 foreach(c2; TupleRange!(0, cols)) {
1840                     ret.matrix[r][c] += matrix[r][c2] * inp.matrix[c2][c];
1841                 }
1842             }
1843         }
1844 
1845         return ret;
1846     }
1847 
1848     Vector!(mt, rows) opBinary(string op : "*", T : Vector!(mt, cols))(T inp) const {
1849         Vector!(mt, rows) ret;
1850         ret.clear(0);
1851 
1852         foreach(c; TupleRange!(0, cols)) {
1853             foreach(r; TupleRange!(0, rows)) {
1854                 ret.vector[r] += matrix[r][c] * inp.vector[c];
1855             }
1856         }
1857 
1858         return ret;
1859     }
1860 
1861     Matrix opBinary(string op : "*")(mt inp) const {
1862         Matrix ret;
1863         mms(inp, ret);
1864         return ret;
1865     }
1866 
1867     Matrix opBinaryRight(string op : "*")(mt inp) const {
1868         return this.opBinary!(op)(inp);
1869     }
1870 
1871     Matrix opBinary(string op)(Matrix inp) const if((op == "+") || (op == "-")) {
1872         Matrix ret;
1873         masm!(op)(inp, ret);
1874         return ret;
1875     }
1876 
1877     void opOpAssign(string op : "*")(mt inp) {
1878         mms(inp, this);
1879     }
1880 
1881     void opOpAssign(string op)(Matrix inp) if((op == "+") || (op == "-")) {
1882         masm!(op)(inp, this);
1883     }
1884 
1885     void opOpAssign(string op)(Matrix inp) if(op == "*") {
1886         this = this * inp;
1887     }
1888 
1889     unittest {
1890         mat2 m2 = mat2(1.0f, 2.0f, 3.0f, 4.0f);
1891         vec2 v2 = vec2(2.0f, 2.0f);
1892         assert((m2*2).matrix == [[2.0f, 4.0f], [6.0f, 8.0f]]);
1893         assert((2*m2).matrix == (m2*2).matrix);
1894         m2 *= 2;
1895         assert(m2.matrix == [[2.0f, 4.0f], [6.0f, 8.0f]]);
1896         assert((m2*v2).vector == [12.0f, 28.0f]);
1897         assert((v2*m2).vector == [16.0f, 24.0f]);
1898         assert((m2*m2).matrix == [[28.0f, 40.0f], [60.0f, 88.0f]]);
1899         assert((m2-m2).matrix == [[0.0f, 0.0f], [0.0f, 0.0f]]);
1900         assert((m2+m2).matrix == [[4.0f, 8.0f], [12.0f, 16.0f]]);
1901         m2 += m2;
1902         assert(m2.matrix == [[4.0f, 8.0f], [12.0f, 16.0f]]);
1903         m2 -= m2;
1904         assert(m2.matrix == [[0.0f, 0.0f], [0.0f, 0.0f]]);
1905 
1906         mat3 m3 = mat3(1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f);
1907         vec3 v3 = vec3(2.0f, 2.0f, 2.0f);
1908         assert((m3*2).matrix == [[2.0f, 4.0f, 6.0f], [8.0f, 10.0f, 12.0f], [14.0f, 16.0f, 18.0f]]);
1909         assert((2*m3).matrix == (m3*2).matrix);
1910         m3 *= 2;
1911         assert(m3.matrix == [[2.0f, 4.0f, 6.0f], [8.0f, 10.0f, 12.0f], [14.0f, 16.0f, 18.0f]]);
1912         assert((m3*v3).vector == [24.0f, 60.0f, 96.0f]);
1913         assert((v3*m3).vector == [48.0f, 60.0f, 72.0f]);
1914         assert((m3*m3).matrix == [[120.0f, 144.0f, 168.0f], [264.0f, 324.0f, 384.0f], [408.0f, 504.0f, 600.0f]]);
1915         assert((m3-m3).matrix == [[0.0f, 0.0f, 0.0f], [0.0f, 0.0f, 0.0f], [0.0f, 0.0f, 0.0f]]);
1916         assert((m3+m3).matrix == [[4.0f, 8.0f, 12.0f], [16.0f, 20.0f, 24.0f], [28.0f, 32.0f, 36.0f]]);
1917         m3 += m3;
1918         assert(m3.matrix == [[4.0f, 8.0f, 12.0f], [16.0f, 20.0f, 24.0f], [28.0f, 32.0f, 36.0f]]);
1919         m3 -= m3;
1920         assert(m3.matrix == [[0.0f, 0.0f, 0.0f], [0.0f, 0.0f, 0.0f], [0.0f, 0.0f, 0.0f]]);
1921 
1922         // test opOpAssign for matrix multiplication
1923         auto m4 = mat4.translation(0,1,2);
1924         m4 *= mat4.translation(0,-1,2);
1925         assert(m4 == mat4.translation(0,0,4));
1926 
1927         //TODO: tests for mat4, mat34
1928     }
1929 
1930     // opEqual => "alias matrix this;"
1931 
1932     bool opCast(T : bool)() const {
1933         return isFinite;
1934     }
1935 
1936     unittest {
1937         assert(mat2(1.0f, 2.0f, 1.0f, 1.0f) == mat2(1.0f, 2.0f, 1.0f, 1.0f));
1938         assert(mat2(1.0f, 2.0f, 1.0f, 1.0f) != mat2(1.0f, 1.0f, 1.0f, 1.0f));
1939 
1940         assert(mat3(1.0f) == mat3(1.0f));
1941         assert(mat3(1.0f) != mat3(2.0f));
1942 
1943         assert(mat4(1.0f) == mat4(1.0f));
1944         assert(mat4(1.0f) != mat4(2.0f));
1945 
1946         assert(!(mat4(float.nan)));
1947         if(mat4(1.0f)) { }
1948         else { assert(false); }
1949     }
1950 
1951 }
1952 
1953 /// Pre-defined matrix types, the first number represents the number of rows
1954 /// and the second the number of columns, if there's just one it's a nxn matrix.
1955 /// All of these matrices are floating-point matrices.
1956 alias Matrix!(float, 2, 2) mat2;
1957 alias Matrix!(float, 3, 3) mat3;
1958 alias Matrix!(float, 3, 4) mat34;
1959 alias Matrix!(float, 4, 4) mat4;
1960 
1961 private unittest {
1962     Matrix!(float,  1, 1) A = 1;
1963     Matrix!(double, 1, 1) B = 1;
1964     Matrix!(real,   1, 1) C = 1;
1965     Matrix!(int,    1, 1) D = 1;
1966     Matrix!(float,  5, 1) E = 1;
1967     Matrix!(double, 5, 1) F = 1;
1968     Matrix!(real,   5, 1) G = 1;
1969     Matrix!(int,    5, 1) H = 1;
1970     Matrix!(float,  1, 5) I = 1;
1971     Matrix!(double, 1, 5) J = 1;
1972     Matrix!(real,   1, 5) K = 1;
1973     Matrix!(int,    1, 5) L = 1;
1974 }
1975 
1976 /// Base template for all quaternion-types.
1977 /// Params:
1978 ///  type = all values get stored as this type
1979 struct Quaternion(type) {
1980     alias type qt; /// Holds the internal type of the quaternion.
1981 
1982     qt[4] quaternion; /// Holds the w, x, y and z coordinates.
1983 
1984     /// Returns a pointer to the quaternion in memory, it starts with the w coordinate.
1985     @property auto value_ptr() const { return quaternion.ptr; }
1986 
1987     /// Returns the current vector formatted as string, useful for printing the quaternion.
1988     @property string as_string() {
1989         return format("%s", quaternion);
1990     }
1991     alias as_string toString;
1992 
1993     @safe pure nothrow:
1994     private @property qt get_(char coord)() const {
1995         return quaternion[coord_to_index!coord];
1996     }
1997     private @property void set_(char coord)(qt value) {
1998         quaternion[coord_to_index!coord] = value;
1999     }
2000 
2001     alias get_!'w' w; /// static properties to access the values.
2002     alias set_!'w' w;
2003     alias get_!'x' x; /// ditto
2004     alias set_!'x' x;
2005     alias get_!'y' y; /// ditto
2006     alias set_!'y' y;
2007     alias get_!'z' z; /// ditto
2008     alias set_!'z' z;
2009 
2010     /// Constructs the quaternion.
2011     /// Takes a 4-dimensional vector, where vector.x = the quaternions w coordinate,
2012     /// or a w coordinate of type $(I qt) and a 3-dimensional vector representing the imaginary part,
2013     /// or 4 values of type $(I qt).
2014     this(qt w_, qt x_, qt y_, qt z_) {
2015         w = w_;
2016         x = x_;
2017         y = y_;
2018         z = z_;
2019     }
2020 
2021     /// ditto
2022     this(qt w_, Vector!(qt, 3) vec) {
2023         w = w_;
2024         quaternion[1..4] = vec.vector[];
2025     }
2026 
2027     /// ditto
2028     this(Vector!(qt, 4) vec) {
2029         quaternion[] = vec.vector[];
2030     }
2031 
2032     /// Returns true if all values are not nan and finite, otherwise false.
2033     @property bool isFinite() const {
2034         foreach(q; quaternion) {
2035             if(isNaN(q) || isInfinity(q)) {
2036                 return false;
2037             }
2038         }
2039         return true;
2040     }
2041     deprecated("Use isFinite instead of ok") alias ok = isFinite;
2042 
2043     unittest {
2044         quat q1 = quat(0.0f, 0.0f, 0.0f, 1.0f);
2045         assert(q1.quaternion == [0.0f, 0.0f, 0.0f, 1.0f]);
2046         assert(q1.quaternion == quat(0.0f, 0.0f, 0.0f, 1.0f).quaternion);
2047         assert(q1.quaternion == quat(0.0f, vec3(0.0f, 0.0f, 1.0f)).quaternion);
2048         assert(q1.quaternion == quat(vec4(0.0f, 0.0f, 0.0f, 1.0f)).quaternion);
2049 
2050         assert(q1.isFinite);
2051         q1.x = float.infinity;
2052         assert(!q1.isFinite);
2053         q1.x = float.nan;
2054         assert(!q1.isFinite);
2055         q1.x = 0.0f;
2056         assert(q1.isFinite);
2057     }
2058 
2059     template coord_to_index(char c) {
2060         static if(c == 'w') {
2061             enum coord_to_index = 0;
2062         } else static if(c == 'x') {
2063             enum coord_to_index = 1;
2064         } else static if(c == 'y') {
2065             enum coord_to_index = 2;
2066         } else static if(c == 'z') {
2067             enum coord_to_index = 3;
2068         } else {
2069             static assert(false, "accepted coordinates are x, y, z and w not " ~ c ~ ".");
2070         }
2071     }
2072 
2073     /// Returns the squared magnitude of the quaternion.
2074     @property real magnitude_squared() const {
2075         return to!real(w^^2 + x^^2 + y^^2 + z^^2);
2076     }
2077 
2078     /// Returns the magnitude of the quaternion.
2079     @property real magnitude() const {
2080         return sqrt(magnitude_squared);
2081     }
2082 
2083     /// Returns an identity quaternion (w=1, x=0, y=0, z=0).
2084     static @property Quaternion identity() {
2085         return Quaternion(1, 0, 0, 0);
2086     }
2087 
2088     /// Makes the current quaternion an identity quaternion.
2089     void make_identity() {
2090         w = 1;
2091         x = 0;
2092         y = 0;
2093         z = 0;
2094     }
2095 
2096     /// Inverts the quaternion.
2097     void invert() {
2098         x = -x;
2099         y = -y;
2100         z = -z;
2101     }
2102     alias invert conjugate; /// ditto
2103 
2104     /// Returns an inverted copy of the current quaternion.
2105     @property Quaternion inverse() const {
2106         return Quaternion(w, -x, -y, -z);
2107     }
2108     alias inverse conjugated; /// ditto
2109 
2110     unittest {
2111         quat q1 = quat(1.0f, 1.0f, 1.0f, 1.0f);
2112 
2113         assert(q1.magnitude == 2.0f);
2114         assert(q1.magnitude_squared == 4.0f);
2115         assert(q1.magnitude == quat(0.0f, 0.0f, 2.0f, 0.0f).magnitude);
2116 
2117         quat q2 = quat.identity;
2118         assert(q2.quaternion == [1.0f, 0.0f, 0.0f, 0.0f]);
2119         assert(q2.x == 0.0f);
2120         assert(q2.y == 0.0f);
2121         assert(q2.z == 0.0f);
2122         assert(q2.w == 1.0f);
2123 
2124         assert(q1.inverse.quaternion == [1.0f, -1.0f, -1.0f, -1.0f]);
2125         q1.invert();
2126         assert(q1.quaternion == [1.0f, -1.0f, -1.0f, -1.0f]);
2127 
2128         q1.make_identity();
2129         assert(q1.quaternion == q2.quaternion);
2130 
2131     }
2132 
2133     /// Creates a quaternion from a 3x3 matrix.
2134     /// Params:
2135     ///  matrix = 3x3 matrix (rotation)
2136     /// Returns: A quaternion representing the rotation (3x3 matrix)
2137     static Quaternion from_matrix(Matrix!(qt, 3, 3) matrix) {
2138         Quaternion ret;
2139 
2140         auto mat = matrix.matrix;
2141         qt trace = mat[0][0] + mat[1][1] + mat[2][2];
2142 
2143         if(trace > 0) {
2144             real s = 0.5 / sqrt(trace + 1.0f);
2145 
2146             ret.w = to!qt(0.25 / s);
2147             ret.x = to!qt((mat[2][1] - mat[1][2]) * s);
2148             ret.y = to!qt((mat[0][2] - mat[2][0]) * s);
2149             ret.z = to!qt((mat[1][0] - mat[0][1]) * s);
2150         } else if((mat[0][0] > mat[1][1]) && (mat[0][0] > mat[2][2])) {
2151             real s = 2.0 * sqrt(1.0 + mat[0][0] - mat[1][1] - mat[2][2]);
2152 
2153             ret.w = to!qt((mat[2][1] - mat[1][2]) / s);
2154             ret.x = to!qt(0.25f * s);
2155             ret.y = to!qt((mat[0][1] + mat[1][0]) / s);
2156             ret.z = to!qt((mat[0][2] + mat[2][0]) / s);
2157         } else if(mat[1][1] > mat[2][2]) {
2158             real s = 2.0 * sqrt(1 + mat[1][1] - mat[0][0] - mat[2][2]);
2159 
2160             ret.w = to!qt((mat[0][2] - mat[2][0]) / s);
2161             ret.x = to!qt((mat[0][1] + mat[1][0]) / s);
2162             ret.y = to!qt(0.25f * s);
2163             ret.z = to!qt((mat[1][2] + mat[2][1]) / s);
2164         } else {
2165             real s = 2.0 * sqrt(1 + mat[2][2] - mat[0][0] - mat[1][1]);
2166 
2167             ret.w = to!qt((mat[1][0] - mat[0][1]) / s);
2168             ret.x = to!qt((mat[0][2] + mat[2][0]) / s);
2169             ret.y = to!qt((mat[1][2] + mat[2][1]) / s);
2170             ret.z = to!qt(0.25f * s);
2171         }
2172 
2173         return ret;
2174     }
2175 
2176     /// Returns the quaternion as matrix.
2177     /// Params:
2178     ///  rows = number of rows of the resulting matrix (min 3)
2179     ///  cols = number of columns of the resulting matrix (min 3)
2180     Matrix!(qt, rows, cols) to_matrix(int rows, int cols)() const if((rows >= 3) && (cols >= 3)) {
2181         static if((rows == 3) && (cols == 3)) {
2182             Matrix!(qt, rows, cols) ret;
2183         } else {
2184             Matrix!(qt, rows, cols) ret = Matrix!(qt, rows, cols).identity;
2185         }
2186 
2187         qt xx = x^^2;
2188         qt xy = x * y;
2189         qt xz = x * z;
2190         qt xw = x * w;
2191         qt yy = y^^2;
2192         qt yz = y * z;
2193         qt yw = y * w;
2194         qt zz = z^^2;
2195         qt zw = z * w;
2196 
2197         ret.matrix[0][0] = 1 - 2 * (yy + zz);
2198         ret.matrix[0][1] = 2 * (xy - zw);
2199         ret.matrix[0][2] = 2 * (xz + yw);
2200 
2201         ret.matrix[1][0] = 2 * (xy + zw);
2202         ret.matrix[1][1] = 1 - 2 * (xx + zz);
2203         ret.matrix[1][2] = 2 * (yz - xw);
2204 
2205         ret.matrix[2][0] = 2 * (xz - yw);
2206         ret.matrix[2][1] = 2 * (yz + xw);
2207         ret.matrix[2][2] = 1 - 2 * (xx + yy);
2208 
2209         return ret;
2210     }
2211 
2212     unittest {
2213         quat q1 = quat(4.0f, 1.0f, 2.0f, 3.0f);
2214 
2215         assert(q1.to_matrix!(3, 3).matrix == [[-25.0f, -20.0f, 22.0f], [28.0f, -19.0f, 4.0f], [-10.0f, 20.0f, -9.0f]]);
2216         assert(q1.to_matrix!(4, 4).matrix == [[-25.0f, -20.0f, 22.0f, 0.0f],
2217                                               [28.0f, -19.0f, 4.0f, 0.0f],
2218                                               [-10.0f, 20.0f, -9.0f, 0.0f],
2219                                               [0.0f, 0.0f, 0.0f, 1.0f]]);
2220         assert(quat.identity.to_matrix!(3, 3).matrix == Matrix!(qt, 3, 3).identity.matrix);
2221         assert(q1.quaternion == quat.from_matrix(q1.to_matrix!(3, 3)).quaternion);
2222 
2223         assert(quat(1.0f, 0.0f, 0.0f, 0.0f).quaternion == quat.from_matrix(mat3.identity).quaternion);
2224 
2225         quat q2 = quat.from_matrix(mat3(1.0f, 3.0f, 2.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f));
2226         assert(q2.x == 0.0f);
2227         assert(almost_equal(q2.y, 0.7071067f));
2228         assert(almost_equal(q2.z, -1.060660));
2229         assert(almost_equal(q2.w, 0.7071067f));
2230     }
2231 
2232     /// Normalizes the current quaternion.
2233     void normalize() {
2234         qt m = to!qt(magnitude);
2235 
2236         if(m != 0) {
2237             w = w / m;
2238             x = x / m;
2239             y = y / m;
2240             z = z / m;
2241         }
2242     }
2243 
2244     /// Returns a normalized copy of the current quaternion.
2245     Quaternion normalized() const {
2246         Quaternion ret;
2247         qt m = to!qt(magnitude);
2248 
2249         if(m != 0) {
2250             ret.w = w / m;
2251             ret.x = x / m;
2252             ret.y = y / m;
2253             ret.z = z / m;
2254         } else {
2255             ret = Quaternion(w, x, y, z);
2256         }
2257 
2258         return ret;
2259     }
2260 
2261     unittest {
2262         quat q1 = quat(1.0f, 2.0f, 3.0f, 4.0f);
2263         quat q2 = quat(1.0f, 2.0f, 3.0f, 4.0f);
2264 
2265         q1.normalize();
2266         assert(q1.quaternion == q2.normalized.quaternion);
2267         //assert(q1.quaternion == q1.normalized.quaternion);
2268         assert(almost_equal(q1.magnitude, 1.0));
2269     }
2270 
2271     /// Returns the yaw.
2272     @property real yaw() const {
2273         return atan2(to!real(2.0*(w*z + x*y)), to!real(1.0 - 2.0*(y*y + z*z)));
2274     }
2275 
2276     /// Returns the pitch.
2277     @property real pitch() const {
2278         return asin(to!real(2.0*(w*y - z*x)));
2279     }
2280 
2281     /// Returns the roll.
2282     @property real roll() const {
2283         return atan2(to!real(2.0*(w*x + y*z)), to!real(1.0 - 2.0*(x*x + y*y)));
2284     }
2285 
2286     unittest {
2287         quat q1 = quat.identity;
2288         assert(q1.pitch == 0.0f);
2289         assert(q1.yaw == 0.0f);
2290         assert(q1.roll == 0.0f);
2291 
2292         quat q2 = quat(1.0f, 1.0f, 1.0f, 1.0f);
2293         assert(almost_equal(q2.yaw, q2.roll));
2294         //assert(almost_equal(q2.yaw, 1.570796f));
2295         assert(q2.pitch == 0.0f);
2296 
2297         quat q3 = quat(0.1f, 1.9f, 2.1f, 1.3f);
2298         //assert(almost_equal(q3.yaw, 2.4382f));
2299         assert(isNaN(q3.pitch));
2300         //assert(almost_equal(q3.roll, 1.67719f));
2301     }
2302 
2303     /// Returns a quaternion with applied rotation around the x-axis.
2304     static Quaternion xrotation(real alpha) {
2305         Quaternion ret;
2306 
2307         alpha /= 2;
2308         ret.w = to!qt(cos(alpha));
2309         ret.x = to!qt(sin(alpha));
2310         ret.y = 0;
2311         ret.z = 0;
2312 
2313         return ret;
2314     }
2315 
2316     /// Returns a quaternion with applied rotation around the y-axis.
2317     static Quaternion yrotation(real alpha) {
2318         Quaternion ret;
2319 
2320         alpha /= 2;
2321         ret.w = to!qt(cos(alpha));
2322         ret.x = 0;
2323         ret.y = to!qt(sin(alpha));
2324         ret.z = 0;
2325 
2326         return ret;
2327     }
2328 
2329     /// Returns a quaternion with applied rotation around the z-axis.
2330     static Quaternion zrotation(real alpha) {
2331         Quaternion ret;
2332 
2333         alpha /= 2;
2334         ret.w = to!qt(cos(alpha));
2335         ret.x = 0;
2336         ret.y = 0;
2337         ret.z = to!qt(sin(alpha));
2338 
2339         return ret;
2340     }
2341 
2342     /// Returns a quaternion with applied rotation around an axis.
2343     static Quaternion axis_rotation(real alpha, Vector!(qt, 3) axis) {
2344         if(alpha == 0) {
2345             return Quaternion.identity;
2346         }
2347         Quaternion ret;
2348 
2349         alpha /= 2;
2350         qt sinaqt = to!qt(sin(alpha));
2351 
2352         ret.w = to!qt(cos(alpha));
2353         ret.x = axis.x * sinaqt;
2354         ret.y = axis.y * sinaqt;
2355         ret.z = axis.z * sinaqt;
2356 
2357         return ret;
2358     }
2359 
2360     /// Creates a quaternion from an euler rotation.
2361     static Quaternion euler_rotation(real roll, real pitch, real yaw) {
2362         Quaternion ret;
2363 
2364         auto cr = cos(roll / 2.0);
2365         auto cp = cos(pitch / 2.0);
2366         auto cy = cos(yaw / 2.0);
2367         auto sr = sin(roll / 2.0);
2368         auto sp = sin(pitch / 2.0);
2369         auto sy = sin(yaw / 2.0);
2370 
2371         ret.w = cr * cp * cy + sr * sp * sy;
2372         ret.x = sr * cp * cy - cr * sp * sy;
2373         ret.y = cr * sp * cy + sr * cp * sy;
2374         ret.z = cr * cp * sy - sr * sp * cy;
2375 
2376         return ret;
2377     }
2378 
2379     unittest {
2380         enum startPitch = 0.1;
2381         enum startYaw = -0.2;
2382         enum startRoll = 0.6;
2383         
2384         auto q = quat.euler_rotation(startRoll,startPitch,startYaw);
2385         
2386         assert(almost_equal(q.pitch,startPitch));
2387         assert(almost_equal(q.yaw,startYaw));
2388         assert(almost_equal(q.roll,startRoll));
2389     }
2390 
2391     /// Rotates the current quaternion around the x-axis and returns $(I this).
2392     Quaternion rotatex(real alpha) {
2393         this = xrotation(alpha) * this;
2394         return this;
2395     }
2396 
2397     /// Rotates the current quaternion around the y-axis and returns $(I this).
2398     Quaternion rotatey(real alpha) {
2399         this = yrotation(alpha) * this;
2400         return this;
2401     }
2402 
2403     /// Rotates the current quaternion around the z-axis and returns $(I this).
2404     Quaternion rotatez(real alpha) {
2405         this = zrotation(alpha) * this;
2406         return this;
2407     }
2408 
2409     /// Rotates the current quaternion around an axis and returns $(I this).
2410     Quaternion rotate_axis(real alpha, Vector!(qt, 3) axis) {
2411         this = axis_rotation(alpha, axis) * this;
2412         return this;
2413     }
2414 
2415     /// Applies an euler rotation to the current quaternion and returns $(I this).
2416     Quaternion rotate_euler(real heading, real attitude, real bank) {
2417         this = euler_rotation(heading, attitude, bank) * this;
2418         return this;
2419     }
2420 
2421     unittest {
2422         assert(quat.xrotation(PI).quaternion[1..4] == [1.0f, 0.0f, 0.0f]);
2423         assert(quat.yrotation(PI).quaternion[1..4] == [0.0f, 1.0f, 0.0f]);
2424         assert(quat.zrotation(PI).quaternion[1..4] == [0.0f, 0.0f, 1.0f]);
2425         assert((quat.xrotation(PI).w == quat.yrotation(PI).w) && (quat.yrotation(PI).w == quat.zrotation(PI).w));
2426         //assert(quat.rotatex(PI).w == to!(quat.qt)(cos(PI)));
2427         assert(quat.xrotation(PI).quaternion == quat.identity.rotatex(PI).quaternion);
2428         assert(quat.yrotation(PI).quaternion == quat.identity.rotatey(PI).quaternion);
2429         assert(quat.zrotation(PI).quaternion == quat.identity.rotatez(PI).quaternion);
2430 
2431         assert(quat.axis_rotation(PI, vec3(1.0f, 1.0f, 1.0f)).quaternion[1..4] == [1.0f, 1.0f, 1.0f]);
2432         assert(quat.axis_rotation(PI, vec3(1.0f, 1.0f, 1.0f)).w == quat.xrotation(PI).w);
2433         assert(quat.axis_rotation(PI, vec3(1.0f, 1.0f, 1.0f)).quaternion ==
2434                quat.identity.rotate_axis(PI, vec3(1.0f, 1.0f, 1.0f)).quaternion);
2435 
2436         quat q1 = quat.euler_rotation(PI, PI, PI);
2437         assert((q1.x > -1e-16) && (q1.x < 1e-16));
2438         assert((q1.y > -1e-16) && (q1.y < 1e-16));
2439         assert((q1.z > -1e-16) && (q1.z < 1e-16));
2440         //assert(q1.w == -1.0f);
2441         assert(quat.euler_rotation(PI, PI, PI).quaternion == quat.identity.rotate_euler(PI, PI, PI).quaternion);
2442     }
2443 
2444     Quaternion opBinary(string op : "*")(Quaternion inp) const {
2445         Quaternion ret;
2446 
2447         ret.w = -x * inp.x - y * inp.y - z * inp.z + w * inp.w;
2448         ret.x = x * inp.w + y * inp.z - z * inp.y + w * inp.x;
2449         ret.y = -x * inp.z + y * inp.w + z * inp.x + w * inp.y;
2450         ret.z = x * inp.y - y * inp.x + z * inp.w + w * inp.z;
2451 
2452         return ret;
2453     }
2454 
2455     auto opBinaryRight(string op, T)(T inp) const if(!is_quaternion!T) {
2456         return this.opBinary!(op)(inp);
2457     }
2458 
2459     Quaternion opBinary(string op)(Quaternion inp) const  if((op == "+") || (op == "-")) {
2460         Quaternion ret;
2461 
2462         mixin("ret.w = w" ~ op ~ "inp.w;");
2463         mixin("ret.x = x" ~ op ~ "inp.x;");
2464         mixin("ret.y = y" ~ op ~ "inp.y;");
2465         mixin("ret.z = z" ~ op ~ "inp.z;");
2466 
2467         return ret;
2468     }
2469 
2470     Vector!(qt, 3) opBinary(string op : "*")(Vector!(qt, 3) inp) const {
2471         Vector!(qt, 3) ret;
2472 
2473         qt ww = w^^2;
2474         qt w2 = w * 2;
2475         qt wx2 = w2 * x;
2476         qt wy2 = w2 * y;
2477         qt wz2 = w2 * z;
2478         qt xx = x^^2;
2479         qt x2 = x * 2;
2480         qt xy2 = x2 * y;
2481         qt xz2 = x2 * z;
2482         qt yy = y^^2;
2483         qt yz2 = 2 * y * z;
2484         qt zz = z * z;
2485 
2486         ret.vector =  [ww * inp.x + wy2 * inp.z - wz2 * inp.y + xx * inp.x +
2487                        xy2 * inp.y + xz2 * inp.z - zz * inp.x - yy * inp.x,
2488                        xy2 * inp.x + yy * inp.y + yz2 * inp.z + wz2 * inp.x -
2489                        zz * inp.y + ww * inp.y - wx2 * inp.z - xx * inp.y,
2490                        xz2 * inp.x + yz2 * inp.y + zz * inp.z - wy2 * inp.x -
2491                        yy * inp.z + wx2 * inp.y - xx * inp.z + ww * inp.z];
2492 
2493        return ret;
2494     }
2495 
2496     Quaternion opBinary(string op : "*")(qt inp) const {
2497         return Quaternion(w*inp, x*inp, y*inp, z*inp);
2498     }
2499 
2500     void opOpAssign(string op : "*")(Quaternion inp) {
2501         qt w2 = -x * inp.x - y * inp.y - z * inp.z + w * inp.w;
2502         qt x2 = x * inp.w + y * inp.z - z * inp.y + w * inp.x;
2503         qt y2 = -x * inp.z + y * inp.w + z * inp.x + w * inp.y;
2504         qt z2 = x * inp.y - y * inp.x + z * inp.w + w * inp.z;
2505         w = w2; x = x2; y = y2; z = z2;
2506     }
2507 
2508     void opOpAssign(string op)(Quaternion inp) if((op == "+") || (op == "-")) {
2509         mixin("w = w" ~ op ~ "inp.w;");
2510         mixin("x = x" ~ op ~ "inp.x;");
2511         mixin("y = y" ~ op ~ "inp.y;");
2512         mixin("z = z" ~ op ~ "inp.z;");
2513     }
2514 
2515     void opOpAssign(string op : "*")(qt inp) {
2516         quaternion[0] *= inp;
2517         quaternion[1] *= inp;
2518         quaternion[2] *= inp;
2519         quaternion[3] *= inp;
2520     }
2521 
2522     unittest {
2523         quat q1 = quat.identity;
2524         quat q2 = quat(3.0f, 0.0f, 1.0f, 2.0f);
2525         quat q3 = quat(3.4f, 0.1f, 1.2f, 2.3f);
2526 
2527         assert((q1 * q1).quaternion == q1.quaternion);
2528         assert((q1 * q2).quaternion == q2.quaternion);
2529         assert((q2 * q1).quaternion == q2.quaternion);
2530         quat q4 = q3 * q2;
2531         assert((q2 * q3).quaternion != q4.quaternion);
2532         q3 *= q2;
2533         assert(q4.quaternion == q3.quaternion);
2534         assert(almost_equal(q4.x, 0.4f));
2535         assert(almost_equal(q4.y, 6.8f));
2536         assert(almost_equal(q4.z, 13.8f));
2537         assert(almost_equal(q4.w, 4.4f));
2538 
2539         quat q5 = quat(1.0f, 2.0f, 3.0f, 4.0f);
2540         quat q6 = quat(3.0f, 1.0f, 6.0f, 2.0f);
2541 
2542         assert((q5 - q6).quaternion == [-2.0f, 1.0f, -3.0f, 2.0f]);
2543         assert((q5 + q6).quaternion == [4.0f, 3.0f, 9.0f, 6.0f]);
2544         assert((q6 - q5).quaternion == [2.0f, -1.0f, 3.0f, -2.0f]);
2545         assert((q6 + q5).quaternion == [4.0f, 3.0f, 9.0f, 6.0f]);
2546         q5 += q6;
2547         assert(q5.quaternion == [4.0f, 3.0f, 9.0f, 6.0f]);
2548         q6 -= q6;
2549         assert(q6.quaternion == [0.0f, 0.0f, 0.0f, 0.0f]);
2550 
2551         quat q7 = quat(2.0f, 2.0f, 2.0f, 2.0f);
2552         assert((q7 * 2).quaternion == [4.0f, 4.0f, 4.0f, 4.0f]);
2553         assert((2 * q7).quaternion == (q7 * 2).quaternion);
2554         q7 *= 2;
2555         assert(q7.quaternion == [4.0f, 4.0f, 4.0f, 4.0f]);
2556 
2557         vec3 v1 = vec3(1.0f, 2.0f, 3.0f);
2558         assert((q1 * v1).vector == v1.vector);
2559         assert((v1 * q1).vector == (q1 * v1).vector);
2560         assert((q2 * v1).vector == [-2.0f, 36.0f, 38.0f]);
2561     }
2562 
2563     int opCmp(ref const Quaternion qua) const {
2564         foreach(i, a; quaternion) {
2565             if(a < qua.quaternion[i]) {
2566                 return -1;
2567             } else if(a > qua.quaternion[i]) {
2568                 return 1;
2569             }
2570         }
2571 
2572         // Quaternions are the same
2573         return 0;
2574     }
2575 
2576     bool opEquals(const Quaternion qu) const {
2577         return quaternion == qu.quaternion;
2578     }
2579 
2580     bool opCast(T : bool)() const  {
2581         return isFinite;
2582     }
2583 
2584     unittest {
2585         assert(quat(1.0f, 2.0f, 3.0f, 4.0f) == quat(1.0f, 2.0f, 3.0f, 4.0f));
2586         assert(quat(1.0f, 2.0f, 3.0f, 4.0f) != quat(1.0f, 2.0f, 3.0f, 3.0f));
2587 
2588         assert(!(quat(float.nan, float.nan, float.nan, float.nan)));
2589         if(quat(1.0f, 1.0f, 1.0f, 1.0f)) { }
2590         else { assert(false); }
2591     }
2592 
2593 }
2594 
2595 /// Pre-defined quaternion of type float.
2596 alias Quaternion!(float) quat;