#!/usr/bin/env julia

# Jacek Poziemski 272389

using LinearAlgebra

include("hilb.jl")
include("matcond.jl")

"""
Errors of the solutions for a system of linear equations
for the Gauss and inverse matrix methods.

# Fields
- `gauss::Float64`: the error for the Gauss method
- `inverse::Float64`: the error for the inverse matrix method
"""
struct Errors
    gauss::Float64
    inverse::Float64
end

"""
    solve(A::Matrix{Float64}, n::Int)::Errors 

Solve the system of linear equations given by the matrix `A`.

# Arguments
- `A::Matrix{Float64}`: the matrix
- `n::Int`: size of the matrix
"""
function solve(A::Matrix{Float64}, n::Int)::Errors
    x::Vector{Float64} = ones(Float64, n)
    b::Vector{Float64} = A * x
    gauss::Float64 = norm((A \ b) - x) / norm(x)
    inverse::Float64 = norm((inv(A) * b) - x) / norm(x)

    return Errors(gauss, inverse)
end

println("hilbert matrices")
println()

for n in 1:30
    A::Matrix{Float64} = hilb(n)
    errors = solve(A, n)
    println("n: $n; cond: $(cond(A)); rank: $(rank(A)); gauss error: $(errors.gauss); inverse error: $(errors.inverse)")
end

println()
println("random matrices")
println()

for n in [5, 10, 20]
    for c in [0, 1, 3, 7, 12, 16]
        A::Matrix{Float64} = matcond(n, 10.0^c)
        errors = solve(A, n)
        println("n: $n; c: $c; cond: $(cond(A)); rank: $(rank(A)); gauss error: $(errors.gauss); inverse error: $(errors.inverse)")
    end
end