from fractions import Fraction as F
import itertools, json
N=[(19,70,-410),(0,289,-186),(-26,-389,-2),(-6,207,-331),(13,-329,101),(1,-167,233),(-5,-263,181),(0,-73,253),(-27,-85,-443),(0,107,201),(0,181,103),(0,287,-27),(32,-394,-134),(0,18,244)]
def dot(a,b): return sum(x*y for x,y in zip(a,b))
def sub(a,b): return tuple(x-y for x,y in zip(a,b))
def crs(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])
c=[dot(n,n) for n in N]
def solve(i,j,k):
    A=[N[i],N[j],N[k]]; b=[c[i],c[j],c[k]]
    d=dot(A[0],crs(A[1],A[2]))
    if d==0: return None
    # Cramer
    x=F(dot((b[0],A[0][1],A[0][2]),crs((b[1],A[1][1],A[1][2]),(b[2],A[2][1],A[2][2]))),d)
    def det(M): return dot(M[0],crs(M[1],M[2]))
    Mx=[(b[r],A[r][1],A[r][2]) for r in range(3)]
    My=[(A[r][0],b[r],A[r][2]) for r in range(3)]
    Mz=[(A[r][0],A[r][1],b[r]) for r in range(3)]
    return (F(det(Mx),d),F(det(My),d),F(det(Mz),d))
V={}
for i,j,k in itertools.combinations(range(14),3):
    p=solve(i,j,k)
    if p is None: continue
    if all(dot(N[m],p)<=c[m] for m in range(14)):
        V.setdefault(p,set()).update((i,j,k))
print('vertices',len(V), 'max planes at a vertex',max(len(s) for s in V.values()))
verts=list(V)
faces=[]
for f in range(14):
    pts=[p for p in verts if f in V[p]]
    # order around centroid in plane
    cx=tuple(sum(p[t] for p in pts)/len(pts) for t in range(3))
    n=N[f]; u=sub(pts[0],cx); w=crs(n,u)
    import math
    pts.sort(key=lambda p: math.atan2(float(dot(sub(p,cx),w)),float(dot(sub(p,cx),u))))
    faces.append(pts)
print('face sizes',[len(f) for f in faces])
# volume centroid via tetrahedra from interior point O
O=tuple(sum(p[t] for p in verts)/len(verts) for t in range(3))
vol=F(0); G=[F(0)]*3
for f in faces:
    for a in range(1,len(f)-1):
        A,B,C=f[0],f[a],f[a+1]
        v=abs(dot(sub(A,O),crs(sub(B,O),sub(C,O))))/6
        vol+=v
        for t in range(3): G[t]+=v*(O[t]+A[t]+B[t]+C[t])/4
G=tuple(g/vol for g in G)
print('G',[float(g) for g in G])
# stable faces: projection of G inside face; margin = min signed distance to edges / scale
stable=[]
R=max(math.sqrt(float(dot(sub(p,G),sub(p,G)))) for p in verts)
for fi,f in enumerate(faces):
    n=N[fi]; nn=dot(n,n)
    t=(c[fi]-dot(n,G))/nn
    P=tuple(G[k]+t*n[k] for k in range(3))
    sg=[]
    for a in range(len(f)):
        A,B=f[a],f[(a+1)%len(f)]
        s=dot(crs(sub(B,A),sub(P,A)),n)
        L=math.sqrt(float(dot(sub(B,A),sub(B,A))))*math.sqrt(float(nn))
        sg.append(float(s)/L)
    inside=all(x>0 for x in sg) or all(x<0 for x in sg)
    m=min(abs(x) for x in sg)
    stable.append((fi,inside,m/R))
print([s for s in stable if s[1]])
print('R',R,'dims',[ (float(min(p[k] for p in verts)),float(max(p[k] for p in verts))) for k in range(3)])
json.dump({'normals':N,'vertices':[[float(x) for x in p] for p in verts],
 'faces':[[verts.index(p) for p in f] for f in faces],'G':[float(g) for g in G],
 'stable':[s[0] for s in stable if s[1]],'margins':{s[0]:s[2] for s in stable}},open('die.json','w'))