#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "l.h"
#define pdm 0
#define mc 1
R min(R x, R y){return x<y?x:y;}
R const pi = 3.1415926535897932385;
R fac(R n){if(n == 0) return 1; if(n==-0.5) return sqrt(pi); return n*fac(n-1);}

int main(){
bp a = g1(), b = g1();
if(0) {void pcm(bp x, char c) {int i, j; for(i=0; i<x.sc; ++i) {
  for(j=0; j<x.d; ++j) printf("%3d ", x.c[i*x.d + j]); printf(" %c[%d]\n", c, i);}}
  pcm(a, 'a'); pcm(b, 'b');}
int const k = a.d + b.d - 1;
R off[k]; {int j=k; while(j--) off[j] = 0;}
R* vog(R x) {off[0] = x; return off;}
mat ob(mat m, bp x, R * vo, mat cen){
   R * a = (R*)malloc(sizeof(R)*x.vc*k);
   mm((mat){a, x.vc, k}, (mat){x.wh, x.vc, x.d}, m);
   {int j=x.vc; while(j--) {int m=k; while(m--) a[j*k+m] += vo[m];}}
   {int i=x.sc; while(i--){int j=k; while(j--){R s=0; int n=x.d; while(n--){
        s += a[k*(x.c[x.d*i + n]/*which vertex*/) + j];}
      cen.a[i*cen.n + j] = s/x.d;}}}
   return (mat){a, x.vc, k};}
R cnta[a.sc][k], cntb[b.sc][k];
void di(R dc){mat at = ({R tr[a.d][k]; {int i=a.d; while(i--){int j=k; while(j--) tr[i][j] = i==j;}}
     ob((mat){(R*)tr, a.d, k}, a, vog(0), (mat){(R*)cnta, a.sc, k});});
mat bt = ({R tr[b.d][k]; {int i=b.d; while(i--){
     int j=k; while(j--) tr[i][j] = (!i && !j) || i && (i+a.d-1==j);}}
     ob((mat){(R*)tr, b.d, k}, b, vog(dc), (mat){(R*)cntb, b.sc, k});});
if(0){pm(at, "at"); pm(bt, "bt");}
if(0){pm((mat){(R*)cnta, a.sc, k}, "cnta"); pm((mat){(R*)cntb, b.sc, k}, "cntb");}

R D[k*k]; R intg = 0, r2m = 1e9;
R pw(R r2){R rk = k&1?sqrt(r2):r2; {int q=(k-1)/2; while (q--) rk *= r2;} return rk;}
{int i=a.sc; while(i--){int j=b.sc; while(j--){R * dp = D;
  {int iv=a.d; while(iv-->1) {
      int y=k; while(y--) dp[y] = at.a[k*a.c[i*a.d+iv] + y] - at.a[k*a.c[i*a.d] + y]; dp += k;}}
  {int iv=b.d; while(iv-->1) {
      int y=k; while(y--) dp[y] = bt.a[k*b.c[j*b.d+iv] + y] - bt.a[k*b.c[j*b.d] + y]; dp += k;}}
  if(mc){{int y=k; while(y--) dp[y] = at.a[k*a.c[i*a.d] + y] - bt.a[k*b.c[j*b.d] + y];}
     int const sz = 1000;
     R ave = 0; {int sn=sz; while(sn--)
        if(0){ave += 1./pw(smp((R*)D, a.d-1, b.d-1));} /// arghh!! .... double printing alternative
        else {R p = smp((R*)D, a.d-1, b.d-1); if(0) printf("vne r=%8.5f\n", sqrt(p)); ave += pw(1./p);}
     if(pdm){pm((mat){D, k, k}, "D"); R d = det(D, k, k); 
       printf("d=%e ave=%e, sz=%d\n", d, ave, sz); intg += d*ave/sz;}
     else intg += det(D, k, k)*ave/sz;}}
  else
    {R r2 = 0; {int y=k; while(y--){R dx = cntb[j][y] - cnta[i][y]; dp[y] = dx;
       if(0) printf("fwhh dx=%e cntb[%d][%d] = %e "
           "cnta[%d][%d] = %e dp-D=%d\n", dx, j, y, cntb[j][y], i, y, cnta[i][y], dp-D);
        r2 += dx*dx;}}
      R rk = pw(r2);
      if(pdm){pm((mat){D, k, k}, "D");
        R d = det(D, k, k); intg += d/rk; printf("(%d %d) d = %e rk=%e r2=%e\n",
           i, j, d, rk, r2);}
      else intg += det(D, k, k)/rk;
      r2m = r2<r2m?r2:r2m;
      if(0) printf("r2=%e rk=%e intg=%e\n", r2, rk, intg);}}}
  {R t = a.d+b.d-1, th = dc<2?t*pow(pi, t/2)/fac(t/2):0;
  printf("%d-brane and %d-brane dc=%4.1f %s:\n"
          "calc = %14.9f, theory = %14.9f, rm = %6.4f\n",
     a.d-1, b.d-1, dc, mc?"mc":"ctr", intg/(fac(a.d-1)*fac(b.d-1)), th, sqrt(r2m));}}}
di(1); di(3); return 0;}