#/***********************************************************
#    gauss.rb -- Gauss (ガウス) 法
#***********************************************************/
require "matutil.rb"

def gauss( n, a)
    for k in 0...(n-1)
        for i in (k+1)...n
            t = a[i][k] / a[k][k]
            for j in (k+1)..n
                a[i][j] -= t * a[k][j]
            end
        end
    end
    (n - 1).downto(0) do |i|
        t = a[i][n]
        for j in (i+1)...n; t -= a[i][j] * a[j][n]; end
        a[i][n] = t / a[i][i]
    end
end

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

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 + 1)
b = new_matrix(n, n + 1)
for i in 0...n
    for j in 0...n
        a[i][j] = b[i][j] = rnd - rnd
    end
end
for i in 0...n
    a[i][n] = b[i][n] = rnd - rnd
end
printf("係数行列 (右辺も含む)\n")
matprint(a, n + 1, 10, "%7.3f")
gauss(n, a)
printf("解と, 解を代入したときの両辺の差\n")
for i in 0...n
    s = b[i][n]
    for j in 0...n; s -= b[i][j] * a[j][n]; end
    printf("%4d: %12.7f %12.7f\n", i, a[i][n], s)
end

exit 0