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.

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

const double pi=3.141592653589;


double Abs(double x)
{
if(x<0) return -x;
return x;
}


double F(double a, double b, double c, double x, int n)
{
double s=0;
double hb=b/2/n, hc=c/2/n;
for(int kb=0; kb<=2*n;++kb)
for(int kc=0; kc<=2*n;++kc) 
{
double y=Abs(x-b/2-c/2+kb*hb+kc*hc);
double g=y;
if(y<1+a) g=( ((y+a)*(y+a)-1)/2.0+(1+a-y)/2.0+(1-(y-a)*(y-a)*(y-a))/6.0 )/2.0/a;
if(y<=1-a) g=(1+y*y+a*a/3)/2.0;

double cckb=(kb==0 || kb==2*n?1/6.0:(kb%2==0?2/6.0:4/6.0)), cckc=(kc==0 || kc==2*n?1/6.0:(kc%2==0?2/6.0:4/6.0));

s+=cckb*cckc*g/n/n;
}
return s;
}

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 t0=pi/2/N;
for(int k=0; k<N;++k)
{
double t=t0+k*pi/N; 
double x1=Abs(X1*cos(t)), y1=Abs(Y1*sin(t)), x2=Abs(X2*cos(t)), y2=Abs(Y2*sin(t)), x=Abs(S1*cos(t)+S2*sin(t));

double U[4]={x1,x2,y1,y2};

for(int kk=0;kk<2;++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;}

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

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

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


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


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,15,2);//D(a1,b1,c1,d1,a2,b2,c2,d2,48,7);
if(r<m) m=r;
if(r>M) M=r;
}

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