§ Computing the smith normal form

#!/usr/bin/env python3.6
# Smith normal form
import numpy as np
from numpy import *
import math

def 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 = 0
def 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 = 2
xs = 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 = 14
xs = 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 = 3
xs = 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 <- INCONSTENT
xs = 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 = 2
xs = 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 = 0
xs = np.asarray([[1, 1], [1, -1]])
ys = np.asarray([1, 0])
print("## INCONSISTENT (SOLVABLE OVER Q NOT Z) ##")
out = smith_normal_form(xs,ys)