import numpy as np
from scipy.linalg import solve


A = np.array(
    [
        [3.0, 1.0, -1.0],
        [2.0, 4.0, 1.0],
        [-1.0, 2.0, 5.0],
    ]
)
b = np.array([4.0, 1.0, 1.0])

x = solve(A, b)
reconstructed = A @ x
residual = reconstructed - b
residual_norm = np.linalg.norm(residual, ord=np.inf)
condition_number = np.linalg.cond(A)

np.set_printoptions(precision=6, suppress=True)
print("matrix A:")
for row in A:
    print(row)
print("right hand side:", b)
print("solution x:", x)
print("A @ x:", reconstructed)
print("residual:", residual)
print(f"residual_inf_norm: {residual_norm:.2e}")
print(f"condition_number: {condition_number:.2f}")
print("verified:", bool(residual_norm < 1e-10))