#include <stdio.h>
#include <math.h>
#include <stdlib.h>
// Units: density = 1; G = 1.
typedef double R;
int sq(int x) {return x*x;}
R p1(R x, R yz){R r2 = x*x+yz; return log((sqrt(r2-2*x+1)-x+1)/(sqrt(r2)-x));}
R p2(R x, R yz){return p1((x+1)/2, yz/4);}
R p3(R L, R x, R yz){R r = 1./L; return p2(x*r, yz*r*r);}
R pste(R L, R B){R B2 = B*B, r = 2*L, r2 = r*r, bd = sqrt(2*B2 + r2);
  return 4*(B2*(log(B*sqrt(2)/(bd-r)) - atan(r/bd)) + B*r*(log((B+bd)/(bd-B)))
         - 0.5*r2*atan(B2/(r*bd)));}
#define lr 20
#define lm (2*lr+1)
const R d2=1/(R)(lr*lr);
int main(){R L[lr+1][lm], P[lr+1][lm];
R pma=0, pmi=1.e10, m0 = 0, m1 = 0, m2 = 0; int ima, jma, imi, jmi;
for(int i=0; i<lr+1; ++i) for(int j=0; j<lm; ++j)
   L[i][j] =({R z = 1.- (i*i + sq(j-lr))*d2; z<0?0:sqrt(z);});
for(int i=0; i<lr+1; ++i) for(int j=0; j<lm; ++j)
    if(L[i][j]>0 || i>0 && L[i-1][j]>0 || i<lr && L[i+1][j]>0 
      || j>0 && L[i][j-1]>0 || j<lm-1 && L[i][j-1]>0) {R p = 0;
  for(int k=0; k<lm; ++k) {const int mir = abs(k-lr);
    for(int l=0; l<lm; ++l) if(L[mir][l]>0) {
      if(i!=k-lr || j!=l) p += p3(L[mir][l], L[i][j], (sq(i+lr-k)+sq(j-l))*d2);
      else p += pste(2*L[i][j], 0.5/lr)/d2;}}
  {R q = p*d2; P[i][j] = q; if(L[i][j]>0){
    if(q>pma) {pma=q; ima=i; jma=j;}
    if(q<pmi) {pmi=q; imi=i; jmi=j;}
    m0 += L[i][j]; m1 += L[i][j]*i; m2 += L[i][j]*(i*i+j*j);}}}
printf("%d by %d grid\n", lr+1, lm);
printf("max %3d %3d %17.14f\n", ima, jma, pma);
printf("min %3d %3d %17.14f\n", imi, jmi, pmi);
printf("xx %17.14f %17.14f\n", pste(0.7, 0.5/lr)/d2, p3(0.35, 0.35, d2));
printf("m0 = %17.14f, m1 = %17.14f, m3 = %17.14e\n", m0*d2, m1*d2, m2*d2);
return 0;}
/*  (4/3)π = 4.1887902047863905
80 by 80 grid
max   4  77  4.18853473663330
min  10   1  4.18799448013306

real	0m6.944s
user	0m6.935s
sys 	0m0.005s
*/