#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
typedef unsigned int nd; // 32 bits
typedef double R;
#define N 100000
#define C const
nd sp[N][8]; // mutable space for blob as octree
C nd E = 0, F = 255;
nd f = 256;

R sq(R x){return x*x;}
// R max(R x, R y){return x<y?y:x;}
R sh(R x){return sq(fabs(x)+1);}
R sg(R x){return sq(fabs(x)-1);} // wrong
nd buy(){nd t = f; f = sp[f][0]; if(t) return t;
  exit(printf("empty\n"));}
void sell(nd n){sp[n][0] = f; f = n;}
nd ball(R x, R y, R z, R r, int d){
  if(fabs(x) > 1+r || fabs(y) > 1+r || fabs(z) > 1+r) return E; // beyond some face
  if(sh(x)+sh(y)+sh(z)<r*r) return F; // cube entirely in ball
  if(fabs(x)>1 || fabs(y)>1 || fabs(z)>1){ // center not in cube
    if(  fabs(x)<1 && sg(y)+sg(z) > r*r // between yz planes but outside rounded box
      || fabs(y)<1 && sg(x)+sg(z) > r*r
      || fabs(z)<1 && sg(x)+sg(y) > r*r) return E; 
    if(sg(x)+sg(y)+sg(z) > r*r) return E;}
  if(d) {nd new = buy(); char pure=1;
    if(0) printf("%11.6f %11.6f %11.6f, r=%11.6f, nd=%d\n", x, y, z, r, new);
    for(int j=0; j<8; ++j) {
      nd n = ball(2*x+(j&4?-1:1), 2*y+(j&2?-1:1), 2*z+(j&1?-1:1), 2*r, d-1);
      pure = pure && (n==F || n==E); sp[new][j] = n;}
    if(pure) {int a=0; for(int j=0; j<8; ++j) a = (a<<1) | (sp[new][j]==F);
      sell(new); return a;}
    return new;}
  return x*x+y*y+z*z<r*r?F:E;}

void pr(int n){for(int i=0; i<n; ++i) {
    for(int j=0; j<8; ++j) printf("%5d", sp[i][j]); printf(" ; %d\n", i);}}

R ex(R x){return x<0?2*x+1:2*x-1;}
int in(R x, R y, R z, nd c) {int v = sp[c][((((x<0)<<1)|(y<0))<<1)|(z<0)];
   if(v==E) return 0; if(v==F) return 1; return in(ex(x), ex(y), ex(z), v);}
void tst(R x, R y, R z, R r){int fn = 6; nd n = ball(x, y, z, r, fn);
    printf("ball at %11.6f %11.6f %11.6f, r=%11.6f, nd=%d\n", x, y, z, r, n);
 {R del = exp2(-fn); int ec=0; R se = 0, me = 0;
   for(R X = -1+del/2; X<1; X += del)
   for(R Y = -1+del/2; Y<1; Y += del)
   for(R Z = -1+del/2; Z<1; Z += del)
     {int w = in(X, Y, Z, n); R s2 = sq(X-x)+sq(Y-y)+sq(Z-z);
     if(!(w == s2 < r*r)) {if(0) {R f = 2/del;
       printf("w%5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %d\n",
       f*X, f*Y, f*Z, f*x, f*y, f*z, w); pr(f);}
       ++ec; {R er = fabs(s2 - r*r); se += er; me = me<er?er:me;}}}
   printf("%d errors, max=%e, sum=%e\n", ec, me, se);}}

  int main(){
  for(nd j = 256; j<N-1; ++j) sp[j][0] = j+1; sp[N-1][0] = 0;
  for(int j=0; j<8; ++j) {sp[0][j] = 0; sp[1][j] = 1;}
// sp[0] is the universal empty cube and sp[1] the full.
  for(int j=0; j<256; ++j) for(int k=0; k<8; ++k)
      sp[j][7-k] = (j>>k)&1?255:0;
// for j to 255 sp[j] is a cube all of whose children are pure.
  if(0) pr(256);
  tst(0, 0, 0, 0.786); pr(f);
  
  return 0;}