Return to Revisions

I tried to implement my proposal in a C-code. That is a mixture of analytic and numeric integration. It does $10^6$ rectangles with half-percent relative precision in about 16 seconds, which is a bit better than the corresponding Iosif's 30 minutes. You can play with parameters to trade speed for precision and vice versa too. The code should be self-explanatory but feel free to ask questions if something is unclear.

Edit: I optimized for speed, so now the guaranteed precision is $0.1\%$ for $10^6$ pairs in 20 seconds but the code is somewhat less readable. The average precision on random pairs of rectangles seems to be much higher but the worst case scenario is about that.

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>

const double pi=3.141592653589;




double gg(double a, double y)
{
double g;
if(y<=1-a) g=(1+y*y+a*a/3)/2.0; 
else
{
g=y;
if(y<1+a) {double z=1+a-y; g+=z*z*z/12/a;}
}
return g;
}

double F(double a, double b, double c, double x, int n)
{
double s1=0, s2=0;
double hb=(b-c)/n, hc=c/n, x1=x-c-b, x2=x-b+c, x3=x+b+c;
for(int k=0; k<=2*n;++k) 
{
double y1=fabs(x1), y2=fabs(x2), y3=fabs(x3);
double g1=gg(a,y2), g2=(gg(a,y1)+gg(a,y3))*k;
double cckc=(k%2==0? (k==0 || k==2*n?1:2):4);
x1+=hc; x2+=hb; x3-=hc;
s1+=cckc*g1; s2+=cckc*g2;
}
return (s1*(b-c)+s2*hc/2)/b/n/6;
}

double D(double a1,double b1, double c1, double d1, double a2,double b2, double c2, double d2, int N, int n)
{
double s=0.0;
double X1=b1-a1, Y1=d1-c1, X2=b2-a2, Y2=d2-c2, S1=(a2+b2-a1-b1)/2, S2=(c2+d2-c1-d1)/2;

double S=sqrt(S1*S1+S2*S2);
double t0=pi/2/N, cs=cos(t0), ss=sin(t0), dcs=2*cs*cs-1, dss=2*cs*ss; 


for(int k=0; k<N;++k)
{ 
double csnew=cs*dcs-ss*dss;
ss=ss*dcs+cs*dss; cs=csnew;
double U[4]={fabs(X1*cs), fabs(Y1*ss), fabs(X2*cs), fabs(Y2*ss)};
double x=fabs(S1*cs+S2*ss);

for(int kk=0;kk<3;++kk)
for(int j=0;j<3;++j)
if(U[j]>U[j+1]) {double u=U[j]; U[j]=U[j+1]; U[j+1]=u;}

for(int kk=0;kk<4;++kk) U[kk]+=0.0000000001*(U[3]+S); 


double S=U[3]/2; x/=S;
for(int j=0;j<4;++j) U[j]/=S; 

s+=S*F(U[2]/2,U[1]/2,U[0]/2,x,n);
}
return pi/2*s/N;
}

double unitrand()
{
return (rand()+0.0)/RAND_MAX;
}


int main()
{
time_t now=time(0);
srand(now); 

int N=50,n=4;

double m=100,M=0;

for(int k=0; k<1000000;++k)
{
if(k%10000==0) {printf("%d %f %f\n",k/10000,m,M);}
double 
a1=unitrand(),b1=a1+unitrand(),
a2=unitrand(),b2=a2+unitrand(),
c1=unitrand(),d1=c1+unitrand(),
c2=unitrand(),d2=c2+unitrand();

double r=D(a1,b1,c1,d1,a2,b2,c2,d2,N,n);

if(k%1000==0)
{
r/=D(a1,b1,c1,d1,a2,b2,c2,d2,300,15);
if(r<m) m=r;
if(r>M) M=r;
}
}

printf("\n%f %f",D(1,2,3,5,4,6,7,8,N,n), D(1,2,3,5,4,6,7,8,2000,40));
printf("\n%f %f",D(0,2,0,2,0,2,0,2,N,n), D(0,2,0,2,0,2,0,2,2000,40));
printf("\n%f %f",D(0,0,0,2,0,0,0,2,N,n), D(0,0,0,2,0,0,0,2,2000,40));
printf("\n%f %f",D(0,2,0,0,0,0,0,2,N,n), D(0,2,0,0,0,0,0,2,2000,40));
printf("\n%f %f",D(0,0,0,0,3,3,4,4,N,n), D(0,0,0,0,3,3,4,4,2000,40));
return 0;
}