From 6595e8a713b96dd82b0efeaae380f635302d1658 Mon Sep 17 00:00:00 2001
From: shahsaurabh0605 <shahsaurabh0605@gmail.com>
Date: Fri, 29 Jan 2016 21:22:30 +0530
Subject: [PATCH] corrected inverse_exact

---
 ext/nmatrix/math.cpp | 30 ++++++++++++++++++++----------
 lib/nmatrix/math.rb  |  9 ++++++++-
 2 files changed, 28 insertions(+), 11 deletions(-)

diff --git a/ext/nmatrix/math.cpp b/ext/nmatrix/math.cpp
index 095f8909..75c55371 100644
--- a/ext/nmatrix/math.cpp
+++ b/ext/nmatrix/math.cpp
@@ -385,10 +385,11 @@ namespace nm {
           rb_raise(nm_eNotInvertibleError, 
               "matrix must have non-zero determinant to be invertible (not getting this error does not mean matrix is invertible if you're dealing with floating points)");
         }
+        DType index= A[0];
         B[0] = A[lda+1] / det;
         B[1] = -A[1] / det;
         B[ldb] = -A[lda] / det;
-        B[ldb+1] = -A[0] / det;
+        B[ldb+1] = index/ det;
 
       } else if (M == 3) {
         // Calculate the exact determinant.
@@ -398,16 +399,24 @@ namespace nm {
           rb_raise(nm_eNotInvertibleError, 
               "matrix must have non-zero determinant to be invertible (not getting this error does not mean matrix is invertible if you're dealing with floating points)");
         }
-
+        DType a1= A[0];
+        DType b1= A[1];
+        DType c1= A[2];
+        DType d1= A[lda];
+        DType e1= A[lda+1];
+        DType f1= A[lda+2];
+        DType g1= A[2*lda];
+        DType h1= A[2*lda+1];
+        DType i1= A[2*lda+2];
         B[0]      = (  A[lda+1] * A[2*lda+2] - A[lda+2] * A[2*lda+1]) / det; // A = ei - fh
-        B[1]      = (- A[1]     * A[2*lda+2] + A[2]     * A[2*lda+1]) / det; // D = -bi + ch
-        B[2]      = (  A[1]     * A[lda+2]   - A[2]     * A[lda+1])   / det; // G = bf - ce
-        B[ldb]    = (- A[lda]   * A[2*lda+2] + A[lda+2] * A[2*lda])   / det; // B = -di + fg
-        B[ldb+1]  = (  A[0]     * A[2*lda+2] - A[2]     * A[2*lda])   / det; // E = ai - cg
-        B[ldb+2]  = (- A[0]     * A[lda+2]   + A[2]     * A[lda])     / det; // H = -af + cd
-        B[2*ldb]  = (  A[lda]   * A[2*lda+1] - A[lda+1] * A[2*lda])   / det; // C = dh - eg
-        B[2*ldb+1]= ( -A[0]     * A[2*lda+1] + A[1]     * A[2*lda])   / det; // F = -ah + bg
-        B[2*ldb+2]= (  A[0]     * A[lda+1]   - A[1]     * A[lda])     / det; // I = ae - bd
+        B[1]      = (- A[1]     * i1 + A[2]     * h1) / det; // D = -bi + ch
+        B[2]      = (  b1     * f1   - c1     * e1)   / det; // G = bf - ce
+        B[ldb]    = (- d1   * i1 + f1 * g1)   / det; // B = -di + fg
+        B[ldb+1]  = (  a1     * i1 - c1     * g1)   / det; // E = ai - cg
+        B[ldb+2]  = (- a1     * f1   + c1     * d1)     / det; // H = -af + cd
+        B[2*ldb]  = (  d1   * h1 - e1 * g1)   / det; // C = dh - eg
+        B[2*ldb+1]= ( -a1     * h1 + b1     * g1)   / det; // F = -ah + bg
+        B[2*ldb+2]= (  a1     * e1   - b1     * d1)     / det; // I = ae - bd
       } else if (M == 1) {
         B[0] = 1 / A[0];
       } else {
@@ -1125,3 +1134,4 @@ void nm_math_transpose_generic(const size_t M, const size_t N, const void* A, co
 
 
 } // end of extern "C" block
+
diff --git a/lib/nmatrix/math.rb b/lib/nmatrix/math.rb
index a8f3ceaf..3a15d5ff 100644
--- a/lib/nmatrix/math.rb
+++ b/lib/nmatrix/math.rb
@@ -80,7 +80,13 @@ def invert!
     raise(DataTypeError, "Cannot invert an integer matrix in-place") if self.integer_dtype?
 
     #No internal implementation of getri, so use this other function
-    __inverse__(self, true)
+    n= self.shape[0]
+    if n<4
+      __inverse_exact__(self, n,n)
+    else
+      __inverse__(self,true)
+      
+    end
   end
 
   #
@@ -1102,3 +1108,4 @@ def dtype_for_floor_or_ceil
     end
   end
 end
+