#!/usr/bin/env python3.6# Smith normal formimport numpy as npfrom numpy import *import mathdef row_for_var(xs, vi: int):    NEQNS, NVARS = xs.shape    for r in range(vi, NVARS):        if xs[r][vi] != 0: return r# return numbers (ax, ay) such that ax*x - ay*y = 0def elim_scale_factor(x: int, y: int): return (y, x)def smith_normal_form(xs, ys):    NEQNS, NVARS = xs.shape    assert(NEQNS == ys.shape[0])    # eliminate variable 'vi' by finding a row 'r' and then using the row    # to eliminate.    for vi in range(NVARS):        ri = row_for_var(xs, vi)        if ri is None:            return (f"unable to find non-zero row for variable: {vi}")        print(f"-eliminating variable({vi}) using row ({ri}:{xs[ri]})")        # eliminate all other rows using this row        for r2 in range(NEQNS):            # skip the eqn ri.            if r2 == ri: continue            # eliminate.            (scale_ri, scale_r2) = elim_scale_factor(xs[ri][vi], xs[r2][vi])            print(f"-computing xs[{r2}] = {scale_r2}*xs[{r2}]:{xs[r2]} - {scale_ri}*xs[{ri}]:{xs[ri]}")            xs[r2] = scale_r2 * xs[r2] - scale_ri * xs[ri]            ys[r2] = scale_r2 * ys[r2] - scale_ri * ys[ri]        print(f"-state after eliminating variable({vi})")        print(f"xs:\n{xs}\n\nys:{ys}")    sols = [None for _ in range(NVARS)]    for vi in range(NVARS):        r = row_for_var(xs, vi)        if r is None:            print(f"unable to find row for variable {vi}")            return None        assert(xs[r][vi] != 0)        if ys[r] % xs[r][vi] != 0:            print(f"unable to solve eqn for variable {vi}: xs:{xs[r]} = y:{ys[r]}")            return None        else:            sols[vi] = ys[r] // xs[r][vi]    # now check solutions if we have more equations than solutions    for r in range(NEQNS):        lhs = 0        for i in range(NVARS): lhs += xs[r][i] * sols[i]        if lhs != ys[r]:            print(f"-solution vector sols:{sols} cannot satisfy row xs[{r}] = ys[{r}]:{xs[r]} = {ys[i]}")            return None    return sols# consistent system# x = 6, y = 4# x + y = 10# x - y = 2xs = np.asarray([[1, 1], [1, -1]])ys = np.asarray([10, 2])print("## CONSISTENT ##")out = smith_normal_form(xs,ys)print("xs:\n%s\n\nys:\n%s\n\nsoln:\n%s" % (xs, ys, out,))# consistent, over-determined system# x = 6, y = 4# x + y = 10# x - y = 2# x + 2y = 14xs = np.asarray([[1, 1], [1, -1], [1, 2]])ys = np.asarray([10, 2, 14])print("## CONSISTENT OVER DETERMINED ##")out = smith_normal_form(xs,ys)print("xs:\n%s\n\nys:\n%s\n\nsoln:\n%s" % (xs, ys, out,))# consistent, under-determined system# x = 6, y = 4, z = 1# x + y + z = 11# x - y + z = 3xs = np.asarray([[1, 1], [1, -1], [1, 2]])ys = np.asarray([10, 2, 14])print("## CONSISTENT UNDER DETERMINED ##")out = smith_normal_form(xs,ys)print("xs:\n%s\n\nys:\n%s\n\nsoln:\n%s" % (xs, ys, out,))# inconsistent system# x = 6, y = 4# x + y = 10# x - y = 2# x + 2y = 1 <- INCONSTENTxs = np.asarray([[1, 1], [1, -1], [1, 2]])ys = np.asarray([10, 2, 1])print("## INCONSISTENT (OVER DETERMINED) ##")out = smith_normal_form(xs,ys)# inconsistent system# x + y = 10# x - y = 2xs = np.asarray([[1, -1], [1, 2], [1, 1]])ys = np.asarray([10, 2, 1])print("## INCONSISTENT (OVER DETERMINED) ##")out = smith_normal_form(xs,ys)# consistent system over Q, not Z# x = y = 0.5# x + y = 1# x - y = 0xs = np.asarray([[1, 1], [1, -1]])ys = np.asarray([1, 0])print("## INCONSISTENT (SOLVABLE OVER Q NOT Z) ##")out = smith_normal_form(xs,ys)