#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "l.h"
R det(R *, int, int);
R min(R x, R y){return x<y?x:y;}
R ig = 0;

int main(){
bp a = gb(), b = gb();
if(1) {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];
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[a.d][k]; {int i=a.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(1), (mat){(R*)cntb, b.sc, k});});
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;
{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[i*b.d+iv] + y] - bt.a[k*b.c[i*b.d] + y]; dp += k;}}
  {R r2 = 0; int y=k; while(y--){R d = cnta[i][y] - cntb[j][y]; dp[y] = d; r2 += d*d;}
    if(0) pm((mat){D, k, k}, "D");
    R rk = k&1?sqrt(r2):r2; {int q=(k-1)/2; while (q--) rk *= r2;}
    if(0){R d = det(D, k, k); intg += d/rk; printf("d = %e (%d %d)\n", d, i, j);}
    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);}}}}
  printf("Integral = %14.9f, rm = %e\n", intg, sqrt(r2m));
  // Integral =   16.679314586, rm = 8.819171e-01 looks like 16/3π = 16.755160819145562
return 0;}