import numpy as np

g = np.array([[1, 2, 1], [2, 1, 1], [1, 1, 1]])
f = np.array([[3, 4, 6], [5, 2, 4], [3, 3, 3]])

matr = np.concatenate([g.T, f.T], axis=1)

# matr = np.array([
#     [4, 2, 1, 2, 1, 1],
#     [1, 1, 1, 3, -1, 3],
#     [3, 2, 1, 3, 1, 2]
# ])

n, m = matr.shape
aug = 3


def swap(i, j):
    matr[i], matr[j] = matr[j].copy(), matr[i].copy()
    return f"\\vec{{ {i + 1} }} \\leftrightarrow \\vec {{ {j + 1} }}"


def multiply(i, k):
    matr[i] *= k
    return f"\\vec{{ {i + 1} }} \\to \\cdot {k}"


def divide(i, k):
    matr[i] //= k
    return f"\\vec{{ {i + 1} }} \\to /{k}"


def first_nonzero(i):
    for k in range(m):
        if matr[i, k] != 0:
            return k
    return -1


def subtract(i, reverse=False):
    r = first_nonzero(i)
    gen = range(i - 1, -1, -1) if reverse else range(i + 1, n)
    for j in gen:
        k = first_nonzero(i) if reverse else first_nonzero(j)
        if k == -1 or k > r:
            continue
        if matr[j, k] % matr[i, r] != 0:
            raise ValueError("Not divisible")
        matr[j] -= (matr[j, k] // matr[i, r]) * matr[i]
    return f"-k \\cdot \\vec {{ {i + 1} }}"


def matrix_to_latex(linebreak):
    latex = ""
    if not linebreak:
        latex += "&"
    else:
        latex += "&&"
    latex += "\\begin{bmatrix}\n\t" if aug == 0 else f"\\begin{{abmatrix}}[{m - aug}]{{{aug}}}\n\t"
    latex += "\\\\\n\t".join([" & ".join([str(k) for k in j]) for j in matr])
    latex += "\n\\end{bmatrix}" if aug == 0 else "\n\\end{abmatrix}"
    if linebreak:
        latex += "\\\\"
    return latex


def print_op(latex, linebreak):
    if linebreak:
        print("&&", end='')
    print(f"\\underset {{ {latex} }} {{\\sim}}\n{matrix_to_latex(linebreak)}")
    return not linebreak


def gauss_step(normalized_lines, linebreak):
    for i in range(n):
        k = first_nonzero(i)
        if k != -1 and matr[i, k] < 0:
            linebreak = print_op(multiply(i, -1), linebreak)

    nonzero_left = min(first_nonzero(i) for i in range(normalized_lines, n))
    line_to_raise = min(((matr[i, nonzero_left], i)
                        for i in range(normalized_lines, n) if matr[i, nonzero_left] != 0), default=(-1, -1))[1]
    if line_to_raise == -1:
        return n
    if line_to_raise != normalized_lines:
        linebreak = print_op(swap(normalized_lines, line_to_raise), linebreak)
    main_line = matr[normalized_lines]
    main_elem = main_line[nonzero_left]
    for i in range(normalized_lines + 1, n):
        cur_elem = matr[i, first_nonzero(i)]
        k = np.lcm(main_elem, cur_elem) // cur_elem
        if k != 1:
            linebreak = print_op(multiply(i, k), linebreak)
    linebreak = print_op(subtract(normalized_lines), linebreak)
    return 1, linebreak


def normalize_step(linebreak):
    for i in range(0, n):
        k = first_nonzero(i)
        gg = matr[i, k]
        for e in range(k, m):
            gg = np.gcd(gg, matr[i, e])
        if gg != 1:
            linebreak = print_op(divide(i, gg), linebreak)
            return True, linebreak
    for i in range(n-1, -1, -1):
        k = first_nonzero(i)

        if any(matr[e, k] != 0 for e in range(i - 1, -1, -1)):
            for e in range(i - 1, -1, -1):
                ll = np.lcm(matr[i, k], matr[e, k])
                if ll != matr[e, k]:
                    linebreak = print_op(multiply(e, ll // matr[e, k]), linebreak)
            linebreak = print_op(subtract(i, True), linebreak)
            return True, linebreak
    return False, linebreak


def gauss():
    print(matrix_to_latex(False))
    linebreak = True
    normalized_lines = 0
    while normalized_lines < n:
        lines_diff, linebreak_new = gauss_step(normalized_lines, linebreak)
        normalized_lines += lines_diff
        linebreak = linebreak_new
    return linebreak


def normalize(linebreak):
    while True:
        diff, linebreak_new = normalize_step(linebreak)
        if not diff:
            break
        linebreak = linebreak_new


def main():
    normalize(gauss())


if __name__ == '__main__':
    main()