#include <math.h>
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include "blob.h"
static const char D = 0;

static R t;
static int qx = 0;
static int Qx = 0;

static void pQ(T i){printf("T: %03d %03d %03d\n", i.a, i.b, i.c);}
static void proc(Pt z){if(0) prP("h", z); assert(qx<VN); q[qx++] = z;}
static void w(R x, R y, R z){proc((Pt){x, y, z});
    assert(fabs(x*x + y*y + z*z - (1+t*t)) < E);}
static void v(R o, R t){w(o, t, 0); w(t, 0, o); w(0, o, t);}
static void x(R o){v(o, t); v(o, -t);}

// static R ip(int i, int j){return q[i].x*q[j].x + q[i].y*q[j].y + q[i].z*q[j].z;}
static R Ip(int i, int j){return ip(q[i], q[j]);}
static int ln(int i, int j){return (fabs(Ip(i, j) - (1 - t*t)) < E);} //sph
static int lx(int i, int j){Pt d = sub(q[i], q[j]); // sph
   return fabs(d.x*d.x + d.y*d.y + d.z*d.z - 4*t*t/(n*n)) < E;}
int sph(){t = 0.5*(sqrt(5) - 1);
  assert(3*SN < 0xfffe);
  printf("wwrr %x %x %x\n", SN, 3*SN, 0xfffe);
  x(1); x(-1); assert (qx==12);
  if(0) for(int j=0; j<12; ++j)
     printf("%2d %10.7f %10.7f %10.7f\n", j, q[j].x, q[j].y, q[j].z);
  {int m=0; int vc[12]; for(int k=0; k<12; ++k) vc[k]=0;
    for(int i=0; i<12; ++i) for(int j=i+1; j<12; ++j)
      if(ln(i, j)) {Pt d = sm(1./n, sub(q[j], q[i]));
        {Pt x = q[i]; for(int m=1; m<n; ++m) {x = add(x, d); proc(x);}}
        for(int k=j+1; k<12; ++k) if(ln(i, k) && ln(j, k)) {++m;
        {Pt e = sm(1./n, sub(q[k], q[i]));
          {Pt x = q[i]; for(int m=1; m<n; ++m) {x = add(x, d);
            Pt y = x; for(int N=m+1; N<n; ++N) {y = add(y, e); proc(y);}}}}
        ++vc[i]; ++vc[j]; ++vc[k];}}
    assert(m==20); for(int y=0; y<12; ++y) assert(vc[y] == 5);}
    assert(qx == VN);
    for(int j=0; j<VN; ++j) h[j] = -1;
    for(int i=0; i<VN; ++i) for(int j=i+1; j<VN; ++j)
      if(lx(i, j)) for(int k=j+1; k<VN; ++k) if(lx(i, k) && lx(j, k)) {
        R v = det(q[i], q[j], q[k]);
        assert(Qx<SN); assert(fabs(fabs(v) - ((sqrt(5)-1)/(n*n))) < E);
        T I = (v<0) ? (T){i, k, j}:(T){i, j, k};
        {ls w = h[I.a]; h[I.a] = 3*Qx  ; L[3*Qx  ] = w;}
        {ls w = h[I.b]; h[I.b] = 3*Qx+1; L[3*Qx+1] = w;}
        {ls w = h[I.c]; h[I.c] = 3*Qx+2; L[3*Qx+2] = w;}
        Q[Qx++] = I;}
if(D) for(int i=0; i<SN; ++i) {
   printf("%3d %3d L: %3d %3d %3d, ", i, 3*i, L[3*i], L[3*i+1], L[3*i+2]); pQ(Q[i]);}
if(D) for(int i=0; i<VN; ++i) {printf("%3d h=%3d ", i, h[i]); prP("vtx", q[i]);}
printf("qx = %d, VN = %d\n", qx, VN);
printf("Qx = %d, SN = %d\n", Qx, SN);
for(int k = 0; k<SN; ++k)
  {R y = det(q[Q[k].a], q[Q[k].b], q[Q[k].c]) - 2.*t/(n*n);
    if(fabs(y) > E) exit(printf("%e %d wqqw\n", y, k));
  assert(fabs(det(q[Q[k].a], q[Q[k].b], q[Q[k].c]) - 2.*t/(n*n)) < E);}
for(int k = 0; k<VN; ++k) for(ls c = h[k]; c<(ls)-1; c = L[c]) {int d = c/3, r = c%3;
  assert(k == Q[d].a && r==0 || k == Q[d].b && r==1 || k == Q[d].c && r == 2);}
for(int j=0; j<VN; ++j) q[j] = sm(1./sqrt(Ip(j, j)), q[j]);
{R tm = 0; for(int j=0; j<SN; ++j) tm += det(q[Q[j].a], q[Q[j].b], q[Q[j].c])/6;
  printf("total mass = %14.11f, 4/3 π = %14.11f\n", tm, (4./3.)*(4.*atan(1)));}
return n;}