def forward(A, b):
    n = len(b)
    x = np.zeros(n)

    for i in range(n):
        x[i] = (b[i] - np.dot(A[i, :i], x[:i])) / A[i, i]

    return x



def gaussseidel(A, b, x0, maxiter=10**3, tol=1e-6):
    x = x0.copy()
    B = np.tril(A)
    d = b - A @ x
    k = 0

    while True:
        x = x + forward(B, d)
        d = b - A @ x
        k += 1

        if np.linalg.norm(d) < tol or k > maxiter:
            break

    return x, k