## ยง 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)