WebSVN – gelsvn – Diff – /trunk/src/CGLA/ArithSqMat4x4Float.cpp


/* ----------------------------------------------------------------------- *
 * This file is part of GEL, http://www.imm.dtu.dk/GEL
 * Copyright (C) the authors and DTU Informatics
 * For license and list of authors, see ../../doc/intro.pdf
 * ----------------------------------------------------------------------- */
 
#include "ArithSqMat4x4Float.h"
#include "Mat4x4f.h"
#include "Mat4x4d.h"
 
namespace CGLA {
	
namespace
{
	
	/* Aux func. computes 3x3 determinants. */
	template<class T>
	inline T d3x3f( T a0, T a1, T a2, 
									T b0, T b1, T b2, 
									T c0, T c1, T c2 )
	{
		return a0*b1*c2 +
			a1*b2*c0 +
			a2*b0*c1 -
			a2*b1*c0 -
			a0*b2*c1 -
			a1*b0*c2 ;
	}
}
 
 
 
	template<class V, class M>
	M invert_affine(const ArithSqMat4x4Float<V,M>& this_mat)
	{
		//   The following com[3]e has been copied from a gem in Graphics Gems II by
		//   Kevin Wu.                      
        //   From the EULA: "Using the code is permitted in any program, product, or 
        //   library, non-commercial or commercial. Giving credit is not required, though is a nice gesture"
		//   The function is very fast, but it can only invert affine matrices. An
		//   exception NotAffine is thrown if the matrix is not affine, and another
		//   exception Singular is thrown if the matrix is singular.
		
		typedef typename M::ScalarType ScalarType;
 
		M new_mat;
		ScalarType det_1;
		ScalarType pos, neg, temp;
      
 
		if (!(is_tiny(this_mat[3][0]) && 
					is_tiny(this_mat[3][1]) && 
					is_tiny(this_mat[3][2]) && 
					is_tiny(this_mat[3][3]-1.0)))
			throw(Mat4x4fNotAffine("Can only invert affine matrices"));
    
#define ACCUMULATE if (temp >= 0.0) pos += temp; else neg += temp
  
		/*
		 * Calculate the determinant of submatrix A and determine if the
		 * the matrix is singular as limited by the float precision
		 * floating-point this_mat representation.
		 */
  
		pos = neg = 0.0;
		temp =  this_mat[0][0] * this_mat[1][1] * this_mat[2][2];
		ACCUMULATE;
		temp =  this_mat[1][0] * this_mat[2][1] * this_mat[0][2];
		ACCUMULATE;
		temp =  this_mat[2][0] * this_mat[0][1] * this_mat[1][2];
		ACCUMULATE;
		temp = -this_mat[2][0] * this_mat[1][1] * this_mat[0][2];
		ACCUMULATE;
		temp = -this_mat[1][0] * this_mat[0][1] * this_mat[2][2];
		ACCUMULATE;
		temp = -this_mat[0][0] * this_mat[2][1] * this_mat[1][2];
		ACCUMULATE;
		det_1 = pos + neg;
  
		/* Is the submatrix A singular? */
		if ((det_1 == 0.0) || (fabs(det_1 / (pos - neg)) < MINUTE)) 
			{
				/* Mat4x4f M has no inverse */
				throw(Mat4x4fSingular("Tried to invert Singular matrix"));
			}
 
		else {
 
			/* Calculate inverse(A) = adj(A) / det(A) */
			det_1 = 1.0 / det_1;
			new_mat[0][0] =   ( this_mat[1][1] * this_mat[2][2] -
													this_mat[2][1] * this_mat[1][2] )
				* det_1;
			new_mat[0][1] = - ( this_mat[0][1] * this_mat[2][2] -
													this_mat[2][1] * this_mat[0][2] )
				* det_1;
			new_mat[0][2] =   ( this_mat[0][1] * this_mat[1][2] -
													this_mat[1][1] * this_mat[0][2] )
				* det_1;
			new_mat[1][0] = - ( this_mat[1][0] * this_mat[2][2] -
													this_mat[2][0] * this_mat[1][2] )
				* det_1;
			new_mat[1][1] =   ( this_mat[0][0] * this_mat[2][2] -
													this_mat[2][0] * this_mat[0][2] )
				* det_1;
			new_mat[1][2] = - ( this_mat[0][0] * this_mat[1][2] -
													this_mat[1][0] * this_mat[0][2] )
				* det_1;
			new_mat[2][0] =   ( this_mat[1][0] * this_mat[2][1] -
													this_mat[2][0] * this_mat[1][1] )
				* det_1;
			new_mat[2][1] = - ( this_mat[0][0] * this_mat[2][1] -
													this_mat[2][0] * this_mat[0][1] )
				* det_1;
			new_mat[2][2] =   ( this_mat[0][0] * this_mat[1][1] -
													this_mat[1][0] * this_mat[0][1] )
				* det_1;
 
			/* Calculate -C * inverse(A) */
			new_mat[0][3] = - ( this_mat[0][3] * new_mat[0][0] +
													this_mat[1][3] * new_mat[0][1] +
													this_mat[2][3] * new_mat[0][2] );
			new_mat[1][3] = - ( this_mat[0][3] * new_mat[1][0] +
													this_mat[1][3] * new_mat[1][1] +
													this_mat[2][3] * new_mat[1][2] );
			new_mat[2][3] = - ( this_mat[0][3] * new_mat[2][0] +
													this_mat[1][3] * new_mat[2][1] +
													this_mat[2][3] * new_mat[2][2] );
    
			/* Fill in last column */
			new_mat[3][0] = new_mat[3][1] = new_mat[3][2] = 0.0;
			new_mat[3][3] = 1.0;
 
			return new_mat;
	
		}
 
#undef ACCUMULATE
	}
 
	template<class V, class M>
	M adjoint(const ArithSqMat4x4Float<V,M>& in)
	{
		double a1, a2, a3, a4, b1, b2, b3, b4;
		double c1, c2, c3, c4, d1, d2, d3, d4;
 
		/* assign to individual variable names to aid  */
		/* selecting correct values  */
	
		a1 = in[0][0]; b1 = in[0][1]; 
		c1 = in[0][2]; d1 = in[0][3];
 
		a2 = in[1][0]; b2 = in[1][1]; 
		c2 = in[1][2]; d2 = in[1][3];
 
		a3 = in[2][0]; b3 = in[2][1];
		c3 = in[2][2]; d3 = in[2][3];
 
		a4 = in[3][0]; b4 = in[3][1]; 
		c4 = in[3][2]; d4 = in[3][3];
 
 
		/* row column labeling reversed since we transpose rows & columns */
	
		M out;
		out[0][0]  =   d3x3f( b2, b3, b4, c2, c3, c4, d2, d3, d4);
		out[1][0]  = - d3x3f( a2, a3, a4, c2, c3, c4, d2, d3, d4);
		out[2][0]  =   d3x3f( a2, a3, a4, b2, b3, b4, d2, d3, d4);
		out[3][0]  = - d3x3f( a2, a3, a4, b2, b3, b4, c2, c3, c4);
	
		out[0][1]  = - d3x3f( b1, b3, b4, c1, c3, c4, d1, d3, d4);
		out[1][1]  =   d3x3f( a1, a3, a4, c1, c3, c4, d1, d3, d4);
		out[2][1]  = - d3x3f( a1, a3, a4, b1, b3, b4, d1, d3, d4);
		out[3][1]  =   d3x3f( a1, a3, a4, b1, b3, b4, c1, c3, c4);
        
		out[0][2]  =   d3x3f( b1, b2, b4, c1, c2, c4, d1, d2, d4);
		out[1][2]  = - d3x3f( a1, a2, a4, c1, c2, c4, d1, d2, d4);
		out[2][2]  =   d3x3f( a1, a2, a4, b1, b2, b4, d1, d2, d4);
		out[3][2]  = - d3x3f( a1, a2, a4, b1, b2, b4, c1, c2, c4);
	
		out[0][3]  = - d3x3f( b1, b2, b3, c1, c2, c3, d1, d2, d3);
		out[1][3]  =   d3x3f( a1, a2, a3, c1, c2, c3, d1, d2, d3);
		out[2][3]  = - d3x3f( a1, a2, a3, b1, b2, b3, d1, d2, d3);
		out[3][3]  =   d3x3f( a1, a2, a3, b1, b2, b3, c1, c2, c3);
 
		return out;
	}
 
 
	template<class V, class M>
	M invert(const ArithSqMat4x4Float<V,M>& in)
	{
		double det = determinant( in );
		if (is_tiny(det)) 
			throw(Mat4x4fSingular("Tried to invert Singular matrix"));
		M out = adjoint(in);
		out/=det;
		return out;
	}
 
 
	template class ArithSqMat4x4Float<Vec4f, Mat4x4f>;
	template Mat4x4f adjoint(const ArithSqMat4x4Float<Vec4f,Mat4x4f>&);
	template Mat4x4f invert(const ArithSqMat4x4Float<Vec4f,Mat4x4f>&);
	template Mat4x4f invert_affine(const ArithSqMat4x4Float<Vec4f,Mat4x4f>&);
 
	template class ArithSqMat4x4Float<Vec4d, Mat4x4d>;
	template Mat4x4d adjoint(const ArithSqMat4x4Float<Vec4d,Mat4x4d>&);
	template Mat4x4d invert(const ArithSqMat4x4Float<Vec4d,Mat4x4d>&);
	template Mat4x4d invert_affine(const ArithSqMat4x4Float<Vec4d,Mat4x4d>&);
 
 
}
 

Subversion Repositories gelsvn

(root)/trunk/src/CGLA/ArithSqMat4x4Float.cpp – Rev 327 → 595

Rev 327		Rev 595
-		1	`/* ----------------------------------------------------------------------- *`
-		2	`* This file is part of GEL, http://www.imm.dtu.dk/GEL`
-		3	`* Copyright (C) the authors and DTU Informatics`
-		4	`* For license and list of authors, see ../../doc/intro.pdf`
-		5	`* ----------------------------------------------------------------------- */`
-		6
1	`#include "ArithSqMat4x4Float.h"`	7	`#include "ArithSqMat4x4Float.h"`
2	`#include "Mat4x4f.h"`	8	`#include "Mat4x4f.h"`
3	`#include "Mat4x4d.h"`	9	`#include "Mat4x4d.h"`
4		10
5	`namespace CGLA {`	11	`namespace CGLA {`
6		12
7	`namespace`	13	`namespace`
8	`{`	14	`{`
9		15
10	`/* Aux func. computes 3x3 determinants. */`	16	`/* Aux func. computes 3x3 determinants. */`
11	`template<class T>`	17	`template<class T>`
12	`inline T d3x3f( T a0, T a1, T a2,`	18	`inline T d3x3f( T a0, T a1, T a2,`
13	`T b0, T b1, T b2,`	19	`T b0, T b1, T b2,`
14	`T c0, T c1, T c2 )`	20	`T c0, T c1, T c2 )`
15	`{`	21	`{`
16	`return a0b1c2 +`	22	`return a0b1c2 +`
17	`a1b2c0 +`	23	`a1b2c0 +`
18	`a2b0c1 -`	24	`a2b0c1 -`
19	`a2b1c0 -`	25	`a2b1c0 -`
20	`a0b2c1 -`	26	`a0b2c1 -`
21	`a1b0c2 ;`	27	`a1b0c2 ;`
22	`}`	28	`}`
23	`}`	29	`}`
24		30
25		31
26		32
27	`template<class V, class M>`	33	`template<class V, class M>`
28	`M invert_affine(const ArithSqMat4x4Float<V,M>& this_mat)`	34	`M invert_affine(const ArithSqMat4x4Float<V,M>& this_mat)`
29	`{`	35	`{`
30	`// The following com[3]e has been copied from a gem in Graphics Gems II by`	36	`// The following com[3]e has been copied from a gem in Graphics Gems II by`
31	`// Kevin Wu.`	37	`// Kevin Wu.`
-		38	`// From the EULA: "Using the code is permitted in any program, product, or`
-		39	`// library, non-commercial or commercial. Giving credit is not required, though is a nice gesture"`
32	`// The function is very fast, but it can only invert affine matrices. An`	40	`// The function is very fast, but it can only invert affine matrices. An`
33	`// exception NotAffine is thrown if the matrix is not affine, and another`	41	`// exception NotAffine is thrown if the matrix is not affine, and another`
34	`// exception Singular is thrown if the matrix is singular.`	42	`// exception Singular is thrown if the matrix is singular.`
35		43
36	`typedef typename M::ScalarType ScalarType;`	44	`typedef typename M::ScalarType ScalarType;`
37		45
38	`M new_mat;`	46	`M new_mat;`
39	`ScalarType det_1;`	47	`ScalarType det_1;`
40	`ScalarType pos, neg, temp;`	48	`ScalarType pos, neg, temp;`
41		49
42		50
43	`if (!(is_tiny(this_mat[3][0]) &&`	51	`if (!(is_tiny(this_mat[3][0]) &&`
44	`is_tiny(this_mat[3][1]) &&`	52	`is_tiny(this_mat[3][1]) &&`
45	`is_tiny(this_mat[3][2]) &&`	53	`is_tiny(this_mat[3][2]) &&`
46	`is_tiny(this_mat[3][3]-1.0)))`	54	`is_tiny(this_mat[3][3]-1.0)))`
47	`throw(Mat4x4fNotAffine("Can only invert affine matrices"));`	55	`throw(Mat4x4fNotAffine("Can only invert affine matrices"));`
48		56
49	`#define ACCUMULATE if (temp >= 0.0) pos += temp; else neg += temp`	57	`#define ACCUMULATE if (temp >= 0.0) pos += temp; else neg += temp`
50		58
51	`/*`	59	`/*`
52	`* Calculate the determinant of submatrix A and determine if the`	60	`* Calculate the determinant of submatrix A and determine if the`
53	`* the matrix is singular as limited by the float precision`	61	`* the matrix is singular as limited by the float precision`
54	`* floating-point this_mat representation.`	62	`* floating-point this_mat representation.`
55	`*/`	63	`*/`
56		64
57	`pos = neg = 0.0;`	65	`pos = neg = 0.0;`
58	`temp = this_mat[0][0] * this_mat[1][1] * this_mat[2][2];`	66	`temp = this_mat[0][0] * this_mat[1][1] * this_mat[2][2];`
59	`ACCUMULATE;`	67	`ACCUMULATE;`
60	`temp = this_mat[1][0] * this_mat[2][1] * this_mat[0][2];`	68	`temp = this_mat[1][0] * this_mat[2][1] * this_mat[0][2];`
61	`ACCUMULATE;`	69	`ACCUMULATE;`
62	`temp = this_mat[2][0] * this_mat[0][1] * this_mat[1][2];`	70	`temp = this_mat[2][0] * this_mat[0][1] * this_mat[1][2];`
63	`ACCUMULATE;`	71	`ACCUMULATE;`
64	`temp = -this_mat[2][0] * this_mat[1][1] * this_mat[0][2];`	72	`temp = -this_mat[2][0] * this_mat[1][1] * this_mat[0][2];`
65	`ACCUMULATE;`	73	`ACCUMULATE;`
66	`temp = -this_mat[1][0] * this_mat[0][1] * this_mat[2][2];`	74	`temp = -this_mat[1][0] * this_mat[0][1] * this_mat[2][2];`
67	`ACCUMULATE;`	75	`ACCUMULATE;`
68	`temp = -this_mat[0][0] * this_mat[2][1] * this_mat[1][2];`	76	`temp = -this_mat[0][0] * this_mat[2][1] * this_mat[1][2];`
69	`ACCUMULATE;`	77	`ACCUMULATE;`
70	`det_1 = pos + neg;`	78	`det_1 = pos + neg;`
71		79
72	`/* Is the submatrix A singular? */`	80	`/* Is the submatrix A singular? */`
73	`if ((det_1 == 0.0) \|\| (fabs(det_1 / (pos - neg)) < MINUTE))`	81	`if ((det_1 == 0.0) \|\| (fabs(det_1 / (pos - neg)) < MINUTE))`
74	`{`	82	`{`
75	`/* Mat4x4f M has no inverse */`	83	`/* Mat4x4f M has no inverse */`
76	`throw(Mat4x4fSingular("Tried to invert Singular matrix"));`	84	`throw(Mat4x4fSingular("Tried to invert Singular matrix"));`
77	`}`	85	`}`
78		86
79	`else {`	87	`else {`
80		88
81	`/* Calculate inverse(A) = adj(A) / det(A) */`	89	`/* Calculate inverse(A) = adj(A) / det(A) */`
82	`det_1 = 1.0 / det_1;`	90	`det_1 = 1.0 / det_1;`
83	`new_mat[0][0] = ( this_mat[1][1] * this_mat[2][2] -`	91	`new_mat[0][0] = ( this_mat[1][1] * this_mat[2][2] -`
84	`this_mat[2][1] * this_mat[1][2] )`	92	`this_mat[2][1] * this_mat[1][2] )`
85	`* det_1;`	93	`* det_1;`
86	`new_mat[0][1] = - ( this_mat[0][1] * this_mat[2][2] -`	94	`new_mat[0][1] = - ( this_mat[0][1] * this_mat[2][2] -`
87	`this_mat[2][1] * this_mat[0][2] )`	95	`this_mat[2][1] * this_mat[0][2] )`
88	`* det_1;`	96	`* det_1;`
89	`new_mat[0][2] = ( this_mat[0][1] * this_mat[1][2] -`	97	`new_mat[0][2] = ( this_mat[0][1] * this_mat[1][2] -`
90	`this_mat[1][1] * this_mat[0][2] )`	98	`this_mat[1][1] * this_mat[0][2] )`
91	`* det_1;`	99	`* det_1;`
92	`new_mat[1][0] = - ( this_mat[1][0] * this_mat[2][2] -`	100	`new_mat[1][0] = - ( this_mat[1][0] * this_mat[2][2] -`
93	`this_mat[2][0] * this_mat[1][2] )`	101	`this_mat[2][0] * this_mat[1][2] )`
94	`* det_1;`	102	`* det_1;`
95	`new_mat[1][1] = ( this_mat[0][0] * this_mat[2][2] -`	103	`new_mat[1][1] = ( this_mat[0][0] * this_mat[2][2] -`
96	`this_mat[2][0] * this_mat[0][2] )`	104	`this_mat[2][0] * this_mat[0][2] )`
97	`* det_1;`	105	`* det_1;`
98	`new_mat[1][2] = - ( this_mat[0][0] * this_mat[1][2] -`	106	`new_mat[1][2] = - ( this_mat[0][0] * this_mat[1][2] -`
99	`this_mat[1][0] * this_mat[0][2] )`	107	`this_mat[1][0] * this_mat[0][2] )`
100	`* det_1;`	108	`* det_1;`
101	`new_mat[2][0] = ( this_mat[1][0] * this_mat[2][1] -`	109	`new_mat[2][0] = ( this_mat[1][0] * this_mat[2][1] -`
102	`this_mat[2][0] * this_mat[1][1] )`	110	`this_mat[2][0] * this_mat[1][1] )`
103	`* det_1;`	111	`* det_1;`
104	`new_mat[2][1] = - ( this_mat[0][0] * this_mat[2][1] -`	112	`new_mat[2][1] = - ( this_mat[0][0] * this_mat[2][1] -`
105	`this_mat[2][0] * this_mat[0][1] )`	113	`this_mat[2][0] * this_mat[0][1] )`
106	`* det_1;`	114	`* det_1;`
107	`new_mat[2][2] = ( this_mat[0][0] * this_mat[1][1] -`	115	`new_mat[2][2] = ( this_mat[0][0] * this_mat[1][1] -`
108	`this_mat[1][0] * this_mat[0][1] )`	116	`this_mat[1][0] * this_mat[0][1] )`
109	`* det_1;`	117	`* det_1;`
110		118
111	`/* Calculate -C * inverse(A) */`	119	`/* Calculate -C * inverse(A) */`
112	`new_mat[0][3] = - ( this_mat[0][3] * new_mat[0][0] +`	120	`new_mat[0][3] = - ( this_mat[0][3] * new_mat[0][0] +`
113	`this_mat[1][3] * new_mat[0][1] +`	121	`this_mat[1][3] * new_mat[0][1] +`
114	`this_mat[2][3] * new_mat[0][2] );`	122	`this_mat[2][3] * new_mat[0][2] );`
115	`new_mat[1][3] = - ( this_mat[0][3] * new_mat[1][0] +`	123	`new_mat[1][3] = - ( this_mat[0][3] * new_mat[1][0] +`
116	`this_mat[1][3] * new_mat[1][1] +`	124	`this_mat[1][3] * new_mat[1][1] +`
117	`this_mat[2][3] * new_mat[1][2] );`	125	`this_mat[2][3] * new_mat[1][2] );`
118	`new_mat[2][3] = - ( this_mat[0][3] * new_mat[2][0] +`	126	`new_mat[2][3] = - ( this_mat[0][3] * new_mat[2][0] +`
119	`this_mat[1][3] * new_mat[2][1] +`	127	`this_mat[1][3] * new_mat[2][1] +`
120	`this_mat[2][3] * new_mat[2][2] );`	128	`this_mat[2][3] * new_mat[2][2] );`
121		129
122	`/* Fill in last column */`	130	`/* Fill in last column */`
123	`new_mat[3][0] = new_mat[3][1] = new_mat[3][2] = 0.0;`	131	`new_mat[3][0] = new_mat[3][1] = new_mat[3][2] = 0.0;`
124	`new_mat[3][3] = 1.0;`	132	`new_mat[3][3] = 1.0;`
125		133
126	`return new_mat;`	134	`return new_mat;`
127		135
128	`}`	136	`}`
129		137
130	`#undef ACCUMULATE`	138	`#undef ACCUMULATE`
131	`}`	139	`}`
132		140
133	`template<class V, class M>`	141	`template<class V, class M>`
134	`M adjoint(const ArithSqMat4x4Float<V,M>& in)`	142	`M adjoint(const ArithSqMat4x4Float<V,M>& in)`
135	`{`	143	`{`
136	`double a1, a2, a3, a4, b1, b2, b3, b4;`	144	`double a1, a2, a3, a4, b1, b2, b3, b4;`
137	`double c1, c2, c3, c4, d1, d2, d3, d4;`	145	`double c1, c2, c3, c4, d1, d2, d3, d4;`
138		146
139	`/* assign to individual variable names to aid */`	147	`/* assign to individual variable names to aid */`
140	`/* selecting correct values */`	148	`/* selecting correct values */`
141		149
142	`a1 = in[0][0]; b1 = in[0][1];`	150	`a1 = in[0][0]; b1 = in[0][1];`
143	`c1 = in[0][2]; d1 = in[0][3];`	151	`c1 = in[0][2]; d1 = in[0][3];`
144		152
145	`a2 = in[1][0]; b2 = in[1][1];`	153	`a2 = in[1][0]; b2 = in[1][1];`
146	`c2 = in[1][2]; d2 = in[1][3];`	154	`c2 = in[1][2]; d2 = in[1][3];`
147		155
148	`a3 = in[2][0]; b3 = in[2][1];`	156	`a3 = in[2][0]; b3 = in[2][1];`
149	`c3 = in[2][2]; d3 = in[2][3];`	157	`c3 = in[2][2]; d3 = in[2][3];`
150		158
151	`a4 = in[3][0]; b4 = in[3][1];`	159	`a4 = in[3][0]; b4 = in[3][1];`
152	`c4 = in[3][2]; d4 = in[3][3];`	160	`c4 = in[3][2]; d4 = in[3][3];`
153		161
154		162
155	`/* row column labeling reversed since we transpose rows & columns */`	163	`/* row column labeling reversed since we transpose rows & columns */`
156		164
157	`M out;`	165	`M out;`
158	`out[0][0] = d3x3f( b2, b3, b4, c2, c3, c4, d2, d3, d4);`	166	`out[0][0] = d3x3f( b2, b3, b4, c2, c3, c4, d2, d3, d4);`
159	`out[1][0] = - d3x3f( a2, a3, a4, c2, c3, c4, d2, d3, d4);`	167	`out[1][0] = - d3x3f( a2, a3, a4, c2, c3, c4, d2, d3, d4);`
160	`out[2][0] = d3x3f( a2, a3, a4, b2, b3, b4, d2, d3, d4);`	168	`out[2][0] = d3x3f( a2, a3, a4, b2, b3, b4, d2, d3, d4);`
161	`out[3][0] = - d3x3f( a2, a3, a4, b2, b3, b4, c2, c3, c4);`	169	`out[3][0] = - d3x3f( a2, a3, a4, b2, b3, b4, c2, c3, c4);`
162		170
163	`out[0][1] = - d3x3f( b1, b3, b4, c1, c3, c4, d1, d3, d4);`	171	`out[0][1] = - d3x3f( b1, b3, b4, c1, c3, c4, d1, d3, d4);`
164	`out[1][1] = d3x3f( a1, a3, a4, c1, c3, c4, d1, d3, d4);`	172	`out[1][1] = d3x3f( a1, a3, a4, c1, c3, c4, d1, d3, d4);`
165	`out[2][1] = - d3x3f( a1, a3, a4, b1, b3, b4, d1, d3, d4);`	173	`out[2][1] = - d3x3f( a1, a3, a4, b1, b3, b4, d1, d3, d4);`
166	`out[3][1] = d3x3f( a1, a3, a4, b1, b3, b4, c1, c3, c4);`	174	`out[3][1] = d3x3f( a1, a3, a4, b1, b3, b4, c1, c3, c4);`
167		175
168	`out[0][2] = d3x3f( b1, b2, b4, c1, c2, c4, d1, d2, d4);`	176	`out[0][2] = d3x3f( b1, b2, b4, c1, c2, c4, d1, d2, d4);`
169	`out[1][2] = - d3x3f( a1, a2, a4, c1, c2, c4, d1, d2, d4);`	177	`out[1][2] = - d3x3f( a1, a2, a4, c1, c2, c4, d1, d2, d4);`
170	`out[2][2] = d3x3f( a1, a2, a4, b1, b2, b4, d1, d2, d4);`	178	`out[2][2] = d3x3f( a1, a2, a4, b1, b2, b4, d1, d2, d4);`
171	`out[3][2] = - d3x3f( a1, a2, a4, b1, b2, b4, c1, c2, c4);`	179	`out[3][2] = - d3x3f( a1, a2, a4, b1, b2, b4, c1, c2, c4);`
172		180
173	`out[0][3] = - d3x3f( b1, b2, b3, c1, c2, c3, d1, d2, d3);`	181	`out[0][3] = - d3x3f( b1, b2, b3, c1, c2, c3, d1, d2, d3);`
174	`out[1][3] = d3x3f( a1, a2, a3, c1, c2, c3, d1, d2, d3);`	182	`out[1][3] = d3x3f( a1, a2, a3, c1, c2, c3, d1, d2, d3);`
175	`out[2][3] = - d3x3f( a1, a2, a3, b1, b2, b3, d1, d2, d3);`	183	`out[2][3] = - d3x3f( a1, a2, a3, b1, b2, b3, d1, d2, d3);`
176	`out[3][3] = d3x3f( a1, a2, a3, b1, b2, b3, c1, c2, c3);`	184	`out[3][3] = d3x3f( a1, a2, a3, b1, b2, b3, c1, c2, c3);`
177		185
178	`return out;`	186	`return out;`
179	`}`	187	`}`
180		188
181		189
182	`template<class V, class M>`	190	`template<class V, class M>`
183	`M invert(const ArithSqMat4x4Float<V,M>& in)`	191	`M invert(const ArithSqMat4x4Float<V,M>& in)`
184	`{`	192	`{`
185	`double det = determinant( in );`	193	`double det = determinant( in );`
186	`if (is_tiny(det))`	194	`if (is_tiny(det))`
187	`throw(Mat4x4fSingular("Tried to invert Singular matrix"));`	195	`throw(Mat4x4fSingular("Tried to invert Singular matrix"));`
188	`M out = adjoint(in);`	196	`M out = adjoint(in);`
189	`out/=det;`	197	`out/=det;`
190	`return out;`	198	`return out;`
191	`}`	199	`}`
192		200
193		201
194	`template class ArithSqMat4x4Float<Vec4f, Mat4x4f>;`	202	`template class ArithSqMat4x4Float<Vec4f, Mat4x4f>;`
195	`template Mat4x4f adjoint(const ArithSqMat4x4Float<Vec4f,Mat4x4f>&);`	203	`template Mat4x4f adjoint(const ArithSqMat4x4Float<Vec4f,Mat4x4f>&);`
196	`template Mat4x4f invert(const ArithSqMat4x4Float<Vec4f,Mat4x4f>&);`	204	`template Mat4x4f invert(const ArithSqMat4x4Float<Vec4f,Mat4x4f>&);`
197	`template Mat4x4f invert_affine(const ArithSqMat4x4Float<Vec4f,Mat4x4f>&);`	205	`template Mat4x4f invert_affine(const ArithSqMat4x4Float<Vec4f,Mat4x4f>&);`
198		206
199	`template class ArithSqMat4x4Float<Vec4d, Mat4x4d>;`	207	`template class ArithSqMat4x4Float<Vec4d, Mat4x4d>;`
200	`template Mat4x4d adjoint(const ArithSqMat4x4Float<Vec4d,Mat4x4d>&);`	208	`template Mat4x4d adjoint(const ArithSqMat4x4Float<Vec4d,Mat4x4d>&);`
201	`template Mat4x4d invert(const ArithSqMat4x4Float<Vec4d,Mat4x4d>&);`	209	`template Mat4x4d invert(const ArithSqMat4x4Float<Vec4d,Mat4x4d>&);`
202	`template Mat4x4d invert_affine(const ArithSqMat4x4Float<Vec4d,Mat4x4d>&);`	210	`template Mat4x4d invert_affine(const ArithSqMat4x4Float<Vec4d,Mat4x4d>&);`
203		211
204		212
205	`}`	213	`}`
206		214