Unofficial Dub Package Documentation 
1 module des.math.linear.matrix;
2 
3 import std.math;
4 import std.traits;
5 import std.algorithm;
6 import std.range;
7 import std.array;
8 import std.exception;
9 
10 import des.util.testsuite;
11 
12 import des.math.util;
13 import des.math.linear.vector;
14 import des.math.linear.quaterni;
15 import des.math.basic.traits;
16 
17 ///
18 template isMatrix(E)
19 {
20     enum isMatrix = is( typeof( impl(E.init) ) );
21     void impl(size_t H,size_t W,T)( Matrix!(H,W,T) ){}
22 }
23 
24 ///
25 template isStaticMatrix(E)
26 {
27     static if( !isMatrix!E )
28         enum isStaticMatrix = false;
29     else enum isStaticMatrix = E.isStatic;
30 }
31 
32 ///
33 template isDynamicMatrix(E)
34 {
35     static if( !isMatrix!E )
36         enum isDynamicMatrix = false;
37     else enum isDynamicMatrix = E.isDynamic;
38 }
39 
40 unittest
41 {
42     static assert( !isStaticMatrix!float );
43     static assert( !isDynamicMatrix!float );
44 }
45 
46 private @property
47 {
48     import std..string;
49 
50     string identityMatrixDataString(size_t S)()
51     { return diagMatrixDataString!(S,S)(1); }
52 
53     string zerosMatrixDataString(size_t H, size_t W)()
54     { return diagMatrixDataString!(H,W)(0); }
55 
56     string diagMatrixDataString(size_t H, size_t W)( size_t val )
57     {
58         string[] ret;
59         foreach( i; 0 .. H )
60         {
61             string[] buf;
62             foreach( j; 0 .. W )
63                 buf ~= format( "%d", i == j ? val : 0 );
64             ret ~= "[" ~ buf.join(",") ~ "]";
65         }
66         return "[" ~ ret.join(",") ~ "]";
67     }
68 
69     unittest
70     {
71         assert( identityMatrixDataString!3 == "[[1,0,0],[0,1,0],[0,0,1]]" );
72         assert( zerosMatrixDataString!(3,3) == "[[0,0,0],[0,0,0],[0,0,0]]" );
73     }
74 
75     string castArrayString(string type, size_t H, size_t W)()
76     { return format( "cast(%s[%d][%d])", type, W, H ); }
77 }
78 
79 /++
80  +/
81 struct Matrix(size_t H, size_t W, E)
82 {
83     ///
84     alias selftype = Matrix!(H,W,E);
85 
86     ///
87     alias datatype = E;
88 
89     /// `H == 0 || W == 0`
90     enum bool isDynamic = H == 0 || W == 0;
91     /// `H != 0 && W != 0`
92     enum bool isStatic = H != 0 && W != 0;
93 
94     /// `H == 0`
95     enum bool isDynamicHeight = H == 0;
96     /// `H != 0`
97     enum bool isStaticHeight = H != 0;
98 
99     /// `W == 0`
100     enum bool isDynamicWidth = W == 0;
101     /// `W != 0`
102     enum bool isStaticWidth = W != 0;
103 
104     /// `isStaticHeight && isDynamicWidth`
105     enum bool isStaticHeightOnly = isStaticHeight && isDynamicWidth;
106     /// `isStaticWidth && isDynamicHeight`
107     enum bool isStaticWidthOnly = isStaticWidth && isDynamicHeight;
108 
109     /// `isDynamicHeight && isDynamicWidth`
110     enum bool isDynamicAll = isDynamicHeight && isDynamicWidth;
111     /// `isStaticWidthOnly || isStaticHeightOnly`
112     enum bool isDynamicOne = isStaticWidthOnly || isStaticHeightOnly;
113 
114     static if( isStatic )
115     {
116         static if( isNumeric!E )
117         {
118             static if( H == W )
119                 /// static data ( if isNumeric!E fills valid numbers: if squred then identity matrix else zeros )
120                 E[W][H] data = mixin( castArrayString!("E",H,H) ~ identityMatrixDataString!H );
121             else 
122                 E[W][H] data = mixin( castArrayString!("E",H,W) ~ zerosMatrixDataString!(H,W) );
123         }
124         else E[W][H] data;
125     }
126     else static if( isStaticWidthOnly )
127         /// static width only
128         E[W][] data;
129     else static if( isStaticHeightOnly )
130         /// static height only
131         E[][H] data;
132     else
133         /// full dynamic
134         E[][] data;
135 
136     ///
137     alias data this;
138 
139 pure:
140 
141     private static bool allowSomeOp(size_t A, size_t B )
142     { return A==B || A==0 || B==0; }
143 
144     static if( isStatic || isDynamicOne )
145     {
146         /++ fill data and set dynamic size for `isDynamicOne` matrix
147 
148             only:
149                 if( isStatic || isDynamicOne )
150          +/
151         this(X...)( in X vals )
152             if( is(typeof(flatData!E(vals))) )
153         {
154             static if( isStatic )
155             {
156                 static if( hasNoDynamic!X )
157                 {
158                     static if( X.length > 1 )
159                     {
160                         static assert( getElemCount!X == W*H, "wrong args count" );
161                         static assert( isConvertable!(E,X), "wrong args type" );
162                         mixin( matrixStaticFill!("E","data","vals",W,E,X) );
163                     }
164                     else static if( X.length == 1 && isStaticMatrix!(X[0]) )
165                     {
166                         static assert( X[0].width == W && X[0].height == H );
167                         foreach( y; 0 .. H )
168                             foreach( x; 0 .. W )
169                                 data[y][x] = vals[0][y][x];
170                     }
171                     else enum __DF=true;
172                 }
173                 else enum __DF=true;
174 
175                 static if( is(typeof(__DF)) )
176                     _fillData( flatData!E(vals) );
177             }
178             else static if( isStaticWidthOnly )
179             {
180                 auto buf = flatData!E(vals);
181                 enforce( !(buf.length%W), "wrong args length" );
182                 resize( buf.length/W, W );
183                 _fillData( buf );
184             }
185             else static if( isStaticHeightOnly )
186             {
187                 auto buf = flatData!E(vals);
188                 enforce( !(buf.length%H), "wrong args length" );
189                 resize( H, buf.length/H );
190                 _fillData( buf );
191             }
192         }
193     }
194     else
195     {
196         /++ set size and fill data
197 
198             only:
199                 if( isDynamicAll )
200          +/
201         this(X...)( size_t iH, size_t iW, in X vals )
202         {
203             auto buf = flatData!E(vals);
204             enforce( buf.length == iH * iW || buf.length == 1, "wrong args length" );
205             resize( iH, iW );
206             _fillData( buf );
207         }
208 
209         /++ set size
210 
211             only:
212                 if( isDynamicAll )
213          +/
214         this( size_t iH, size_t iW ) { resize( iH, iW ); }
215     }
216 
217     ///
218     this(size_t oH, size_t oW, X)( in Matrix!(oH,oW,X) mtr )
219         if( is( typeof( E(X.init) ) ) )
220     {
221         static if( isDynamic )
222             resize( mtr.height, mtr.width );
223         foreach( i; 0 .. mtr.height )
224             foreach( j; 0 .. mtr.width )
225                 data[i][j] = E(mtr[i][j]);
226     }
227 
228     static if( isDynamic )
229     {
230         /++
231             only:
232                 if( isDynamic )
233          +/
234         this(this)
235         {
236             data = data.dup;
237             foreach( ref row; data )
238                 row = row.dup;
239         }
240 
241         /++
242             only:
243                 if( isDynamic )
244          +/
245         ref typeof(this) opAssign(size_t bH, size_t bW, X)( in Matrix!(bH,bW,X) b )
246             if( allowSomeOp(H,bH) && allowSomeOp(W,bW) && is(typeof(E(X.init))) )
247         {
248             static if( isStaticHeight && b.isDynamicHeight ) enforce( height == b.height, "wrong height" );
249             static if( isStaticWidth && b.isDynamicWidth ) enforce( width == b.width, "wrong width" );
250             static if( isDynamic ) resize(b.height,b.width);
251             _fillData(flatData!E(b.data));
252             return this;
253         }
254     }
255 
256     private void _fillData( in E[] vals... )
257     {
258         enforce( vals.length > 0, "no vals to fill" );
259         if( vals.length > 1 )
260         {
261             enforce( vals.length == height * width, "wrong vals length" );
262             size_t k = 0;
263             foreach( ref row; data )
264                 foreach( ref v; row )
265                     v = vals[k++];
266         }
267         else foreach( ref row; data ) row[] = vals[0];
268     }
269 
270     /// fill data
271     ref typeof(this) fill( in E[] vals... )
272     {
273         _fillData( vals );
274         return this;
275     }
276 
277     static if( (W==H && isStatic) || isDynamicOne )
278     {
279         /++ get diagonal matrix, diagonal elements fills cyclically
280 
281             only:
282                 if( (W==H && isStatic) || isDynamicOne )
283          +/
284         static auto diag(X...)( in X vals )
285         if( X.length > 0 && is(typeof(E(0))) && is(typeof(flatData!E(vals))) )
286         {
287             selftype ret;
288             static if( isStaticHeight ) auto L = H;
289             else auto L = W;
290 
291             static if( ret.isDynamic )
292                 ret.resize(L,L);
293 
294             ret._fillData( E(0) );
295             ret.fillDiag( flatData!E(vals) );
296 
297             return ret;
298         }
299     }
300 
301     ///
302     ref typeof(this) fillDiag( in E[] vals... )
303     {
304         enforce( vals.length > 0, "no vals to fill" );
305         enforce( height == width, "not squared" );
306         size_t k = 0;
307         foreach( i; 0 .. height )
308             data[i][i] = vals[k++%$];
309         return this;
310     }
311 
312     static if( isStaticHeight )
313         /++
314             only:
315                 if( isStaticHeight )
316          +/
317         enum height = H;
318     else
319     {
320         /// if isStaticHeight it's enum (not available to set)
321         @property size_t height() const { return data.length; }
322         /// if isStaticHeight it's enum (not available to set)
323         @property size_t height( size_t nh )
324         {
325             resize( nh, width );
326             return nh;
327         }
328     }
329 
330     static if( isStaticWidth )
331         /++
332             only:
333                 if( isStaticWidth )
334          +/
335         enum width = W;
336     else
337     {
338         private size_t _width = W;
339         /// if isStaticWidth it's enum (not available to set)
340         @property size_t width() const { return _width; }
341         /// if isStaticWidth it's enum (not available to set)
342         @property size_t width( size_t nw )
343         {
344             resize( height, nw );
345             return nw;
346         }
347     }
348 
349     static if( isDynamic )
350     {
351         /++ resize dynamic matrix, new height/width must equals static height/width
352 
353             only:
354                 if( isDynamic )
355          +/
356         ref typeof(this) resize( size_t nh, size_t nw )
357         {
358             static if( isStaticHeight ) enforce( nh == H, "height is static" );
359             static if( isStaticWidth ) enforce( nw == W, "width is static" );
360 
361             static if( isDynamicHeight ) data.length = nh;
362             static if( isDynamicWidth )
363             {
364                 _width = nw;
365                 foreach( i; 0 .. height )
366                     data[i].length = nw;
367             }
368             return this;
369         }
370     }
371 
372     ///
373     auto expandHeight(size_t bH, size_t bW, X)( in Matrix!(bH,bW,X) mtr ) const
374         if( (bW==W||W==0||bW==0) && is(typeof(E(X.init))) )
375     {
376         static if( isDynamicWidth || mtr.isDynamicWidth )
377             enforce( mtr.width == width, "wrong width" );
378 
379         auto ret = Matrix!(0,W,E)(this);
380         auto last_height = height;
381         ret.resize( ret.height + mtr.height, width );
382         foreach( i; 0 .. mtr.height )
383             foreach( j; 0 .. width )
384                 ret.data[last_height+i][j] = E(mtr.data[i][j]);
385         return ret;
386     }
387 
388     ///
389     auto expandHeight(X...)( in X vals ) const
390         if( is(typeof(Matrix!(1,W,E)(vals))) )
391     { return expandHeight( Matrix!(1,W,E)(vals) ); }
392 
393     ///
394     auto expandWidth(size_t bH, size_t bW, X)( in Matrix!(bH,bW,X) mtr ) const
395         if( (bH==H||H==0||bH==0) && is(typeof(E(X.init))) )
396     {
397         static if( isDynamicHeight || mtr.isDynamicHeight )
398             enforce( mtr.height == height, "wrong height" );
399         auto ret = Matrix!(H,0,E)(this);
400         auto last_width = width;
401         ret.resize( height, ret.width + mtr.width );
402         foreach( i; 0 .. height )
403             foreach( j; 0 .. mtr.width )
404                 ret.data[i][last_width+j] = E(mtr.data[i][j]);
405         return ret;
406     }
407 
408     ///
409     auto expandWidth(X...)( in X vals ) const
410         if( is(typeof(Matrix!(H,1,E)(vals))) )
411     { return expandWidth( Matrix!(H,1,E)(vals) ); }
412 
413     ///
414     auto asArray() const @property
415     {
416         auto ret = new E[](width*height);
417         foreach( i; 0 .. height )
418             foreach( j; 0 .. width )
419                 ret[i*width+j] = data[i][j];
420         return ret;
421     }
422 
423     ///
424     auto sliceHeight( size_t start, size_t count=0 ) const
425     {
426         enforce( start < height );
427         count = count ? count : height - start;
428         auto ret = Matrix!(0,W,E)(this);
429         ret.resize( count, width );
430         foreach( i; 0 .. count )
431             ret.data[i][] = data[start+i][];
432         return ret;
433     }
434 
435     ///
436     auto sliceWidth( size_t start, size_t count=0 ) const
437     {
438         enforce( start < width );
439         count = count ? count : width - start;
440         auto ret = Matrix!(H,0,E)(this);
441         ret.resize( height, count );
442         foreach( i; 0 .. height )
443             foreach( j; 0 .. count )
444             ret.data[i][j] = data[i][start+j];
445         return ret;
446     }
447 
448     ///
449     auto opUnary(string op)() const
450         if( op == "-" && is( typeof( E.init * (-1) ) ) )
451     {
452         auto ret = selftype(this);
453         foreach( ref row; ret.data )
454             foreach( ref v; row )
455                 v = v * -1;
456         return ret;
457     }
458 
459     private void checkCompatible(size_t bH, size_t bW, X)( in Matrix!(bH,bW,X) mtr ) const
460     {
461         static if( isDynamicHeight || mtr.isDynamicHeight )
462             enforce( height == mtr.height, "wrong height" );
463         static if( isDynamicWidth || mtr.isDynamicWidth )
464             enforce( width == mtr.width, "wrong width" );
465     }
466 
467     ///
468     auto opBinary(string op, size_t bH, size_t bW, X)( in Matrix!(bH,bW,X) mtr ) const
469         if( (op=="+"||op=="-") && allowSomeOp(H,bH) && allowSomeOp(W,bW) )
470     {
471         checkCompatible( mtr );
472 
473         auto ret = selftype(this);
474 
475         foreach( i; 0 .. height )
476             foreach( j; 0 .. width )
477                 mixin( `ret[i][j] = E(data[i][j] ` ~ op ~ ` mtr[i][j]);` );
478         return ret;
479     }
480 
481     ///
482     auto opBinary(string op,X)( in X b ) const
483         if( (op!="+" && op!="-") && isValidOp!(op,E,X) )
484     {
485         auto ret = selftype(this);
486         foreach( i; 0 .. height )
487             foreach( j; 0 .. width )
488                 mixin( `ret[i][j] = E(data[i][j] ` ~ op ~ ` b);` );
489         return ret;
490     }
491 
492     ///
493     auto opBinary(string op,size_t bH, size_t bW, X)( in Matrix!(bH,bW,X) mtr ) const
494         if( op=="*" && allowSomeOp(W,bH) && isValidOp!("*",E,X) )
495     {
496         static if( isDynamic || mtr.isDynamic )
497             enforce( width == mtr.height, "incompatible sizes for mul" );
498         Matrix!(H,bW,E) ret;
499         static if( ret.isDynamic ) ret.resize(height,mtr.width);
500 
501         foreach( i; 0 .. height )
502             foreach( j; 0 .. mtr.width )
503             {
504                 ret[i][j] = E( data[i][0] * mtr.data[0][j] );
505                 foreach( k; 1 .. width )
506                     ret[i][j] = ret[i][j] + E( data[i][k] * mtr.data[k][j] );
507             }
508 
509         return ret;
510     }
511 
512     ///
513     ref typeof(this) opOpAssign(string op, E)( in E b )
514         if( mixin( `is( typeof( selftype.init ` ~ op ~ ` E.init ) == selftype )` ) )
515     { mixin( `return this = this ` ~ op ~ ` b;` ); }
516 
517     ///
518     bool opCast(E)() const if( is( E == bool ) )
519     { 
520         foreach( v; asArray ) if( !isFinite(v) ) return false;
521         return true;
522     }
523 
524     /// transponate
525     auto T() const @property
526     {
527         Matrix!(W,H,E) ret;
528         static if( isDynamic ) ret.resize(width,height);
529         foreach( i; 0 .. height)
530             foreach( j; 0 .. width )
531                 ret[j][i] = data[i][j];
532         return ret;
533     }
534 
535     ///
536     ref typeof(this) setCol(X...)( size_t no, in X vals )
537         if( is(typeof(flatData!E(vals))) )
538     {
539         enforce( no < width, "bad col index" );
540         auto buf = flatData!E(vals);
541         enforce( buf.length == height, "bad data length" );
542         foreach( i; 0 .. height )
543             data[i][no] = buf[i];
544         return this;
545     }
546 
547     ///
548     ref typeof(this) setRow(X...)( size_t no, in X vals )
549         if( is(typeof(flatData!E(vals))) )
550     {
551         enforce( no < height, "bad row index" );
552         auto buf = flatData!E(vals);
553         enforce( buf.length == width, "bad data length" );
554         data[no][] = buf[];
555         return this;
556     }
557 
558     ///
559     ref typeof(this) setRect(size_t bH, size_t bW, X)( size_t pos_row, size_t pos_col, in Matrix!(bH,bW,X) mtr )
560         if( is(typeof(E(X.init))) )
561     {
562         enforce( pos_row < height, "bad row index" );
563         enforce( pos_col < width, "bad col index" );
564         enforce( pos_row + mtr.height <= height, "bad height size" );
565         enforce( pos_col + mtr.width <= width, "bad width size" );
566 
567         foreach( i; 0 .. mtr.height )
568             foreach( j; 0 .. mtr.width )
569                 data[i+pos_row][j+pos_col] = E(mtr[i][j]);
570 
571         return this;
572     }
573 
574     ///
575     auto col( size_t no ) const
576     {
577         enforce( no < width, "bad col index" );
578         Matrix!(H,1,E) ret;
579         static if( ret.isDynamic )
580             ret.resize(height,1);
581         foreach( i; 0 .. height )
582             ret[i][0] = data[i][no];
583         return ret;
584     }
585 
586     ///
587     auto row( size_t no ) const
588     {
589         enforce( no < height, "bad row index" );
590         Matrix!(1,W,E) ret;
591         static if( ret.isDynamic )
592             ret.resize(1,width);
593         ret[0][] = data[no][];
594         return ret;
595     }
596 
597     ///
598     auto opBinary(string op,size_t K,X)( in Vector!(K,X) v ) const
599         if( op=="*" && allowSomeOp(W,K) && isValidOp!("*",E,X) && isValidOp!("+",E,E) )
600     {
601         static if( isDynamic || v.isDynamic )
602             enforce( width == v.length, "wrong vector length" );
603 
604         Vector!(H,E) ret;
605 
606         static if( ret.isDynamic )
607             ret.length = height;
608 
609         foreach( i; 0 .. height )
610         {
611             ret[i] = data[i][0] * v[0];
612 
613             foreach( j; 1 .. width )
614                 ret[i] = ret[i] + E(data[i][j] * v[j]);
615         }
616 
617         return ret;
618     }
619 
620     ///
621     auto opBinaryRight(string op,size_t K,X)( in Vector!(K,X) v ) const
622         if( op=="*" && isVector!(typeof(selftype.init.T * typeof(v).init)) )
623     { return this.T * v; }
624 
625     static if( W==4 && ( H==4 || H==3 ) && isFloatingPoint!E )
626     {
627         /++ transform vector, equals `( m * vec4( v, fc ) ).xyz`
628 
629             only:
630             W==4 && ( H==4 || H==3 ) && isFloatingPoint!E
631          +/
632         Vector!(3,E) tr(X)( in Vector!(3,X) v, E fc ) const
633         { return ( this * Vector!(4,E)( v, fc ) ).xyz; }
634 
635         /++ project vector
636          +
637          + equals:
638          + ---
639          + auto r = m * vec4( v, fc );
640          + return r.xyz / r.w;
641          + ---
642          + only:
643          +    W==4 && ( H==4 || H==3 ) && isFloatingPoint!E
644          +/
645         Vector!(3,E) prj(X)( in Vector!(3,X) v, E fc ) const
646         {
647             auto r = this * Vector!(4,E)( v, fc );
648             return r.xyz / r.w;
649         }
650        
651         @property
652         {
653             ///
654             Vector!(3,E) offset() const
655             { return Vector!(3,E)( cast(E[3])( col(3).data[0..3] ) ); }
656 
657             ///
658             Vector!(3,X) offset(X)( in Vector!(3,X) v )
659             { setCol(3, Vector!(4,E)( v, data[3][3] ) ); return v; }
660         }
661     }
662 
663     static private size_t[] getIndexesWithout(size_t max, in size_t[] arr)
664     {
665         size_t[] ret;
666         foreach( i; 0 .. max ) if( !canFind(arr,i) ) ret ~= i;
667         return ret;
668     }
669 
670     ///
671     auto subWithout( size_t[] without_rows=[], size_t[] without_cols=[] ) const
672     {
673         auto with_rows = getIndexesWithout(height,without_rows);
674         auto with_cols = getIndexesWithout(width,without_cols);
675 
676         return sub( with_rows,with_cols );
677     }
678 
679     /// get sub matrix
680     auto sub( size_t[] with_rows, size_t[] with_cols ) const
681     {
682         auto wrows = array( uniq(with_rows) );
683         auto wcols = array( uniq(with_cols) );
684 
685         enforce( all!(a=>a<height)( wrows ), "bad row index" );
686         enforce( all!(a=>a<width)( wcols ), "bad col index" );
687 
688         Matrix!(0,0,E) ret;
689         ret.resize( wrows.length, wcols.length );
690 
691         foreach( i, orig_i; wrows )
692             foreach( j, orig_j; wcols )
693                 ret[i][j] = data[orig_i][orig_j];
694 
695         return ret;
696     }
697 
698     static if( isFloatingPoint!E && ( isDynamicAll ||
699                 ( isStaticHeightOnly && H > 1 ) ||
700                 ( isStaticWidthOnly && W > 1 ) ||
701                 ( isStatic && W == H ) ) )
702     {
703         ///
704         E cofactor( size_t i, size_t j ) const
705         { return subWithout([i],[j]).det * coef(i,j); }
706 
707         private static nothrow @trusted auto coef( size_t i, size_t j )
708         { return ((i+j)%2?-1:1); }
709         
710         ///
711         auto det() const @property 
712         {
713             static if( isDynamic )
714                 enforce( width == height, "not square matrix" );
715 
716             static if( isDynamic )
717             {
718                 if( width == 1 )
719                     return data[0][0];
720                 else if( width == 2 )
721                     return data[0][0] * data[1][1] -
722                            data[0][1] * data[1][0];
723                 else return classicDet;
724             }
725             else
726             {
727                 static if( W == 1 )
728                     return data[0][0];
729                 else static if( W == 2 )
730                     return data[0][0] * data[1][1] -
731                            data[0][1] * data[1][0];
732                 else return classicDet;
733             }
734         }
735 
736         private @property classicDet() const
737         {
738             auto i = 0UL; // TODO: find max zeros line
739             auto ret = data[i][0] * cofactor(i,0);
740 
741             foreach( j; 1 .. width )
742                 ret = ret + data[i][j] * cofactor(i,j);
743 
744             return ret;
745         }
746 
747         ///
748         auto inv() const @property
749         {
750             static if( isDynamic )
751                 enforce( width == height, "not square matrix" );
752             else static assert( W==H, "not square matrix" );
753 
754             selftype buf;
755 
756             static if( isDynamic )
757                 buf.resize(height,width);
758 
759             foreach( i; 0 .. height )
760                 foreach( j; 0 .. width )
761                     buf[i][j] = cofactor(i,j);
762 
763             auto i = 0UL; // TODO: find max zeros line
764             auto d = data[i][0] * buf[i][0];
765 
766             foreach( j; 1 .. width )
767                 d = d + data[i][j] * buf[i][j];
768 
769             return buf.T / d;
770         }
771 
772         static if( (isStaticHeightOnly && H==4) ||
773                    (isStaticWidthOnly && W==4) ||
774                    (isStatic && H==W && H==4) ||
775                    isDynamicAll )
776         {
777             /++ only for transform matrix +/
778             @property auto speedTransformInv() const
779             {
780                 static if( isDynamic )
781                 {
782                     enforce( width == height, "not square matrix" );
783                     enforce( width == 4, "matrix must be 4x4" );
784                 }
785 
786                 selftype ret;
787                 static if( isDynamic )
788                     ret.resize( height, width );
789 
790                 foreach( i; 0 .. 3 )
791                     foreach( j; 0 .. 3 )
792                         ret[i][j] = this[j][i];
793 
794                 auto a22k = 1.0 / this[3][3];
795 
796                 ret[0][3] = -( ret[0][0] * this[0][3] + ret[0][1] * this[1][3] + ret[0][2] * this[2][3] ) * a22k;
797                 ret[1][3] = -( ret[1][0] * this[0][3] + ret[1][1] * this[1][3] + ret[1][2] * this[2][3] ) * a22k;
798                 ret[2][3] = -( ret[2][0] * this[0][3] + ret[2][1] * this[1][3] + ret[2][2] * this[2][3] ) * a22k;
799 
800                 ret[3][0] = -( this[3][0] * ret[0][0] + this[3][1] * ret[1][0] + this[3][2] * ret[2][0] ) * a22k;
801                 ret[3][1] = -( this[3][0] * ret[0][1] + this[3][1] * ret[1][1] + this[3][2] * ret[2][1] ) * a22k;
802                 ret[3][2] = -( this[3][0] * ret[0][2] + this[3][1] * ret[1][2] + this[3][2] * ret[2][2] ) * a22k;
803                 
804                 ret[3][3] = a22k * ( 1.0 - ( this[3][0] * ret[0][3] + this[3][1] * ret[1][3] + this[3][2] * ret[2][3] ) );
805 
806                 return ret;
807             }
808         }
809 
810         ///
811         auto rowReduceInv() const @property
812         {
813             static if( isDynamic )
814                 enforce( width == height, "not square matrix" );
815 
816             auto ln = height;
817 
818             auto orig = selftype(this);
819             selftype invt;
820             static if( isDynamic )
821             {
822                 invt.resize(ln,ln);
823                 foreach( i; 0 .. ln )
824                     foreach( j; 0 .. ln )
825                         invt[i][j] = E(i==j);
826             }
827 
828             foreach( r; 0 .. ln-1 )
829             {
830                 auto k = E(1) / orig[r][r];
831                 foreach( c; 0 .. ln )
832                 {
833                     orig[r][c] *= k;
834                     invt[r][c] *= k;
835                 }
836                 foreach( rr; r+1 .. ln )
837                 {
838                     auto v = orig[rr][r];
839                     foreach( c; 0 .. ln )
840                     {
841                         orig[rr][c] -= orig[r][c] * v;
842                         invt[rr][c] -= invt[r][c] * v;
843                     }
844                 }
845             }
846 
847             foreach_reverse( r; 1 .. ln )
848             {
849                 auto k = E(1) / orig[r][r];
850                 foreach( c; 0 .. ln )
851                 {
852                     orig[r][c] *= k;
853                     invt[r][c] *= k;
854                 }
855                 foreach_reverse( rr; 0 .. r )
856                 {
857                     auto v = orig[rr][r];
858                     foreach( c; 0 .. ln )
859                     {
860                         orig[rr][c] -= orig[r][c] * v;
861                         invt[rr][c] -= invt[r][c] * v;
862                     }
863                 }
864             }
865 
866             return invt;
867         }
868     }
869 }
870 
871 alias Matrix2(T)   = Matrix!(2,2,T); ///
872 alias Matrix2x3(T) = Matrix!(2,3,T); ///
873 alias Matrix2x4(T) = Matrix!(2,4,T); ///
874 alias Matrix2xD(T) = Matrix!(2,0,T); ///
875 alias Matrix3x2(T) = Matrix!(3,2,T); ///
876 alias Matrix3(T)   = Matrix!(3,3,T); ///
877 alias Matrix3x4(T) = Matrix!(3,4,T); ///
878 alias Matrix3xD(T) = Matrix!(3,0,T); ///
879 alias Matrix4x2(T) = Matrix!(4,2,T); ///
880 alias Matrix4x3(T) = Matrix!(4,3,T); ///
881 alias Matrix4(T)   = Matrix!(4,4,T); ///
882 alias Matrix4xD(T) = Matrix!(4,0,T); ///
883 alias MatrixDx2(T) = Matrix!(0,2,T); ///
884 alias MatrixDx3(T) = Matrix!(0,3,T); ///
885 alias MatrixDx4(T) = Matrix!(0,4,T); ///
886 alias MatrixDxD(T) = Matrix!(0,0,T); ///
887 
888 alias Matrix!(2,2,float) mat2;   ///
889 alias Matrix!(2,3,float) mat2x3; ///
890 alias Matrix!(2,4,float) mat2x4; ///
891 alias Matrix!(2,0,float) mat2xD; ///
892 alias Matrix!(3,2,float) mat3x2; ///
893 alias Matrix!(3,3,float) mat3;   ///
894 alias Matrix!(3,4,float) mat3x4; ///
895 alias Matrix!(3,0,float) mat3xD; ///
896 alias Matrix!(4,2,float) mat4x2; ///
897 alias Matrix!(4,3,float) mat4x3; ///
898 alias Matrix!(4,4,float) mat4;   ///
899 alias Matrix!(4,0,float) mat4xD; ///
900 alias Matrix!(0,2,float) matDx2; ///
901 alias Matrix!(0,3,float) matDx3; ///
902 alias Matrix!(0,4,float) matDx4; ///
903 alias Matrix!(0,0,float) matD;   ///
904 
905 alias Matrix!(2,2,double) dmat2;   ///
906 alias Matrix!(2,3,double) dmat2x3; ///
907 alias Matrix!(2,4,double) dmat2x4; ///
908 alias Matrix!(2,0,double) dmat2xD; ///
909 alias Matrix!(3,2,double) dmat3x2; ///
910 alias Matrix!(3,3,double) dmat3;   ///
911 alias Matrix!(3,4,double) dmat3x4; ///
912 alias Matrix!(3,0,double) dmat3xD; ///
913 alias Matrix!(4,2,double) dmat4x2; ///
914 alias Matrix!(4,3,double) dmat4x3; ///
915 alias Matrix!(4,4,double) dmat4;   ///
916 alias Matrix!(4,0,double) dmat4xD; ///
917 alias Matrix!(0,2,double) dmatDx2; ///
918 alias Matrix!(0,3,double) dmatDx3; ///
919 alias Matrix!(0,4,double) dmatDx4; ///
920 alias Matrix!(0,0,double) dmatD;   ///
921 
922 alias Matrix!(2,2,real) rmat2;   ///
923 alias Matrix!(2,3,real) rmat2x3; ///
924 alias Matrix!(2,4,real) rmat2x4; ///
925 alias Matrix!(2,0,real) rmat2xD; ///
926 alias Matrix!(3,2,real) rmat3x2; ///
927 alias Matrix!(3,3,real) rmat3;   ///
928 alias Matrix!(3,4,real) rmat3x4; ///
929 alias Matrix!(3,0,real) rmat3xD; ///
930 alias Matrix!(4,2,real) rmat4x2; ///
931 alias Matrix!(4,3,real) rmat4x3; ///
932 alias Matrix!(4,4,real) rmat4;   ///
933 alias Matrix!(4,0,real) rmat4xD; ///
934 alias Matrix!(0,2,real) rmatDx2; ///
935 alias Matrix!(0,3,real) rmatDx3; ///
936 alias Matrix!(0,4,real) rmatDx4; ///
937 alias Matrix!(0,0,real) rmatD;   ///
938 
939 unittest
940 {
941     static assert( Matrix!(3,3,float).sizeof == float.sizeof * 9 );
942     static assert( Matrix!(3,3,float).isStatic );
943     static assert( Matrix!(0,3,float).isDynamic );
944     static assert( Matrix!(0,3,float).isDynamicHeight );
945     static assert( Matrix!(0,3,float).isStaticWidth );
946     static assert( Matrix!(3,0,float).isDynamic );
947     static assert( Matrix!(3,0,float).isDynamicWidth );
948     static assert( Matrix!(3,0,float).isStaticHeight );
949     static assert( Matrix!(0,0,float).isDynamic );
950     static assert( Matrix!(0,0,float).isDynamicHeight );
951     static assert( Matrix!(0,0,float).isDynamicWidth );
952 }
953 
954 unittest
955 {
956     static assert( isStaticMatrix!(mat3) );
957     static assert( !isStaticMatrix!(matD) );
958     static assert( !isStaticMatrix!(int[]) );
959 }
960 
961 unittest
962 {
963     assert( eq( mat3.init, [[1,0,0],[0,1,0],[0,0,1]] ) );
964     assert( eq( mat2.init, [[1,0],[0,1]] ) );
965     assert( eq( mat2x3.init, [[0,0,0],[0,0,0]] ) );
966 }
967 
968 unittest
969 {
970     auto a = Matrix!(3,2,double)( 36, 0, 3, 3, 0, 0 );
971     auto b = Matrix!(3,2,double)( [ 36, 0, 3, 3, 0, 0 ] );
972     assert( eq( a.asArray, b.asArray ) );
973     assert( eq( a.asArray, [ 36, 0, 3, 3, 0, 0 ] ) );
974 }
975 
976 unittest
977 {
978     auto a = mat3( 1,2,3,4,5,6,7,8,9 );
979     assert( a.height == 3 );
980     assert( a.width == 3 );
981     assert( eq( a, [[1,2,3],[4,5,6],[7,8,9]] ) );
982     static assert( !__traits(compiles,a.resize(1,1)) );
983     a.fill(1);
984     assert( eq( a, [[1,1,1], [1,1,1], [1,1,1]] ) );
985     a[0][0] *= 3;
986     assert( eq( a, [[3,1,1], [1,1,1], [1,1,1]] ) );
987     a[1][0] += 3;
988     assert( eq( a, [[3,1,1], [4,1,1], [1,1,1]] ) );
989 }
990 
991 ///
992 unittest
993 {
994     static struct Test
995     {
996         union
997         {
998             mat3 um;
999             float[9] uf;
1000         }
1001     }
1002 
1003     Test tt;
1004 
1005     foreach( i, ref v; tt.uf ) v = i+1;
1006     assert( eq( tt.um, [[1,2,3],[4,5,6],[7,8,9]] ) );
1007 }
1008 
1009 ///
1010 unittest
1011 {
1012     auto a = matDx3( 1,2,3,4,5,6,7,8,9 );
1013     assert( mustExcept( {matDx3(1,2,3,4);} ) );
1014     assert( a.height == 3 );
1015     assert( a.width == 3 );
1016     assert( eq( a, [[1,2,3],[4,5,6],[7,8,9]] ) );
1017     assert( mustExcept({ a.resize(2,2); }) );
1018     a.resize(2,3);
1019     assert( eq( a, [[1,2,3],[4,5,6]] ) );
1020     assert( a.height == 2 );
1021     a.fill(1);
1022     assert( eq( a, [[1,1,1],[1,1,1]] ) );
1023     a.resize(0,3);
1024     assert( a.width == 3 );
1025     a.resize(2,3);
1026     a.fill(1);
1027 
1028     auto b = a;
1029     assert( eq( b, [[1,1,1],[1,1,1]] ) );
1030 }
1031 
1032 unittest
1033 {
1034     auto a = mat3xD( 1,2,3,4,5,6,7,8,9 );
1035     assert( mustExcept( {mat3xD(1,2,3,4);} ) );
1036     assert( a.height == 3 );
1037     assert( a.width == 3 );
1038     assert( eq( a, [[1,2,3],[4,5,6],[7,8,9]] ) );
1039     assert( mustExcept({ a.resize(2,2); }) );
1040     a.resize(3,2);
1041     assert( eq( a, [[1,2],[4,5],[7,8]] ) );
1042     assert( a.width == 2 );
1043     a.fill(1);
1044     assert( eq( a, [[1,1],[1,1],[1,1]] ) );
1045 
1046     auto b = a;
1047     assert( eq( b, [[1,1],[1,1],[1,1]] ) );
1048 }
1049 
1050 unittest
1051 {
1052     auto a = matD( 3,3, 1,2,3,4,5,6,7,8,9 );
1053     assert( mustExcept( {matD(1,2,3,4,5);} ) );
1054     assert( a.height == 3 );
1055     assert( a.width == 3 );
1056     assert( eq( a, [[1,2,3],[4,5,6],[7,8,9]] ) );
1057     a.resize(2,2);
1058     assert( eq( a, [[1,2],[4,5]] ) );
1059     auto b = matD(2,2);
1060     b.fill(1);
1061     assert( eq( b, [[1,1],[1,1]] ) );
1062     b = a;
1063     assert( eq( a, b ) );
1064     auto c = matD( Matrix!(2,4,float)(1,2,3,4,5,6,7,8) );
1065     assert( eq( c, [[1,2,3,4],[5,6,7,8]] ) );
1066     assert( c.height == 2 );
1067     assert( c.width == 4 );
1068     assert( b.height == 2 );
1069     assert( b.width == 2 );
1070     b = c;
1071     assert( b.height == 2 );
1072     assert( b.width == 4 );
1073     assert( eq( b, c ) );
1074     b[0][0] = 666;
1075     assert( !eq( b, c ) );
1076 }
1077 
1078 unittest
1079 {
1080     auto a = mat3xD( 1,2,3,4,5,6,7,8,9 );
1081     matD b;
1082     assert( b.height == 0 );
1083     assert( b.width == 0 );
1084     b = a;
1085     assert( b.height == 3 );
1086     assert( b.width == 3 );
1087     assert( eq( a, b ) );
1088     a[0][0] = 22;
1089     assert( !eq( a, b ) );
1090     a = b;
1091     assert( eq( a, b ) );
1092     b.height = 4;
1093     assert(  mustExcept({ a = b; }) );
1094     assert( !mustExcept({ b = a; }) );
1095 }
1096 
1097 ///
1098 unittest
1099 {
1100     auto a = mat3( 1,2,3,4,5,6,7,8,9 );
1101     auto b = matD( 3,3, 1,2,3,4,5,6,7,8,9 );
1102     auto c = mat3xD( 1,2,3,4,5,6,7,8,9 );
1103     auto d = matDx3( 1,2,3,4,5,6,7,8,9 );
1104 
1105     auto eha = a.expandHeight( 8,8,8 );
1106     auto ehb = b.expandHeight( 8,8,8 );
1107     auto ehc = c.expandHeight( 8,8,8 );
1108     auto ehd = d.expandHeight( 8,8,8 );
1109 
1110     assert( eq( eha, [[1,2,3],[4,5,6],[7,8,9],[8,8,8]] ));
1111     assert( eha.height == 4 );
1112     assert( ehb.height == 4 );
1113     assert( ehc.height == 4 );
1114     assert( ehd.height == 4 );
1115     assert( eq( eha, ehb ) );
1116     assert( eq( eha, ehc ) );
1117     assert( eq( eha, ehd ) );
1118 
1119     static assert( is(typeof(eha) == Matrix!(0,3,float)) );
1120     static assert( is(typeof(ehd) == Matrix!(0,3,float)) );
1121 
1122     static assert( is(typeof(ehb) == Matrix!(0,0,float)) );
1123     static assert( is(typeof(ehc) == Matrix!(0,0,float)) );
1124 
1125     auto ewa = a.expandWidth( 8,8,8 );
1126     auto ewb = b.expandWidth( 8,8,8 );
1127     auto ewc = c.expandWidth( 8,8,8 );
1128     auto ewd = d.expandWidth( 8,8,8 );
1129 
1130     assert( eq( ewa, [[1,2,3,8],[4,5,6,8],[7,8,9,8]] ));
1131     assert( ewa.width == 4 );
1132     assert( ewb.width == 4 );
1133     assert( ewc.width == 4 );
1134     assert( ewd.width == 4 );
1135     assert( eq( ewa, ewb ) );
1136     assert( eq( ewa, ewc ) );
1137     assert( eq( ewa, ewd ) );
1138 
1139     static assert( is(typeof(ewa) == Matrix!(3,0,float)) );
1140     static assert( is(typeof(ewc) == Matrix!(3,0,float)) );
1141 
1142     static assert( is(typeof(ewb) == Matrix!(0,0,float)) );
1143     static assert( is(typeof(ewd) == Matrix!(0,0,float)) );
1144 
1145     auto aa = a.expandHeight(a);
1146     assert( eq( aa, [[1,2,3],[4,5,6],[7,8,9],[1,2,3],[4,5,6],[7,8,9]] ));
1147     assert( aa.height == 6 );
1148     static assert( is(typeof(aa) == Matrix!(0,3,float)) );
1149 }
1150 
1151 unittest
1152 {
1153     auto a = mat3( 1,2,3,4,5,6,7,8,9 );
1154     assert( a.asArray == [1.0f,2,3,4,5,6,7,8,9] );
1155 }
1156 
1157 ///
1158 unittest
1159 {
1160     auto a = matD(4,4,0).fillDiag(1);
1161     assert( eq( a, mat4() ) );
1162     assert( eq( a.inv, a ) );
1163 }
1164 
1165 ///
1166 unittest
1167 {
1168     auto a = mat4x2( 1,2,3,4,5,6,7,8 );
1169     auto b = mat4( 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16 );
1170     auto c = a.T * b * a;
1171     static assert( c.height == 2 && c.width == 2 );
1172 }
1173 
1174 unittest
1175 {
1176     matD a;
1177     matD b;
1178     a.resize( 10, 4 );
1179     b.resize( 10, 10 );
1180     auto c = a.T * b * a;
1181     assert( c.height == 4 && c.width == 4 );
1182 }
1183 
1184 ///
1185 unittest
1186 {
1187     auto a = mat3.diag(1);
1188     assert( eq(a,[[1,0,0],[0,1,0],[0,0,1]]) );
1189     auto b = mat3xD.diag(1,2);
1190     assert( eq(b,[[1,0,0],[0,2,0],[0,0,1]]) );
1191     auto c = mat3xD.diag(1,2,3);
1192     assert( eq(c,[[1,0,0],[0,2,0],[0,0,3]]) );
1193     static assert( !__traits(compiles,matD.diag(1)) );
1194     auto d = matD(3,3).fill(0).fillDiag(1);
1195     assert( eq(d,[[1,0,0],[0,1,0],[0,0,1]]) );
1196 }
1197 
1198 ///
1199 unittest
1200 {
1201     auto a = mat3( 1,2,3,4,5,6,7,8,9 );
1202     auto sha = a.sliceHeight(1);
1203     assert( eq(sha,[[4,5,6],[7,8,9]]) );
1204     auto swsha = sha.sliceWidth(0,1);
1205     assert( eq(swsha,[[4],[7]]) );
1206 }
1207 
1208 ///
1209 unittest
1210 {
1211     auto a = mat3.diag(1);
1212     assert( eq( -a,[[-1,0,0],[0,-1,0],[0,0,-1]]) );
1213 }
1214 
1215 unittest
1216 {
1217     auto a = mat3.diag(1);
1218     auto b = a*3-a;
1219     assert( eq( b,a*2 ) );
1220     b /= 2;
1221     assert( eq( b,a ) );
1222 }
1223 
1224 ///
1225 unittest
1226 {
1227     auto a = rmat3( 3,2,2, 1,3,1, 5,3,4 );
1228     auto ainv = rmat3xD( 9,-2,-4,  1,2,-1, -12,1,7 ) / 5;
1229 
1230     auto b = a * ainv;
1231     static assert( is( typeof(b) == Matrix!(3,0,real) ) );
1232     assert( eq( b, mat3.diag(1) ) );
1233 
1234     static assert( !__traits(compiles,a *= ainv ) );
1235     a *= rmat3(ainv);
1236     assert( eq( a, mat3.diag(1) ) );
1237 }
1238 
1239 ///
1240 unittest
1241 {
1242     alias Matrix!(4,5,float) mat4x5;
1243 
1244     auto a = mat2x4.init * mat4xD.init;
1245     assert( mustExcept({ mat4xD.init * mat2x4.init; }) );
1246     auto b = mat2x4.init * mat4x5.init;
1247 
1248     static assert( is( typeof(a) == Matrix!(2,0,float) ) );
1249     static assert( is( typeof(b) == Matrix!(2,5,float) ) );
1250     static assert( is( typeof(mat4xD.init * mat2x4.init) == Matrix!(4,4,float) ) );
1251 }
1252 
1253 unittest
1254 {
1255 
1256     auto a = mat2x3( 1,2,3, 4,5,6 );
1257     auto b = mat4xD( 1,2,3,4 );
1258 
1259     auto aT = a.T;
1260     auto bT = b.T;
1261 
1262     static assert( is( typeof(aT) == Matrix!(3,2,float) ) );
1263     static assert( is( typeof(bT) == Matrix!(0,4,float) ) );
1264 
1265     assert( bT.height == b.width );
1266 }
1267 
1268 unittest
1269 {
1270     auto a = mat3.diag(1);
1271 
1272     a.setRow(2,4,5,6);
1273     assert( eq( a, [[1,0,0],[0,1,0],[4,5,6]] ) );
1274     a.setCol(1,8,4,2);
1275     assert( eq( a, [[1,8,0],[0,4,0],[4,2,6]] ) );
1276 
1277     assert( eq( a.row(0), [[1,8,0]] ) );
1278     assert( eq( a.col(0), [[1],[0],[4]] ) );
1279 }
1280 
1281 ///
1282 unittest
1283 {
1284     auto b = mat3.diag(2) * vec3(2,3,4);
1285     assert( is( typeof(b) == vec3 ) );
1286     assert( eq( b, [4,6,8] ) );
1287 
1288     auto c = vec3(1,1,1) * mat3.diag(1,2,3);
1289     assert( is( typeof(c) == vec3 ) );
1290     assert( eq( c, [1,2,3] ) );
1291 }
1292 
1293 ///
1294 unittest
1295 {
1296     auto mtr = matD(2,3).fill( 1,2,3, 
1297                                3,4,7 );
1298 
1299     auto a = mtr * vec3(1,1,1);
1300 
1301     assert( is( typeof(a) == Vector!(0,float) ) );
1302     assert( a.length == 2 );
1303     assert( eq( a, [ 6, 14 ] ) );
1304 
1305     auto b = vec3(1,1,1) * mtr.T;
1306     assert( eq( a, b ) );
1307 }
1308 
1309 ///
1310 unittest
1311 {
1312     void check(E)( E mtr ) if( isMatrix!E )
1313     {
1314         mtr.fill( 1,2,3,4,
1315                   5,6,7,8,
1316                   9,10,11,12,
1317                   13,14,15,16 );
1318         auto sm = mtr.subWithout( [0,3,3,3], [1,2] );
1319         assert( is( typeof(sm) == matD ) );
1320         assert( sm.width == 2 );
1321         assert( sm.height == 2 );
1322         assert( eq( sm, [ [5,8], [9,12] ] ) );
1323         assert( mustExcept({ mtr.sub( [0,4], [1,2] ); }) );
1324         auto sm2 = mtr.subWithout( [], [1,2] );
1325         assert( sm2.width == 2 );
1326         assert( sm2.height == 4 );
1327         assert( eq( sm2, [ [1,4],[5,8],[9,12],[13,16] ] ) );
1328     }
1329 
1330     check( matD(4,4) );
1331     check( mat4() );
1332 }
1333 
1334 unittest
1335 {
1336     assert( eq( matD(4,4,0).fillDiag(1,2,3,4).det, 24 ) );
1337 }
1338 
1339 unittest
1340 {
1341     auto mtr = rmatD(4,4).fill( 1,2,5,2,
1342                                 5,6,1,4,
1343                                 9,1,3,0,
1344                                 9,2,4,2 );
1345     auto xx = mtr * mtr.inv;
1346     assert( eq( xx, matD(4,4,0).fillDiag(1) ) );
1347 }
1348 
1349 unittest
1350 {
1351     auto mtr = matD(4,4).fill( 0,1,0,2,
1352                                1,0,0,4,
1353                                0,0,1,1,
1354                                0,0,0,1 );
1355     auto vv = vec4(4,2,1,1);
1356     auto rv = mtr.speedTransformInv * (mtr*vv);
1357     assert( eq( rv, vv ) );
1358 }
1359 
1360 unittest
1361 {
1362     auto mtr = matD(4,4).fillDiag(1);
1363     auto vv = vec4(4,2,1,1);
1364     auto rv = vv * mtr;
1365     auto vr = mtr.T * vv * mtr;
1366 }
1367 
1368 unittest
1369 {
1370     auto mtr = rmat4.diag(1);
1371     auto vv = vec4(4,2,1,1);
1372     auto rv = vv * mtr;
1373     auto vr = mtr.T * vv * mtr;
1374 
1375     auto xx = mtr * mtr.inv;
1376     assert( eq( xx, mat4.diag(1) ) );
1377 }
1378 
1379 unittest
1380 {
1381     auto mtr = rmat4( 1,2,5,2,
1382                       5,6,1,4,
1383                       9,1,3,0,
1384                       9,2,4,2 );
1385 
1386     auto xx = mtr * mtr.rowReduceInv;
1387     assert( eq( xx, mat4() ) );
1388 }
1389 
1390 ///
1391 unittest
1392 {
1393     auto mtr = mat4().setRect(0,0,mat3.diag(1)).setCol(3,1,2,3,4);
1394     assert( eq( mtr, [[1,0,0,1],
1395                       [0,1,0,2],
1396                       [0,0,1,3],
1397                       [0,0,0,4]] ) );
1398 }
1399 
1400 unittest
1401 {
1402     auto stm = mat4();
1403     assert( stm );
1404     auto dnm = matD();
1405     assert( dnm );
1406 }
1407 
1408 unittest
1409 {
1410     auto t = mat4.diag(1,2,3,4);
1411     assertEq( t.offset, vec3(0) );
1412     t.setCol( 3, vec4(3,2,1,1) );
1413     assertEq( t.offset, vec3(3,2,1) );
1414     t.offset = vec3(1,2,3);
1415     assertEq( t.offset, vec3(1,2,3) );
1416 }
1417 
1418 ///
1419 auto quatToMatrix(E)( in Quaterni!E iq )
1420 {
1421     auto q = iq / iq.len2;
1422 
1423     E wx, wy, wz, xx, yy, yz, xy, xz, zz, x2, y2, z2;
1424 
1425     x2 = q.i + q.i;
1426     y2 = q.j + q.j;
1427     z2 = q.k + q.k;
1428     xx = q.i * x2;   xy = q.i * y2;   xz = q.i * z2;
1429     yy = q.j * y2;   yz = q.j * z2;   zz = q.k * z2;
1430     wx = q.a * x2;   wy = q.a * y2;   wz = q.a * z2;
1431 
1432     return Matrix!(3,3,E)( 1.0-(yy+zz),  xy-wz,        xz+wy,
1433                            xy+wz,        1.0-(xx+zz),  yz-wx,
1434                            xz-wy,        yz+wx,        1.0-(xx+yy) );
1435 }
1436 
1437 ///
1438 unittest
1439 {
1440     auto q = rquat.fromAngle( PI_2, vec3(0,0,1) );
1441 
1442     auto m = quatToMatrix(q);
1443     auto r = mat3(0,-1, 0, 1, 0, 0, 0, 0, 1);
1444     assert( eq_approx( m, r, 1e-7 ) );
1445 }
1446 
1447 ///
1448 auto quatAndPosToMatrix(E,V)( in Quaterni!E iq, in V pos )
1449     if( isSpecVector!(3,E,V) )
1450 {
1451     auto q = iq / iq.len2;
1452 
1453     E wx, wy, wz, xx, yy, yz, xy, xz, zz, x2, y2, z2;
1454 
1455     x2 = q.i + q.i;
1456     y2 = q.j + q.j;
1457     z2 = q.k + q.k;
1458     xx = q.i * x2;   xy = q.i * y2;   xz = q.i * z2;
1459     yy = q.j * y2;   yz = q.j * z2;   zz = q.k * z2;
1460     wx = q.a * x2;   wy = q.a * y2;   wz = q.a * z2;
1461 
1462     return Matrix!(4,4,E)( 1.0-(yy+zz),  xy-wz,        xz+wy,       pos.x,
1463                            xy+wz,        1.0-(xx+zz),  yz-wx,       pos.y,
1464                            xz-wy,        yz+wx,        1.0-(xx+yy), pos.z,
1465                            0,            0,            0,           1     );
1466 }
1467 
1468 ///
1469 unittest
1470 {
1471     auto q = rquat.fromAngle( PI_2, vec3(0,0,1) );
1472     auto m = quatAndPosToMatrix(q, vec3(1,2,3) );
1473     assert( eq_approx( m, [[0,-1,0,1],
1474                            [1, 0,0,2],
1475                            [0, 0,1,3],
1476                            [0, 0,0,1]], 1e-7 ) );
1477 
1478 }
Page generated by adrdox