typedef double R;
typedef unsigned short ls;
// float over double makes a difference with -O2
static R E = sizeof(R)<8?1.e-7:1.e-15;
typedef struct{R x; R y; R z;} Pt;
typedef struct{ls a; ls b; ls c;} T;
static int const n=2;
static int const VN = 10*n*n+2; // VN vertices
static int const SN = 20*n*n; // SN slivers
ls h[VN]; // heads of queues of slivers
Pt add(Pt, Pt);
Pt sub(Pt, Pt);
Pt sm(R, Pt);
R len(Pt);
void prP(char *, Pt); // print a Pt.
T Q[SN]; // indices for the three bounding vertices
Pt q[VN]; // coordinates of vertices
ls L[3*SN]; // Links for inverted Q.
int sph(void); // Initialize shape and writes h, Q, , q, tm, L.
R det(Pt, Pt, Pt); // Determinant
Pt cp(Pt, Pt); // Cross product
Pt nv(int); // normal vector of factet j.
R ip(Pt, Pt); // inner product
Pt puff(int); // vector sum of normals of facets neighboring vertex j.
Pt centr(int); // centroid of facet j relative to vertices p.
R moments(R m1[3], R m2[3][3]); // Moments:
void pR(R x);
void rm(char * c, R m0, R m1[3], R m2[3][3]);
void eigen(R m[3][3], R l[3], R vc[3][3]);