From: John Tsiombikas <nuclear@member.fsf.org>
Date: Sat, 13 Oct 2018 07:36:27 +0000 (+0300)
Subject: fixed inverse and submatrix, now testing identical to gmath
X-Git-Url: http://git.mutantstargoat.com/user/nuclear/?p=gph-cmath;a=commitdiff_plain;h=28ac94545d140161b5a772b7436f176bbbde3802

fixed inverse and submatrix, now testing identical to gmath
---

diff --git a/src/cgmath.h b/src/cgmath.h
index e6648ec..675daba 100644
--- a/src/cgmath.h
+++ b/src/cgmath.h
@@ -134,9 +134,9 @@ static inline void cgm_mmul(float *a, const float *b);
 
 static inline void cgm_msubmatrix(float *m, int row, int col);
 static inline void cgm_mupper3(float *m);
-static inline float cgm_msubdet(float *m, int row, int col);
-static inline float cgm_mcofactor(float *m, int row, int col);
-static inline float cgm_mdet(float *m);
+static inline float cgm_msubdet(const float *m, int row, int col);
+static inline float cgm_mcofactor(const float *m, int row, int col);
+static inline float cgm_mdet(const float *m);
 static inline void cgm_mtranspose(float *m);
 static inline void cgm_mcofmatrix(float *m);
 static inline int cgm_minverse(float *m);	/* returns 0 on success, -1 for singular */
diff --git a/src/cgmmat.inl b/src/cgmmat.inl
index c1b44fa..f64b259 100644
--- a/src/cgmmat.inl
+++ b/src/cgmmat.inl
@@ -41,11 +41,11 @@ static inline void cgm_msubmatrix(float *m, int row, int col)
 
 	subi = 0;
 	for(i=0; i<4; i++) {
-		if(i == col) continue;
+		if(i == row) continue;
 
 		subj = 0;
 		for(j=0; j<4; j++) {
-			if(j == row) continue;
+			if(j == col) continue;
 
 			m[subi * 4 + subj++] = orig[i * 4 + j];
 		}
@@ -61,19 +61,28 @@ static inline void cgm_mupper3(float *m)
 	m[15] = 1.0f;
 }
 
-static inline float cgm_msubdet(float *m, int row, int col)
+static inline float cgm_msubdet(const float *m, int row, int col)
 {
-	cgm_msubmatrix(m, row, col);
-	return cgm_mdet(m);
+	float tmp[16];
+	float subdet00, subdet01, subdet02;
+
+	cgm_mcopy(tmp, m);
+	cgm_msubmatrix(tmp, row, col);
+
+	subdet00 = tmp[5] * tmp[10] - tmp[6] * tmp[9];
+	subdet01 = tmp[4] * tmp[10] - tmp[6] * tmp[8];
+	subdet02 = tmp[4] * tmp[9] - tmp[5] * tmp[8];
+
+	return tmp[0] * subdet00 - tmp[1] * subdet01 + tmp[2] * subdet02;
 }
 
-static inline float cgm_mcofactor(float *m, int row, int col)
+static inline float cgm_mcofactor(const float *m, int row, int col)
 {
 	float min = cgm_msubdet(m, row, col);
 	return (row + col) & 1 ? -min : min;
 }
 
-static inline float cgm_mdet(float *m)
+static inline float cgm_mdet(const float *m)
 {
 	return m[0] * cgm_msubdet(m, 0, 0) - m[1] * cgm_msubdet(m, 0, 1) +
 		m[2] * cgm_msubdet(m, 0, 2) - m[3] * cgm_msubdet(m, 0, 3);
@@ -123,6 +132,7 @@ static inline int cgm_minverse(float *m)
 			m[i * 4 + j] = cgm_mcofactor(tmp, j, i) * inv_det;	/* transposed */
 		}
 	}
+	return 0;
 }
 
 static inline void cgm_mtranslation(float *m, float x, float y, float z)
diff --git a/test/Makefile b/test/Makefile
index c23e3f1..4c8e22a 100644
--- a/test/Makefile
+++ b/test/Makefile
@@ -1,4 +1,5 @@
 obj = test.o
+dep = $(obj:.o=.d)
 bin = test
 
 CXXFLAGS = -pedantic -Wall -g -I../src
@@ -7,6 +8,15 @@ LDFLAGS = -lgmath -lm
 $(bin): $(obj)
 	$(CXX) -o $@ $(obj) $(LDFLAGS)
 
+-include $(dep)
+
+%.d: %.cc
+	@$(CXX) $(CXXFLAGS) $< -MM -MT $(@:.d=.o) >$@
+
 .PHONY: clean
 clean:
 	rm -f $(obj) $(bin)
+
+.PHONY: cleandep
+cleandep:
+	rm -f $(dep)
diff --git a/test/test.cc b/test/test.cc
index 640cd18..162c82d 100644
--- a/test/test.cc
+++ b/test/test.cc
@@ -80,6 +80,30 @@ bool test_mat()
 	a.pre_rotate_y(trans.x);
 	ASSERT_EQ(a, ca);
 
+	cgm_mcopy(cb, ca);
+	cgm_mtranspose(cb);
+	b = a;
+	b.transpose();
+	ASSERT_EQ(b, cb);
+
+	for(int i=0; i<4; i++) {
+		for(int j=0; j<4; j++) {
+			cgm_mcopy(cb, ca);
+			b = a;
+
+			cgm_msubmatrix(cb, i, j);
+			Mat3 sub = b.submatrix(i, j);
+			b = Mat4::identity;
+			float *bp = b[0];
+			float *sp = sub[0];
+			bp[0] = sp[0]; bp[1] = sp[1]; bp[2] = sp[2];
+			bp[4] = sp[3]; bp[5] = sp[4]; bp[6] = sp[5];
+			bp[8] = sp[6]; bp[9] = sp[7]; bp[10] = sp[8];
+
+			ASSERT_EQ(b, cb);
+		}
+	}
+
 	cgm_minverse(ca);
 	a.inverse();
 	ASSERT_EQ(a, ca);