#/***********************************************************
#    matinv.rb -- 逆行列
#***********************************************************/
require "matutil.rb"  # 行列操作用小道具集 

def lu( n, a, ip)  # LU分解 
    weight = new_vector(n);    # weight[0..n-1] の記憶領域確保 
    det = 0;                   # 行列式 
    for k in 0...n
        ip[k] = k              # 行交換情報の初期値 
        u = 0                  # その行の絶対値最大の要素を求める 
        for j in 0...n
            t = a[k][j].abs;  if (t > u); u = t; end
        end
        if (u == 0); return det; end  # 0 なら行列はLU分解できない 
        weight[k] = 1 / u     # 最大絶対値の逆数 
    end
    det = 1                   # 行列式の初期値 
    for k in 0...n
        u = -1
        for i in k...n
            ii = ip[i]            # 重み×絶対値 が最大の行を見つける 
            t = a[ii][k].abs * weight[ii]
            if (t > u);  u = t;  j = i;  end
        end
        ik = ip[j]
        if (j != k) 
            ip[j] = ip[k];  ip[k] = ik;  # 行番号を交換 
            det = -det;  # 行を交換すれば行列式の符号が変わる 
        end
        u = a[ik][k];  det *= u;  # 対角成分 
        if (u == 0); return det; end    # 0 なら行列はLU分解できない 
        for i in (k + 1)...n
            ii = ip[i]
            t = (a[ii][k] /= u)
            for j in (k + 1)...n
                a[ii][j] -= t * a[ik][j]
            end
        end
    end
    return det;           # 戻り値は行列式 
end

def matinv( n, a, a_inv)
    ip = []
    det = lu(n, a, ip)
    if (det != 0)
        for k in 0...n
            for i in 0...n
                ii = ip[i];  t = (ii == k) ? 1 : 0
                for j in 0...i
                    t -= a[ii][j] * a_inv[j][k]
                end
                a_inv[i][k] = t
            end
            (n - 1).downto(0) do |i|
                t = a_inv[i][k];  ii = ip[i]
                for j in (i + 1)...n
                    t -= a[ii][j] * a_inv[j][k]
                end
                a_inv[i][k] = t / a[ii][i]
            end
        end
    end
    return det
end

#************ 以下はテスト用 ***************

ULONG_MAX = 4294967295
ULONG_MAXplus = 4294967296
$seed = 123456789

def rnd  # 乱数  0 < rnd() < 1 
    $seed = $seed * 69069 % ULONG_MAXplus
    return $seed / (ULONG_MAX + 1.0)
end

printf("n = ");  n = gets.to_i
a = new_matrix(n, n);  asave = new_matrix(n, n)
a_inv = new_matrix(n, n)
for i in 0...n
    for j in 0...n
        a[i][j] = asave[i][j] = rnd() - rnd()
    end
end
printf("A\n")
matprint(a, n, 7, "%10.6f")
s = matinv(n, a, a_inv)
printf("行列式 = %g\n", s)
printf("A^-1end\n")
matprint(a_inv, n, 7, "%10.6f")
t = 0
for i in 0...n
    for j in 0...n
        s = (i == j) ? 1 : 0
        for k in 0...n
            s -= asave[i][k] * a_inv[k][j]
        end
        t += s * s
    end
end
printf("A A^-1end の成分の二乗平均誤差 %g\n", \
    Math::sqrt(t / (n * n)))
t = 0
for i in 0...n
    for j in 0...n
        s = (i == j) ? 1 : 0
        for k in 0...n
            s -= a_inv[i][k] * asave[k][j]
        end
        t += s * s
    end
end
printf("A^-1end A の成分の二乗平均誤差 %g\n", \
    Math::sqrt(t / (n * n)))

exit 0