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

def invr( n, x )  # 上三角行列の逆行列 
    for k in 0...n
        x[k][k] = 1 / x[k][k]
        (k-1).downto(0) do |j|
            s = 0
            for i in j..k
                s -= x[i][j] * x[i][k]
            end
            x[j][k] = s * x[j][j]
        end
    end
    x
end

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

$seed = 123456789 # 奇数
ULONG_MAX = 4294967295
ULONG_MAXplus = 4294967296

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

printf("n = ");  n = gets.to_i
r = new_matrix(n, n);  
r_inv = new_matrix(n, n)
rr = new_matrix(n, n)
for i in 0...n
    for j in 0...n
        r[i][j] = r_inv[i][j] = 0
    end
end
for i in 0...n
    for j in i...n
        r[i][j] = r_inv[j][i] = rnd - rnd
    end
end
printf("$R$\n")
matprint(r, n, 5, "% -14g")

r_inv = invr(n, r_inv)
for i in 0...n
    for j in 0...i
        r_inv[i][j] = 0
    end
end
printf("$R^-1end$\n")
matprint(r_inv, n, 5, "% -14g")
for i in 0...n
    for j in 0...n
        s = 0
        for k in 0...n
            s += r_inv[i][k] * r[k][j]
        end
        rr[i][j] = s
    end
end

printf("$R^-1end * R$\n")
matprint(rr, n, 5, "% -14g")

for i in 0...n
    for j in 0...i
        s = 0
        for k in 0...n
            s += r[i][k] * r_inv[k][j]
        end
        rr[i][j] = s
    end
end
printf("$R * R^-1end$\n")
matprint(rr, n, 5, "% -14g")

t = 0
for i in 0...n
    for j in 0...n
        if (i == j); s = 1; else; s = 0; end
        for k in 0...n
            s -= r[i][k] * r_inv[k][j]
        end
        t += s * s
    end
end
printf("$R R^-1end$ の成分の二乗平均誤差 %g\n", Math::sqrt(t / (n * n)))

t = 0
for i in 0...n
    for j in 0...n
        if (i == j); s = 1; else; s = 0; end
        for k in 0...n
            s -= r_inv[i][k] * r[k][j]
        end
        t += s * s
    end
end
printf("$R^-1end R$ の成分の二乗平均誤差 %g\n", Math::sqrt(t / (n * n)))
exit 0