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: This is the best and the fastest version. $n$ is gone now and the guaranteed relative precision is $1/N^2$ (the constant $1$ is correct, so if you want $10^{-3}$ accuracy (to compare with Mathematica time), just set $N=34$ and get $10^6$ pairs in under 10 seconds. The time is essentially proportional to $N$ For $10^{-5}$ accuracy $N=340$ and 83 seconds suffice. I'll explain the algorithm a bit later; now it makes sense :-)

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

const double pi=3.141592653589, ppi=pi/57.6, pi2=pi/2, dl=sqrt(0.6)/2;

double gghh(double a, double b, double c, double d, double x, double t)
{
double y=fabs(t), g=y*a, h=2*d;
if(y<=a-b) g=(a*a+y*y+b*b/3)*0.5; 
else if(y<a+b) {double z=a+b-y; g+=z*z*z/(12*b);}

y=fabs(t-x); 
if(y>c-d) h-=(y-c+d);
return g*h;
}

double F(double a, double b, double c, double d, double x, double aa, double bb)
{
double t2=(aa+bb)*0.5, bbaa=bb-aa, dt=dl*bbaa;
return bbaa*(gghh(a,b,c,d,x,t2-dt)+gghh(a,b,c,d,x,t2+dt)+1.6*gghh(a,b,c,d,x,t2))/(a*c*d);
}

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

double t0=pi2/N, cs=cos(t0), ss=sin(t0), dcs=2*cs*cs-1, dss=2*cs*ss; 
double SS=fabs(S1)+fabs(S2)+fabs(X1)+fabs(X2)+fabs(Y1)+fabs(Y2);
SS*=0.00000001;
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)+SS, fabs(Y1*ss)+SS, fabs(X2*cs)+SS, fabs(Y2*ss)+SS};
double x=-fabs(S1*cs+S2*ss);


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


double U0=U[0], U1=U[1], U2=U[2], U3=U[3];

double V[4]={-U3-U2,-U3+U2,U3-U2,U3+U2}, 
VV[4]={x-U1-U0,x-U1+U0,x+U1-U0,x+U1+U0};

double W[8]; 
int i=0, ii=0, kstart=-1, kfinish=-1;
while(ii<4)
{
++kfinish; 
if(V[i]<VV[ii]) {W[kfinish]=V[i]; ++i;}
else {W[kfinish]=VV[ii]; if(ii==0) kstart=kfinish; ++ii;} 
}

for(int kk=kstart;kk<kfinish;++kk)
s+=F(U3,U2,U1,U0,x,W[kk],W[kk+1]);
}
return ppi*s/N;
}



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


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

int N=1000;

double m=100,M=0;

for(int k=0; k<1000000;++k)
{
if(k%10000==0) {printf("%d %.12f %.12f\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);

if(k%1000==0)
{
r/=D(a1,b1,c1,d1,a2,b2,c2,d2,600);
if(r<m) m=r;
if(r>M) M=r;
}
}
printf("\n%.12f",D(1,2,3,5,4,6,7,8,4000));
printf("\n%.12f",D(1,2,3,5,4,6,7,8,N)/D(1,2,3,5,4,6,7,8,2000)-1);
printf("\n%.12f",D(0,2,0,2,0,2,0,2,N)/D(0,2,0,2,0,2,0,2,2000)-1);
printf("\n%.12f",D(0,3,0,0.0001,0,3,0,0.0001,N)-1);
printf("\n%.12f",D(0,2,0,0,0,0,0,2,N)/D(0,2,0,0,0,0,0,2,2000)-1);
printf("\n%.12f",D(0,0,0,0,3,3,4,4,N)/5-1);
return 0;
}