#/***********************************************************
#    gjmatinv.c -- 逆行列
#***********************************************************/
require "matutil.rb"

def matinv( n, a)
    det = 1
    for k in 0...n
        t = a[k][k];  det *= t
        for i in 0...n; a[k][i] /= t; end
        a[k][k] = 1 / t
        for j in 0...n
            if (j != k) 
                u = a[j][k]
                for i in 0...n
                    if (i != k); a[j][i] -= a[k][i] * u
                    else;        a[j][i] = -u / t; end
                end
            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);  b = new_matrix(n, n)
for i in 0...n
    for j in 0...n
        a[i][j] = b[i][j] = rnd() - rnd()
    end
end
printf("A\n")
matprint(a, n, 7, "%10.6f")
s = matinv(n, b)
printf("行列式 = %g\n", s)
printf("A^-1end\n")
matprint(b, 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 -= a[i][k] * b[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 -= b[i][k] * a[k][j]
        end
        t += s * s
    end
end
printf("A^-1end A の成分の二乗平均誤差 %g\n", \
    Math::sqrt(t / (n * n)))

exit 0