program Bol;
{
This software has been designed and it is
CopyLefted by Han de Bruijn:

                          (===)
                         @-O^O-@
                          #/_\#
                           ###
}
Uses SysUtils;

type
  punt = record
    x,y,z : double;
  end;
var
  basis : integer; { Refinement factor }
  diepte : integer; { Depth of refinement }
  Vol,Uit : double; { Inner / Outer volume }
  teller : array of integer; { # of cubes }
  xmin,xmax,ymin,ymax,zmin,zmax : double;

procedure Read_Parameters(var basis,diepte : integer);
{
  Syntax: [program] [number base] [refinement depth]
}
const
  maximum : array[2..3] of integer = (10,6);
var
  OK : boolean;
  tel : integer;
  woord : string;
begin
  OK := true;
  tel := ParamCount;
  if not (tel < 3) then OK := false;
  if not OK then Halt;
  basis := 2; diepte := maximum[2];
  if tel = 0 then Exit;

  woord := ParamStr(1);
  if (Length(woord) > 1)
  or not (woord[1] in ['2','3']) then
    Writeln(woord,' out of [2..3] range; default = 2')
  else basis := StrToInt(woord);
  if basis = 3 then diepte := maximum[3];
  if tel = 1 then Exit;

  woord := ParamStr(2);
  OK := not (Length(woord) > 2);
  OK := OK and (woord[1] in ['0'..'9']);
  if Length(woord) = 2 then
    OK := OK and (woord[2] in ['0'..'9']);
  if OK then
    diepte := StrToInt(woord)
  else
    Writeln(woord,' out of range; default chosen');
  if diepte > maximum[basis] then
  begin
    diepte := maximum[basis];
    Writeln(woord,' out of range; default chosen');
  end;
end;

procedure coordinaat(diep,getal : integer; var c0,c1 : double);
{
  Calculate normed coordinates (c) of cubes
  with depth (diep) and given coding (getal)
}
var
  k,macht,g,m : integer;
  c : double;
  rij : array of integer;
begin
  SetLength(rij,diep);
  g := abs(getal);
  for k := 1 to diep do
  begin
    m := g mod basis;
    rij[k-1] := m;
    g := g div basis;
  end;
  macht := 1;
  c := 0;
  for k := diep downto 1 do
  begin
    macht := macht*basis;
    m := rij[k-1];
    c := c + m/macht;
  end;
  c0 := c; c1 := c0 + 1/macht;
end;

function F(x,y,z : double) : double;
{
  Unit Sphere
}
begin
  F := sqr(x)+sqr(y)+sqr(z)-1;
end;

procedure brutaal(d,i,j,k : integer);
{
  Recursive Refinement
}
var
  h : array[0..7] of punt;
  binnen,buiten : boolean;
  m,p,q,r : integer;
  c0,c1,x0,y0,x1,y1,z0,z1 : double;
begin
  if d > diepte then Exit; 
  teller[d] := teller[d] + 1;

{ Calculate cubes coordinates}
  coordinaat(d,i,c0,c1);
  x0 := xmin + (xmax-xmin)*c0;
  x1 := xmin + (xmax-xmin)*c1;
  coordinaat(d,j,c0,c1);
  y0 := ymin + (ymax-ymin)*c0;
  y1 := ymin + (ymax-ymin)*c1;
  coordinaat(d,k,c0,c1);
  z0 := zmin + (zmax-zmin)*c0;
  z1 := zmin + (zmax-zmin)*c1;
  h[0].x := x0; h[0].y := y0; h[0].z := z0;
  h[1].x := x1; h[1].y := y0; h[1].z := z0;
  h[2].x := x0; h[2].y := y1; h[2].z := z0;
  h[3].x := x1; h[3].y := y1; h[3].z := z0;
  h[4].x := x0; h[4].y := y0; h[4].z := z1;
  h[5].x := x1; h[5].y := y0; h[5].z := z1;
  h[6].x := x0; h[6].y := y1; h[6].z := z1;
  h[7].x := x1; h[7].y := y1; h[7].z := z1;

{ Volume inside }
  binnen := true;
  for m := 0 to 7 do
    binnen := binnen and (F(h[m].x,h[m].y,h[m].z) <= 0);
  if binnen then
  begin
    Vol := Vol + (x1-x0)*(y1-y0)*(z1-z0);
    teller[d] := teller[d] - 1;
    Exit;
  end;
{ Area outside }
  buiten := true;
  for m := 0 to 7 do
    buiten := buiten and (F(h[m].x,h[m].y,h[m].z) > 0);
  if buiten then
  begin
    Uit := Uit + (x1-x0)*(y1-y0)*(z1-z0);
    teller[d] := teller[d] - 1;
    Exit;
  end;

{ Recursion }
  for p := 0 to basis-1 do
    for q := 0 to basis-1 do
      for r := 0 to basis-1 do
        brutaal(d+1,basis*i+p,basis*j+q,basis*k+r);
end;

function macht(n : integer) : integer;
{
  Used for (basis^n)
}
var
  p,k : integer;
begin
  p := 1;
  for k := 1 to n do
    p := p*basis;
  macht := p;
end;

procedure just_doit;
{
  Just do it
}
var
  k,p,q,r : integer;
  Rand,dx,dy,dz : double;
begin
  Vol := 0; Uit := 0;
  SetLength(teller,diepte+1);
  for k := 0 to diepte do
    teller[k] := 0;
  for p := 0 to basis-1 do
    for q := 0 to basis-1 do
      for r := 0 to basis-1 do
        brutaal(1,p,q,r);
 
  for k := 0 to diepte do
    Writeln('# cubes at depth = ',k,' : ',teller[k]);
  dx := (xmax-xmin)/macht(diepte);
  dy := (ymax-ymin)/macht(diepte);
  dz := (zmax-zmin)/macht(diepte);
  Rand := teller[diepte]*dx*dy*dz;
  Writeln('Boundary + Inner Volume + Outer Volume = Volume of Region:');
  Writeln(Rand,' +',Vol,' +',Uit,' =');
  Writeln(Rand + Vol + Uit,' =',(xmax-xmin)*(ymax-ymin)*(zmax-zmin));
  Writeln('Final result:');
  Uit := (xmax-xmin)*(ymax-ymin)*(zmax-zmin)-Uit;
  Writeln(Vol:9:6,' < Volume = ',4*Pi/3:9:6,' <',Uit:9:6);
end;

begin
  Read_Parameters(basis,diepte);
{ xmin := 0; xmax := 2;
  ymin := 0; ymax := 2;
  zmin := 0; zmax := 2; }
  xmin := -1.2; xmax := 1.2;
  ymin := -1.2; ymax := 1.2;
  zmin := -1.2; zmax := 1.2;
  just_doit;
end.