# Tetrahedrons and Spheres

Tetrahedrons and Spheres

There are a tetrahedrons and b spheres in the 3D-splace, you’re asked to calculate the volume occupied by at least one of them (i.e. volume of the union of the objects).

# Input

There will be at most 20 test cases. Each case begins with two integers a, b, the number of tetrahedrons and the number of spheres (1<=a,b<=5). The next a lines each contains 12 integers: x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, the coordinates (xi, yi, zi)(1<=i<=4) of the four vertices  of a tetrahedron. The next b lines each contains 4 integers x, y, z, r, the coordinates of the center (x, y, z) and the radius r (r<=3). All the coordinate values are integers with absolute values no more than 5. The input is terminated by a=b=0.

# Output

For each test case, print a single line, the volume occupied by at least one of them, rounded to three decimal points.

``````1 1
0 0 4 1 0 4 0 1 4 0 0 5
0 0 0 1
0 0
``````

``````4.356
``````

``````#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#include <vector>
using namespace std;
#define sqr(x) ((x)*(x))

const double eps=1e-6;
const double inf=8;

const double pi=3.14159265358979323846;

inline bool cmpDouble(const double &a,const double &b)
{return fabs(a-b)<eps;}
struct Tpoint
{
double x,y;
Tpoint(){}
Tpoint(double a,double b){x=a;y=b;}
inline double norm(){ return sqrt( sqr(x)+sqr(y) ); }
};

inline Tpoint operator +(const Tpoint &a,const Tpoint &b){return Tpoint(a.x+b.x,a.y+b.y);}
inline Tpoint operator -(const Tpoint &a,const Tpoint &b){return Tpoint(a.x-b.x,a.y-b.y);}
inline Tpoint operator *(const double &a,const Tpoint &b){return Tpoint(a*b.x,a*b.y);}
inline Tpoint operator *(const Tpoint &b,const double &a){return Tpoint(a*b.x,a*b.y);}
inline Tpoint operator /(const Tpoint &a,const double &b){return Tpoint(a.x/b,a.y/b);}
inline bool operator <(const Tpoint &a,const Tpoint &b){return a.x+eps<b.x || fabs(a.x-b.x)<eps && a.y+eps<b.y;}
inline bool operator ==(const Tpoint &a,const Tpoint &b){return fabs(a.x-b.x)<eps && fabs(a.y-b.y)<eps;}
inline double det(const Tpoint &a,const Tpoint &b){return a.x*b.y-a.y*b.x;}
inline double dot(const Tpoint &a,const Tpoint &b){return a.x*b.x+a.y*b.y;}

struct P3
{
double x,y,z;

P3(){}
P3(double a,double b,double c){x=a;y=b;z=c;}
};
inline P3 operator +(const P3 &a,const P3 &b){return P3(a.x+b.x,a.y+b.y,a.z+b.z);}
inline P3 operator -(const P3 &a,const P3 &b){return P3(a.x-b.x,a.y-b.y,a.z-b.z);}
inline P3 operator *(const double &a,const P3 &b){return P3(a*b.x,a*b.y,a*b.z);}
inline P3 operator *(const P3 &b,const double &a){return P3(a*b.x,a*b.y,a*b.z);}
inline P3 operator /(const P3 &a,const double &b){return P3(a.x/b,a.y/b,a.z/b);}

struct Tcir
{
double r;
Tpoint o;

Tcir(){}
Tcir(Tpoint x,double y){o=x,r=y;}
};
vector <Tcir> Circle;
typedef vector <Tpoint> TP;
vector <TP> Poly;

struct Tinter
{
double x,y,Area,mid;
int delta;
Tinter(){}
Tinter(double xx,double yy,double mm,int dd,double aa)
{
x=xx;y=yy;mid=mm;
delta=dd;Area=aa;
}
};
inline bool operator <(const Tinter &a,const Tinter &b){return a.mid>b.mid+eps;}
inline bool operator ==(const Tinter &a,const Tinter &b){return fabs(a.mid-b.mid)<eps;}
vector <Tinter> inter;
vector <double> bak;

inline double dist(const Tpoint &a,const Tpoint &b)
{
return sqr(a.x-b.x)+sqr(a.y-b.y);
}

{
bak.push_back(x);
}

inline void CircleIntersectCircle(const Tcir &a,const Tcir &b)
{
double l=dist(a.o,b.o);
double s=((a.r-b.r)*(a.r+b.r)/l+1)/2;
double t=sqrt(-(l-sqr(a.r+b.r))*(l-sqr(a.r-b.r))/(l*l*4));
double ux=b.o.x-a.o.x,uy=b.o.y-a.o.y;
double ix=a.o.x+s*ux+t*uy,iy=a.o.y+s*uy-t*ux;
double jx=a.o.x+s*ux-t*uy,jy=a.o.y+s*uy+t*ux;
}

inline bool InLine(const Tpoint &a,const Tpoint &b,const Tpoint &c)
{
return fabs(det(b-a,c-a))<eps && dot(a-c,b-c)<eps;
}

inline void LineToLine(const Tpoint &a,const Tpoint &b,const Tpoint &c,const Tpoint &d)
{
double s1=det(c-a,b-a),s2=det(b-a,d-a);
if (s1*s2<-eps) return;
Tpoint e=c+(d-c)*s1/(s1+s2);

if (InLine(a,b,e) && InLine(c,d,e))
{
}
}

inline void LineToCircle(const Tpoint &a,const Tpoint &b,const Tcir &c)
{
double h=fabs(det(c.o-a,b-a))/(b-a).norm();
if (h>c.r+eps) return;
double lamda=dot(c.o-a,b-a);
lamda/=dist(a,b);
Tpoint x=a+(b-a)*lamda;
double d=sqrt( sqr(c.r)-sqr(h) );
d/=(b-a).norm();
Tpoint e=x+(b-a)*d;
Tpoint f=x-(b-a)*d;
if (InLine(a,b,e))
if (InLine(a,b,f))
return;
}

inline void CircleToCircle()
{
for (int i=0;i<Circle.size();++i)
{
for (int j=i+1;j<Circle.size();++j)
if (dist(Circle[i].o,Circle[j].o)<=sqr(Circle[i].r+Circle[j].r))
if (dist(Circle[i].o,Circle[j].o)>=sqr(Circle[i].r-Circle[j].r))
CircleIntersectCircle(Circle[i],Circle[j]);
}
}

inline void CircleToPoly()
{
for (int i=0;i<Circle.size();++i)
for (int j=0;j<Poly.size();++j)
for (int v=0;v<Poly[j].size();++v)
LineToCircle(Poly[j][v],Poly[j][(v+1)%Poly[j].size()],Circle[i]);
}

inline void PolyToPoly()
{
for (int i=0;i<Poly.size();++i)
for (int u=0;u<Poly[i].size();++u)
for (int i=0;i<Poly.size();++i)
for (int j=i+1;j<Poly.size();++j)
for (int u=0;u<Poly[i].size();++u)
for (int v=0;v<Poly[j].size();++v)
LineToLine(Poly[i][u],Poly[i][(u+1)%Poly[i].size()],Poly[j][v],Poly[j][(v+1)%Poly[j].size()]);
}

inline void Get(const Tcir &c,double x,double &l,double &r)
{
double y=fabs(c.o.x-x);
double d=sqrt(fabs( sqr(c.r)-sqr(y) ));
l=c.o.y+d;
r=c.o.y-d;
}

inline double arcArea(const Tcir &a,double l,double x,double r,double y)
{
double len=sqrt(sqr(l-r) + sqr(x-y));
double d=sqrt(sqr(a.r)-sqr(len)/4.0);
double angle=atan(len/2.0/d);
return fabs(angle*sqr(a.r)-d*len/2.0);
}

inline void Get_Interval(const Tcir &a,double l,double r)
{
double L1,L2,R1,R2,M1,M2;
Get(a,l,L1,L2);
Get(a,r,R1,R2);
Get(a,(l+r)/2.0,M1,M2);
int D1=1,D2=-1;
double A1=arcArea(a,l,L1,r,R1),A2=arcArea(a,l,L2,r,R2);
inter.push_back( Tinter(L1,R1,M1,D1,A1) );
inter.push_back( Tinter(L2,R2,M2,D2,A2) );
}

inline double calcSlice(double xl,double xr)
{
inter.clear();
double lmost=-inf,rmost=inf;
for (int i=0;i<Poly.size();++i)
{
int cc=0;
Tinter I[5];
for (int u=0;u<Poly[i].size();++u)
{
Tpoint x=Poly[i][u];
Tpoint y=Poly[i][(u+1)%Poly[i].size()];
double l=min(x.x,y.x),r=max(x.x,y.x);
if (l<=xl+eps && xr<=r+eps)
{
if (fabs(l-r)<eps) continue;
Tpoint d=y-x;
Tpoint Left=x+d/d.x*(xl-x.x);
Tpoint Right=x+d/d.x*(xr-x.x);
Tpoint Mid=(Left+Right)/2;
I[cc++]=Tinter(Left.y,Right.y,Mid.y,1,0);
}
}
sort(I,I+cc);
if (cc==2)
{
I[1].delta=-1;
inter.push_back(I[0]);
inter.push_back(I[1]);
lmost=max(lmost,I[1].mid);
rmost=min(rmost,I[0].mid);
}
}

for (int i=0;i<Circle.size();++i)
if (fabs(Circle[i].o.x-xl)<Circle[i].r+eps && fabs(Circle[i].o.x-xr)<Circle[i].r+eps)
Get_Interval(Circle[i],xl,xr);

if (!inter.size()) return 0;
double ans=0;
sort(inter.begin(),inter.end());
int cnt=0;
for (int i=0;i<inter.size();++i)
{
if (cnt>0)
{
ans+=(fabs(inter[i-1].x-inter[i].x)+fabs(inter[i-1].y-inter[i].y))*(xr-xl)/2.0;
ans+=inter[i-1].delta*inter[i-1].Area;
ans-=inter[i].delta*inter[i].Area;
}
cnt+=inter[i].delta;
}
return ans;
}

#define maxn 105
int n,m;
struct Poly4
{
P3 p[4];
}a[maxn];
struct Sphere
{
P3 o;
double r;
}b[maxn];

inline bool equal(double a,double b)
{
return fabs(a-b)<eps;
}

inline bool InterSect(const Tpoint &a,const Tpoint &b,const Tpoint &c,const Tpoint &d)
{
double s1=det(b-a,c-a),s2=det(b-a,d-a);
if (s1*s2>-eps) return false;
s1=det(d-c,a-c),s2=det(d-c,b-c);
if (s1*s2>-eps) return false;
return true;
}

inline void ToHull(vector <Tpoint> &a)
{
sort(a.begin(),a.end());
int hull[10],len,limit=1;
hull[len=1]=0;
for (int i=1;i<4;++i)
{
while (len>limit && det(a[hull[len]]-a[hull[len-1]],a[i]-a[hull[len]])>=0) --len;
hull[++len]=i;
}
limit=len;
for (int i=2;i>=0;--i)
{
while (len>limit && det(a[hull[len]]-a[hull[len-1]],a[i]-a[hull[len]])>=0) --len;
hull[++len]=i;
}
vector <Tpoint> b=a;
a.resize(len-1);
for (int i=0;i<len-1;++i)
a[i]=b[hull[i+1]];
}

inline double calcArea(double z)
{
Poly.clear();
Circle.clear();
bak.clear();
for (int i=0;i<n;++i)
{
vector <Tpoint> cross;
for (int j=0;j<4;++j)
for (int k=j+1;k<4;++k)
if (!equal(a[i].p[j].z,a[i].p[k].z))
{
double l=min(a[i].p[j].z,a[i].p[k].z),r=max(a[i].p[j].z,a[i].p[k].z);
if (l<=z+eps && z<=r+eps)
{
P3 d=a[i].p[k]-a[i].p[j];
d=d/d.z;
d=d*(z-a[i].p[j].z);
d=d+a[i].p[j];
Tpoint x(d.x,d.y);
cross.push_back(x);
}
}
sort(cross.begin(),cross.end());
cross.erase(unique(cross.begin(),cross.end()),cross.end());
if (cross.size()>2)
{
if (cross.size()>4){puts("Too Many Points!!!");while (1);}
if (cross.size()==4)
ToHull(cross);
Poly.push_back(cross);
}
}

for (int i=0;i<m;++i)
if (fabs(z-b[i].o.z)+eps<b[i].r)
{
Tpoint o(b[i].o.x,b[i].o.y);
double r=sqrt( sqr(b[i].r)-sqr(z-b[i].o.z) );
Circle.push_back(Tcir(o,r));
}

CircleToCircle();
CircleToPoly();
PolyToPoly();

sort(bak.begin(),bak.end());
double res=0;
for (int i=0;i+1<bak.size();++i)
if (fabs(bak[i+1]-bak[i])>eps)
res+=calcSlice(bak[i],bak[i+1]);
return res;
}

int main()
{
for (;scanf("%d%d",&n,&m)!=EOF && (n+m);)
{
for (int i=0;i<n;++i)
for (int j=0;j<4;++j)

for (int i=0;i<m;++i)

double last=0,Ans=calcArea(-inf)+calcArea(inf);
const int Block=4000;
double h=(inf+inf)/(double)Block;
for (int i=1;i<Block;i+=2)
Ans+=4*calcArea(-inf+i*h);
for (int i=2;i<Block;i+=2)
Ans+=2*calcArea(-inf+i*h);
Ans=Ans*h/3.0;
printf("%.3f\n",Ans);
}
return 0;
}
``````

``````#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#include <vector>
using namespace std;
#define sqr(x) ((x)*(x))

const double eps=1e-6;
const double inf=8;

const double pi=3.14159265358979323846;

inline bool cmpDouble(const double &a,const double &b)
{return fabs(a-b)<eps;}
struct Tpoint
{
double x,y;
Tpoint(){}
Tpoint(double a,double b){x=a;y=b;}
inline double norm(){ return sqrt( sqr(x)+sqr(y) ); }
};

inline Tpoint operator +(const Tpoint &a,const Tpoint &b){return Tpoint(a.x+b.x,a.y+b.y);}
inline Tpoint operator -(const Tpoint &a,const Tpoint &b){return Tpoint(a.x-b.x,a.y-b.y);}
inline Tpoint operator *(const double &a,const Tpoint &b){return Tpoint(a*b.x,a*b.y);}
inline Tpoint operator *(const Tpoint &b,const double &a){return Tpoint(a*b.x,a*b.y);}
inline Tpoint operator /(const Tpoint &a,const double &b){return Tpoint(a.x/b,a.y/b);}
inline bool operator <(const Tpoint &a,const Tpoint &b){return a.x+eps<b.x || fabs(a.x-b.x)<eps && a.y+eps<b.y;}
inline bool operator ==(const Tpoint &a,const Tpoint &b){return fabs(a.x-b.x)<eps && fabs(a.y-b.y)<eps;}
inline double det(const Tpoint &a,const Tpoint &b){return a.x*b.y-a.y*b.x;}
inline double dot(const Tpoint &a,const Tpoint &b){return a.x*b.x+a.y*b.y;}

struct P3
{
double x,y,z;

P3(){}
P3(double a,double b,double c){x=a;y=b;z=c;}
};
inline P3 operator +(const P3 &a,const P3 &b){return P3(a.x+b.x,a.y+b.y,a.z+b.z);}
inline P3 operator -(const P3 &a,const P3 &b){return P3(a.x-b.x,a.y-b.y,a.z-b.z);}
inline P3 operator *(const double &a,const P3 &b){return P3(a*b.x,a*b.y,a*b.z);}
inline P3 operator *(const P3 &b,const double &a){return P3(a*b.x,a*b.y,a*b.z);}
inline P3 operator /(const P3 &a,const double &b){return P3(a.x/b,a.y/b,a.z/b);}

struct Tcir
{
double r;
Tpoint o;

Tcir(){}
Tcir(Tpoint x,double y){o=x,r=y;}
};
vector <Tcir> Circle;
typedef vector <Tpoint> TP;
vector <TP> Poly;

struct Tinter
{
double x,y,Area,mid;
int delta;
Tinter(){}
Tinter(double xx,double yy,double mm,int dd,double aa)
{
x=xx;y=yy;mid=mm;
delta=dd;Area=aa;
}
};
inline bool operator <(const Tinter &a,const Tinter &b){return a.mid>b.mid+eps;}
inline bool operator ==(const Tinter &a,const Tinter &b){return fabs(a.mid-b.mid)<eps;}
vector <Tinter> inter;
vector <double> bak;

inline double dist(const Tpoint &a,const Tpoint &b)
{
return sqr(a.x-b.x)+sqr(a.y-b.y);
}

{
bak.push_back(x);
}

inline void CircleIntersectCircle(const Tcir &a,const Tcir &b)
{
double l=dist(a.o,b.o);
double s=((a.r-b.r)*(a.r+b.r)/l+1)/2;
double t=sqrt(-(l-sqr(a.r+b.r))*(l-sqr(a.r-b.r))/(l*l*4));
double ux=b.o.x-a.o.x,uy=b.o.y-a.o.y;
double ix=a.o.x+s*ux+t*uy,iy=a.o.y+s*uy-t*ux;
double jx=a.o.x+s*ux-t*uy,jy=a.o.y+s*uy+t*ux;
}

inline bool InLine(const Tpoint &a,const Tpoint &b,const Tpoint &c)
{
return fabs(det(b-a,c-a))<eps && dot(a-c,b-c)<eps;
}

inline void LineToLine(const Tpoint &a,const Tpoint &b,const Tpoint &c,const Tpoint &d)
{
double s1=det(c-a,b-a),s2=det(b-a,d-a);
if (s1*s2<-eps) return;
Tpoint e=c+(d-c)*s1/(s1+s2);

if (InLine(a,b,e) && InLine(c,d,e))
{
}
}

inline void LineToCircle(const Tpoint &a,const Tpoint &b,const Tcir &c)
{
double h=fabs(det(c.o-a,b-a))/(b-a).norm();
if (h>c.r+eps) return;
double lamda=dot(c.o-a,b-a);
lamda/=dist(a,b);
Tpoint x=a+(b-a)*lamda;
double d=sqrt( sqr(c.r)-sqr(h) );
d/=(b-a).norm();
Tpoint e=x+(b-a)*d;
Tpoint f=x-(b-a)*d;
if (InLine(a,b,e))
if (InLine(a,b,f))
return;
}

inline void CircleToCircle()
{
for (int i=0;i<Circle.size();++i)
{
for (int j=i+1;j<Circle.size();++j)
if (dist(Circle[i].o,Circle[j].o)<=sqr(Circle[i].r+Circle[j].r))
if (dist(Circle[i].o,Circle[j].o)>=sqr(Circle[i].r-Circle[j].r))
CircleIntersectCircle(Circle[i],Circle[j]);
}
}

inline void CircleToPoly()
{
for (int i=0;i<Circle.size();++i)
for (int j=0;j<Poly.size();++j)
for (int v=0;v<Poly[j].size();++v)
LineToCircle(Poly[j][v],Poly[j][(v+1)%Poly[j].size()],Circle[i]);
}

inline void PolyToPoly()
{
for (int i=0;i<Poly.size();++i)
for (int u=0;u<Poly[i].size();++u)
for (int i=0;i<Poly.size();++i)
for (int j=i+1;j<Poly.size();++j)
for (int u=0;u<Poly[i].size();++u)
for (int v=0;v<Poly[j].size();++v)
LineToLine(Poly[i][u],Poly[i][(u+1)%Poly[i].size()],Poly[j][v],Poly[j][(v+1)%Poly[j].size()]);
}

inline void Get(const Tcir &c,double x,double &l,double &r)
{
double y=fabs(c.o.x-x);
double d=sqrt(fabs( sqr(c.r)-sqr(y) ));
l=c.o.y+d;
r=c.o.y-d;
}

inline double arcArea(const Tcir &a,double l,double x,double r,double y)
{
double len=sqrt(sqr(l-r) + sqr(x-y));
double d=sqrt(sqr(a.r)-sqr(len)/4.0);
double angle=atan(len/2.0/d);
return fabs(angle*sqr(a.r)-d*len/2.0);
}

inline void Get_Interval(const Tcir &a,double l,double r)
{
double L1,L2,R1,R2,M1,M2;
Get(a,l,L1,L2);
Get(a,r,R1,R2);
Get(a,(l+r)/2.0,M1,M2);
int D1=1,D2=-1;
double A1=arcArea(a,l,L1,r,R1),A2=arcArea(a,l,L2,r,R2);
inter.push_back( Tinter(L1,R1,M1,D1,A1) );
inter.push_back( Tinter(L2,R2,M2,D2,A2) );
}

inline double calcSlice(double xl,double xr)
{
inter.clear();
double lmost=-inf,rmost=inf;
for (int i=0;i<Poly.size();++i)
{
int cc=0;
Tinter I[5];
for (int u=0;u<Poly[i].size();++u)
{
Tpoint x=Poly[i][u];
Tpoint y=Poly[i][(u+1)%Poly[i].size()];
double l=min(x.x,y.x),r=max(x.x,y.x);
if (l<=xl+eps && xr<=r+eps)
{
if (fabs(l-r)<eps) continue;
Tpoint d=y-x;
Tpoint Left=x+d/d.x*(xl-x.x);
Tpoint Right=x+d/d.x*(xr-x.x);
Tpoint Mid=(Left+Right)/2;
I[cc++]=Tinter(Left.y,Right.y,Mid.y,1,0);
}
}
sort(I,I+cc);
if (cc==2)
{
I[1].delta=-1;
inter.push_back(I[0]);
inter.push_back(I[1]);
lmost=max(lmost,I[1].mid);
rmost=min(rmost,I[0].mid);
}
}

for (int i=0;i<Circle.size();++i)
if (fabs(Circle[i].o.x-xl)<Circle[i].r+eps && fabs(Circle[i].o.x-xr)<Circle[i].r+eps)
Get_Interval(Circle[i],xl,xr);

if (!inter.size()) return 0;
double ans=0;
sort(inter.begin(),inter.end());
int cnt=0;
for (int i=0;i<inter.size();++i)
{
if (cnt>0)
{
ans+=(fabs(inter[i-1].x-inter[i].x)+fabs(inter[i-1].y-inter[i].y))*(xr-xl)/2.0;
ans+=inter[i-1].delta*inter[i-1].Area;
ans-=inter[i].delta*inter[i].Area;
}
cnt+=inter[i].delta;
}
return ans;
}

#define maxn 105
int n,m;
struct Poly4
{
P3 p[4];
}a[maxn];
struct Sphere
{
P3 o;
double r;
}b[maxn];

inline bool equal(double a,double b)
{
return fabs(a-b)<eps;
}

inline bool InterSect(const Tpoint &a,const Tpoint &b,const Tpoint &c,const Tpoint &d)
{
double s1=det(b-a,c-a),s2=det(b-a,d-a);
if (s1*s2>-eps) return false;
s1=det(d-c,a-c),s2=det(d-c,b-c);
if (s1*s2>-eps) return false;
return true;
}

inline void ToHull(vector <Tpoint> &a)
{
sort(a.begin(),a.end());
int hull[10],len,limit=1;
hull[len=1]=0;
for (int i=1;i<4;++i)
{
while (len>limit && det(a[hull[len]]-a[hull[len-1]],a[i]-a[hull[len]])>=0) --len;
hull[++len]=i;
}
limit=len;
for (int i=2;i>=0;--i)
{
while (len>limit && det(a[hull[len]]-a[hull[len-1]],a[i]-a[hull[len]])>=0) --len;
hull[++len]=i;
}
vector <Tpoint> b=a;
a.resize(len-1);
for (int i=0;i<len-1;++i)
a[i]=b[hull[i+1]];
}

inline double calcArea(double z)
{
Poly.clear();
Circle.clear();
bak.clear();
for (int i=0;i<n;++i)
{
vector <Tpoint> cross;
for (int j=0;j<4;++j)
for (int k=j+1;k<4;++k)
if (!equal(a[i].p[j].z,a[i].p[k].z))
{
double l=min(a[i].p[j].z,a[i].p[k].z),r=max(a[i].p[j].z,a[i].p[k].z);
if (l<=z+eps && z<=r+eps)
{
P3 d=a[i].p[k]-a[i].p[j];
d=d/d.z;
d=d*(z-a[i].p[j].z);
d=d+a[i].p[j];
Tpoint x(d.x,d.y);
cross.push_back(x);
}
}
sort(cross.begin(),cross.end());
cross.erase(unique(cross.begin(),cross.end()),cross.end());
if (cross.size()>2)
{
if (cross.size()>4){puts("Too Many Points!!!");while (1);}
if (cross.size()==4)
ToHull(cross);
Poly.push_back(cross);
}
}

for (int i=0;i<m;++i)
if (fabs(z-b[i].o.z)+eps<b[i].r)
{
Tpoint o(b[i].o.x,b[i].o.y);
double r=sqrt( sqr(b[i].r)-sqr(z-b[i].o.z) );
Circle.push_back(Tcir(o,r));
}

CircleToCircle();
CircleToPoly();
PolyToPoly();

sort(bak.begin(),bak.end());
double res=0;
for (int i=0;i+1<bak.size();++i)
if (fabs(bak[i+1]-bak[i])>eps)
res+=calcSlice(bak[i],bak[i+1]);
return res;
}

int main()
{
for (;scanf("%d%d",&n,&m)!=EOF && (n+m);)
{
for (int i=0;i<n;++i)
for (int j=0;j<4;++j)

for (int i=0;i<m;++i)

double last=0,Ans=calcArea(-inf)+calcArea(inf);
const int Block=4000;
double h=(inf+inf)/(double)Block;
for (int i=1;i<Block;i+=2)
Ans+=4*calcArea(-inf+i*h);
for (int i=2;i<Block;i+=2)
Ans+=2*calcArea(-inf+i*h);
Ans=Ans*h/3.0;
printf("%.3f\n",Ans);
}
return 0;
}
``````