Sweet Snippet 之矩阵求逆

-- matrix store is table in row order
-- e.g. 2 x 2 matrix is stored as table { m11, m12, m21, m22 }

-- determinant 2 with scalar elements
function det2s(e11, e12, e21, e22)
    return e11 * e22 - e21 * e12
end

-- determinant 3 with scalar elements
function det3s(e11, e12, e13, e21, e22, e23, e31, e32, e33)
    return e11 * det2s(e22, e23, e32, e33) - e12 * det2s(e21, e23, e31, e33) + e13 * det2s(e21, e22, e31, e32)
end
 
-- determinant 2
function det2(m2)
    return det2s(table.unpack(m2))
end
 
-- determinant 3
function det3(m3)
    return det3s(table.unpack(m3))
end
 
-- multiply 2 
function mul2(lm2, rm2)
    if lm2 and rm2 and #lm2 == #rm2 and #lm2 == 4 then
        local m2 = {}
        
        m2[1] = lm2[1] * rm2[1] + lm2[2] * rm2[3]
        m2[2] = lm2[1] * rm2[2] + lm2[2] * rm2[4]
        m2[3] = lm2[3] * rm2[1] + lm2[4] * rm2[3]
        m2[4] = lm2[3] * rm2[2] + lm2[4] * rm2[4]
        
        return m2
    end
end

-- multiply 3
function mul3(lm3, rm3)
    if lm3 and rm3 and #lm3 == #rm3 and #lm3 == 9 then
        local m3 = {}
        
        m3[1] = lm3[1] * rm3[1] + lm3[2] * rm3[4] + lm3[3] * rm3[7]
        m3[2] = lm3[1] * rm3[2] + lm3[2] * rm3[5] + lm3[3] * rm3[8]
        m3[3] = lm3[1] * rm3[3] + lm3[2] * rm3[6] + lm3[3] * rm3[9]
        m3[4] = lm3[4] * rm3[1] + lm3[5] * rm3[4] + lm3[6] * rm3[7]
        m3[5] = lm3[4] * rm3[2] + lm3[5] * rm3[5] + lm3[6] * rm3[8]
        m3[6] = lm3[4] * rm3[3] + lm3[5] * rm3[6] + lm3[6] * rm3[9]
        m3[7] = lm3[7] * rm3[1] + lm3[8] * rm3[4] + lm3[9] * rm3[7]
        m3[8] = lm3[7] * rm3[2] + lm3[8] * rm3[5] + lm3[9] * rm3[8]
        m3[9] = lm3[7] * rm3[3] + lm3[8] * rm3[6] + lm3[9] * rm3[9]
        
        return m3
    end
end

-- transpose 2
function trp2(m2)
    if m2 and #m2 == 4 then
        local trp_m2 = {}
        
        trp_m2[1] = m2[1]
        trp_m2[2] = m2[3]
        trp_m2[3] = m2[2]
        trp_m2[4] = m2[4]
        
        return trp_m2
    end
end

-- transpose 3
function trp3(m3)
    if m3 and #m3 == 9 then
        local trp_m3 = {}
        
        trp_m3[1] = m3[1]
        trp_m3[2] = m3[4]
        trp_m3[3] = m3[7]
        trp_m3[4] = m3[2]
        trp_m3[5] = m3[5]
        trp_m3[6] = m3[8]
        trp_m3[7] = m3[3]
        trp_m3[8] = m3[6]
        trp_m3[9] = m3[9]
        
        return trp_m3
    end
end

-- inverse 2
function inv2(m2)
    local det = det2(m2)
    if det and det ~= 0 then
        local inv_m2 = {}
        
        local adjm = 
        { 
            m2[4], 
            m2[3], 
            m2[2], 
            m2[1] 
        }
        adjm = trp2(adjm)
        
        inv_m2[1] = adjm[1] / det
        inv_m2[2] = -adjm[2] / det
        inv_m2[3] = -adjm[3] / det
        inv_m2[4] = adjm[4] / det
        
        return inv_m2
    end
end

-- inverse 3
function inv3(m3)
    local det = det3(m3)
    if det and det ~= 0 then
        local inv_m3 = {}
        
        local adjm = 
        { 
            det2s(m3[5], m3[6], m3[8], m3[9]),
            det2s(m3[4], m3[6], m3[7], m3[9]),
            det2s(m3[4], m3[5], m3[7], m3[8]),
            det2s(m3[2], m3[3], m3[8], m3[9]),
            det2s(m3[1], m3[3], m3[7], m3[9]),
            det2s(m3[1], m3[2], m3[7], m3[8]),
            det2s(m3[2], m3[3], m3[5], m3[6]),
            det2s(m3[1], m3[3], m3[4], m3[6]),
            det2s(m3[1], m3[2], m3[4], m3[5]),
        }
        adjm = trp3(adjm)
        
        inv_m3[1] = adjm[1] / det
        inv_m3[2] = -adjm[2] / det
        inv_m3[3] = adjm[3] / det
        inv_m3[4] = -adjm[4] / det
        inv_m3[5] = adjm[5] / det
        inv_m3[6] = -adjm[6] / det
        inv_m3[7] = adjm[7] / det
        inv_m3[8] = -adjm[8] / det
        inv_m3[9] = adjm[9] / det
        
        return inv_m3
    end
end

local m3 = 
{ 
    1, 2, 3, 
    4, 5, 6, 
    7, 8, 9 
}