#include <stdio.h>
#include <math.h>
typedef double R;
#define lc 3000
R pi = 3.14159265358979323846;
typedef struct{R x; R y; R z;} v;
R sq(R x){return x*x;}
v sp(R s, v vv){return (v){s*vv.x, s*vv.y, s*vv.z};}
R det(v x, v y, v z){return
   x.x * (y.y * z.z - y.z * z.y)
 + x.y * (y.z * z.x - y.x * z.z)
 + x.z * (y.x * z.y - y.y * z.x);}
v vd(v * a, v * b){return (v){a->x - b->x, a->y - b->y, a->z - b->z};}
v va(v a, v b){return (v){a.x + b.x, a.y + b.y, a.z + b.z};}
void dcir(v ans[lc+1], v c, v r1, v r2){
R da = 2.*pi/lc;
for(int k=0; k<lc; ++k)
   ans[k] = va(c, va(sp(sin(da*k), r1), sp(cos(da*k), r2)));
ans[lc] = ans[0];}
v l1[lc+1], l2[lc+1]; // places on the respective curves.
int main(){
  dcir(l1, (v){0,   0, 0}, (v){1, 0, 0}, (v){0, 1, 0});
  dcir(l2, (v){1,   0, 0}, (v){1, 0, 0}, (v){0, 0, 1});

// printf("det=%e\n", det((v){0,2,0}, (v){0,0,1}, (v){1,0,0}));
{R s = 0, rm = 1000.;
for(int j=0; j<lc; ++j) for(int k=0; k<lc; ++k) {
   v r = vd(&l1[j], &l2[k]);
   R v2 = sq(r.x) + sq(r.y)+ sq(r.z);
   s += det(vd(&l1[j+1], &l1[j]), vd(&l2[k+1], &l2[k]), r)/(sqrt(v2)*v2);
   if(v2<rm) rm=v2;}
printf("s = %11.7f, 4π = %11.7f, rm=%e\n", s, 4.*pi, sqrt(rm));}
return 0;}
// s = -12.5663522, 4π =  12.5663706, rm=1.000000e+00